Original Research ARTICLE
Cerebral Autoregulation Evidenced by Synchronized Low Frequency Oscillations in Blood Pressure and Resting-State fMRI
- 1Cardiff University Brain Research Imaging Centre (CUBRIC), School of Physics and Astronomy, Cardiff University, Cardiff, United Kingdom
- 2CUBRIC, School of Psychology, Cardiff University, Cardiff, United Kingdom
- 3Department of Physical Therapy and Human Movement Sciences, Feinberg School of Medicine, Northwestern University, Chicago, IL, United States
Resting-state functional magnetic resonance imaging (rs-fMRI) is a widely used technique for mapping the brain’s functional architecture, so delineating the main sources of variance comprising the signal is crucial. Low frequency oscillations (LFO) that are not of neural origin, but which are driven by mechanisms related to cerebral autoregulation (CA), are present in the blood-oxygenation-level-dependent (BOLD) signal within the rs-fMRI frequency band. In this study we use a MR compatible device (Caretaker, Biopac) to obtain a non-invasive estimate of beat-to-beat mean arterial pressure (MAP) fluctuations concurrently with rs-fMRI at 3T. Healthy adult subjects (n = 9; 5 male) completed two 20-min rs-fMRI scans. MAP fluctuations were decomposed into different frequency scales using a discrete wavelet transform, and oscillations at approximately 0.1 Hz show a high degree of spatially structured correlations with matched frequency fMRI fluctuations. On average across subjects, MAP fluctuations at this scale of the wavelet decomposition explain ∼2.2% of matched frequency fMRI signal variance. Additionally, a simultaneous multi-slice multi-echo acquisition was used to collect 10-min rs-fMRI at three echo times at 7T in a separate group of healthy adults (n = 5; 5 male). Multiple echo times were used to estimate the R2∗ decay at every time point, and MAP was shown to strongly correlate with this signal, which suggests a purely BOLD (i.e., blood flow related) origin. This study demonstrates that there is a significant component of the BOLD signal that has a systemic physiological origin, and highlights the fact that not all localized BOLD signal changes necessarily reflect blood flow supporting local neural activity. Instead, these data show that a proportion of BOLD signal fluctuations in rs-fMRI are due to localized control of blood flow that is independent of local neural activity, most likely reflecting more general systemic autoregulatory processes. Thus, fMRI is a promising tool for studying flow changes associated with cerebral autoregulation with high spatial resolution.
Functional connectivity in the brain can be assessed with blood-oxygenation-level-dependent (BOLD) functional magnetic resonance imaging (fMRI). The source of BOLD contrast is the difference in magnetic susceptibility between oxy- and deoxyhemoglobin, which has an effect on apparent transverse relaxation (R2∗), and thus imparts sensitivity to blood oxygenation in the MR signal (Buxton, 2013). Neurovascular coupling (NVC) allows brain activity to be mapped using BOLD fMRI, because localized increases in cerebral blood flow (CBF), which are proportionally larger than changes in oxygen metabolism, cause increases in local venous blood oxygenation. An implicit assumption that predicates BOLD fMRI as a tool for mapping neural activity in the brain is that NVC related changes in CBF are the predominant source of signal variance. There are, however, other mechanisms besides NVC that regulate CBF, such as arterial blood gas concentration, particularly carbon dioxide (CO2), which is a potent vasodilator with a strong effect on CBF (Battisti-Charbonney et al., 2011). Furthermore, systemic control of the brain’s blood supply is governed by numerous homeostatic mechanisms that are broadly defined as cerebral autoregulation (CA) (Willie et al., 2014), the theoretical process that modulates cerebrovascular resistance to ensure CBF is kept at a sufficient level in the face of transient changes in systemic haemodynamics (e.g., blood pressure and cardiac output).
Understanding non-neuronal sources of variance in CBF fluctuations is especially important with regard to resting state fMRI (rs-fMRI) paradigms for two reasons. Firstly, unlike traditional task based paradigms for which the timing and duration of evoked BOLD signal changes is known a priori, the timing of spontaneous neural fluctuations can’t be assumed, meaning non-neuronal effects can’t be mitigated as they are in task based designs by averaging over trials. Secondly, the low frequency range (∼0.01–0.1 Hz) over which functional connectivity is observed overlaps with the spectrum at which other systemic physiological effects occur (Murphy et al., 2013). Spontaneous fluctuations in breathing cause cerebrovascular reactivity (CVR) to CO2 to manifest as low frequency (<0.05 Hz) oscillations in the BOLD signal (Wise et al., 2004), and endogenous fluctuations (<0.1 Hz) in vascular tone have been reported in various different vascular beds across multiple species (Nilsson and Aalkjaer, 2003). Recently, low frequency oscillations (LFO) of a systemic origin have been observed by correlating fMRI signals with functional near-infrared spectroscopy (fNIRS) measures of peripheral haemodynamics (Tong and Frederick, 2010; Tong et al., 2012). Intriguingly, these systemic LFOs appear to propagate throughout the brain with spatially structured temporal delays (Erdogan et al., 2016).
Arterial blood pressure (ABP) is dynamic over multiple time scales, including at a beat-to-beat level, and so is likely to contribute significantly to fluctuations in CBF. Transcranial Doppler ultrasound (TCD) studies have consistently demonstrated how ABP fluctuations modulate cerebral blood flow velocity (CBFV) in large intracranial arteries (Aaslid et al., 1989; Zhang et al., 1998), and that they account for a considerable proportion of the total variance, approximately 60% of the total predictive power of CBFV fluctuations in right middle cerebral artery (MCA) (Mitsis et al., 2004). Evidence for how blood pressure dynamics affect fMRI fluctuations is scarce, mostly limited to animal studies on the relationship between transient changes and evoked neural responses (Wang et al., 2006; Qiao et al., 2007; Uchida et al., 2017). Similar to fMRI, ABP time series have a 1/f power spectrum, but also show distinct oscillations (∼0.1 Hz in humans) known as Mayer waves (Mayer, 1876), which are independent of respiration and tightly coupled to efferent sympathetic nervous activity (SNA) (Julien, 2006). Oscillations at this frequency have also been observed in cerebral haemodynamics measured with fNIRS (Obrig et al., 2000; Yucel et al., 2016) and intraoperative multispectral optical intrinsic signal imaging (Rayshubskiy et al., 2014). However, the origins of such signals are unclear, and separating the effects of ABP fluctuations from vasomotion (which is usually regarded as distinct) on cerebral haemodynamics is an open challenge.
Nevertheless, the TCD literature provides compelling reason to believe that ABP fluctuations should contribute to the BOLD fMRI signal. Measurement of the coupling between fluctuations in ABP and CBFV in intracranial arteries has found widespread use as a clinically useful means of assessing CA (Zhang et al., 1998; Panerai et al., 2002; Panerai, 2008), and nonlinear modeling estimates that ABP accounts for 60% of the predictive power of CBFV fluctuations in the MCA (Mitsis et al., 2004). Although TCD has been widely used to assess both CVR and CA in research and clinical practice (Willie et al., 2011), more recently fMRI has emerged as a powerful tool for mapping CVR across the brain (Pillai and Mikulis, 2015), and the feasibility of obtaining BOLD fMRI based measures of CVR from spontaneous CO2 fluctuations has also been demonstrated (Golestani et al., 2016). So far, these advances in measuring cerebrovascular function with fMRI have not yet extended into the domain of CA. However, the TCD literature proves that blood pressure related spontaneous CBFV fluctuations provide an effective means of characterizing CA, which is promising for the development of an equivalent whole-brain fMRI method.
In this study we explore the relationship between systemic fluctuations in blood pressure and the resting-state fMRI signal. We measure beat-to-beat blood pressure fluctuations concurrently with single-echo fMRI at 3T and multi-echo fMRI at 7T, and show that widespread patterns of correlations exist in low frequency BOLD signals, which we posit are due to fluctuations in CBF associated with CA.
Materials and Methods
Magnetic Resonance Imaging Acquisition
The study consisted of two separate experiments conducted on two different scanners. Nine healthy volunteers (age 22–37 years) were recruited for a 3T session to collect single-echo fMRI data (3T) and five additional healthy volunteers (age 30–41 years) were recruited for a 7T session to collect multi-echo fMRI data (7T-ME). All participants gave written informed consent, and the School of Psychology Cardiff University Ethics Committee approved the study in accordance with the guidelines stated in the Cardiff University Research Framework (version 4.0, 2010).
The 7T scan protocol was added to enable us to address the origin of MAP correlated fMRI signal changes. Primarily this was achieved by using a multi-echo acquisition, which allows us to separate the pure BOLD component from the fMRI signal. Furthermore, the session consisted of two scans with differing acquisition parameters (see section “Multi-Echo Fit” below for theoretical details).
Two twenty-minute rs-fMRI runs were acquired on a 3T GE HDx scanner (GE Healthcare, Milwaukee, WI, United States) with an eight-channel receive head-coil using a gradient-echo EPI readout with a single echo time (TR = 2000 ms; TE = 35 ms; flip angle (α) = 90°; FOV = 224 mm; 3.5 mm2 in-plane resolution; 33 slices (3.5 + 0.5 mm gap), SENSE (GE ASSET) acceleration factor = 2). Whole-brain T1-weighted anatomical images were acquired using an FSPGR sequence (FOV = 256 mm, TR = 7900 ms, TE = 3 ms, 172 contiguous sagittal slices, 1 mm3 isotropic).
Two ten-minute eyes-closed rs-fMRI runs were acquired on a 7T Siemens MAGNETOM scanner (Siemens Healthcare GmbH, Erlangen, Germany) equipped with a single-channel transmit/32-channel receive head coil (Nova Medical, Wilmington, MA, United States). The CMRR SMS-EPI sequence (R015) (Moeller et al., 2010) was used to acquire multi-echo multiband EPI data with three echoes using the following parameters: Scan 1 – [TR = 1000 ms; TE1/2/3 = 8.14/21.47/34.8 ms; flip angle (α) = 35°; FOV = 220 mm; 2.4 mm2 in-plane resolution; 36 slices (2.5 mm thick); multiband factor = 4; GRAPPA acceleration factor = 2] and Scan 2 – [TR = 500 ms; TE1/2/3 = 8.14/21.47/34.8 ms; flip angle (α) = 90°; FOV = 220 mm; 2.4 mm2 in-plane resolution; six slices (2.5 mm thick); multiband factor = 1; GRAPPA acceleration factor = 2]. Whole brain T1-weighted anatomical images were acquired using an MPRAGE sequence (FOV = 220 mm, TR = 2200 ms, TE = 3 ms, TI = 1050 ms, 224 contiguous sagittal slices, 0.7 mm3 isotropic).
Concurrent physiological traces were recorded for all runs and sampled at 500 Hz (CED, Cambridge, United Kingdom). This included using photoplethysmography (PPG) to measure pulse waveforms for deriving cardiac information, a pneumatic respiratory belt for timing and relative respiration volume measures, capnography for measuring expired partial pressure of end-tidal carbon dioxide (PETCO2). The CareTaker system (Biopac) was used to measure beat-to-beat blood pressure with a cuff attached to the first digit of the hand (thumb). The system uses the cuff to pneumatically sensor the arterial pressure wave, and estimates beat-to-beat systolic and diastolic blood pressure via analysis of the timing between different components of the pulse waveform (Baruch et al., 2011), and has been validated against gold standard arterial line measurements (Baruch et al., 2014).
Data were preprocessed and registered to a standard space using a pipeline created with AFNI, FSL, and in-house code. Preprocessing of 3T and 7T-ME data consisted of the same following steps: (1) De-spiking; (2) Motion correction by registering all volumes to the first one. For 7T-ME scans steps 1–2 were performed separately for each echo time dataset, then a nonlinear fit was performed to create S0 and R2∗ datasets (see section “Multi-Echo Fit” below). Subsequent steps were performed separately for S0 and R2∗ datasets. (3) Nuisance regression with pre-whitening (Bright et al., 2017) to remove cardiac and respiratory related noise (Glover et al., 2000; Birn et al., 2008; Chang et al., 2009), end-tidal CO2 fluctuations (convolved with HRF), and six estimated motion parameters; (4) Slice time correction; (5) Non-linear registration to 2 mm MNI space; (6) De-trending and motion censoring in a single step, with the top 5% of volumes most severely corrupted by motion (according to framewise displacement) being censored. Censored time points were replaced with interpolated values calculated from neighboring (non-censored) time points (NTRP option in 3dTproject) in order to keep to the data temporally consistent for subsequent wavelet decomposition. A discrete wavelet transform was then performed on the preprocessed data (see section “Maximum Overlap Discrete Wavelet Transform” below). Note that physiological noise correction was performed on unfiltered data.
For each subject gray matter (GM) masks were created from segmented T1 images, with GM voxels defined as those with a partial volume estimate greater than 66%. GM masks were used in subsequent parts of the analysis, and GM mask averaged time series were calculated for 3T and 7T-ME data for estimating the global lag with blood pressure (see section “Blood Pressure Correlation” below).
Assuming a mono-exponential decay, the signal across multiple echo times can be described according to Eq. 1.
Where S0 is spin density weighted signal intensity at zero echo time and R2∗ is the apparent transverse relaxation rate (inverse of relaxation time T2∗). S0 is modulated by changes in apparent T1 (e.g., due to inflow) and bulk motion and related spin history effects, whereas R2∗ reflects magnetic field homogeneity and thus, due to blood oxygenation induced changes in microscopic susceptibility, the source of the BOLD effect (Buxton, 2013). NVC related functional responses in gradient-echo fMRI are considered to be driven almost entirely by R2∗ changes, which has motivated the use of multi-echo acquisitions for separating neuronal from non-neuronal signal components (Posse, 2012; Kundu et al., 2017), based on the rationale that non-neuronal components (i.e., not flow related) are mostly restricted to changes in S0. It should be noted that inflow effects on S0 are driven by changes in flow velocity through an imaging slice, and so may partially reflect changes in CBF, however, this effect is generally considered to be small in multi-slice acquisitions with standard TRs (Gao and Liu, 2012).
However, as discussed above, there are non-neuronal contributors to flow, whose activity in principle should manifest primarily in the form of R2∗ changes. Ignoring the effects of macroscopic field inhomogeneity, more generally R2∗ changes are driven by CBF dynamics (via the changes in blood oxygenation they produce), whether they are neuronally driven or not. The motivation for collecting 7T-ME data as part of this study is that it allows us to determine to what extent MAP-fMRI correlations are driven by changes in R2∗ and S0, and thus whether or not they are related to changes in CBF. Moreover, as S0 is determined by the steady-state longitudinal magnetization it is intrinsically dependent on acquisition parameters like TR and flip angle. Thus, any significant S0 fluctuations are likely to be modulated between 7T-ME scans 1 and 2, whereas R2∗ fluctuations will not.
For each voxel and TR, 7T-ME data were fit to the mono-exponential signal model with a nonlinear least-squares approach using the Levenberg-Marquardt algorithm (Galassi, 2009), creating S0 and R2∗ datasets. In multi-echo fMRI studies numerous physical limitations restrict the number of echo times that can be achieved to a small number, in this case three, which presents a challenge for accurate estimation of relaxation rates (Gowland and Bowtell, 2007), and sample-by-sample parameter estimates are considered to be noisy (Kundu et al., 2012). Compared with previous studies, here we benefit from the higher SNR afforded by 7T to improve sample-by-sample parameter estimates, and we performed simulations to better understand the precision with which R2∗ and S0 can be measured using our nonlinear fit approach (details included in Supplementary Material). These simulations demonstrate that across the expected range of R2∗ values, our choice of echo times allows us to estimate both R2∗ and S0 parameters without bias.
Maximum Overlap Discrete Wavelet Transform
Wavelet transforms provide a way of decomposing the total variance within a time series into different frequency scales [see Bullmore et al. (2004) for an fMRI focused review]. They are conceptually similar to Fourier transforms, but because wavelets are compactly supported (i.e., transient, not extending infinitely like sine waves), they provide sensitivity to non-stationary features within the scales of the decomposition. Thus, wavelet coefficients (WC) provide a “time-frequency” representation of data, including both temporal and spectral information, analogous to moving-window Fourier transforms. This makes the discrete wavelet transform useful for analyzing real-world physiological data, which are expected to be non-stationary, and the multi-resolution analysis allows signal energy to be decomposed into distinct frequency bands. Given that fluctuations in MAP are expected to show more power in certain frequency bands, likely reflecting different underlying mechanisms, in this study we used a wavelet transform to identify the frequency scale of interest for the relationship between MAP and fMRI signals.
For each voxel time series the maximum overlap discrete wavelet transform (MODWT) was used to decompose the signal into six scales (Witcher, 2015). The MODWT used a fourth-order Daubechies wavelet filter as has been used previously for fMRI applications (Bullmore et al., 2001; Patel and Bullmore, 2016). The central frequency (fc) and band between lowest (flow) and highest (fhigh) frequencies contained within each scale (j) depends on the sampling rate (TR), and is given by Eqs 2 and 3.
In summary, each fMRI time series consisting of N volumes was decomposed into six frequency scales, each composed of N WCs. For reference the scale frequency bands for the TRs used in this study are given in Table 1 and the frequency response of the filters at each scale are shown in Supplementary Figure 1.
Blood Pressure Correlation
Beat-to-beat systolic (SBP) and diastolic (DBP) blood pressure time series were processed in-house with a robust outlier removal procedure and resampled to the relevant TR. Mean arterial pressure (MAP) was estimated according to Eq. 4 (Brzezinski, 1990).
MAP time series were decomposed with a MODWT into the same six frequency scales as the fMRI data. Figure 1 shows a MAP time series and its WCs for a representative subject. For each run the cross-correlation between the average gray matter signal and temporally displaced MAP signals was calculated. The raw cross-correlation functions were then fit to a set of Legendre polynomial functions of increasing order (until R2 > 0.95, or up to a maximum order of 10) to obtain a smooth function from which the lag between MAP and fMRI fluctuations was obtained. Voxelwise correlations between MAP and scale matched fMRI WCs for all frequency scales were calculated with the globally optimized lag. Additionally, although more susceptible to noise, we calculated a voxelwise lag using the same procedure and then calculated corresponding correlations between MAP and scale matched fMRI WCs.
Figure 1. A (i) Shows the mean gray matter signal from 3TSE data in a representative subject along with corresponding mean arterial pressure trace (iii). Wavelet coefficients at 6 scales are shown for both GM (ii) and MAP (iv). (B) Group mean cross-correlation function between MAP and GM for each scale (i–vi), with shaded area representing SEM.
It should be noted that in contrast with the discrete wavelet transform (DWT), which is orthogonal (like the discrete Fourier transform), the MODWT is an oversampled transform in which there is some redundancy (i.e., coefficients are not completely independent). This means that the effective degrees of freedom are reduced, which impacts on statistical parameters. Additionally, the practice of filtering time series and allowing temporal shifts changes the null distribution of Pearson’s correlations (Bright et al., 2017). To address these potential confounds on statistical inference we performed a permutation analysis to determine the correct null hypothesis. The approach, which has been used previously (Bright et al., 2017), estimates the null distribution with phase-randomized versions of the MAP time series. For each subject 1000 phase-randomized MAP traces were correlated with fMRI to get a voxel-wise estimate of the correlation null distribution from which p-values were calculated.
Group Level Analysis
For each frequency scale, subject level MAP - fMRI WC correlation maps were entered into an independent two-tailed t-test, from which group-level correlation and Z-score maps were derived. The mean GM correlation and Z-score values were calculated as a means of identifying the scale with the strongest MAP vs. fMRI coupling. Subsequent analyses of 7T-ME data are restricted to the scale with the matched frequency. Test–retest repeatability was assessed using spatial correlation between Scan 1 and Scan 2 for 3T. For 7T-ME data, Scan 2 parameters were chosen differently from Scan 1 from a theoretical perspective to maximize any potential inflow effect in the S0 signal (Gao and Liu, 2012). For the different sources of image contrast, R2∗ and S0, to assess their relative contributions to MAP – fMRI correlations, the absolute value of average GM vs. MAP correlations were entered into a two way repeated measures ANOVA.
Note, that although a wavelet transform was used to decompose MAP and fMRI signals into different frequency scales, we also looked at the effect of MAP on unfiltered data by regressing unfiltered optimal lag MAP traces onto voxel-wise unfiltered fMRI data (see Supplementary Figure 4).
Figure 1A, shows an example taken from a representative subject of GM and MAP traces, and their respective WCs at the six frequency scales listed in Table 1. Figure 1B shows the group mean cross-correlation functions between MAP and GM WC time series for each scale of the MODWT. Scales 1–3 all show clear maxima with a similar degree of lag (5.75, 5.50, and 5.25 s for scales 1–3, respectively).
MAP vs. fMRI Correlation
Voxelwise group average correlations (at optimal lag) are shown for scale 2 WCs only in Figure 2C, along with corresponding Z-scores in Figure 2D. Figure 2A shows the mean correlation values and Z-scores within the voxelwise group maps for all scales (at optimal lag), and it can be seen the scale 2 has the highest mean correlation value and highest mean Z-score value. It can be seen by the standard deviation error bars that the majority of voxels at scale 2 have a Z-score >3.1 (p = 0.001). Figure 2B demonstrates that the spatial correlation between scans (i.e., within subject agreement) is also highest for scale 2 compared with the other scales. As documented in Table 1, the frequency band for scale 2 for 3TSE data (TR = 2s) is 0.063–0.125 Hz, which corresponds to a central frequency () of ∼0.1 Hz.
Figure 2. 3TSE data. (A) Mean correlations and Z-scores (across voxels) from group mean statistical parameter maps for each scale. Dotted lines indicated the p-value for give z-score vales. (B) Spatial correlation scans 1 and 2 MAP-fMRI correlation maps for each scale (∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001, and ∗∗∗∗p < 0.0001). (C) Group mean MAP – fMRI correlation map for scale 2 WC and corresponding Z-scores. (D) Correlations are those at the optimal lag values as shown in Figure 1.
As stated in Section “Maximum Overlap Discrete Wavelet Transform”, permutation tests based on phase-randomized MAP time series were used to estimate the correlation null distribution for each subject on a voxel-wise basis. Supplementary Figure 3A shows the correlation null distributions estimated for each subject, and the mean GM correlation value. For each subject the null distribution is non-zero, as expected due to the effect of temporal shifting, but in every case the true mean GM correlation with MAP is more than three standard deviations removed from the null correlation. Supplementary Figure 2 shows subject level correlation maps along with associated threshold permutation test p-values.
For reference, Supplementary Figure 4 shows the estimated effect size of unfiltered MAP fluctuations on fMRI data. The average GM effect size across subjects (±SD) is 0.01% BOLD/mm Hg (±0.006). Across subjects the average absolute maximum deviation in MAP is ∼12 mm Hg, which suggests fairly modest total BOLD signal changes on the order of ∼0.12% are expected, but given the heterogeneity evident in Supplementary Figure 4, it is clear that in some regions this may be as large as ∼0.5%.
Figure 3A shows group level MAP – fMRI correlations as a function of lag with respect to MAP. It shows how a spatially structured pattern of fMRI signal changes evolves over time in response to MAP fluctuations. Figures 3B,C shows the maximum correlation and the lag at which it is seen, respectively. The lag time in cortical gray matter appears relatively uniformly distributed at ∼5 s, in good agreement with the GM signal global lags shown in Figure 1. Interestingly, there is a correlation pattern that emerges earlier (∼2–4 s), which appears in deep white matter structures and in the areas bordering the lateral ventricles. Figure 3B shows that there are widespread correlations with MAP, albeit with different lags, extending across the whole brain.
Figure 3. 3TSE data. (A) MAP – fMRI as a function of lag with respect to MAP, in 2 s intervals. (B) The maximum correlation (i.e., arg max of cross-correlation function) and the associated lag time (C).
The 7T-ME data allows us to tease apart the different sources of contrast underlying the MAP – fMRI correlations. Figure 4A shows group level voxelwise MAP – R2∗ and MAP – S0 correlations for Scan 1 and Scan 2. R2∗ – MAP correlations show a similar pattern to 3T fMRI – MAP correlations, with well defined gray/white matter contrast and matching areas of high magnitude correlations (e.g., in occipital cortex). Note that negative R2∗ – MAP correlations are equivalent to positive MAP – fMRI (3T) correlations, as a decrease in R2∗ corresponds to a lengthening of T2∗ and a positive increase in BOLD signal. Figure 4B shows the spatial correlations between R2∗ – MAP and S0 – MAP correlation maps, and 3T BOLD – MAP correlation maps. Compared with S0 – MAP correlations, R2∗ – MAP correlation maps are more spatially similar to 3T BOLD – MAP correlation maps, with a Pearson’s correlation of −0.68 vs. 0.33 of S0 – MAP, which amounts to ∼4 times as much variance explained. Note that the negative correlation in Figure 4Bi is due to the inverse relationship between R2∗ and the BOLD signal (a BOLD signal increase results from less dephasing, i.e., a decrease in relaxation rate).
Figure 4. 7TME data. (A) group level MAP – fMRI correlation maps for R2∗ and S0 and scans 1 and 2, at the frequency scale corresponding to 0.063 – 0.125 Hz (scales 3 and 4 for scans 1 and 2 respectively). (B) (i) Spatial correlation between 7TME MAP – R2∗ correlation maps, and 3TSE MAP – fMRI correlation maps. (ii) Spatial correlation between 7TME MAP – S0 correlation maps, and 3TSE MAP – fMRI correlation maps. (C) Bar chart showing group mean GM correlations (absolute value) for R2∗/S0 and scans 1 and 2 (∗p < 0.05, ∗∗p < 0.01).
Figure 4C shows the group mean GM absolute correlation values for R2∗ and S0 scans 1 and 2. A two-way repeated measures ANOVA revealed a significant effect of contrast (R2∗ > S0), but no effect of scan. Following the rationale outlined in Section “Multi-Echo Fit,” this would suggest that the MAP correlated fMRI signal has a BOLD origin related to changes in CBF. Furthermore, there appears to be minimal non-BOLD contribution, as scan number did not significantly modulate the MAP-S0 correlation values.
Blood Pressure Correlation
To our knowledge, this study is first to demonstrate that MAP LFOs are positively correlated with fMRI LFOs within the frequency band between 0.063 and 0.125 Hz. These correlations appear highly spatially structured, with strong gray/white matter contrast, and are repeatable between subjects with a spatial correlation of ∼0.42. Results from the 7T-ME data suggest that fluctuations in MAP lead to gray matter signal fluctuations in BOLD fMRI that are primarily related to CBF, given that they are related to changes in R2∗ and relatively independent of acquisition parameters. This is consistent with a large TCD literature that shows beat-to-beat fluctuations in blood pressure result in measurable changes in CBFV in large intracranial arteries (Aaslid et al., 1989; Diehl et al., 1991; Blaber et al., 1997; Kuo et al., 1998; Zhang et al., 1998). The same studies also make it clear that there is a lag in CBFV fluctuations with respect to MAP (i.e., the CBFV response occurs later), which has been estimated at ∼2 s (Kuo et al., 1998). The delay time ∼5.5 s estimated from these data is markedly longer, but this is perhaps not surprising given that BOLD fMRI is sensitive to flow in a different vascular compartment (i.e., deoxygenated/venous blood flow). An increased lag time would be expected to account for the fact the flow changes will take time to propagate along the vascular tree. Delayed fMRI responses to hypercapnia challenges are frequently observed on the order of 8–15 s (Blockley et al., 2011; Murphy et al., 2011), although potentially longer in patient groups (Duffin et al., 2015; Donahue et al., 2016), and are presumed to contain both gas bolus transit time and vascular reactivity information. Further research into the causal relationship underlying MAP correlated fMRI signals is needed in order to better interpret what this delay time reflects.
Somewhat more difficult to explain is the apparent shorter lag time in white matter, in particular bordering the ventricles, as seen in the voxel-wise lag time in Figure 3. Assuming a purely flow based explanation, one would expect a longer delay time in white matter, as blood arrival time is extended compared to gray matter (van Gelderen et al., 2008). However, this finding is consistent with the anti-correlated signal changes observed during in the same areas during respiratory challenges (Bright et al., 2014). It suggests that in regions located at the edge of CSF spaces, signal changes may be dominated by fluctuations in CSF partial volume as opposed to flow related (BOLD) changes. This interpretation is supported by Figure 3A, in which it can be seen that that the earliest signal changes occur in the regions proximal to the ventricles. This is expected according to the Monroe-Kellie doctrine, i.e., within a rigid skull volume changes in different fluid compartments must be balanced. Thus, it is possible that changes in arterial blood volume preceding the BOLD signal change are associated with matched changes in CSF volume, although this is only speculative and requires further evidence.
From these data we cannot know the exact origin of these BP correlated fMRI fluctuations. Previous fMRI studies have found signal fluctuations that are correlated with peripheral measures of vascular tone, such as NIRS in fingers/toes (Tong and Frederick, 2010; Tong et al., 2012, 2013) or amplitude of photoplethysmography (PPG) (van Houdt et al., 2010; Ozbay et al., 2018). These observations support the existence of endogenous systemic LFOs, which propagate throughout the entire cardiovascular system, appearing as synchronized, but out of phase, oscillations at different vascular sites. A common systemic source is one explanation for the MAP correlated fMRI signals we have measured, and a potential candidate for this systemic origin is SNA. The MODWT reveals that MAP coupled fMRI LFOs are strongest in the frequency band centered at ∼0.1 Hz, the frequency of Mayer waves, which are defined in terms of their coherence with SNA (Julien, 2006). SNA regulates blood pressure via modulation of peripheral vascular tone (Fisher and Paton, 2012), but also potentially influences cerebrovascular tone (Brassard et al., 2017), and increases in SNA elicited by post exercise induced ischemia have been shown to decrease compliance in the brain’s major arteries (Warnert et al., 2016). Orthostatic challenges such as lower body negative pressure (LBNP), which are associated with increases in SNA, lead to considerable reductions in MCA CBFV (Levine et al., 1994; Serrador et al., 2000; Zhang and Levine, 2007), and reductions in blood volume indicative of vasoconstriction in the brain’s largest arteries (Whittaker et al., 2017). Furthermore, studies have shown that ganglion blockade designed to dampen SNA, significantly alter the dynamics between MAP and CBFV (Zhang et al., 2002; Mitsis et al., 2009), suggesting autonomic neural control cerebrovascular tone likely plays a role in beat-to-beat CBF regulation.
The time scale of the BP correlated LFO and its basis predominantly being changes in apparent transverse relaxation is strongly indicative of a CBF related cause. Compared with respiratory challenges for which the CBF response primarily probes CVR, the observed flow response associated with BP is likely to be related to the process of CA. Thus, whereas CVR is a measure of localized vascular integrity, i.e., the ability of arterial vessels to change their resistance, measures of CA relate to the systemic orchestrated vascular mechanisms that regulate CBF (Carrera et al., 2009). The LFO fluctuations we have observed in this study are correlated with BP measured in the periphery, and so are more related to CA than CVR. Impairments in CA associated with adverse cerebrovascular events such as ischaemic stroke and severe head injury have been well studied (Panerai, 2008), and are increasingly thought to play a role in the development of vascular dementia (Toth et al., 2017) and Alzheimer’s disease (Claassen and Zhang, 2011; den Abeelen et al., 2014). Despite the widespread clinical implications of pathological CA, its underlying mechanisms are still relatively poorly understood.
Transcranial Doppler ultrasound is the most widely used modality for measuring CA, which despite having excellent temporal resolution and high suitability for clinical settings, is ultimately of limited value since the measurements are restricted to only the largest intracranial arteries. In contrast, fMRI has whole-brain sensitivity with millimeter resolution and so is a desirable tool for better understanding CA, and has the potential to deliver more predictive clinical measures. For example, CA is critical for keeping stable CBF in the penumbra region following ischaemic stroke (Xiong et al., 2017), so a method such as fMRI, which has the spatial resolution to resolve localized alterations, is promising as a more informative prognostic tool. In the TCD literature the transfer function between BP and CBFV is used to characterize CA, primarily through gain and phase shift. It is commonly assumed that a phase shift and low gain constitutes good CA (i.e., CBFV fluctuations are delayed with respect to BP and are dampened) (van Beek et al., 2008). In this study we observed a lag in fMRI with respect to MAP, which may be related to the phase shifts measured in TCD. Furthermore, although the effect-size of MAP on fMRI measured here appears small (Supplementary Figure 4), this may be due to the young healthy subject group. In patient groups with less effective CA both effect-size and lag may be modulated.
Effect on Resting-State fMRI
This study provides the first step in characterizing the relationship between MAP and the fMRI signal, but further work is needed to address the degree to which MAP impacts functional connectivity measures. The estimated effect size on unfiltered data is relatively modest (∼0.01% BOLD/mm Hg), which across subjects on average amounts to total BOLD signal changes on the order of ∼0.1% across all gray matter, although total signal changes as large ∼0.5% are possible, depending on individual subject response variability. The results of the wavelet transformed data show that MAP fluctuations effect fMRI within a particular frequency band, and as in practice resting-state fMRI analyses never use raw unfiltered data, it is likely that the effect of MAP on functional connectivity metrics will depend on a variety of analysis choices, such as filter passband or window length in dynamic connectivity studies. These data also serve as a reminder that not all sources of BOLD contrast are neuronal in origin, and so multi-echo based approaches like ME-ICA (Kundu et al., 2012) are likely to be less effective. An interesting avenue of future research would be look at the effect of different echo combination and de-noising schemes to determine their impact.
In this study we have shown that beat-to-beat fluctuations in BP are correlated with fluctuations in the resting-state fMRI signal that are delayed by approximately 5.5 s, and which are strongest at the frequency band centered at ∼0.1 Hz. Using a multi-echo acquisition we were able to isolate the pure BOLD (R2∗) component of the BP correlated fMRI signal and have shown that it is the main source of contrast. This would indicate that it is changes in CBF that mediate this low frequency BP correlated signal, which we hypothesize is related to the process of CA. We propose that resting-state fMRI is a promising new tool for assessment of dynamic CA with high spatial resolution, which may prove to be a useful biomarker in a range of cerebrovascular and neurological conditions.
The School of Psychology Cardiff University Ethics Committee approved the study in accordance with the guidelines stated in the Cardiff University Research Framework (version 4.0, 2010).
JW and KM conceived of the presented idea. JW, ID, MV, and MB collected the data. JW analyzed the data. All the authors provided critical feedback and helped shape the research, analysis, and manuscript.
This work was supported by the Wellcome Trust [WT200804].
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.
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2019.00433/full#supplementary-material
Baruch, M. C., Kalantari, K., Gerdt, D. W., and Adkins, C. M. (2014). Validation of the pulse decomposition analysis algorithm using central arterial blood pressure. Biomed. Eng. Online 13:96. doi: 10.1186/1475-925X-13-96
Baruch, M. C., Warburton, D. E., Bredin, S. S., Cote, A., Gerdt, D. W., and Adkins, C. M. (2011). Pulse decomposition analysis of the digital arterial pulse during hemorrhage simulation. Nonlinear Biomed. Phys. 5:1. doi: 10.1186/1753-4631-5-1
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
Blaber, A. P., Bondar, R. L., Stein, F., Dunphy, P. T., Moradshahi, P., Kassam, M. S., et al. (1997). Transfer function analysis of cerebral autoregulation dynamics in autonomic failure patients. Stroke 28, 1686–1692. doi: 10.1161/01.str.28.9.1686
Blockley, N. P., Driver, I. D., Francis, S. T., Fisher, J. A., and Gowland, P. A. (2011). An improved method for acquiring cerebrovascular reactivity maps. Magn. Reson. Med. 65, 1278–1286. doi: 10.1002/mrm.22719
Brassard, P., Tymko, M. M., and Ainslie, P. N. (2017). Sympathetic control of the brain circulation: appreciating the complexities to better understand the controversy. Auton. Neurosci. 207, 37–47. doi: 10.1016/j.autneu.2017.05.003
Bright, M. G., Bianciardi, M., de Zwart, J. A., Murphy, K., and Duyn, J. H. (2014). Early anti-correlated BOLD signal changes of physiologic origin. Neuroimage 87, 287–296. doi: 10.1016/j.neuroimage.2013.10.055
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
Carrera, E., Lee, L. K., Giannopoulos, S., and Marshall, R. S. (2009). Cerebrovascular reactivity and cerebral autoregulation in normal subjects. J. Neurol. Sci. 285, 191–194. doi: 10.1016/j.jns.2009.06.041
den Abeelen, A. S., Lagro, J., van Beek, A. H., and Claassen, J. A. (2014). Impaired cerebral autoregulation and vasomotor reactivity in sporadic Alzheimer’s disease. Curr. Alzheimer Res. 11, 11–17. doi: 10.2174/1567205010666131119234845
Diehl, R. R., Diehl, B., Sitzer, M., and Hennerici, M. (1991). Spontaneous oscillations in cerebral blood flow velocity in normal humans and in patients with carotid artery disease. Neurosci. Lett. 127, 5–8. doi: 10.1016/0304-3940(91)90880-3
Donahue, M. J., Strother, M. K., Lindsey, K. P., Hocke, L. M., Tong, Y., and Frederick, B. D. (2016). Time delay processing of hypercapnic fMRI allows quantitative parameterization of cerebrovascular reactivity and blood flow delays. J. Cereb. Blood Flow Metab. 36, 1767–1779. doi: 10.1177/0271678x15608643
Duffin, J., Sobczyk, O., Crawley, A. P., Poublanc, J., Mikulis, D. J., and Fisher, J. A. (2015). The dynamics of cerebrovascular reactivity shown with transfer function analysis. Neuroimage 114, 207–216. doi: 10.1016/j.neuroimage.2015.04.029
Erdogan, S. B., Tong, Y., Hocke, L. M., Lindsey, K. P., and de, B. F. B. (2016). Correcting for blood arrival time in global mean regression enhances functional connectivity analysis of resting state fMRI-BOLD Signals. Front. Hum. Neurosci. 10:311. doi: 10.3389/fnhum.2016.00311
Glover, G. H., Li, T. Q., and Ress, D. (2000). Image-based method for retrospective correction of physiological motion effects in fMRI: RETROICOR. Magn. Reson. Med. 44, 162–167. doi: 10.1002/1522-2594(200007)44:1<162::aid-mrm23>3.3.co;2-5
Golestani, A. M., Wei, L. L., and Chen, J. J. (2016). Quantitative mapping of cerebrovascular reactivity using resting-state BOLD fMRI: validation in healthy adults. Neuroimage 138, 147–163. doi: 10.1016/j.neuroimage.2016.05.025
Kundu, P., Inati, S. J., Evans, J. W., Luh, W. M., and Bandettini, P. A. (2012). Differentiating BOLD and non-BOLD signals in fMRI time series using multi-echo EPI. Neuroimage 60, 1759–1770. doi: 10.1016/j.neuroimage.2011.12.028
Kundu, P., Voon, V., Balchandani, P., Lombardo, M. V., Poser, B. A., and Bandettini, P. A. (2017). Multi-echo fMRI: a review of applications in fMRI denoising and analysis of BOLD signals. Neuroimage 154, 59–80. doi: 10.1016/j.neuroimage.2017.03.033
Kuo, T. B., Chern, C. M., Sheng, W. Y., Wong, W. J., and Hu, H. H. (1998). Frequency domain analysis of cerebral blood flow velocity and its correlation with arterial blood pressure. J. Cereb. Blood Flow Metab. 18, 311–318. doi: 10.1097/00004647-199803000-00010
Levine, B. D., Giller, C. A., Lane, L. D., Buckey, J. C., and Blomqvist, C. G. (1994). Cerebral versus systemic hemodynamics during graded orthostatic stress in humans. Circulation 90, 298–306. doi: 10.1161/01.cir.90.1.298
Mitsis, G. D., Poulin, M. J., Robbins, P. A., and Marmarelis, V. Z. (2004). Nonlinear modeling of the dynamic effects of arterial pressure and CO2 variations on cerebral blood flow in healthy humans. IEEE Trans. Biomed. Eng. 51,1932–1943. doi: 10.1109/tbme.2004.834272
Mitsis, G. D., Zhang, R., Levine, B. D., Tzanalaridou, E., Katritsis, D. G., and Marmarelis, V. Z. (2009). Autonomic neural control of cerebral hemodynamics. IEEE Eng. Med. Biol. Mag. 28, 54–62. doi: 10.1109/MEMB.2009.934908
Moeller, S., Yacoub, E., Olman, C. A., Auerbach, E., Strupp, J., Harel, N., et al. (2010). Multiband multislice GE-EPI at 7 tesla, with 16-fold acceleration using partial parallel imaging with application to high spatial and temporal whole-brain fMRI. Magn. Reson. Med. 63, 1144–1153. doi: 10.1002/mrm.22361
Murphy, K., Harris, A. D., and Wise, R. G. (2011). Robustly measuring vascular reactivity differences with breath-hold: normalising stimulus-evoked and resting state BOLD fMRI data. Neuroimage 54, 369–379. doi: 10.1016/j.neuroimage.2010.07.059
Obrig, H., Neufang, M., Wenzel, R., Kohl, M., Steinbrink, J., Einhaupl, K., et al. (2000). Spontaneous low frequency oscillations of cerebral hemodynamics and metabolism in human adults. Neuroimage 12, 623–639. doi: 10.1006/nimg.2000.0657
Ozbay, P. S., Chang, C., Picchioni, D., Mandelkow, H., Moehlman, T. M., Chappel-Farley, M. G., et al. (2018). Contribution of systemic vascular effects to fMRI activity in white matter. Neuroimage 176, 541–549. doi: 10.1016/j.neuroimage.2018.04.045
Panerai, R. B., Hudson, V., Fan, L., Mahony, P., Yeoman, P. M., Hope, T., et al. (2002). Assessment of dynamic cerebral autoregulation based on spontaneous fluctuations in arterial blood pressure and intracranial pressure. Physiol. Meas. 23, 59–72. doi: 10.1088/0967-3334/23/1/306
Patel, A. X., and Bullmore, E. T. (2016). A wavelet-based estimator of the degrees of freedom in denoised fMRI time series for probabilistic testing of functional connectivity and brain graphs. Neuroimage 142, 14–26. doi: 10.1016/j.neuroimage.2015.04.052
Qiao, M., Rushforth, D., Wang, R., Shaw, R. A., Tomanek, B., Dunn, J. F., et al. (2007). Blood-oxygen-level-dependent magnetic resonance signal and cerebral oxygenation responses to brain activation are enhanced by concurrent transient hypertension in rats. J. Cereb. Blood Flow Metab. 27, 1280–1289. doi: 10.1038/sj.jcbfm.9600436
Rayshubskiy, A., Wojtasiewicz, T. J., Mikell, C. B., Bouchard, M. B., Timerman, D., Youngerman, B. E., et al. (2014). Direct, intraoperative observation of ∼0.1 Hz hemodynamic oscillations in awake human cortex: implications for fMRI. Neuroimage 87, 323–331. doi: 10.1016/j.neuroimage.2013.10.044
Serrador, J. M., Picot, P. A., Rutt, B. K., Shoemaker, J. K., and Bondar, R. L. (2000). MRI measures of middle cerebral artery diameter in conscious humans during simulated orthostasis. Stroke 31, 1672–1678. doi: 10.1161/01.str.31.7.1672
Tong, Y., and Frederick, B. D. (2010). Time lag dependent multimodal processing of concurrent fMRI and near-infrared spectroscopy (NIRS) data suggests a global circulatory origin for low-frequency oscillation signals in human brain. Neuroimage 53, 553–564. doi: 10.1016/j.neuroimage.2010.06.049
Tong, Y., Hocke, L. M., Licata, S. C., and Frederick, B. (2012). Low-frequency oscillations measured in the periphery with near-infrared spectroscopy are strongly correlated with blood oxygen level-dependent functional magnetic resonance imaging signals. J. Biomed. Opt. 17:106004. doi: 10.1117/1.JBO.17.10.106004
Tong, Y., Hocke, L. M., Nickerson, L. D., Licata, S. C., Lindsey, K. P., and Frederick, B. (2013). Evaluating the effects of systemic low frequency oscillations measured in the periphery on the independent component analysis results of resting state networks. Neuroimage 76, 202–215. doi: 10.1016/j.neuroimage.2013.03.019
Toth, P., Tarantini, S., Csiszar, A., and Ungvari, Z. (2017). Functional vascular contributions to cognitive impairment and dementia: mechanisms and consequences of cerebral autoregulatory dysfunction, endothelial impairment, and neurovascular uncoupling in aging. Am. J. Physiol. Heart Circ. Physiol. 312, H1–H20. doi: 10.1152/ajpheart.00581.2016
Uchida, S., Bois, S., Guillemot, J. P., Leblond, H., and Piche, M. (2017). Systemic blood pressure alters cortical blood flow and neurovascular coupling during nociceptive processing in the primary somatosensory cortex of the rat. Neuroscience 343, 250–259. doi: 10.1016/j.neuroscience.2016.12.014
van Beek, A. H., Claassen, J. A., Rikkert, M. G., and Jansen, R. W. (2008). Cerebral autoregulation: an overview of current concepts and methodology with special focus on the elderly. J. Cereb. Blood Flow Metab. 28, 1071–1085. doi: 10.1038/jcbfm.2008.13
van Gelderen, P., de Zwart, J. A., and Duyn, J. H. (2008). Pittfalls of MRI measurement of white matter perfusion based on arterial spin labeling. Magn. Reson. Med. 59, 788–795. doi: 10.1002/mrm.21515
van Houdt, P. J., Ossenblok, P. P., Boon, P. A., Leijten, F. S., Velis, D. N., Stam, C. J., et al. (2010). Correction for pulse height variability reduces physiological noise in functional MRI when studying spontaneous brain activity. Hum. Brain Mapp. 31, 311–325. doi: 10.1002/hbm.20866
Wang, R., Foniok, T., Wamsteeker, J. I., Qiao, M., Tomanek, B., Vivanco, R. A., et al. (2006). Transient blood pressure changes affect the functional magnetic resonance imaging detection of cerebral activation. Neuroimage 31, 1–11. doi: 10.1016/j.neuroimage.2005.12.004
Warnert, E. A., Hart, E. C., Hall, J. E., Murphy, K., and Wise, R. G. (2016). The major cerebral arteries proximal to the circle of willis contribute to cerebrovascular resistance in humans. J. Cereb. Blood Flow Metab. 36,1384–1395. doi: 10.1177/0271678X15617952
Whittaker, J. R., Bright, M. G., Driver, I. D., Babic, A., Khot, S., and Murphy, K. (2017). Changes in arterial cerebral blood volume during lower body negative pressure measured with MRI. Neuroimage 187, 166–175. doi: 10.1016/j.neuroimage.2017.06.041
Willie, C. K., Colino, F. L., Bailey, D. M., Tzeng, Y. C., Binsted, G., Jones, L. W., et al. (2011). Utility of transcranial doppler ultrasound for the integrative assessment of cerebrovascular function. J. Neurosci. Methods 196, 221–237. doi: 10.1016/j.jneumeth.2011.01.011
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
Xiong, L., Liu, X., Shang, T., Smielewski, P., Donnelly, J., Guo, Z. N., et al. (2017). Impaired cerebral autoregulation: measurement and application to stroke. J. Neurol. Neurosurg. Psychiatry 88, 520–531.
Yucel, M. A., Selb, J., Aasted, C. M., Lin, P. Y., Borsook, D., Becerra, L., et al. (2016). Mayer waves reduce the accuracy of estimated hemodynamic response functions in functional near-infrared spectroscopy. Biomed. Opt. Express 7, 3078–3088. doi: 10.1364/BOE.7.003078
Zhang, R., and Levine, B. D. (2007). Autonomic ganglionic blockade does not prevent reduction in cerebral blood flow velocity during orthostasis in humans. Stroke 38, 1238–1244. doi: 10.1161/01.str.0000260095.94175.d0
Keywords: cerebral autoregulation, resting-state fMRI, blood pressure, cerebral physiology, LFO, BOLD, CBF
Citation: Whittaker JR, Driver ID, Venzi M, Bright MG and Murphy K (2019) Cerebral Autoregulation Evidenced by Synchronized Low Frequency Oscillations in Blood Pressure and Resting-State fMRI. Front. Neurosci. 13:433. doi: 10.3389/fnins.2019.00433
Received: 31 October 2018; Accepted: 15 April 2019;
Published: 07 May 2019.
Edited by:Jean Chen, University of Toronto, Canada
Reviewed by:Changwei Wesley Wu, Taipei Medical University, Taiwan
Thomas T. Liu, University of California, San Diego, United States
Copyright © 2019 Whittaker, Driver, Venzi, Bright and Murphy. 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: Joseph R. Whittaker, email@example.com