Axial variation of deoxyhemoglobin density as a source of the low-frequency time lag structure in blood oxygenation level-dependent signals

Perfusion-related information is reportedly embedded in the low-frequency component of a blood oxygen level-dependent (BOLD) functional magnetic resonance imaging (fMRI) signal. The blood-propagation pattern through the cerebral vascular tree is detected as an interregional lag variation of spontaneous low-frequency oscillations (sLFOs). Mapping of this lag, or phase, has been implicitly treated as a projection of the vascular tree structure onto real space. While accumulating evidence supports the biological significance of this signal component, the physiological basis of the “perfusion lag structure,” a requirement for an integrative resting-state fMRI-signal model, is lacking. In this study, we conducted analyses furthering the hypothesis that the sLFO is not only largely of systemic origin, but also essentially intrinsic to blood, and hence behaves as a virtual tracer. By summing the small fluctuations of instantaneous phase differences between adjacent vascular regions, a velocity response to respiratory challenges was detected. Regarding the relationship to neurovascular coupling, the removal of the whole lag structure, which can be considered as an optimized global-signal regression, resulted in a reduction of inter-individual variance while preserving the fMRI response. Examination of the T2* and S0, or non-BOLD, components of the fMRI signal revealed that the lag structure is deoxyhemoglobin dependent, while paradoxically presenting a signal-magnitude reduction in the venous side of the cerebral vasculature. These findings provide insight into the origin of BOLD sLFOs, suggesting that they are highly intrinsic to the circulating blood.


1
Introduction 47 In functional magnetic resonance imaging (fMRI), there are 2 established physiological 48 bases of signal change: neurovascular coupling (NVC) and autoregulation. The former 49 involves a local blood flow increase of 30-70% which gives rise to a 0.5-2% blood 50 oxygenation level-dependent (BOLD) signal increase with around a 5-s delay (Buxton, 51 2013). This is the target phenomenon of fMRI as a tool for brain mapping, due to its limited 52 spatial extent mainly involving local arterioles (Hillman, 2014). In contrast, autoregulatory 53 responses in vessel diameter are found in a wide range of arteries including the internal 54 carotid or middle cerebral arteries (Hoiland et al., 2016). Detection of a compromised 55 response in vascular disorders has proven useful for clinical purposes (Murphy et al., 2011). 56 Importantly, there is no clear physiological distinction between these 2 phenomena as each 57 involves multiple pathways (Willie et al., 2014). Their traces in the fMRI signal are also 58 uniformly postulated to reflect the increased cerebral blood flow and eventual dilution of 59 deoxy-hemoglobin (Hb) in the postcapillary part of the vasculature, with the additional 60 effect of a local blood volume increase (Kim & Ress, 2016). 61 In efforts to improve the efficiency of fMRI, studies have revealed various systemic 62 physiological components in BOLD signal fluctuations. Physiological parameters, such as 63 cardiac pulsation (Chang et al., 2009), blood pressure, and end-tidal carbon dioxide (CO 2 ) 64 (Wise et al., 2004;Murphy et al., 2013), are considered artifact sources. The contamination 65 is expected to be emphasized in resting-state fMRI (rs-fMRI), where signals are evaluated 66 without trial averaging (Winder et al., 2017). However, discrimination between neuronal 67 and non-neuronal components has been a major challenge due to the lack of validation 68 techniques with spatial and temporal precision comparable to that of fMRI. Another source 69 of difficulty is the fact that many neural and non-neural parameters are intercorrelated in 70 this low-frequency range (Murphy et al., 2013). 71 A focus of recent studies exploring this matter is the spontaneous BOLD low-72 frequency oscillation (sLFO, < 0.1 Hz) possibly encompassing multiple artifact sources 73 (Zhu et al., 2015;Tong et al., 2017). One well-known sLFO source is the respiratory 74 volume fluctuation involving a chemoreflex loop (Birn et al., 2008). An sLFO in systemic 75 blood pressure, known as the Traube-Hering-Mayer wave, has been shown to originate 76 from another autonomic loop (Guyton & Harris, 1951;Julien, 2006). Moreover, associated 77 sLFOs in blood flow and velocity have been found (Killip, 1962;Fagrell et al., 1977) and, 78 later, transcranial Doppler ultrasonography and optical methods confirmed their traces 79 (volume fraction of erythrocytes to the blood volume, typically around 40%). The oxygen 136 extraction fraction (OEF) represents the only source of deoxy-Hb under the assumption of 137 100% oxygen saturation (SaO 2 ) in the inflow. 138 Variations of this base susceptibility can occur due to the local Hct and SaO 2 changes 139 and, in fact, have been shown to cause intersession variabilities (Cohen et al., 2002;140 Tuunanen & Kauppinen, 2006); however, within-session fluctuations have rarely been 141 considered (Thomas et al., 2000). For example, even at a constant Hb density and oxygen 142 partial pressure, CO 2 fluctuation, a driving factor of vasomotion, alone can modify SaO 2 143 through pH changes (Collins et al., 2015a). Although the assumption of constant base 144 deoxy-Hb concentration may be sufficient for modeling its dilution by NVC (Ogawa et al., 145 1998), other effects might not be negligible in non-neuronal fluctuations (Fig. 1b). 146   In this exploratory study, we conducted 3 investigations to further advance our 178 knowledge on the BOLD lag structure and its underlying physiology. We first evaluated if 179 changes in blood transit velocity are embedded in the BOLD low-frequency phase to 180 confirm its behavior as a virtual tracer [a preliminary analysis of these data was presented 181 previously as a poster (Aso et al., 2017b)]. Next, we investigated the effect of eliminating 182 the lag structure from task-based fMRI, which had not been tested previously. Finally, we 183 used multi-echo imaging to assess the components of the BOLD signal that determine the 184 lag structure. One of the recent approaches toward fMRI denoising has focused on S 0 185 fluctuations (signal at TE = 0, which is the baseline MR signal from the fluid 186 compartments), as it is weakly associated with neural activity (Posse et al., 1999;Wu et al., 187 2012;Kundu et al., 2012;Yen et al., 2017). In contrast, the total-Hb sLFO, detected by 188 NIRS, is interpreted as a local CBV change under the assumption of a constant Hct, which 189 should also affect the non-BOLD component via changes in plasma volume and inflow 190 (Rostrup et al., 2005). Notably, the contributions of the T2* and S 0 components may differ 191 from NVC contributions to the BOLD lag structure; hence their impact may have been 192 overlooked in studies based on trial averaging. From the influences of 3 different fMRI task 193 paradigms, including a simple reaction-time visuomotor task, a short breath-holding task, 194 and a hyperventilation task on the neural and non-neural components of the fMRI signal, 195 we sought further validation of our hypothetical model of the BOLD lag structure. 196 2 Materials and methods 197

Participants and experimental procedures 198
Twenty-one healthy participants (8 women, 19-26 years of age) participated in Experiment 199 1; only 1 person was excluded from the analysis because of an abrupt head motion, which 200 prevented BOLD lag mapping (see Data processing). The remaining 20 participants 201 performed the sparse visuomotor and 10-s breath-holding tasks, but the hyperventilation 202 task was only performed by 18 participants, as it was introduced after the first 2 individuals 203 had concluded their participation. Another 21 participants were recruited for Experiment 2, 204 involving multi-echo acquisition, all of whom performed the above 3 tasks but with an 18-s 205 version of the breath-holding task. To avoid vigilance level fluctuations, all MRI sessions 206 were scheduled in the morning and the participants were encouraged to sleep well the 207 previous night. 208 The protocol for this study was approved by the internal ethics review board of Kyoto 209 university. The participants provided written informed consent in advance, according to the 210 Declaration of Helsinki, for the analysis of anonymized MRI scans and simultaneously 211 acquired physiological data. 212

Image acquisition 213
A Tim-Trio 3-Tesla scanner (Siemens, Erlangen, Germany) with a 32-channel phased-array 214 head coil was used for MRI acquisition. For Experiment 1, T2*-weighted echo-planar 215 images were acquired using multiband gradient-echo echo-planar imaging (EPI) ( (1,080 volumes) were acquired for each of the 3 task conditions. The same pulse sequence 220 program was used in Experiment 2 but with multi-echo settings: TE1 = 7.76 ms and TE2 = 221 25.82 ms for the first 6 participants and TE1 = 11.2 ms and TE2 = 32.78 ms for the 222 remaining 15 participants. A smaller multiband factor of 2 was selected to allow for the 223 short TE in combination with parallel imaging using GeneRalized Autocalibrating Partial 224 Parallel Acquisition. Other acquisition parameters were: TR = 1,300 ms, flip angle = 65°; 225 36-slice interleave, FOV = 256 × 192 mm 2 , 64 × 48 matrix, and 3.5-mm slice thickness. 226 Three 7-min (323 TR) runs were acquired. Seven participants in Experiment 2 underwent 2 227 additional runs with a shorter TR of 700 ms and a flip angle of 45° to examine the 228 sensitivity of the respiration-related signal component to inflow modulation. At the end of 229 every experimental session, a 3-dimensional (3D) magnetization-prepared rapid acquisition 230 with gradient echo (MPRAGE) T1-weighted image was acquired for obtaining anatomical 231 information (Aso et al., 2017a). A dual-echo gradient-echo dataset for B 0 field mapping 232 was also acquired after the BOLD scan in the same orientation. 233

Task conditions 234
Throughout the experimental session, task instructions were presented via a liquid crystal 235 display (LCD) monitor inside the scanner room, viewed through a mirror. Beat-to-beat 236 fluctuations in the mean arterial pressure and heart rate were obtained via a non-invasive 237 MR-compatible device (Caretaker, BIOPAC Systems, Inc., Goleta, CA, USA). Careful 238 instructions were provided to the participants on how to avoid motion, especially during the 239 respiratory challenges. 240

Sparse visuomotor task 241
A simple visuomotor task with a varying intertrial interval of 6 to 24 s was performed 242 during the first run. Participants were instructed to press a button with their right index 243 finger as soon as the computer screen changed from "Please hold still" to "Press the button." 244 The screen returned to "Please hold still" at the button press or after 3 s, if the participant 245 had not pressed the button. 246

Breath-holding task 247
To minimize head motion induced by the tasks, the 2 respiratory challenges were adapted 248 to be less strenuous than in earlier studies. In Experiment 1, the breath-holding task was 249 cued by a "Hold your breath" instruction on the screen, at which point the participants were 250 asked to immediately hold their breath, irrespective of the respiration phase. The holding 251 periods lasted 10 s and were separated by 90-s intervals. This short duration was selected to 252 minimize strain that can cause body movements, while evoking a detectable autoregulatory 253 response (Murphy et al., 2011). In Experiment 2, a longer holding period of 18 s after a 254 brief inhalation for 2 s was used to evoke a more pronounced vasodilation. 255

Hyperventilation task 256
The hyperventilation task involved paced breathing at 0.2 Hz for 25 s, separated by a 30-s 257 rest. Each 5-s cycle began with the screen presenting "Please inhale" for 1.5 s, followed by 258 "Please exhale slowly," lasting 3.5 s. Participants were instructed to breathe as deeply as 259 possible, while avoiding head movement. A short inspiration period was selected to 260 suppress motion by minimizing movements in the thoracic cage and spine. The global mean signal was used to select the initial seed that defined the reference 290 phase (lag = 0). First, voxels presenting a cross-correlogram peak at 0 with the global signal 291 were determined. The time course averaged over this set of voxels served as the initial 292 reference. In each step of the recursive procedure toward up-and downstream, a cross-293 correlogram was calculated between the time series obtained from the previous seed voxels 294 and every undetermined voxel to find a set of voxels with a peak at ± 0.5 s, which then 295 served as the new seed. This tracking part retained some voxels without any lag values 296 because cross-correlogram peaks < 0.3 were not used following earlier works (Tong et al., 297 2017). These voxels were later filled 1 by 1 with average phases from voxels with similar 298 time courses and correlation coefficients > 0.3. There were single isolated holes even after 299 this procedure, which were filled by linear interpolation, using the 6 neighbors. There is a 300 concern that this correlation coefficient threshold is too low to accurately claim a 301 significant correlation. However, recursive tracking involves finding the cross-correlogram 302 peak precisely at ± 0.5 s, which conveys different information from its height. Besides, 303 most earlier works empirically supporting the biological significance of this phenomenon 304 involved no thresholding. When the threshold correlation coefficient was increased to 0.6, 305 most brain voxels required the hole-filling procedure, but we still obtained lag maps (by 306 between-voxel intra-class correlation (2,1) > 0.4) in 16 out of 20 participants during resting 307 state and 17 during 10-s breath holding. 308

fMRI analysis on "cleaned" BOLD datasets 309
Individual BOLD data from Experiment 1, after the above pre-processing steps including 310 the motion parameter regression, served as the reference or "raw" dataset. GSR with a low-311 pass filtered global signal, a normal GSR, global scaling implemented in SPM12, and 312 without the perfusion lag structure ("deperfusioned") were compared with the raw dataset. The temporal profile of the instantaneous velocity (green curves) was roughly in 445 phase with the global BOLD response but had different onset and peak timings, indicating 446 different physiological bases. To evaluate the relationship of this phenomenon with the 447 vascular structure, the inlet, center, and outlet regions were separately analyzed (Fig. 4) the motion-related variances were removed (Fig. 5). However, after the deperfusioning 469 procedure, primary motor cortex activation was successfully detected by decreased 470 interindividual variances (black circles). The effect of these procedures was not uniform 471 between the respiratory challenges (Fig. 6), but some interpretable clusters were selectively 472 captured by the "deperfusioned" signals despite the reduced cluster number.  (Fig. 7). The correlation of the 499 neuronal response with the global signal (or the extracted sLFO) was near 0, as described 500 above, but it varied across participants and tended to be positive. This trace of NVC may 501 have created the spurious deactivation after regression. Notably, this effect was very weak 502 after deperfusioning across all 3 conditions (green plots), despite a large amount of variance 503 removed by the procedure. 504  In lag maps from the 3 T2*-weighted images, there were some effects of the 530 respiratory challenges, but the gross cerebral vascular structure was preserved across tasks; 531 the periventricular regions and major venous sinuses were uniformly found downstream 532 (i.e., with negative arrival time) of the global signal phase, while the cortical territory of the 533 middle cerebral arteries exhibited earlier arrival (Fig. 9A). Only the S 0 image presented a 534 different lag structure, according to the image similarity (Fig. 9A, right panels). 535 Temporal analysis of the signal components revealed significant main effects of both 536 region [F (2,532) = 16.877, p < 10 -6 ] and T2* weighting [F (2,532) = 280.786, p < 10 -6 ] 537 (Supp . Fig. 1A). The S 0 time series failed to show a correlation with the T2*-weighted 538 signals, but the z-value in the inlet (i.e., the arterial side) differed significantly from that in 539 other regions (p < 0.05, Tukey's HSD). A region effect was also found for phase 540 relationships (Supp. Fig. 1B). Similarly, the main effects for both region [F (2,532) = 541 27.548, p < 10 -6 ] and T2* weighting [F (2,532) = 44.679, p < 10 -6 ], as well as their 542 interaction [F (4,532) = 24.902, p < 10 -6 ], were significant. The phase of the TE1 signal, 543 which is less T2*-weighted than that of the BOLD signal, gradually advanced, finally 544 showing a phase lead in the outlet region, further supporting an interaction between the 545 signal components and vascular regions. These signal phase dissociations within regions 546 are displayed in Supp. Fig. 1C. Significant differences among the 3 T2*-weighted signals 547 were also found after the post-hoc test (p < 0.05). The signal-region interaction was evident 548 in the signal response to the respiratory challenges shown in Fig. 9B. Note that the traces 549 contain higher frequency components that were eliminated prior to lag mapping. In contrast 550 to the changes in T2* responses for both phase and magnitude, the S 0 component was 551 stable across vascular regions. involvement. Importantly, this spatial pattern was not that of typical motion artifacts that 584 can accompany volitional respiration in spite of the careful instruction. 585 Finally, Supp. Fig. 2C shows the data from a subset of participants, obtained using a 586 different TR/flip-angle setting to manipulate the inflow effect. The T2* response was 587 smaller than that shown in Fig. 2A, presumably due to the short TR. The different TR also 588 contributed to the rich high-frequency components by the fast sampling rate. The slow S 0 589 change was also diminished, but the respiration-related fast component was relatively 590 preserved suggesting the absence of a strong inflow effect. 591

4
Discussion 592 The principal findings of this study are summarized as follows. First, based on the 593 instantaneous phase difference within the BOLD lag structure, we observed a small blood 594 flow velocity change selectively in the inlet region of the vasculature. Next, the complete 595 elimination of the lag structure reduced interindividual variance and spurious deactivation, 596 supporting our hypothesis that NVC could be observed more specifically by this 597 deperfusioning procedure. This finding is in agreement with the results of earlier work on 598 resting-state fMRI (Erdoğan et al., 2016). Finally, the lag structures in the S 0 (or non-599 BOLD) component did not correlate with that from T2*, either spatially or temporally. We 600 also found a vascular region-dependent change in the T2* sLFO, with a decreased 601 amplitude in the outlet part close to major veins, in contrast to the S 0 response that 602 remained constant; this finding replicates a previous observation in the raw BOLD signal 603 (Aso et al., 2017a). The S 0 component exhibited a unique brain region-dependent response 604 to the respiratory phase, suggesting that certain perfusion parameters specifically contribute 605 to this component, but not the perfusion lag. Overall, the BOLD low-frequency phase 606 behaved as a deoxy-Hb-based virtual contrast agent in the present data, leaving a global 607 noise component for the fMRI analysis. 608 The observation that the velocity on the arterial side exhibits changes alongside 609 respiratory variations is consistent with the findings of previous reports using transcranial 610 Doppler ultrasonography (Malatino et al., 1992;Poulin et al., 1996). This information was 611 extracted from the BOLD lag structure, which itself presented autoregulatory response 612 consistent with earlier work (Murphy et al., 2011). During the initial whole-brain CBV 613 increase in response to autoregulatory vasodilation, a sole increase in the inflow should first 614 occur to meet the volume demand. It is therefore reasonable that this effect is absent in the 615 outlet (i.e., the venous side of the gross vasculature). This observation seems to support the 616 model in which the BOLD lag structure is derived from an axial non-uniformity in the 617 vessels, already present in the inflow (Tong et al., 2018). A distinct mechanism of the lag 618 structure was also suggested by the diminished magnitude observed in the outlet side of the 619 vasculature, since a CBF increase should evoke larger response in the downstream (Krings 620 et al., 1999). 621 It is unclear what proportion of this axial variation is systemic, i.e., originates from 622 the autonomic loops, mediated by peripheral baro-or chemoreceptors. However, even when 623 the neural activity is contributing to the sLFO time course as demonstrated previously (Aso 624 et al, 2017), the resulting lag structure largely reflects the vasculature. In this work, we 625 focused on non-neuronal mechanisms to account for the BOLD lag structure as much as 626 possible, in the hope that it may ultimately help achieve a better understanding and provide 627 improved modeling approaches of the fMRI signal. 628

Source of BOLD low-frequency oscillation signals 629
Previous studies on sLFOs have reported that both Hb species fluctuate, but with varying 630 phase differences that are selectively found in the brain (Obrig et  such as oxygen saturation (Fantini, 2014). In support of the conventional theory, 636 Rayshubskiy and colleagues reported, in their human intraoperative study, that slow Hb 637 oscillations correlated with vasomotion in the superficial arteries (Rayshubskiy et al., 2014). 638 However, it remains unclear whether an equivalent vasomotion exists in the non-arterial 639 vessels to fully account for the observed lag structure. Hence, it is worth considering other 640 sources of deoxy-Hb variation. 641 The concept of vasomotion stems from an active diameter change in the precapillary 642 vessels, driving local velocity fluctuations, termed "flowmotion" (Intaglietta, 1990). This 643 flowmotion can reportedly accompany the fluctuation of local Hct that should affect deoxy-644 Hb density (Fagrell et al., 1980;Hudetz et al., 1999). Another possible source for the 645 deoxy-Hb fluctuations is a change in SaO 2 that ranges from 94-98% in the artery 646 (Intaglietta et al., 1996;Collins et al., 2015). For example, the respiration-related BOLD 647 signal component is supposed to be mediated by the blood CO 2 level and pH, which can 648 shift the oxygen dissociation curve (Birn et al., 2006;Chang & Glover, 2009a). These 649 parameters are considered to fluctuate in the blood as part of the autonomic loop, possibly 650 driving local vasomotion, which persists after denervation (Morita et al., 1995). As 651 mentioned above, the phase difference between the 2 Hb species remains to be elucidated 652  (Fagrell et al., 1980;Mayhew et al., 1996). 660 Furthermore, the reduction in T2* LFO amplitude in the outlet side can be explained by the 661 high tissue deoxy-Hb density, which likely diminishes the proportional effect of intrinsic 662 deoxy-Hb fluctuations, unless the OEF is completely coupled to this variation. Importantly, 663 temporal dispersion alone would not fully account for the amplitude reduction, as it was 664 only found in the venous side, despite the fact that the lag structure was tracked both up-665 and downstream from the global phase. These results provide good contrast with the stable 666 S 0 response, reflecting its insensitivity to oxygen saturation. Although it would be too 667 challenging to incorporate complex rheological parameters, a reconsideration of the 668 constant deoxy-Hb assumption may help improve BOLD signal modelling. 669 In the hyperventilation condition, we observed a fast response to each ventilation 670 cycle, accompanying blood pressure changes. This is consistent with reports using 671 optimized acquisition techniques (Dresel et al., 2005;Pattinson et al., 2009), supported by 672 anatomical (Simonyan & Jürgens, 2003), as well as electrophysiological (Radna & 673 MacLean, 1981) studies. The premotor peak at coordinates [+56, 0, 40] was also very close 674 to the reported activation site for volitional respiration (McKay et al., 2008). Although the 675 effect of respiratory movement cannot be fully excluded, the spatial pattern in Supp. Fig.  676 2B is not that of a typical motion artifact centered on the brain surface (Krings et al., 2001). 677 In healthy participants, inhalation increases systemic venous return through decreased 678 intrathoracic pressure, causing an elevation of cardiac output with some delay. In contrast, 679 exhalation is considered to cause CBV increases through elevated cerebral venous pressure. 680 To our knowledge, the timing order of these events has not been studied at the precision of 681 the current data; further studies are needed to determine the source of this S 0 fluctuation 682 (Yen et al., 2017). The only available clue in our results is the spatial distribution, such as 683 the interesting anterior-posterior segmentation (resembling the unique "S 0 lag structure" in 684 has become a de facto standard. It is indeed computationally closer to our deperfusioning in 696 that the regional phase difference is somehow tolerated. However, because this practice 697 lacks a strong theoretical background (Chang et al., 2009;Liu et al., 2017), currently, the 698 identification and elimination of bodily movements and physiological noise are more 699 widely recommended. There are various approaches to this end, such as simultaneous 700 physiological measurements (Chang & Glover, 2009b), as well as data-driven methods that 701 only use fMRI data (Smith et al., 2013). To date, however, objective criteria for 702 distinguishing neural activity from noise components remains an issue (Salimi-Khorshidi et 703 al., 2014). 704 In the present study, the lag structure was treated as a broadly distributed, structured 705 noise for fMRI. Indeed, it can be partly eliminated by sophisticated denoising techniques 706 (Aso et al., 2017a). However, the specificity of lag mapping in isolating information on a 707 purely vascular origin remains unclear. For example, measurements of velocity changes 708 critically depend on a recursive lag-tracking method that incorporates the gradual change in 709 LFO over regions (Tong & Frederick, 2014). Adaptation for changes in the waveform that 710 may arise from different paths of the blood was demonstrated to increase the 711 reproducibility of the lag map (Aso et al., 2017a). However, as the changes in waveform 712 can also reflect NVC, removing the whole lag structure may lead to type II errors in the 713 fMRI results. Hence, the favorable impacts of the deperfusioning procedure that we 714 observed on the fMRI results are clearly insufficient to prove the advantage of this 715 technique and require further confirmation. 716 Importantly, the detection of the lag structure itself largely depends on the data 717 quality, especially in terms of head movement. When a head movement results in a 718 synchronized deflection that exceeds the LFO amplitude, it would obscure the phase 719 variation. However, it can be also questioned if the correlational structure of neural activity 720 is reliably detected from such motion-contaminated data. In general, NVC should have a 721 limited spatial extent and signal magnitude without time-locked averaging (Power et al., 722 2012). In turn, successful tracking of a lag structure may even be considered as evidence of 723 "clean" data. The elimination of this identified lag structure can be a relatively 724 straightforward approach to reduce structured physiological noise (Caballero-Gaudes & 725 Reynolds, 2017). 726 In conclusion, by investigating various aspects of the BOLD sLFO, we compiled 727 supporting evidence for a component intrinsic to flowing blood that has been a focus of 728 interest in earlier works (Tong et al., 2017). To establish a framework by which the fMRI 729 signal can be fully modeled, more detailed characterization of the lag structure as part of 730 the "global noise" is needed (Glasser et al., 2018). 731