Low Frequency Systemic Hemodynamic “Noise” in Resting State BOLD fMRI: Characteristics, Causes, Implications, Mitigation Strategies, and Applications

Advances in functional magnetic resonance imaging (fMRI) acquisition have improved signal to noise to the point where the physiology of the subject is the dominant noise source in resting state fMRI data (rsfMRI). Among these systemic, non-neuronal physiological signals, respiration and to some degree cardiac fluctuations can be removed through modeling, or in the case of newer, faster acquisitions such as simultaneous multislice acquisition, simple spectral filtering. However, significant low frequency physiological oscillation (∼0.01–0.15 Hz) remains in the signal. This is problematic, as it is the precise frequency band occupied by the neuronally modulated hemodynamic responses used to study brain connectivity, precluding its removal by spectral filtering. The source of this signal, and its method of production and propagation in the body, have not been conclusively determined. Here, we summarize the defining characteristics of the systemic low frequency noise signal, and review some current theories about the signal source and the evidence supporting them. The strength and distribution of the systemic LFO signal make characterizing and removing it essential for accurate quantification, especially for resting state connectivity, when no stimulation can be compared with the signal. Widespread correlated non-neuronal signals obscure and distort the more localized patterns of neuronal correlations between interacting brain regions; they may even cause apparent connectivity between regions with no neuronal interaction. Here, we discuss a simple method we have developed to parse the global, moving, blood-borne signal from the stationary, neuronal connectivity signals, substantially reducing the negative correlations that result from global signal regression. Finally, we will discuss some of the uses to which the moving systemic low frequency oscillation can be put if we consider it a “signal” carrying information, rather than simply “noise” complicating the interpretation of resting state connectivity. Properly utilizing this signal may offer insights into subtle hemodynamic alterations that can be used as early indicators of circulatory dysfunction in a number of neuropsychiatric conditions, such as prodromal stroke, moyamoya, and Alzheimer’s disease.


INTRODUCTION
Resting state functional magnetic resonance imaging (rsfMRI) seeks to elucidate neuronal connectivity throughout the brain by examining of fMRI signal correlations between regions during scans when the brain is not performing any specific task (it is "at rest"). However, the BOLD signal does not measure neuronal activation directly. It is a composite measurement of hemodynamic properties of blood flow, volume, and oxygenation changes in response to neuronal activity (Buxton et al., 1998) (i.e., neurovascular coupling). In short, increased neuronal firing leads to an increase in regional blood flow, which brings an oversupply of oxygenated blood (Fox and Raichle, 1986). The increase in blood flow and oxygenation leads to an elevation of the BOLD signal. As a result, the observable BOLD signal (blood-related) is much slower (∼s) than neuronal firing (∼ms) (Logothetis et al., 2001) and cannot reflect fast changes of the field potential of neuronal firings. In mathematical terms, the BOLD signal is the result of the convolution of the fast neuronal signals with the slow hemodynamic functions. Therefore, the frequencies of "neuronal" BOLD signals are generally below 0.15 Hz (Josephs and Henson, 1999).
However, neuronal activations are not the only contributors to BOLD signals in the low frequency band. Advances in fMRI acquisition techniques and hardware have improved signal to noise to the point where the physiology of the subject being studied is the dominant noise source in rsfMRI data. In addition to neuronal BOLD, there are systemic, non-neuronal fluctuations in brain hemodynamics due to heartbeat, respiration, and so-called "low frequency oscillations" (LFOs). These signals are unavoidable, and taken together can account for 20-70% of the BOLD signal variance (see Figure 1), depending on acquisition, and locations of the voxels (Liu, 2017). Numerous processing strategies have been devised to mitigate them. Recent improvements in hardware and pulse sequences (specifically simultaneous multislice acquisition protocols) have pushed the temporal resolution of fMRI high enough that respiration, and to some degree cardiac fluctuations (see Figure 1) can be removed through simple spectral filtering, or through more advanced modeling methods when these signals are aliased (Glover et al., 2000;Birn et al., 2006Birn et al., , 2008Behzadi et al., 2007;Chang et al., 2009Chang et al., , 2013. These methods are well described elsewhere and are not the focus of this manuscript. However, significant signal power remains in the "low frequency oscillation band, " a loosely defined region from ∼0.01-0.15 Hz. Nonneuronal signal in this frequency band accounts for at least 30% of the signal variance in gray matter (Frederick et al., 2012a). This is problematic as it is the precise frequency band occupied by the neuronally modulated hemodynamic responses used to study brain connectivity, precluding its removal by filtering.
The strength and distribution of the systemic LFO signal make characterizing and removing it essential for accurate quantification of neuronal connectivity. Moreover, these signals might not be "noise" after all. Understanding their origins and characteristics will help to develop novel methods to assess brain physiology which could greatly compliment the functional findings.

CHARACTERISTICS OF LOW FREQUENCY OSCILLATIONS
LFOs in BOLD fMRI have been found and studied extensively, but as noted above, there are numerous potential explanations for this signal, and even some variation in what in particular should be considered a "low frequency oscillation." For the remainder of this paper, we will apply two criteria to our discussions of LFOs; the frequency band of the signal, and whether the signal is stationary, or moves with the blood.

Frequency Content
The first criterion is simply a definition. LFOs are signals that occur in the brain (and in some cases throughout the body) that have frequencies between ∼0.009 and 0.2 Hz. The exact endpoints of this band are extremely variable in the literature. Biswal's original paper on resting state connectivity used 0.01-0.1 Hz (Biswal et al., 1995), but later papers expanded this range; for the purpose of this discussion, we will use the range of 0.01-0.15. 0.15 Hz has a particular significance as the top of the range, as this is the highest expected frequency in neuronally generated hemodynamic signals (based on the shape of the canonical hemodynamic response function) (Josephs and Henson, 1999), so this defines the frequency range where spectral filtering cannot be used to remove non-neuronal signal. The frequency content of a typical BOLD signal is shown in Figure 1.

Dynamic Versus Stationary Noise Signals
The second defining characteristic of LFOs is less commonly considered. Our research into low frequency physiological noise in fMRI has established that a significant fraction of the low frequency variance in fMRI data can be modeled quite effectively as a single low frequency signal with varying delay times across the brain. Moreover, the pattern of relative delay times in different regions of the brain is consistent with the delays that would be expected if the signal were moving through the brain with blood as it flowed through the vasculature. We refer to this dynamic signal as "systemic low frequency oscillations" (sLFOs) (Tong et al., 2015). The realization that a significant fraction of the low frequency "noise" in fMRI appears to be moving has important implications for how to identify, remove, or even utilize this signal , 2014bTong et al., 2011bTong et al., ,c, 2013Tong et al., , 2014Tong et al., , 2017Tong et al., , 2018Frederick et al., 2012b).

Temporal Pattern -sLFOs Propagate on a Hemodynamic Timescale
Our research on sLFOs strongly suggest that the underlying oscillations propagate on hemodynamic, physiological, rather than on neuronal, timescales, taking several seconds to fully transit the brain rather than milliseconds. Circulatory measurements of the traversal of a Tc99 tracer through the brain vasculature (from the carotid, to the internal brain arteries, through the parenchyma, to the superior sagittal sinus) showed that the transit time of blood from the anterior and middle cerebral arteries to the superior sagittal sinus takes ∼6.7 s in healthy middle-aged controls (Crandell et al., 1973), timing consistent with results found by other imaging methods. For example, using echo contrast-enhanced ultrasound, delays of 7.5 ± 1.8 s from the carotids to the jugulars were found in 64 healthy subjects (Schreiber et al., 2002). Similar, but smaller delays (4.9-6.4 s) were found in other ultrasound studies (Schreiber et al., 2005). Direct evidence can also be found in digital subtraction angiography (DSA), where an x-ray contrast bolus was injected directly into the ICA and followed as it passed through the brain (Monti et al., 2015;Jann et al., 2016).

Spatial Pattern -sLFOs Travel Along the Vasculature
The patterns of the delays of the sLFO signal clearly suggest that this signal is related to blood flow, evolving in a pattern that reflects the vasculature -the sLFO signal appears first in the center of the brain, propagates out through the parenchyma, and ends up in the superior sagittal sinus, with a range of delays of ∼6.5 s , a pattern we have seen consistently in many subsequent studies over several years (Tong et al., 2011a(Tong et al., ,b,c, 2013(Tong et al., , 2014(Tong et al., , 2015Frederick et al., 2012a;Frederick, 2012, 2014a). In our most recent studies we have directly confirmed the association between the sLFO delay pattern and blood flow, first by performing sLFO analysis and time resolved dynamic susceptibility imaging in the same scan session , and more recently by following the sLFO signal all the way from the internal carotid arteries through the draining veins .

Origin -sLFOs Seem to Originate Outside the Brain
The sLFO BOLD signal identified from the carotids preceded the signal found in any voxel of the brain. Indeed, the same sLFO signal can be found throughout the body (measured in the periphery using NIRS) Tong and Frederick, 2010;Li et al., 2018). The delays in the periphery are symmetric across the midline of the body, and the arrival time of the sLFO signal found in the fingers and toes precede the arrival in many brain voxels . While it is possible that some process in the brain is the ultimate source of this moving signal, there is no evidence whatsoever in the fMRI data that this is the case -to the contrary the implication is that the signal does not originate in the brain.

Summary
From the growing body of evidence from our group and others we can summarize that the sLFO BOLD signal: (1) is a spontaneous physiological oscillation, (2) travels with the blood, and (3) has an extracerebral origin. Given the large amount of LFO signal variance that is clearly attributable to the moving component (at least 30% of the low frequency signal variance in gray matter, Frederick et al., 2012a), we believe this constitutes the majority of the physiological LFO signal power. In addition, these qualities provide the key to isolating the signal from the neuronal signals of interest.
In contrast, there may also be non-neuronal LFOs which do not propagate. However, because of the difficulty in separating these signals from putative neuronal signals, they are far harder to characterize. Certain mechanisms (detailed below) such as the Mayer wave, are thought to be synchronous throughout the body, and therefore stationary. Isolating the contribution of stationary LFOs to the resting state signal would require as yet undeveloped processing strategies, which is why we will discuss these possible sources, however we will focus primarily on the dynamic portion of the signal during the remainder of the paper.

POTENTIAL CAUSES OF THE LOW FREQUENCY OSCILLATION
The source of LFO signal, and its mechanism of production and propagation in the body, have not been conclusively determined -LFOs have been variously attributed to alterations in sympathetic nervous system tone, partial pressure of carbon dioxide (paCO 2 ) fluctuations modulated by respiration, blood pressure regulation, low frequency neuronal "waves, " and even gastric motility. It may in fact be a combination of multiple, independent signals with distinct sources. We will review these current theories about the signal source and the evidence supporting them. This section is summarized in Table 1.

Variations in Heart Rate and Respiratory Volumes
One major potential source of LFOs comes from variations of the heart rate and respiration. For example respiration volume per time (RVT) (Birn et al., 2006) and respiration variation (convolved with the respiratory response function, Birn et al., 2008;Chang et al., 2009) are two methods to model fluctuations from CO 2 concentration. In specific, the former model describes the depth of the respiration, whereas the latter reflects the variation in respiration. Previously, these methods have been shown to explain additive variation in data with longer TRs between 9 and 11% (Birn et al., 2006(Birn et al., , 2008Chang et al., 2009) in the voxels affected. From the same group a model taking the variation in heart rate into account was developed (Chang et al., 2009), explaining 3% more variance in the affected voxels. The underlying mechanism of this change is still not well understood, but was proposed to relate to neuronal activity linked with changes in levels of arousal (Chang et al., 2009). These models typically incorporate a delay of several seconds to best match the modeled noise waveforms with the fMRI data, suggesting that these signals are in the moving category (although the standard implementations of these methods do not account for regional delays within the brain).

Carbon Dioxide
Carbon dioxide changes are closely related to the previous topic, and they partially share the same mechanism, namely induced changes in cerebral blood flow and volume due to CO 2 induced vasodilation. However, here we discuss the direct effects of paCO 2 by comparing BOLD fMRI to measured fluctuations in the partial pressure of end-tidal carbon dioxide, especially in the lower range of the LF band, namely 0-0.05 Hz (Wise et al., 2004;Sassaroli et al., 2012), rather than the indirect estimation of this effect from the respiratory and/or cardiac waveforms. Highly synchronous within the brain -it seems unlikely that it would appear to move.

Yes? Not likely
Vasomotion Vasomotion is a spontaneous oscillation (0.01∼0.3 Hz) in the vascular tone, which is independent of respiration, pulsation and neuronal activity.
The induced vascular variations could propagate with the blood, resulting in a mix of stationary and moving signal.
Unclear Possible CO 2 Carbon dioxide is a potent vassal dilator. It can travel with the blood and induce changes in cerebral blood flow and volume.
Clearly a moving signal -the CO 2 travels in the blood.

No Possible
Variations in heart rate and respiratory volumes The variations of the heart rate and respiration, including the depth of the respiration, are in the low frequency range.
Depth of respiration and heart rate variation can alter blood volume (through CO 2 and pressure changes), and the effects should move with the blood.
No Possible

Gastric oscillations
The electrogastrogram signal (i.e., synchronized gut motions at ∼0.05 Hz) significant correlates with BOLD fMRI data with time delays.
May be controlled neuronally, but the effects seem to originate in the gut and move with blood.

Mostly no Possible
Aliased signals of cardiac and respiration Aliased signals of cardiac and respiration due to long TR are in the low frequency range.
These signals are in the right frequency range and can travel with blood, but are not highly correlated with the sLFO signal when tested by fMRI data with very short TR or fNIRS.

No
Not likely  Wise et al. (2004) found that paCO 2 levels measured with end tidal CO 2 in the 0-0.05 Hz band were significantly correlated with both increased middle cerebral artery blood velocity and increased BOLD fMRI signal in gray and white matter. This is attributed to the vasodilatory effect of CO 2 -increased CO 2 leads to increased arterial diameter and blood volume. As a consequence, this signal can be clearly placed in the group of sLFOs that move with the blood. Wise found up to 28% of the low frequency signal variance in the BOLD signal was attributable to the paCO 2 . This is likely the lower limit, because while the peak correlation delays between the BOLD and paCO 2 timecourses were between 3 and 13 s, only a single delay of 6.3 s was used for all voxels, which would reduce the apparent correlation.

Mayer Waves
Mayer waves are spontaneous LFOs at around 0.1 Hz (Julien, 2006) and are associated with variations in mean arterial pressure, and have been associated with a sympathetic autoregulation mechanism (Tsuji et al., 2000), in particular below 0.1 Hz (Obrig et al., 2000;van Beek et al., 2008;Sassaroli et al., 2012). These studies connected mean arterial blood pressure to cerebral intravascular oxygenation when monitoring cerebral circulation and blood flow with NIRS. However, a characteristic of Mayer waves is their largely synchronous nature within the brain (Sassaroli et al., 2012), differing from the sLFOs, which have differing delays across the brain Tong et al., 2011b). We would therefore classify Mayer waves as stationary LFO signals.

Vasomotion From Oscillations in the Vascular Tone
Vasomotion is a spontaneous oscillation (0.01∼0.3 Hz) in the vascular tone, which is independent of respiration, pulsation and neuronal activity (Hundley et al., 1988;Mayhew et al., 1996;Rivadulla et al., 2011). These oscillations of the lumen diameter modify blood flow in a corresponding fashion resulting in periodic oscillations in the blood flow (Aalkjaer et al., 2011). Another LF signal is attributed to vasomotion, referring to the oscillations in the vascular tone, thought to be generated movement within the vessel walls (Zhang et al., 1998;Aalkjaer et al., 2011;Sassaroli et al., 2012;Müller and Österreich, 2014). These changes are highly localized and have been linked to oscillatory intracellular calcium (Aalkjaer et al., 2011). This would give the signals a local, stationary nature, however the vascular variations induced would then propagate with the blood, resulting in a mix of stationary and moving signal.

Aliased Signals of Cardiac and Respiration
Because fMRI is usually not sampled fast enough to resolve cardiac or respiratory waveforms, some fraction of the energy in these signals will be aliased to the low frequency band. In order to determine the significance of this component of the signal, we evaluated sLFOs in a dataset with high temporal resolution in which the respiration and cardiac waveforms are fully sampled . We found that even with fully sampled data, in which the respiration and cardiac bands can be isolated with spectral filters, the purely non-neuronal sLFOs (as determined by time-delayed correlation with a peripheral NIRS signal) still account for over 13% of the total BOLD signal variance across all frequency bands.

Gastric Oscillations
A somewhat more recent theory for the cause of sLFOs is proposed by Rebollo et al. (2018) who found significant correlations between electrogastrogram signals (which measure synchronized gut motions at ∼0.05 Hz) and BOLD fMRI data. Unfortunately, the technique used was unable to determine directionality. However, the authors found delays between the earliest (somatosensory cortices) and the latest (dorsal precuneous and extrastriate body area) nodes of the proposed "gastric network" were about 3.3 s. These later regions lie in close proximity to the superior sagittal and transverse sinuses, respectively, vessels at the end of the vascular path through the brain. While Rebollo's analysis could not ascertain directionality -this observation is consistent with a hemodynamic perturbation generated in the stomach which then propagated through the cerebral vasculature, which suggests that gastric signals likely contribute to the "moving signal" category.
This finding is in good agreement with previous work by Yacin et al. (2011) which showed a strong relationship between gastric activity and systemic LFOs in the periphery. Yacin was able to reconstruct the gastric slow wave signal from a fingertip photoplethysmogram, using a deep learning approach. The reconstructed signal correlated with the measured electrogastrogram slow wave with R >= 0.9, clearly establishing that the gastric signal contributed a significant portion of the sLFO variance observable in the periphery.

IMPLICATIONS FOR RESTING STATE ANALYSIS
As mentioned previously, BOLD fMRI infers neuronal activation indirectly through neurovascular coupling. As a result, the neuronal activation will appear in the low frequency range (∼0.01-0.15 Hz) of BOLD signal in both resting state and task fMRI studies. This frequency range significantly overlaps with that of sLFO. Therefore, the presence of sLFO in the BOLD signal will confound the results of fMRI analyses, especially in resting state studies, as the neuronal firing is also spontaneous and of unknown timecourse (like sLFO), unlike task activation which can be modeled.

Pure Physiological sLFOs in Resting State Networks
In the following, we describe a previously published study demonstrating the confounding effect of sLFO to the resting state network analysis in both simulation as well as using concurrent fMRI/fNIRS data . The signals measured by NIRS (concentration changes in oxyhemoglobin and deoxyhemoglobin: [HbO] and [Hb]) are, like BOLD fMRI signal, blood-related measures. High consistency between these NIRS and fMRI signals has been demonstrated in concurrent studies (Strangman et al., 2002;Sassaroli et al., 2006;Cui et al., 2011). In this particular concurrent study, [HbO] and [Hb] were measured in the periphery (i.e., on finger and toes) using NIRS, instead of the brain. We found that the LFO band component of [HbO] in the periphery was highly correlated with the sLFOs of the BOLD signal in the brain, with a time delay. Most importantly, since the data was recorded in the periphery, the sLFOs measured here represented "pure" physiological fluctuations that were clearly not contaminated by the neuronal LFO as they would be in the brain.
We used this pure signal to assess the physiological influence of sLFOs on each resting state network. In detail, each subject's connectivity networks were identified from rsfMRI data using ICA (Melodic and FSL), from which, the signature BOLD signals of each network were extracted. To assess the sLFO signal's influence, we correlated each network's signature BOLD signal with the subject's own concurrent peripheral data (i.e., [HbO] of the fingertip). The networks with high correlations are the ones being significantly influenced by the sLFOs. The results showed that in addition to clearly vascular components, such as the superior sagittal sinus, several sensory networks (i.e., visual, auditory, etc.) are strongly affected by the sLFOs (Figure 2).

The Dynamic Nature of sLFOs in Resting State Networks
With this in mind, we sought to determine if the dynamic nature of the sLFO signal in BOLD could, on its own, lead to artifactual connectivity in "Can apparent network patterns be generated by the moving physiological sLFO signal alone?" (Tong et al., 2015). First, we performed a simple test on simulated data, which consisted of a sinusoidal wave with gradually increased time delay along the traveling direction (Figure 3a, the direction of the arrows), representing the traveling sLFO BOLD signal in the brain, additive noise, and a constant offset (see Figure 3). We then applied a standard resting state analysis methodology, namely ICA. The result showed that multiple "networks" along the traveling direction were identified, even though the only difference between the time series of these networks is the time delay (see Figure 3c). This simulation demonstrated that methods like ICA are prone to being confounded by time delayed versions of identical signals in different voxels. We then performed this test on data much closer to real rsfMRI data. Initially, we calculated the delay between every brain voxel and the peripheral NIRS signal using the subjects' real sLFO BOLD data as described above . Then, a sinusoidal signal, adjusted with these delays (the real delays of the subjects resting state data) replaced the real BOLD signal at each voxel (i.e., the voxel-specific time delay was decided by the delay value of that voxel). After that, we applied ICA on these simulated fMRI data with identical time series at each voxel differing only by the time delay. As a result, several "resting state networks" (RSNs) were identified (see Figure 4), some of which closely matched standard networks described in the literature. Finally, we applied seed analysis on the same simulated data and were able to identify "RSNs" as well (it is known that seed analysis is sensitive to time delays). This study demonstrated that physiological noise signals, depending only on vascular time delays can generate network patterns similar to well-known RSNs through common analytical procedures.
While clearly worrisome, these results should not be interpreted as suggesting that RSNs are nothing but a vascular artifact. There is extensive evidence for the existence of neuronal RSNs, from both animal and human studies, with a range of imaging technologies (Martinez-Montes et al., 2004;Hillman et al., 2007;Brookes et al., 2011;Ma et al., 2016). Moreover, while we have shown that up to 30% of the low frequency gray matter variance (13% of the total variance across all frequency bands, Hocke et al., 2016) is due to non-neuronal sLFOs (Frederick et al., 2012a), this means, necessarily, that 70% of the variance is not due to sLFOs, and likely represents neuronal signal. However, it is clear that there are both vascular and neuronal "connectivity" networks, with significant spatial overlap.
In order to support the metabolic demands of neurons, vascular networks are formed according to various factors, such as neuron density and metabolic demand. Areas of the brain which routinely coactivate likely develop similar vascular supplies, at similar times. This could lead to overlapping of both networks. It is also possible that, despite being a map of sLFO arrival time, the delay map may nevertheless contain some neuronal information. We hypothesize that elevated neuronal activation within a network will increase the blood flow locally, minimizing the time delays within the network. These subtle differences could be identified by ICA. Nevertheless, we demonstrated that sLFO signals will confound the quantification of some RSNs, in both the spatial and temporal domains, unless some steps are taken to disentangle these signals.

MITIGATION STRATEGIES
As discussed, these widespread, correlated, non-neuronal sLFOs obscure and distort the more localized patterns of neuronal correlations between interacting brain regions, and may even cause apparent connectivity between regions with no neuronal interaction -numerous vascular "networks" are commonly seen in data driven connectivity analyses. Standard mitigation methods such as global signal regression (GSR) have serious drawbacks, and may in fact induce artifactual negative

B A
FIGURE 4 | Results from group ICA on 11 subjects' real BOLD data were shown in panel (A). Results from group ICA on 11 subjects' synthetic data were shown in panel (B). The value in each result showed the spatial correlation coefficient calculated between that component and the corresponding RSN from the template (Beckmann et al., 2005). The two components in the red block are the same (Figure adapted from Tong et al., 2015).
correlations between brain regions. Newer methods such as CompCor (Behzadi et al., 2007) avoid some of these problems, but may not fully remove the sLFO signal. By examining the structure of temporal cross-correlations with non-zero time delays throughout the brain (or by making simultaneous, independent measurements in the periphery), it is possible to parse the global, moving, blood-borne signal from the stationary, neuronal connectivity signals. In the following we describe a simple method to remove this signal, leaving the neuronal connectivity intact, while substantially reducing the negative (potentially spurious) correlations that result from global signal regression.
We have conducted continuous research on isolating, characterizing, and separating the neuronal LFO and sLFO in the resting state BOLD fMRI during the past 8 years , 2014aTong et al., 2011aTong et al., ,b,c, 2012Tong et al., , 2013Tong et al., , 2014Tong et al., , 2015Tong et al., , 2016Tong et al., , 2017Tong et al., , 2018Frederick et al., 2012a;Erdogan et al., 2016;Hocke et al., 2016). The methods developed utilized the key differences between these two oscillations, part of which are discussed in "Characteristics of LFOs": (1) neuronal LFO is regional, while sLFO is global; (2) the neuronal LFO signal does not "propagate" in space, while sLFO does (sLFOs are dynamic), traversing the brain on the time scale of seconds; (3) while neuronal LFO should be found largely in the capillary bed (it is known to biased toward veins), sLFO BOLD is also detected near/in the large blood vessels, especially veins as well as capillaries; (4) neuronal LFO originates in the brain, while sLFO has extracerebral origins and propagates into and through the brain with the blood. These differences do allow us to effectively parse the sLFO and neuronal components of the low frequency BOLD signal, and examine them separately using a technique we call Regressor Interpolation at Progressive Time Delays (RIPTiDe) which is described in more detail below.

Regressor Interpolation at Progressive Time Delays (RIPTiDe)
The RIPTiDe procedure has been described in detail previously (Frederick et al., 2012a). Initially, we developed a method to determine relative blood arrival times in the voxels of a resting state fMRI (rsfMRI) dataset by using simultaneous NIRS to estimate the non-neuronal, systemic signal. We use a crosscorrelation procedure to determine the delay time between this peripheral signal and the timecourse of each voxel in the fMRI dataset. Determination of the precise arrival time of blood-borne signal at every brain voxel also allows correct determination of the fraction of that voxel's signal that is accounted for by the moving blood signal, which depends on the voxel's relative cerebral blood volume (rCBV) and the oxygenation of blood in the voxel . We have used this NIRS based method to remove non-neuronal signal from BOLD data (Frederick et al., 2012a), as a probe to study physiological signal partitioning in brain Frederick, 2012, 2014a;Tong et al., , 2013, and to measure cerebrovascular reserve (Tong et al., 2011a).
Our previous time-delay image analysis work focused closely on two physiological inputs -endogenous hemodynamic fluctuations in normal controls measured directly and concurrently with scanning using NIRS, and exogenous hemodynamic fluctuations (caused by a carbogen gas challenge) extracted post hoc from BOLD imaging data collected in symptomatic IC stenosis patients undergoing CVR experiments (Donahue et al., 2016). However, we have determined that in many cases, the signal can be extracted from the fMRI data itself, either from a region rich in venous blood (such as in the superior sagittal sinus) , or more simply from the global mean average of the data .
Using the global mean data has a number of advantages relative to other methods, the most obvious being that it requires nothing other than the fMRI data itself -no external recordingsand processing is extremely simple, as there is no need to define anatomic regions a priori. Moreover, we showed in our recent study  that the global mean signal is highly correlated with the BOLD signal extracted from the SSS (i.e., maximum cross-correlation values are 0.81 ± 0.1), which indicates that essential components of global mean overlap with those of large veins (with little neuronal contamination). The drawback is that the global mean signal is essentially a temporally "blurred" version of the physiological regressor, because it contains contributions from voxels over a range of delay values. Moreover, each voxel contains fluctuations caused by local neuronal activity (which in this case are considered noise). To overcome these drawbacks, we have developed bootstrap sharpening method to recover the source signal, which is diagrammed in Figure 5. We have released an open source software package, "rapidtide" 1 , to perform this fitting procedure 1 https://github.com/bbfrederick/rapidtide and isolate and remove the sLFO signal from resting state (or task) fMRI data.
One of the most widespread methods of preprocessing rsfMRI to remove low frequency physiological noise is GSR. In this procedure, the mean signal of all voxels over time is regressed out of all of the BOLD time series prior to resting state analysis. While this does remove a significant amount of physiological signal, it has a serious drawback -the creation of spurious, negative correlations between brain regions (Carbonell et al., 2011). This is an unavoidable consequence of simple regression. The global mean signal is a summation of many copies of the sLFO signal with a range of delays reflecting the blood arrival time throughout the brain. Because the sLFO signal is low frequency (below 0.1-0.15 Hz), when copies of the signal over the range of delays found in the brain are summed, the signal strongly resembles the driving sLFO signal. However, this signal is not properly aligned in time in the vast majority of the voxels of the brain -it is shifted forward or backward relative to each voxels' signal, but the correlation with each voxel is generally high. Regressing out a delayed version of the driving signal at the wrong time delay unavoidably results in a lower amplitude, inverted version of the global signal being added to the voxel at the correct time delay. This will necessarily create artifactual negative correlations in GSR processed data, as shown by Erdogan in both real and simulated data (see Figure 6, reproduced from Erdogan et al., 2016). In a detailed comparison of static global signal regresson (sGSR) with the dynamic global signal regression (dGSR) performed by rapidtide, we found that by regressing the sLFO signal out of each voxel at the proper time delays, the efficacy of noise removal was improved. More importantly, we demonstrated that by removing the sLFO dynamically, negative correlations, which were present in the results of sGSR processing, were substantially attenuated . We would argue that many of these negative correlations are potentially spurious, being generated by the removal process itself, but this has yet to be confirmed. Aso et al. (2017) have found that an analogous noise removal procedure yield similar improvements in task based analyses, and also showed significant increases in reproducibility of analyses over time.
We also compared sLFOs derived from peripheral NIRS [HbO] with other LFO models , namely the model-based methods for respiration and cardiac listed above in "Potential causes of the low frequency oscillation." With high temporal resolution (TR of 400 ms), we found only small contributions (1-5%) of explained variance by the models considering respiration and cardiac variation. We also found that sLFO explained significantly more variance (up to 16%) when aliased respiration or cardiac signals do not play a critical role when fully sampled and filtered. In addition, sLFO was also substantially different from the variation models with very little temporal and spatial overlap , even though NIRS [HbO] is closely related to modulations in CO 2 concentration. This study showed that sLFO is not an artificial signal created by suboptimal acquisition parameters, but a real and pervasive physiological signal accounting for a substantial amount of the variance in the BOLD LFO, and which is not accounted of by previous methods. It should be noted that this procedure cannot remove stationary non-neuronal signal. However, this type of noise, should it exist, could be removed using ICA techniques, which are well suited to detecting (and removing) spatial patterns of synchronized signals.

Significance Determination
Because RIPTiDe analysis relies on the cross-correlation of low pass filtered signals, there is some concern that the correlations derived by the method may be spurious. Until recently, determination of the statistical significance of RIPTiDe metrics has been problematic. As the method is based on the cross-correlation of time series, it is tempting to use standard formulae, which determine significance based on number of degrees of freedom and the correlation coefficient. However, this greatly overestimates the statistical significance of the data obtained through our procedure, for two reasons. The first is that both our test regressor and the fMRI data are bandpass filtered to select the LFO component prior to correlation, which effectively reduces the degrees of freedom in the correlation. One method proposed to correct for this, specifically for the case of fMRI data, is to apply a correction factor to the degrees of freedom based on the portion of the spectrum retained by the filtering procedure (Davey et al., 2013). While this improves the estimation of significance in filtered correlations, to be strictly correct, both the exact transfer function of the filter function and the spectrum of the data being filtered must be known a priori and included in the calculation of the correction factor. In practice, the power spectra of fMRI data in general and the systemic low frequency oscillation signal in particular are not white (see section "Limitations" for further discussion on this topic), even over the limited frequency band of the sLFO, and vary in space. This second condition makes this procedure cumbersome.
The second, and more difficult aspect of the analysis to address, is that the peak correlation value within a range is selected to determine the "optimal" time lag, which necessarily serves to inflate the correlation value, and bias it toward more positive values. Proper application of correction for multiple comparisons requires accounting for the smoothness of the correlation function, which in turn is determined by the factors listed above. There are analytical methods for doing this (Olden and Neff, 2001), however they too are somewhat intractable for fMRI data.
While analytical calculation of the significance is difficult, there are two straightforward numerical methods to achieve this goal. The first, and simplest, method is to perform mismatched correlations. In this case the voxel timecourses are correlated with an sLFO signal from a different subject (or from the same subject at an extremely large time delay of several minutes). In this case any correlations between the timecourses are known to be spurious; a distribution of spurious correlations can be calculated to find various significance thresholds. We have employed this method in many of our analyses where the data permits (most recently here, Yao et al., 2019).
However, in cases where datasets are small, or time delayed sLFO signals are not available, there is a more general method which is also quite straightforward -we can estimate the distribution of null correlations using a Monte Carlo approach . The probe regressor is permuted by randomizing the time indices, preserving the distribution of intensity values but destroying any temporal correlations, and the RIPTiDe procedure (filtering, cross-correlation, and peak finding) is performed on this timecourse with the unpermuted regressor. The procedure is repeated a sufficient number of times that the null distribution of correlation coefficients can be estimated, so the p-Values of different correlation coefficients can be directly determined. This procedure is rapid -our analysis software estimates this distribution from 10000 iterations at the beginning of each refinement cycle for each subject's data in under 25 s. By default, the results presented by rapidtide are thresholded to the p < 0.05% level. This is probably too stringent in general, as there will be many voxels with true, lower correlation due to low rCBV. In practice for high quality data though, this is not too much of a problem.
It is important to note that for autocorrelated data, such as fMRI, one should only permute samples within exchangeability blocks to maintain the autocorrelation -time index shuffling is not generally correct. A better (but significantly slower) procedure is to randomize the phase of the Fourier transform of the data, which preserves the autocorrelation structure (Handwerker et al., 2012). However, in the case of RIPTiDe processing, both the regressor and shuffled data are filtered to the LFO band after shuffling, and the effect of this filter dominates the autocorrelation properties of the inputs to the correlation, so time shuffling performs well. This was verified in the HCP data, where we estimated the p < 0.05 correlation level both by the shuffling procedure described, and with the mismatched correlation method. The results were in close agreement, with an average spurious correlation threshold of 0.2.

APPLICATIONS OF THE sLFO SIGNAL
Up to now we have discussed the moving sLFO as a noise signal contaminating rsfMRI. But whether the moving signal is noise or signal is simply a matter of perspective. If we reframe the moving, blood-borne variance as a "signal" carrying information, rather than simply a nuisance complicating the interpretation of resting state connectivity, we can use it as a sensitive measure of hemodynamic function. The relative arrival time and strength of this signal as it propagates through the cerebral vasculature carries information regarding the distribution and timing of blood flow in the brain. We propose that the moving hemodynamic signal is a unique contrast mechanism in its own right, which provides information not currently available to other techniques. It may offer insights into subtle hemodynamic alterations that can be used as early indicators of circulatory dysfunction in a number of neuropsychiatric conditions, such as prodromal stroke, moyamoya, and Alzheimer's disease. Characterization of the sLFO signal throughout the brain allows for continuous monitoring of blood arrival time delay without a dedicated acquisition, with high sensitivity, and over a wide range of delay times, without any externally administered contrast.

Cerebrovascular Reactivity (CVR) Mapping
The moving sLFO signal can be used to determine cerebrovascular reactivity (CVR) to CO 2 changes both for clinical evaluation and to calibrate the BOLD response. CVR is typically measured using a hypercapnic challenge (either exogenously applied gas or breathhold) (Heyn et al., 2010;Bright and Murphy, 2013;Hare et al., 2013;Donahue et al., 2014). Accounting for the particular dynamics of the sLFO signal can give more accurate estimates of the regional response (Tong et al., 2011a). Golestani et al. (2016a,b) fully exploited the effect of the moving paCO 2 sLFO waveform on voxel wise BOLD to perform quantitative CVR mapping throughout the brain using only resting state data. By correcting for time delay and correlating the end-tidal CO 2 (a proxy for paCO 2 ) with the BOLD signal (with additional noise removal and modeling), they were able to determine the local BOLD response per percent change in paCO 2 simply from the resting state signal fluctuations.

Quantitative Blood Flow Imaging
In addition to assessing CVR, sLFO signals can be used to track blood throughout the brain to reveal both normal circulation patterns, and circulatory alterations in response to tasks, pharmacological challenges, or pathology. The temporal resolution of delay measured using cross-correlation depends on the length of the input signals, rather than the signal repetition time, so very fine delay distinctions are possible with normal fMRI data. Furthermore, because the delay measurement relies on a pattern of pseudorandom signal fluctuations over the entire timecourse, rather than a single tag, such as that used in arterial spin labeling (ASL), extremely long delay times can be measured. Newer ASL techniques, such as multidelay ASL offer some information on arrival time, but are still limited by the short lifetime of the ASL tag (under 3 s). Velocity-selective ASL does remove the restrictions on tag lifetime, and has been used for delays over 6 s (Qiu et al., 2012), however, the ability to quantify delay over a range of 10 to 100 s of seconds has not been demonstrated, and this does require another scan in addition to the resting state.

Healthy Circulation
We have used RIPTiDe analysis extensively in healthy subjects to quantify typical blood flow patterns, have validated these measurements against gold standard dynamic susceptibility contrast (DSC) imaging data (see Figure 7) collected in the same session , as shown in Figure 7, and have followed signal through the head from the inflowing carotid arteries to the exiting jugular veins  to establish that the sLFO signal does indeed move with the blood.
Because the bootstrap RIPTiDe analysis requires only fMRI data to assess hemodynamic parameters, it can be applied retrospectively to existing data [as in the myconnectome analysis described above , and a wider analysis of myconnectome, ABCD, and HCP data (Yao et al., 2019)]. We have also applied this technique to the resting state data from the Human Connectome Project , and were able to produce very detailed maps of average time delay and correlation strength throughout the brain (shown in Figure 8). This can provide a standard comparison dataset for young healthy controls.
It is important to note that the time delay measurement is always relative. Our standard is to set "zero time" at the peak of the histogram of delay values. As a result, there are necessarily always positive and negative delay times -positive delays tend to be in the later parts of the vascular tree -large draining veins, etc., while negative delays correspond to parts of the brain close to the source arteries (so they get blood before most other regions). This is a somewhat arbitrary choice -we could (and have) used an anatomical reference region such as the cerebellum as the delay time origin (Donahue et al., 2016) depending on the application.
We note that we consistently observe higher positive delay values in white matter (by a few seconds) than in the cortex; presumably this is because gray matter circulation is prioritized. This observation has been confirmed with DSC imaging, the gold standard MR technique for perfusion measurement.

Circulatory Pathology
This technique is not limited to studying healthy circulation, however. Compromised circulation due to arterial occlusion from to atherosclerosis, moyamoya disease, or stroke, can lead to extremely long delay times (up to 10 s of seconds) which are not quantifiable using conventional methods such as ASL or CT angiography. We demonstrated that time delay imaging with exogenous CO 2 manipulation can reveal and accurately quantitate extremely long delays in vivo Donahue et al., 2016). However, the exogenous CO 2 manipulation is not necessary -the hemodynamic delays arising from moyamoya disease can also be quantified using the endogenous sLFO signal from rsfMRI alone as a probe (Christen et al., 2014). Others have used the endogenous sLFO in a number of clinical applications.

Stroke
There have been several studies using global signal correlation delays to investigate hypoperfusion acutely and longterm after strokes. Lv et al. (2013) showed that "significant delay in BOLD signal corresponded to areas of hypoperfusion identified by contrast-based perfusion MRI" in 11 subjects acutely after ischemic stroke. Amemiya et al. (2014) found similar agreement in the results of BOLD delay and contrast perfusion measurements in five patients with chronic hypoperfusion and six with acute stroke. In 2016, Siegel assessed 130 subjects 2 weeks, and 3 and 12 months post-stroke onset (and 30 controls), and showed that increased BOLD lag delay was strongly correlated with decreased blood flow assessed with ASL. Furthermore, they found that removing the delayed hemodynamic signal somewhat normalized functional connectivity measurements (which were distorted by delayed hemodynamics) (Siegel et al., 2016). Khalil et al. (2017) assessed delay changes of sLFO BOLD signals among acute stroke patients in two separate studies. In one study, they found the delay maps were highly correlated with the time-to-peak maps derived from DSC-MRI in ischemic stroke. In another longitudinal study, they found the sizes of the extended-delay regions and the corresponding delay values changed according to the vessel conditions (Khalil et al., 2018). Recently, Nishida et al. (2018) has found the delay maps from patients with the arterial occlusive disease were correlated with the CVR maps from SPECT. These multimodal studies validated the method of BOLD delay and showed its great potential in perfusion assessment of cerebrovascular disease. Christen et al. (2014) were able to quantify delayed blood flow in moyamoya disease using a cross-correlation technique. In addition, Yan et al. (2018) have successfully detected perfusion deficits (from BOLD delay) in patients with Alzheimer's disease and mild cognitive impairment.

LIMITATIONS
There are two rather significant factors that can complicate the crosscorrelation method for tracking sLFO's -unfortunate spectral characteristics of the sLFO signal, and initializing the proper regressor to perform the calculation. These factors are discussed below.

Spectral Characteristics
The first problem is the random nature of the sLFO signal itself. Usually this benefits us in that it allows us to find the signal throughout the brain despite whatever local signal variations are present, and to determine its time delay -the cross-correlation of the sLFO with any other random variation (such as neuronal signals) is, in general, low, and the correlation with delayed copies of itself is high and strongly peaked at the appropriate time delay. The spurious correlation threshold can be determined through permutation or comparison of correlations between subjects, as discussed above, so we can decide when correlation is "real." However, all of this rests on the assumption that that sLFO signal is (1) truly random, or at least not determined by factors that will influence other noise signals found in the brain, and that (2) the signal is sufficiently "white" within the band of interest. While the first condition seems to be satisfied in general, the second frequently is not. The random nature of the signal means that by chance it sometimes has undesirable spectral properties which make it less suitable for our purposes. This leads to two problems -pseudoperiodicity and non-uniform spectra. The first is a special case of the second, but is common enough to be discussed on its own.

Pseudoperiodicity
The first uncontrolled quantity is pseudoperiodicity. From time to time, signal energy in the 0.09-0.15 Hz band will be strongly concentrated in one or more spectral peaks. Whether this is completely random, or due to some pathological or congenital condition that affects circulation is not known -it seems for the most part to be purely by chance, as it is occasionally seen when looking at multiple runs in the same subject, where one run is pseudoperiodic while the rest are not.
The effect of this is to cause the crosscorrelation between the probe signal and voxel timecourses to have more than one strong correlation peak. This means that in the presence of noise, or extreme spectral concentration of the sLFO, the wrong crosscorrelation peak can appear larger, leading to an incorrect delay estimation. This is particularly problematic if the pseudoperiod is shorter than the reciprocal of the search window (for example, if the search window for correlation peaks is between −5 and +5 s, and the sLFO has a strong spectral component at 0.1 Hz or higher, more than one correlation peak will occur within the search window). As the width of the search range increases, the spectral range of potentially confounding spectral peaks covers more of the sLFO frequency band.

Implications of pseudoperiodicity
The extent to which pseudoperiodicity is a problem depends on the application. In the case of noise removal, where the goal is to remove the global sLFO signal, and leave the local or networked neuronal signal variance, it turns out not to be much of a problem at all. If the sLFO signal in given voxel is sufficiently periodic that that the correctly delayed signal is indistinguishable from the signal one or more periods away, then it doesn't matter which signal is removed -the resulting denoised signal is the same.

Mitigation of pseudoperiodicity
While we continue to work on fully resolving this issue, we have a number of ways of dealing with this. First of all, spectral analysis of the sLFO signal allows us to determine if the signal may be problematic. Rapidtide checks the autocorrelation function of the sLFO signal for large sidelobes with periods within the delay search window and issues a warning when these signals are present. Then after delay maps are calculated, they are processed with an iterative despeckling process analogous to phase unwrapping. The delay of each voxel is compared to the median delay of its neighbors. If the voxel delay differs by the period of an identified problematic sidelobe, the delay is constrained to "correct" value and refit. This procedure greatly attenuates, but does not completely solve, the problem of bad sidelobes. A more general solution to the problem of nonuniform spectra will likely improve the correction.

Non-uniform Spectra
As noted before, the pseudoperiodicity is a special case of non-uniformity within the sLFO spectral region. In addition to peaks in the power spectrum, there can be gaps, which are also problematic.

Implications of non-uniform spectra due to spectral gaps
Non-uniform spectra will tend to distort the crosscorrelation between the sLFO and any given voxel signal. In addition to sidelobes which lead to the periodic correlation functions discussed above, gaps in the spectrum, especially in the higher frequency regions, can lead to blurry correlation functions. Our method relies on identifying peaks in the crosscorrelation waveform -the sharpness of these peaks depends on the higher frequency portion of the sLFO power spectrum. If the energy of the sLFO is concentrated in the lower frequency portion of the LFO band, the crosscorrelation peaks become broad, which makes the estimation of the peak location less accurate in the presence of noise. This makes delay maps less accurate, and will tend to lower the correlation values closer to the spurious correlation threshold.

Mitigation of non-uniform spectra
The most straightforward solution to non-uniform spectra is to prewhiten the sLFO and voxel signals (effectively flattening the peaks and troughs of the magnitude spectrum, while preserving phase) prior to performing the correlation; there are numerous variants of this procedure, known as generalized crosscorrelation (Knapp and Carter, 1976;Liang et al., 2015). The methods require some tuning to determine thresholds for magnitude recovery to avoid inflating noise, but are included as options in rapidtide. However, when using the generalize crosscorrelation, the resulting maximum correlation magnitudes are difficult to interpret, as they no longer directly represent the amount of variance explained by the sLFO regressor.

Inhomogeneous Time Delays
Finally, there is the problem of obtaining the sLFO regressor to begin with. In subjects with healthy circulation, the majority of the brain has delays which are relatively tightly clustered over a range of a few seconds; in this case the global mean is a good starting regressor. After multiple refinement passes, the regressor will converge on a stable candidate sLFO signal. However, pathology can lead to significant volumes of brain tissue with a large delay relative to the rest of the brain, but with a small spread of values around that delay. This can result in a global mean regressor which includes one or more "echoes"strong, delayed copies of the true driving regressor.

Implications of Inhomogenous Delays and Mitigation
Having a regressor with multiple delayed copies of the driving sLFO signal will lead to ambiguous delay values, and will keep the refinement process from converging. While it may be possible to clean the signal using a technique analogous to echo cancelation, it is generally easier to avoid the situation to begin with by starting from a region of homogenous delay values. A number of the cited studies in pathology have used regressors derived from the superior sagittal sinus, which is easy to locate and clearly homogenous. However, its location at the end of the vascular tree means that blood in that region may have traversed multiple distinct paths to get there, which could result in multiple delay components. Our current thinking is to use a cerebellar ROI to derive the starting sLFO regressor, as circulation in that region is undisturbed in a wide range of pathologies (Donahue et al., 2016).

Stationarity
One assumption that has been made throughout these discussions is that the time delay in a region is constant over time, or "stationary." While convenient, this is clearly a simplification. We know that these delays are not, in fact, completely constant over time. There are slight variations in delay time within HCP subjects between runs on the same day, and larger variations between days. There is no reason to suspect that there are not variations within runs as well -the correlation delays presented are averaged over the entire run. Moreover, when we calculate the delay maps from the HCP motor task data, we consistently see a regional, average 0.5 s decrease in the blood arrival time in the motor cortex relative to the values from the resting state scans, consistent with increased blood flow due to neuronal load -this is presumably an average change between the active and non-active periods leading to a shorter average delay time. This is a potential interesting area of research that as of yet, does not seem to have been explored. We have experimented some with a windowed version of RIPITiDe analysis analogous to dynamic connectivity studies, but have not worked out whether there is sufficient SNR to do this routinely.

CONCLUSION
Low frequency oscillations contribute significantly to the rsfMRI signal. The signal is defined by its characteristics rather than by its origins -in specific, a low frequency range of ∼0.01-0.15 Hz. In fact, there are likely many sources of signals in this spectral region. We have discussed a number of theories for the origin of LFOs; it is important to note that these theories do not conflictpower in this frequency band is likely due to some combination of the sources we describe. Because of this, it is more useful to talk about how the low frequency contamination in the fMRI signal behaves, and what can be done about and with it.
Unlike other physiological signals, such as respiration and cardiac contamination, LFOs cannot be separated from neuronal signals through spectral filtering -they must be modeled. We have observed that up to 30% of the low frequency signal power in the gray matter moves through the vasculature, and is carried with the blood into and through the brain. This portion of the signal is therefore amenable to detection, quantitation, and removal using cross-correlation techniques, which perform well at the task of noise removal, without introducing the significant artifacts seen with other methods, such as global signal regression.
Finally, we and others have demonstrated that this moving signal can be used as a probe to quantitate cerebral hemodynamic parameters, over a wide range of conditions, without the use of contrast or specialized imaging techniques, making this an ideal method for inferring hemodynamic information both in new studies and in retrospective analysis of existing datasets. Half of these studies have been in the last 2 years, demonstrating the increasing interest in the broad clinical application of the method. As large-scale public databases such as the Human Connectome Project and the UK Biobank become available, we expect the use of these techniques to continue to expand into new research areas.