The impact of real-time fMRI denoising on online evaluation of brain activity and functional connectivity

Objective. Comprehensive denoising is imperative in functional magnetic resonance imaging (fMRI) analysis to reliably evaluate neural activity from the blood oxygenation level dependent signal. In real-time fMRI, however, only a minimal denoising process has been applied and the impact of insufficient denoising on online brain activity estimation has not been assessed comprehensively. This study evaluated the noise reduction performance of online fMRI processes in a real-time estimation of regional brain activity and functional connectivity. Approach. We performed a series of real-time processing simulations of online fMRI processing, including slice-timing correction, motion correction, spatial smoothing, signal scaling, and noise regression with high-pass filtering, motion parameters, motion derivatives, global signal, white matter/ventricle average signals, and physiological noise models with image-based retrospective correction of physiological motion effects (RETROICOR) and respiration volume per time (RVT). Main results. All the processing was completed in less than 400 ms for whole-brain voxels. Most processing had a benefit for noise reduction except for RVT that did not work due to the limitation of the online peak detection. The global signal regression, white matter/ventricle signal regression, and RETROICOR had a distinctive noise reduction effect, depending on the target signal, and could not substitute for each other. Global signal regression could eliminate the noise-associated bias in the mean dynamic functional connectivity across time. Significance. The results indicate that extensive real-time denoising is possible and highly recommended for real-time fMRI applications.


Introduction
Real-time functional magnetic resonance imaging (rtfMRI) is a system for evaluating online brain activity as soon as an imaging volume is acquired [1,2]. Applications of rtfMRI have expanded from online inspection of signal quality [1] to brain-computer interface for manipulating equipment with brain signal [3] and, most significantly, to neurofeedback for self-regulating brain activity [4,5]. The primary challenge of rtfMRI is performing image processing in a short time. The processing includes image reconstruction, data processing for denoising and extracting the target signal, and displaying the processed signal. All these processes need to be done in less than the repetition time (TR) of volume acquisition, usually a few seconds or less than 1 s, to keep the pace of real-time data presentation without accumulating a delay. Thanks to the advancement of computational hardware, image reconstruction and signal presentation times are becoming negligibly short. In contrast, the time of data processing for signal denoising is still not negligible and its impact on signal quality has not yet been well-examined [6].
The blood-oxygen-level-dependent (BOLD) signal in fMRI contains various noise and requires many denoising steps to get a reliable estimate of neural activity. Although comprehensive denoising is critical to obtain a reliable estimate of brain activity, its execution in a rtfMRI environment is so cumbersome that most rtfMRI applications have used minimal denoising steps compared to offline analyses. The real-time estimate of neural activity with a limited denoising process cannot have comparable quality to that available in offline processing. This could hinder the reproducibility of the result of rtfMRI applications like neurofeedback [7]. This limitation has been being compensated for through the implementation of a fast computation algorithm, the development of a processing method adapted for real-time computation, and the improvement of the computational capacity of personal computers and graphics processing units (GPU) [8][9][10]. Several packaged frameworks implementing advanced real-time processing (RTP) have been released as well [11][12][13]. These improvements in realtime signal processing for fMRI have been reducing the gap between offline and real-time evaluations of neural activity, which should improve the efficacy of rtfMRI applications [14].
Notwithstanding that these advanced online processing systems could boost the reliability of neural activity estimates in rtfMRI, their use has not been considered obligatory in rtfMRI applications. Moreover, the impact of real-time denoising steps on the quality of the neurofeedback signal has not been well-recognized, demonstrated by the fact that most neurofeedback studies do not report real-time denoising steps in detail [6]. The CRED-nf checklist [15] is one of the efforts to resolve this undesired state of the field by encouraging researchers to report details of the RTP steps, while it does not indicate what processing is recommended for which rtfMRI application.
We have not known which type of RTP has a benefit in which rtfMRI application due to the dearth of studies quantifying the denoising benefit in improving neurofeedback signal quality. While a few reports are available about the benefit of fMRI RTP [8,[16][17][18], these reports focus on specific parts of the process and do not cover comprehensive denoising steps as employed in offline analysis. Heunis, Lamerichs, Zinger, Caballero-Gaudes, Jansen, Aldenkamp and Breeuwer [6] has described this state of research and highlighted the need for a methodological study that quantifies the effect of the denoising steps on rtfMRI signal quality.
This study aimed to comprehensively evaluate the benefit of fMRI RTP in improving the quality of online neural activity estimates. Specifically, we performed RTP simulation analyses to evaluate how much each stage of RTP reduced the noise effect on the online estimate of voxel-wise signal and functional connectivity (FC). The evaluated RTPs included slice-timing correction, motion correction, spatial smoothing, signal scaling, and regression of low-frequency fluctuations (high-pass filtering (HPF)), motion parameters, motion derivatives, global signal, white matter/ventricle mean signals, and physiological noises modeled by the image-based retrospective correction of physiological motion effects (RETROICOR) [19] and the respiration volume per time (RVT) [20]. The examined noises included motion and physiological noises, which are the most prevalent noises in the fMRI signal and FC [21][22][23][24][25].
At evaluating the RTP performance, we should also consider the cost of computation time. If the benefit of an added process is negligible compared to its computational cost, there would be no need to use it in the limited time of rtfMRI. Also, since using more regressors in noise regression requires more scans before obtaining a usable signal to avoid overfitting [10], we should not include ineffective regressors in the RTP. Thus, we evaluated the computation time of each processing and discussed the number of regressors as the cost of RTP.
In addition, we examined the advantages and disadvantages of different approaches of GLM for real-time analysis, such as cumulative GLM (cGLM) and incremental GLM (iGLM) [8]. The cGLM performs an ordinary regression at each volume using the whole history of available online data, while iGLM updates only the latest estimation without keeping data history. Most rtfMRI applications have used iGLM with its short computation time and small memory consumption. However, iGLM cannot update regressor values retrospectively. This limitation could be critical when a regressor is created online and can be improved by a retrospective update, such as HPF and physiological noise regressors. We demonstrated that this limitation of iGLM hinders real-time physiological noise regression performance so that cGLM is preferred in RTP when using physiological noise regression.
This study performed a series of simulation analyses structured as follows. Firstly, we compared the iGLM and cGLM regarding their quality of onlinemade regressors, such as HPF and physiological noise models. Since the result indicated that cGLM is preferred for RTP, the simulated RTP system implemented cGLM. We used this system for the RTP simulation with resting-state fMRI data and rtfMRI neurofeedback task data targeting the left amygdala activation (LA-NF) [26]. To evaluate RTP's benefit in noise reduction, we calculated the amounts of variance explained by the motion and physiological noises in the real-time estimate. For the resting-state data, evaluation was conducted for voxel-wise signals and dynamic FC in the whole-brain regions. For the LA-NF task data, we conducted the noise analysis on the neurofeedback target signal and calculated the signals' correlation between the RTP and offline process to measure the integrity of the neurofeedback signal.

Human data
We performed the simulation analyses on fMRI resting-state data from 87 healthy participants and rtfMRI neurofeedback task data from 22 participants (14 major depressive disorder patients). The supplementary material (available online at stacks.iop.org/ JNE/18/046092/mmedia), 'fMRI samples for the simulation analyses' describes the details about the participants. All human neuroimaging data were acquired in the studies [26] and [27] approved by the Western Institutional Review Board, Puyallup, WA. Human research in this study was conducted according to the principles expressed in Declaration of Helsinki. All subjects gave written informed consent to participate in the study and received financial compensation.

Comparisons between the incremental and cumulative GLM
Most rtfMRI applications with regression use iGLM [8], which calculates the latest estimates (beta and residual) without keeping data history. Although this has a notable advantage in its short computation time and small memory requirements, iGLM cannot update past regressor values because the current estimation is made based on a previous estimate. This could be a disadvantage in RTP when the regressor values are formed online and the regressor model could be improved with increasing sample points, such as HPF and physiological noise regressors. The cGLM, in contrast, performs an ordinary regression at each sampling, using the whole history of available online data [6]. While the cGLM may take a longer computation time than iGLM, the online update of the past regressor values could improve the accuracy of the online-made regressors.
The aim of the present simulation analysis was to compare the quality of online-made regressors between the iGLM and cGLM approaches for HPF and physiological noise models. Specifically, we evaluated the frequency filtering performance of the HPF regressors and the accuracy of online-made physiological noise regressors of RETROICOR and RVT.

High-pass filtering regression
The Legendre polynomial regressors and discrete cosine transform (DCT) basis sets have been used for HPF in AFNI (https://afni.nimh.nih.gov/) and SPM (www.fil.ion.ucl.acuk/spm/), respectively. We evaluated the frequency filtering performance of these regressors in iGLM and cGLM. We generated a random white-noise signal with a length of 250 samples, at a 2 s sampling interval and unit absolute power for all frequencies. GLM was applied to this signal with the Legendre polynomial or DCT regressors at 50, 100, 150, 200, and 250 lengths to simulate the real-time evaluation with a limited number of samples. Regressors in iGLM were made for the length of 250 TRs, and then the first part of them was extracted. This simulated that iGLM cannot update past regressor values-the same approach used in Bagarinao, Matsuo, Nakai and Sato [8]. In contrast, regressors in cGLM were updated at each length with the order adjustment (e.g. order of the Legendre polynomials or the cycles of DCT) for a designed pass frequency in each sampling length. The order of Legendre polynomials was calculated as 1 + floor(d/150), where d is the scan duration in seconds (default in AFNI). The threshold frequency of DCT was 1/128 Hz (default in SPM). The simulation was repeated 1 000 times with different random white noise signals. The frequency filtering performance was examined with the absolute power spectra of the residual signal.
We also performed the same simulation for resting-state fMRI signals to examine the difference of GLM approaches in frequency-filtering performance. The resting-state images of 87 participants (see supplementary material 'fMRI samples for the simulation analyses') were applied slice-timing correction, motion alignment, and spatial smoothing using the RTP simulation system (supplementary material 'RTP simulation system implementation'). Then, the HPF regression was applied to the signal timecourse of voxels within the intersection of the brain mask and a signal mask (3dAutomask in AFNI was applied to the functional images to make the signal mask). The signal in each voxel was normalized to z-score before applying the GLM. The average result of the whole-brain voxels for all participants is reported.

Physiological noise model
RETROICOR and RVT regressors were evaluated in the online GLM comparison. We used AFNI's RetroTS.py with the default parameters to make the regressors. RETROICOR regressors were Fourier basis sets with the phase synchronized to cardiac or respiration signal time-course [19]. A basis set of four regressors was made for each of the cardiac and respiration signals. RVT is a respiration volume per time calculated by detecting the top and bottom peaks of the respiration time-course to measure respiration volume [20]. The five time-shifted regressors were made.
The simulation examined how the online calculation of the RETROICOR and RVT regressors diverged from that in the offline calculation. We used respiration and cardiac signals recorded with resting-state fMRI (supplementary material 'fMRI samples for the simulation analyses'). The simulation was initiated at the time equal to 70 s in order to wait to obtain enough samples to create the regressors. In both GLM approaches, the regressors were recalculated at every TR (=2 s) using the accumulated history of cardiac and respiration signals from the start of the scan until the current TR. There was no difference in the regressor calculation Overview of the fMRI RTP simulation system. The system monitors a directory where an fMRI data volume file is created in real-time. The WATCH module detects new file creation and reads the data to send it to a process sequence of TSHIFT (slice-timing correction), VOLREG (motion alignment), SMOOTH (spatial smoothing) and REGRESS (scaling and general linear model analysis). The output (residual of the noise regression as a noise-removed fMRI signal) can be sent to an independent application of neurofeedback or a brain-computer interface. The simulation was performed by copying a volume image of a pre-scanned fMRI data into a monitored directory. between the two GLM approaches. The only difference was the update method of regressor values at GLM; in the cumulative approach (cGLM), all values were replaced in the updated regressor values, while in the incremental approach (iGLM), only the latest TR value was updated, which is an inevitable limitation of iGLM. We used Pearson correlation between the online and offline regressors to evaluate the quality of the online regressor. Mean correlation across participants was calculated by applying Fisher's z-transformation, averaging the transformed values, and then transforming back to the correlation coefficient. Figure 1 shows an overview of the simulation system. The WATCH module monitors a directory where a new fMRI volume file is created in real-time and loads it to send the data to the processing modules. In the simulation, we copied a pre-scanned image file volume-by-volume to the monitored directory. The RTP sequence included slice-timing correction with a temporal shifting of sampling time (TSHIFT), motion correction with volume registration (VOLREG), spatial smoothing by convolving Gaussian kernel within a mask (SMOOTH), and noise regression (REGRESS) with cGLM. Signal scaling is done within the REGRESS module. Details of each process are in the supplementary material (RTP simulation system implementation). The code for the simulation system is available on GitHub (https:// github.com/mamisaki/fMRI_RTP_Simulation).

Real-time processing simulation and offline analysis
The simulation was performed on a Linux computer (Ubuntu 16.04 LTS) with dual Intel Xeon CPUs (Gold 6126, 2.6 GHz, 12 cores for each), 256 GB RAM, and an NVIDIA TITAN V GPU (5120 CUDA cores with 12 GB memory). The simulation was started by copying an fMRI volume file in the directory monitored by the WATCH module. The resting-state fMRI data of 87 participants and the rtfMRI neurofeedback task data of 22 participants were used in the simulation. Their details are described in the supplementary material, 'fMRI samples for the simulation analyses. ' We tested ten types of RTP pipelines to evaluate the benefit of each processing step and regressor. Table 1 shows the process and regressors included in each pipeline. RTP0 included only the VOLREG for motion correction, a conventional minimum RTP, used as a baseline of the evaluation. RTP1 added slicetiming correction (TSHIFT) before VOLREG, and RTP2 added spatial smoothing (SMOOTH) to RTP1. From RTP3 to RTP9, noise models were regressed out. The regressors added at each RTP were HPF with Legendre polynomials in RTP3, motion parameters (Mot) of three shifts and three rotations in RTP4, temporal derivatives of the motion parameters (dMot) in RTP5, global signal (GS, mean signal within a brain mask) in RTP6, white matter/ventricle mean signals (WM, Vent) in RTP7, RETRO-ICOR (four cardiac and four respiration basis sets) in RTP8, and RVT (respiration volume per time with five time-shifted regressors) in RTP9. Details of each process and regressor are described in the supplementary material (RTP simulation system implementation). The regression was done with cGLM. Residual volumes of the regression were used as the processed signal [17]. We also ran offline processing using AFNI (details are shown in the supplementary material, 'Offline fMRI processing').
The computation time of each RTP was measured by the time from receiving a volume from the previous step to sending the processed volume to the next step in the simulation system. The start time of the WATCH module was the file creation time.  the voxel-wise signals of RTP0, RTP1, and RTP2. Dynamic FC was calculated for the detrended signal. For the RTP3 and the later pipelines, detrending was included in the HPF regression. Dynamic FC was calculated for all combinations of the 264 functional ROIs (6 mm-radius spheres) defined by Power, Cohen, Nelson, Wig, Barnes, Church, Vogel, Laumann, Miezin, Schlaggar and Petersen [28]. The ROIs in the MNI template space were warped into the participant's native image space using the Advanced Normalization Tools (ANTs, http://stnava.github.io/ ANTs/). We used two methods to calculate the dynamic FC; sliding-window correlation [29] and two-point algorithm [30]. The sliding-window correlation is a z-transformed Pearson correlation between ROIs within a time window [31]. A rectangular window of 5-TR width was used in the simulation. The window was moved at each TR to calculate the timecourse of dynamic FC. Although the usefulness of a narrow-window dynamic FC in neurofeedback has been demonstrated [32][33][34], a wide window-dynamic FC provides a more robust measure of connectivity [35] and has been used for neurofeedback [18]. Therefore, we also performed simulations with a 30-TR wide sliding window, as shown in the supplementary materials (p 17).

Noise amount in the real
The two-point algorithm assesses the agreement of change directions between ROIs. The feedback signal was given a binary value with positive reinforcement feedback (+1) when the two regions had the same change direction and no reinforcement feedback (0) when the change directions were different (if we wanted to train a participant to increase the connectivity). While the original introduction of the two-point method [30] used another ROI as a control to cancel a signal change unspecific to the target connectivity, we did not use a control ROI in the current simulation because another study [36] showed a control ROI is unnecessary when comprehensive noise reduction is applied.
For the neurofeedback task data, the mean signal in the left amygdala (LA) region, which was the target signal in the self-regulation task [27], was extracted for the noise evaluation. We applied the same procedures as for the voxel-wise signals of the resting-state data. The LA mask was defined in the MNI template brain with FreeSurfer (https://surfer.nmr.mgh.harvard.edu/) segmentation and warped into an individual brain space using ANTs.

Noise variance ratio
The motion and physiological noise amounts in the RTP signal were measured by a variance ratio explained by motion, cardiac, and respiration noise time-courses. Specifically, we performed a multiple regression analysis for the RTP voxel-wise and dynamic FC signals with noise regressors to calculate the coefficient of determination. The coefficient of determination with noise regressors, R n 2 , indicates the variance ratio explained by noises. We used this measure as the noise amount in the RTP signal. The first 65 TRs (62 TRs plus three TRs before the steady state) were excluded from the regression analysis because the output of real-time regression with a small number of samples is unreliable due to overfitting [10].
The effect of the motion, cardiac, and respiration noises were evaluated independently. Noise regressors were calculated in an offline process. The motion noise regressor was a time-course of framewise displacement (FD), calculated by the root sum of squared motion temporal derivatives. The cardiac noise regressors included heart rate (HR) [37] and the RETROICOR model [19]. The HR time-course was calculated using BioSPPy library in python (https:// biosppy.readthedocs.io/en/stable/) from a pulse oximetry measurement, resampled at each TR, and convolved with the cardiac response function presented in Chang, Cunningham and Glover [37]. The RETROICOR regressors were calculated offline using AFNI's RetroTS.py. The four basis regressors for the cardiac signal were extracted.
The respiration noise regressors included the respiration variance (RV) [37,38], the RETROICOR model, and RVT. The RV is a variance of respiration signal within a 3 s-window centered at each TR [37]. The RV time-course was convolved with the respiration response function shown in Chang, Cunningham and Glover [37] and provided by Power, Lynch, Dubin, Silver, Martin and Jones [38]. The RETRO-ICOR and RVT regressors were calculated offline using AFNI's RetroTS.py. The four basis regressors for the respiration signal and the five RVT regressors were used.
In the noise regression for dynamic FC, the windowed calculation could obscure the temporal association between the noise regressors and the dynamic FC. For example, in the sliding-window dynamic FC, any motions within the time window could affect the current neurofeedback signal. However, if we evaluated only the immediate association between the noise and RTP signal, we could miss a delayed effect of the motion within a window. Therefore, we added windowed regressors, the mean and standard deviation of the noise regressors within the connectivity calculation window, to the noise regressors. For example, the motion noise regressors for the dynamic FC were composed of the FD at each TR and the sliding-window mean and standard deviation of the FD within the FC calculation window (five TRs for the sliding-window and two TRs for the two-point method). If the node ROI signals were both affected by the same noise, their dynamic FC could increase when the standard deviation of the noise within the window was high. When the noise regressor is a variance measure (e.g. RV), dynamic FC could associate with the mean noise regressor value within the window. We included these possible delayed noise effects in the present noise analysis.
We also investigated the group-level association for participant-wise mean connectivity with the motion and physiological noises. Weiss, Zamoscik, Schmidt, Halli, Kirsch and Gerchen [18] indicated that subject-wise mean FC feedback signal could be correlated with breath rate. To investigate how RTP could eliminate this confounding effect, we calculated the mean dynamic FC across time and ROIs, and examined its association with each participant's mean framewise displacement, standard deviation of HR, and standard deviation of respiration rate. HR and respiration rate were calculated for each TR using the BioSPPy library, and then their standard deviations in a scan run were calculated for each participant.

Integrity of the neurofeedback signal
For the neurofeedback task data, we evaluated the benefit of RTP in improving the integrity of the neurofeedback signal. Integrity, in this case, refers to how well the signal reflects neural activation, free from noises. Although the signal amplitude relative to a rest block or the temporal contrast to noise ratio (TCNR) could be considered a measure of signal quality, we argue that these do not necessarily indicate the RTP benefit. For example, if the motion or physiological noise correlated with a self-regulation task, the task-related signal change could be larger for an RTP without noise corrections. TCNR could also be higher for an RTP with fewer noise corrections if a noise changed with the task time course. Such noise could reduce random signal fluctuation and increase the TCNR.
Since the size of task-related signal change cannot be a reliable measure of neurofeedback signal integrity, we used the signal correlation between RTP and offline processing to evaluate the RTP benefit. Assuming that the offline processed signal has the best available quality, the RTP signal with a higher correlation with the offline processed one should have better neurofeedback signal integrity. We calculated the Pearson correlation between the RTP and offline processed signals in the LA for each participant, applied Fisher's z-transformation and averaged across the participants for each RTP pipeline.

Statistical analysis
Participant-wise average R n 2 across the brain voxels or connectivity was compared between the RTP pipelines. Linear mixed-effect model analysis was applied to the average R n 2 with a fixed effect of the pipeline and a random effect of the participant on the intercept. We used the lme4 [39] with the lmerTest package [40] in R language and statistical computing [41]. The ad-hoc comparison between pipelines was performed in a sequential way to test the individual effect of the real-time process and regressor (e.g. RTP1-RTP0 for evaluating the TSHIFT effect, RTP2-RTP1 for evaluating the SMOOTH effect, RTP3-RTP2 for evaluating the HPF regressor effect). The comparison was performed with the lsmeans package [42] on R with multiple testing correction through the multivariate t-distribution method.
We also tested the noise effect in each voxel and connectivity. Since the BOLD signal, as well as the physiological signals, includes an autocorrelation component, a statistical test with a null hypothesis assuming the Gaussian noise could underestimate the false positive rate [43]. Thus, we used the phase randomization test [43,44] to evaluate the statistical significance of the R n 2 value. Phaserandomization was applied to the noise regressors. For the dynamic FC, the windowed average and standard deviation regressors were calculated after the phase-randomization. The randomized regressors were made 5 000 times for each participant.
All the noise calculation was done in the participant's native space because resampling in spatial normalization could affect the signal time series and change the noise contribution [23]. For the group analysis of the voxel-wise signal, the individual participant's R n 2 values were warped into the MNI template brain to calculate the mean R n 2 across participants in each voxel. R n 2 maps of phase-randomized regressors were also warped into the MNI to make a null distribution of the mean R n 2 for calculating p-value. The false discovery rate (FDR) correction with the Benjamini-Hochberg procedure was applied to the p-values. The maps were thresholded with FDR < 0.05.
For the neurofeedback task data, we used a permutation test to evaluate the statistical significance of the R n 2 value for each pipeline and comparisons between the pipelines, and FDR correction was applied to p-values. We used the same linear mixedeffect model as for the resting-state data to test the correlation between RTP and offline processed signals. Figure 2 shows the HPF performance of Legendre polynomial and DCT regressors for the white noise signal (figure 2(A)) and the resting-state fMRI signal (figure 2(B)). Figure 2(A) indicates that frequencies higher than the designed threshold (vertical dotted line) decreased more with iGLM than with cGLM, especially when the number of TRs was short, demonstrating that iGLM had a less accurate HPF property. Figure 2(B) shows a similar tendency for the BOLD signal, while the difference between the GLM approaches was less prominent than for white noise.  Figure 3 shows the correlation between the online and offline calculations of the RETROICOR and RVT regressors for the incremental (iGLM) and cumulative (cGLM) approaches. No difference between iGLM and cGLM was observed for the cardiac regressors (Card). For the respiration regressors (Resp), cGLM had a higher correlation with the offline calculation. For the RVT regressors, iGLM had considerably lower correlations with the offline calculation. This indicates that online-made physiological noise regressors were inaccurate in iGLM without retrospective correction, especially for the respiration and RVT regressors. These results suggest that cGLM is preferred to iGLM in RTP, as far as computation time permits. Figure 4 shows the average computation time across participants for each processing module at each TR. The processing was applied to whole-brain voxels (149 082 voxels on average). The REGRESS process used the cumulative approach and included all the regressors implemented in the system. The REGRESS process started at the 42nd TR to allow for the acquisition of enough samples for the regression (the initial three volumes were excluded from the process to ensure the fMRI signal reached the steady state; thus, 39 samples were available at the 42nd TR).

Computation time
The REGRESS's initial processing TR took a long time due to the initialization process (see supplementary material, 'RTP simulation system implementation'). Other than the initial TR of REGRESS, the most time-consuming processing was the VOLREG, volume registration for motion correction. We found resampling the volume image in the aligned grid took a long time in VOLREG. The time for REGRESS with CPU increased with TR, which was because cumulative calculation used more data in a later TR. The slope of the increase was less steep with GPU computation. Other processes took less than 100 ms each, and the total process time was less than 400 ms with GPU and less than 500 ms with CPU, which would be short enough for a majority of rtfMRI applications. Figure 5(A) shows distributions of the mean noise variance ratio (R n 2 ) explained by the motion and physiological noises for the voxel-wise signal. The plot shows the distribution of the participant mean values in the brain. The voxels in the white matter and the ventricle areas were excluded from the mean calculation (the masks for the white matter/ventricles signal regression in RTP were used). Significant decreases or increases of the mean noise variance ratio by adding an RTP component were summarized in table 2 (statistical values for each comparison are shown in supplementary table S1).

RTP effects on the noise variance ratio (R n 2 ) in the voxel-wise signal
The cardiac and respiration noise ratios increased with spatial smoothing (SMOOTH, △ in table 2). Voxel-wise analysis results are shown in figures 6-8. Figure 6 shows maps of the voxel-wise significant motion noise R n 2 . Adding the motion parameter regression removed most of the significant voxels (RTP4), and adding the motion derivative regression (RTP5) eliminated a significant motion noise ratio in all voxels. Supplementary figure S1 shows significant differences in the motion noise R n 2 at each sequential contrast. TSHIFT and SMOOTH increased the motion noise ratio in many voxels, especially around the white matter and the ventricle areas.  Figure 7 shows maps of the voxel-wise significant cardiac noise R n 2 . Without RETROICOR, a significant cardiac noise ratio was seen in the whole-brain area and was especially high in the medial to anterior temporal and the anterior cingulate areas, overlapped with the large cerebral blood vessels. Adding the RETROICOR regression (RTP8) removed most of the significant voxels, although a significant effect remained in the large blood vessel areas. Supplementary figure S2 shows significant differences in the cardiac noise R n 2 (FDR < 0.05) at each sequential contrast. Spatial smoothing (SMOOTH), HPF regression (REG[HPF]), and global signal regression (REG[GS]) significantly increased the cardiac noise ratio in many voxels. In contrast, significant decreases were seen by adding the regressors of motion derivatives (REG[dMot]), white matter/ventricle signals (REG[WM, Vent]), and the RETROICOR (REG[RICOR]). The RETROICOR regression had the most substantial effect among them. Figure 8 shows maps of the voxel-wise significant respiration noise R n 2 . Without RETROICOR, a significant respiration noise ratio was seen in the whole-brain region, and most of the significant voxels disappeared by adding RETROICOR. Supplementary figure S3 shows significant differences of the respiration noise R n 2 at each sequential contrast. Spatial smoothing (SMOOTH) and HPF regres- ----, no significant effect; △, significant increase of the mean noise variance ratio; O, significant decrease of the mean noise variance ratio. * * , p < 0.01; * * * , p < 0.001. Detailed statistical values are shown in supplementary table S1.  Connectivity-wise analysis results are shown in figures 9 and 10. Figure 9 shows maps of the connectivity-wise significant motion noise R n 2 . A significant motion noise was seen across whole-brain connectivity before adding the motion derivative regressor (REG[dMot]). After adding REG[dMot] (RTP5), no significant motion noise ratio was seen in any connectivity. No connectivity-wise significant difference of the motion noise R n 2 was seen in any connectivity at any sequential contrast. Figure 10 shows maps of the connectivity-wise significant cardiac noise R n 2 . A significant cardiac noise was seen in the connectivity across the anterior to medial temporal and the anterior cingulate areas. Connectivity in the inferior occipital area also showed a significant cardiac noise ratio. The significant cardiac noise ratio gradually decreased with adding the regressors of motion parameters (REG  no connectivity showed a significant cardiac noise ratio. No connectivity-wise significant difference of the cardiac noise R n 2 was seen in any connectivity at any sequential contrast. For the respiration noise, the connectivity-wise analysis showed no significant R n 2 and R n 2 differences between pipelines in any connectivity.

RTP effects on the noise variance ratio (R n 2 ) in the sliding-window connectivity
The results for the 30-TR sliding-window connectivity are shown in the supplementary material (p. 17). Connectivity-wise analysis results are shown in figures 11 and 12. For the motion noise, the connectivity-wise analysis showed no significant R n 2 and R n 2 differences between pipelines in any connectivity. Figure 11 shows maps of the connectivity-wise significant cardiac noise R n 2 .   Without RETROICOR, a significant cardiac noise ratio was seen in the whole-brain connectivity, and the effect was large in the connectivity between the ROIs with a significant cardiac noise ratio in the voxel-wise analysis ( figure 7). The significant cardiac noise effect disappeared after adding the RETRO-ICOR regressor (REG[RICOR] at RTP8). Supplementary figure S4 shows connectivity-wise significant     differences in the cardiac noise R n 2 at each sequential contrast. A significant increase was seen when adding the spatial smoothing (SMOOTH) in the areas affected by the cardiac noise in the voxel-wise analysis (figure 7). A significant decrease was seen when adding the RETROICOR (REG[RICOR] at RTP8) in the whole brain with a substantial effect in the areas affected by the cardiac noise in the voxel-wise analysis. Figure 12 shows maps of the connectivity-wise significant respiration noise R n 2 . A significant respiration noise ratio was seen only for the RTP0 (only the motion correction) in the whole-brain area. No connectivity-wise significant difference of the respiration noise R n 2 was seen in any connectivity at any sequential contrast.  indicate that the mean connectivity was biased to a positive value for many participants and the size of the bias was correlated with the participant's amplitude of motion, HR, and respiration measures for RTP0. Adding the global signal regression (RTP6) eliminated this bias. However, the correlations (Spearman's rank-order correlation) were still statistically significant even after the global signal regression, although the effect size (slope of the regression) became small. Figure 13(B) shows the same association for the two-point connectivity. Figures for all pipelines are shown in supplementary figures S8-S10. Like the sliding-window connectivity, the mean two-point connectivity was correlated with the participant's amplitude of motion, HR, and respiration measures for RTP0. With the global signal regression (RTP6), the mean connectivity became nearly equal to 0.5 for all participants, while the associations remained statistically significant with a small effect size. Figure 14(A) shows distributions of the noise variance ratio (R n 2 ) explained by the motion and physiological noises for the LA signal in the neurofeedback task run. The RTP labels with boldface indicate R n 2 was significantly high with permutation test. The result was similar to the voxel-wise analysis of the resting-state data in the amygdala region. No significant motion R n 2 was observed. The cardiac R n 2 was significant in all pipelines and was significantly decreased with RETROICOR regression (RTP8). The respiration R n 2 was significant until RTP3, and motion regression (RTP4) significantly reduced it. Figure 14(B) shows the mean and 95% confidence interval of the signal correlation (z-transformed) between the RTP and offline-processed signal. The RTP labels with boldface indicate that the correlation was significantly higher than RTP0. For the neurofeedback task run, the LA signal with RTP8 had the highest mean correlation with the offline-processed signal.

Discussion
The simulation results indicated that (a) cGLM compared to iGLM could substantially improve the accuracy of online-formed regressors with a retrospective correction, (b) extensive RTP with the cGLM could be completed in a short enough time for RTP, and (c) significant noise reduction was seen for many tested RTP steps, though its benefit differs among the voxelwise and online dynamic FC signals.

Cumulative GLM is preferred in RTP
The HPF regression with iGLM filtered higher frequency than the designed threshold at early TRs ( figure 2). This deficit of filtering property is because a piece of HPF regressors made for a full-length scan could fit higher frequency components than the designed threshold. cGLM avoids this problem by updating the regressors to adjust to the designed pass frequency at each length of data. However, the difference between the GLM approaches was not large for the BOLD signal compared to the white noise signal. This is because most of the BOLD signal power is in a low-frequency range, resulting in relatively small differences in the highfrequency range. This result suggests that the iGLM drawback of high-pass-filtering may not be a substantial problem in a practical rtfMRI application. Also, the drawback of HPF regression in iGLM may be avoided by using another online frequency filtering method, such as the FIR filter [13]. However, we should remember that we must apply the same filter to the regressors if we use noise regression other than the HPF; otherwise, an artifact in the residual signal could result [45].
In an online creation of physiological noise regressors, respiration and RVT noise models were not accurate with the incremental approach (figure 3). To detect a low-frequency fluctuation and its phase we need a certain length of observation (at least one cycle of fluctuation). Since the cycle of respiration is usually longer than a TR, the respiration phase estimation can be made only retrospectively. Also, the RVT regressor depends on the peak detection in the respiration signal time-series, but the peak can be identified only retrospectively. Thus, the phase estimation and peak detection cannot be accurate in an online process. As the defect of an inaccurate noise regressor was accumulated with iGLM, the correlation between the online and offline regressors continued decreasing in iGLM. In contrast, cGLM could compensate for this problem by updating the regressors retrospectively at every TR.
The iGLM has been popular in rtfMRI applications thanks to its low computer memory consumption and short computation time. However, the present simulation indicates that the cumulative approach has a notable advantage when an online refinement of the regressor is required. The simulation also demonstrated that the cost of computation time for cGLM is not so high. Consequently, we should use the cGLM in RTP, especially when using online-made physiological noise regressors.

Computation time does not limit the extensive RTP
The computation time was short enough for RTP in our RTP simulation system. The total processing time was less than 400 ms with GPU and less than 500 ms with CPU only (figure 4). The most time-consuming part of the present system was motion correction (VOLREG). We found that resampling a data volume in a registered space took long in this process. While we used a compiled C library of AFNI implementation using CPU only, a GPU implementation, like Scheinost, Hampson, Qiu, Bhawnani, Constable and Papademetris [9], might further shorten the computation time of this process.
The REGRESS took more time in later TRs than in earlier TRs since the number of samples included in the cGLM increased with time. The slope of the processing time increase was less steep with GPU than CPU only, indicating that a cGLM computation with a more extended scan is possible with GPU, as far as the memory space allows. In the current simulation, we used a relatively large matrix size (128 × 128 × 34) as fMRI and the number of voxels was 149 082 (about 600 kB with single-precision float) per volume on average after applying the brain mask. This demonstrates that the data size would not limit the cumulative computation with a recent GPU's large gigabyte memory. Even if the computational capacity is limited compared to the present simulation, we can reduce the burden by limiting the processing regions (e.g. in the gray matter or in a target region). Thus, the computation time should not limit the application of comprehensive real-time fMRI processing.
The simulation result shows that we can almost ignore the cost of computation time in RTP with current computer hardware. However, we still need to consider another cost of RTP, the number of regressors. The more regressors we use, the longer burn-in time we need before utilizing the regression in order to wait for a sufficient number of samples. Hinds, Ghosh, Thompson, Yoo, Whitfield-Gabrieli, Triantafyllou and Gabrieli [17] showed that 25 TRs were required to make a reliable regression result. Our previous investigation [10] indicated that the required TRs depended on the number of regressors and the output of real-time regression with many regressors was unreliable in the early TRs because of overfitting. The long burn-in time could cost scan time and limit experimental design in some applications of rtfMRI. Hence, we should refrain from including a regressor without a significant effect on reducing noise.

RTP noise reduction for the voxel-wise signal
Spatial smoothing (SMOOTH) increased the mean noise variance ratio for the cardiac and respiration noises, and HPF (REG[HPF]) increased the mean noise variance ratio for the cardiac noise. The increased noise ratio with smoothing was partly due to the spread of the noise effect in the voxels neighboring the noise-contaminated areas (e.g. supplementary figures S2 and S3). The increased noise variance ratio is also attributable to the decreased total variance. Since the smoothing did not specifically reduce the physiological noises, the decrease of the total variance by smoothing could increase the relative ratio of the physiological noise variance. The same logic could be applied to HPF. As the BOLD signal is dominated by low-frequency components ( figure 2(B)), HPF could reduce a large amount of total variance. If the cardiac noise was not included in the removed low-frequency range, the ratio of cardiac noise variance to the total variance could increase.
We consider that these increases in noise variance ratio are not problematic. These are the results of removing different kinds of noise variance, such as high-spatial-frequency noise and low-temporalfrequency signal fluctuation. Indeed, HPF regression significantly reduced the mean motion and respiration noise ratios, indicating that this regression is beneficial in removing a low-frequency noise. Also, the slice-timing correction has indicated a significant benefit in a reliable estimation of brain activity [46], and a noticeable effect of spatial smoothing on activation detection has been demonstrated [47]. Regardless of the motion and physiological noise variance ratios, we consider that TSHIFT, VOLREG, SMOOTH, and REG[HPF] should always be included in RTP because these remove known noise components other than the motion and physiological noises.
A significant reduction of the mean noise variance ratio in any of the motion and physiological noises for the voxel-wise signal was seen with regressions of HPF (REG[HPF]), motion parameters (REG[Mot]), white matter/ventricle signals (REG[WM, Vent]), and RETROICOR (REG[RICOR]). In the voxel-wise evaluation, the motion derivatives also had a significant effect on reducing the noise in many voxels (supplementary figures S1-S3) while not significant in the whole-brain average.
The global signal regression significantly increased the cardiac noise ratio (supplementary figure S2) and the RVT regression increased the respiration noise ratio (supplementary figure S3) in many voxels in the voxel-wise evaluation. The increase with global signal regression might be attributable to a reduction of the total signal variance. The increase with RVT regression could be due to the inaccuracy of the online-made RVT regressor. Although the retrospective update with cGLM was better than iGLM, the online-made regressor still showed a reduced correlation with the offline one ( figure 3). The error in the online-made RVT regressor might have synchronized with the respiration signal, which could recall the respiration-correlated signal fluctuation in the residual of the RTP regression. Similar recall of noise has been demonstrated by Hallquist, Hwang and Luna [45] for the HPF regression; when a HPF signal was regressed with a regressor without HPF, the residual could include the filtered lowfrequency component. This harmful effect indicates that RVT regression should not be used in RTP.
Adding the white matter/ventricle signal regression (REG[WM, Vent]) after the global signal regression (REG[GS]) showed a significant reduction of the cardiac noise (RTP7-RT6 in table S1 and figure S2 in supplementary materials), indicating that the mean white matter/ventricle signals were distinct from the global signal. This result is consistent with the evidence that gray matter had the largest effect on the global signal and the effect of the white matter and cerebrospinal fluid was small except for the partialvolume effects [48,49]. Also, a significant cardiac noise reduction was seen when adding RETROICOR (REG[RICOR]) after adding the global signal and white matter/ventricle signal regressors (RTP8-RTP7 in table S1 and figure S2 in supplementary materials), indicating that the global signal and the tissuewise average signals cannot substitute for the RETRO-ICOR regressor. Since the cardiac noise effect was concentrated on the voxels around the large blood vessels, the global and tissue-wise averages could not match such a spatially localized noise. This indicates that each of the global signal regression, white matter/ventricle signal regression, and RETROICOR regression have a distinctive contribution to noise reduction, and they cannot substitute for each other.
A potential risk of physiological noise effect on fMRI results has been indicated in the regions prone to cardiac noise [50]. The regions with significant cardiac noise (figure 7) also overlap with the salience network area [51] and are a clinically meaningful target for intervening emotion-regulation function [52,53]. RETROICOR regression should be vital for getting a robust neurofeedback signal from these areas. Global signal regression and mean white matter/ventricle signal regression cannot substitute for it. We also note that the cardiac noise ratio remained significant even after RETROICOR in the voxels around the large cerebral blood vessels (figure 7, RTP8). Although the remaining noise ratio was small (<0.1), we should be cautious of the cardiac noise when using the signal in these areas for neurofeedback.
In summary, for the online voxel-wise signal evaluation, regression with motion parameters (REG[Mot]), mean white matter/ventricle signals (REG[WM, Vent]), and RETROICOR (REG[RICOR]) significantly reduced the motion or physiological noises on average across the brain, while the effect of RTP was different across voxels. Motion derivative regressor (REG[dMot]) also reduced the noise at several voxels. While TSHIFT, VOLREG, SMOOTH, and REG[HPF] had no significant effect or increased the noise ratios due to a reduced total variance, they should be applied in RTP as they could reduce other kinds of noise. No significant contribution was found with the global signal regression (REG[GS]), and the online-made RVT regressor introduced a respiration noise artifact; thus, these regressors should not be used in RTP for the voxelwise signals.

RTP noise reduction for the sliding-window connectivity
A significant reduction of the mean noise variance ratio in any of the motion and physiological noises was seen with regressions of motion derivatives (REG[dMot]), global signal (REG[GS]), and RETROICOR (REG[RICOR]) for the slidingwindow dynamic FC. Global signal regression had a significant effect on reducing the motion noise variance ratio. This is consistent with the report by Parkes, Fulcher, Yucel and Fornito [54], showing the benefit of global signal regression in reducing motion noise for an offline resting-state fMRI analysis. The significant reduction of the cardiac noise with RETROICOR indicates that the global signal and the tissue-wise average signals cannot substitute for the RETRO-ICOR regressor, which was the same as in the voxelwise signal.
With the connectivity-wise evaluation, significant motion noise was seen in the whole-brain area, that was disappeared after adding the motion derivative regression (RTP5 in figure 9). A significant cardiac noise ratio was localized in the deep brain regions, including the inferior occipital, medial to anterior temporal, and the anterior cingulate areas (figure 10), which overlap with large cerebral blood vessels. The connectivity in these areas was often implicated in emotion-regulation functions [52,53]. The present result demonstrates the importance of extensive realtime denoising when we target these connectivities in a neurofeedback application.
Regarding the group-level association (figure 13(A) and supplementary figures S5-S7), only the global signal regression could eliminate the noiseassociated bias of the mean connectivity, consistent with the observation by Weiss, Zamoscik, Schmidt, Halli, Kirsch and Gerchen [18]. We will discuss the effect of global signal regression in RTP in a later independent section.
In summary, for the online sliding-window dynamic FC evaluation, regression with motion derivative (REG[dMot]), global signal (REG[GS]), and RETROICOR (REG[RICOR]) significantly reduced the motion or physiological noises on average across the brain, while the effect was different across regions. REG[GS] also removed the bias of average connectivity across time and brain regions. While TSHIFT, VOLREG, SMOOTH, and REG[HPF] had no significant effect, we think that they should be applied in RTP as they could reduce other kinds of noise, as we discussed above. Although no significant contribution was found with the motion parameter regression (REG[Mot]), when we removed this regressor from the pipeline, we observed a significant increase in the mean motion noise variance ratio compared to the best-performed pipeline (RTP8; t = 3.114, p = 0.002). This suggests that global signal regression cannot substitute for the motion parameter regression and that these regressors had distinctive contributions to motion noise reduction. Mean white matter/ventricle signals (REG[WM, Vent]) and RVT regressions had no significant contribution to the noise reduction.

RTP noise reduction for the two-point connectivity
Overall, the mean noise variance ratio for the two-point connectivity was less than that for the 5-TR sliding-window connectivity (the vertical axis scale in figure 5(C) is smaller than that in figure 5(B)), which could be due to the narrow window width preventing from spreading the noise effect. A significant reduction of the mean noise variance ratio in any of the motion and physiological noises was seen with TSHIFT and regressions of motion parameters (REG[Mot]), motion derivatives (REG[dMot]), white matter/ventricle signal (REG[WM, Vent]), and RETROICOR (REG[RICOR]). With the connectivity-wise evaluation, no significant effect of motion was seen in any connectivity. This is consistent with our previous investigation, demonstrating that the two-point method was less prone to motion than the slidingwindow methods [36].
The two-point connectivity showed more statistically significant cardiac noise than in the slidingwindow connectivity. The connectivity-wise evaluation showed a significant cardiac noise variance ratio across the whole brain ( figure 11). Because the two-point connectivity evaluates the change direction between the two consecutive time points, it could be more sensitive to a high-temporal-frequency noise than the sliding-window method. The cardiac noise effect was still significant after REG [GS] and REG[WM, Vent], but disappeared after adding the RETROICOR regression ( figure 11, RTP8). This demonstrates that the RETROICOR benefit cannot be replaced by global signal and tissue-wise average signals. Significant respiration noise was removed by adding slice-timing correction (TSHIFT, figure 12), suggesting that the two-point connectivity is sensitive to a small temporal shift by the slice-timing difference. The slice-timing difference has been known to make a deviation from true connectivity [55].
Regarding the group-level association ( figure 13(B), supplementary figures S8, S9, S10), only the global signal regression could eliminate the noise-associated bias of the mean connectivity. This was the same as for the sliding-window connectivity and consistent with the observation by Weiss et al [18].
In summary, for the online two-point dynamic FC evaluation, slice-timing correction (TSHIFT) and regression with motion parameters (REG[Mot]), motion derivatives (REG[dMot]), white matter/ventricle signal (REG[WM, Vent]), and RETROICOR (REG[RICOR]) significantly reduced the physiological noises on average across the brain, while the effect was different across regions. Global signal regression (REG[GS]) removed the bias of average connectivity across time and brain regions. While VOLREG, SMOOTH, and REG[HPF] had no significant effect or increased the noise ratios due to a reduced total variance, we think that they should be applied in RTP as they could reduce other kinds of noise. Regression with mean white matter/ventricle signals (REG[WM, Vent]) and RVT had no significant contribution to the noise reduction.

Global signal regression in RTP
The global signal regression had a significant effect on reducing the motion noise in the 5-TR slidingwindow dynamic FC (table 3). It also had a notable effect on reducing the average connectivity bias across time and brain regions in the grouplevel analysis for both the sliding-window and the two-point connectivity ( figure 13 and supplementary figures S5-S10). This result is consistent with the reports for offline resting-state FC analysis [22] and the neurofeedback study with dynamic FC [18]. The neurofeedback study showed that global signal regression could remove the correlation between the participant-wise average dynamic FC and respiration measures. They also indicated that regression with the physiological noise models, including RETROICOR, could not remove this bias.
However, the present results showed that the RETROICOR regression had an additional benefit in reducing the cardiac noise after the global signal regression in the individual participant's signal time-course. Also, even with the global signal regression, excluding the motion parameters significantly increased the motion noise variance ratio for the sliding-window connectivity (section 4.4), indicating that the distinctive effect of motion parameter regression from the global signal regression. Thus, the global signal regression is not the all-around solution for noise reduction in FC, and the motion parameter and RETROICOR regressors significantly contributed to noise reduction even when accompanied by the global signal regression.
The discrepancy between the group-level association and the individual noise regression analysis suggests that the global signal regression removed a signal component distinct from the motion and RETROICOR regressors. The noise component inducing the average connectivity bias and being removed by the global signal regression might be a BOLD signal fluctuation having a loose temporal association with the noise time-course. Indeed, Power, Schlaggar and Petersen [22] reported that the increase of the FC by motion could be seen even in 8-10 s delays, which might be due to a spin-history noise [21].
The global signal regression is the most controversial processing used in resting-state fMRI. Many studies have investigated the source of the global signal and its association with noise and functional neural activations [49,[56][57][58][59][60]. The side effects of the global signal regression, such as an artifactual negative correlation, have also been demonstrated [61][62][63]. Although we refrain from further considering the source of the global signal, which is out of the scope of the present study, we consider a practical suggestion in using the global signal regression for a rtfMRI application.
In a rtfMRI application, a signal fluctuation with a tight temporal association with a motion or physiological noise should be more apparent for a participant, and such a noise could increase the risk of inducing spurious training with noise regulation. Thus, regardless of using global signal regression, we should use the motion and RETROICOR regressions in RTP. When we use a voxel-wise or an ROI-based neurofeedback signal, the present result indicates that the global signal regression is not needed in RTP. The global average signal was not useful in removing a voxel-wise noise.
When we use dynamic FC as a neurofeedback signal, we should include the global signal regression as it has a notable benefit in reducing the bias in the average connectivity across time. However, we should also be careful about a possible artifact due to the global signal regression. Once we determine the target connectivity, we should perform an RTP simulation analysis, like Ramot and Gonzalez-Castillo [64] and Misaki, Tsuchiyagaito, Al Zoubi, Paulus, Bodurka and Tulsa [36], to confirm that the global signal regression does not make a spurious shift of the connectivity neurofeedback signal. If a notable artifact is seen, we might have to refrain from using the global signal regression in RTP and admit the reduced signalto-noise ratio in the neurofeedback signal or consider another target.
The group-level association for the mean connectivity with the motion and physiological noises remained statistically significant even after the global signal regression (figures 8 and 9). However, we consider that this statistical significance would not matter in the practice of neurofeedback training because the effect size of the noise (slope of the regression line) was minimal after the global signal regression. In the neurofeedback context, we can ignore the noise effect if its effect size was much smaller than that by self-regulation.

RTP noise reduction for the neurofeedback target signal
To confirm that the simulation results for the restingstate data could be applied to a neurofeedback task run, we performed the same simulation analysis for the neurofeedback task data [26]. The result was parallel to the resting-state data's voxel-wise analysis. In this neurofeedback experiment, participants were trained to increase the LA activation by recalling a positive autobiographical memory. Hence, the LA signal should include a self-regulated change so that the noise variance ratios could be less dominant than in a resting state. However, we still observed significant noise variance ratios. While the motion noise was not significant in any pipelines, the respiration noise was significant before adding the motion parameter regression, and the cardiac noise was significant in all RTPs with a significant reduction by RETRO-ICOR regression ( figure 14(A)). This was consistent with the resting-state data analysis result in figure 7, showing that the cardiac noise ratio was significantly reduced with RETROICOR, while it was still significant in the amygdala region. The signal integrity analysis ( figure 14(B)) also demonstrated that the signal correlation between RTP and the offline process improved with the inclusion of additional noise regressors.
We should note that the noise amount and the RTP effect in a neurofeedback signal could vary depending on the target region, the applied regulation strategy, and the individual's self-regulation performance. Thus, we used resting-state data to evaluate the general effect of noise and RTPs. Nevertheless, the results for the neurofeedback task data demonstrated that the noise effect could be significant even with a task-related signal change and RTP noise reduction is beneficial in the neurofeedback application.

Limitations
Several limitations of the current study should be acknowledged. Since we used similar regressors in RTP and the noise regression analysis, it may look trivial that including these regressors in RTP could reduce the noises. However, we should point out that the RTP regressors and offline-made noise regressors used in the noise regression analysis cannot be identical. Online RTP can use limited samples of the current and past time points. Indeed, we found the deficit of the RVT regression in RTP. Furthermore, the physiological noise variance extracted in the noise regression analysis was not specific to the RETRO-ICOR model. The RETROICOR regressor has a general form of the Fourier basis set that can fit any temporal fluctuation synchronized with the cardiac and respiration fluctuations. We also included the HR and RV convolved with their hemodynamic response functions in the noise regression analysis. Thus, the benefit of RETROICOR regressors in RTP is never a trivial result.
The present result is limited to the noise measures used for the evaluation. We used the motion (FD), HR, RV, RETROICOR, and RVT in the noise evaluation. Although these are major noise sources for the BOLD signal, there are yet other sources that could have a significant effect on the signal change. For example, Power, Lynch, Dubin, Silver, Martin and Jones [38] indicated that deep breath had a large effect on the fMRI signal changes and that that cannot be captured by the respiration envelope, RV, and RVT measures. In a rtfMRI application, complete noise-cleaning may be too costly since using more RTP regressors increases the initial burn-in time further. However, at the offline analysis in evaluating the neurofeedback effect, we should examine the noise effect that has not been removed at RTP. The apparent training effect could result from a change in motion or physiological noise pattern [18].
Using resting-state data in the simulation might also limit the implication of the present results. When there was a task-induced signal change, the noise variance ratio may be less significant than in the present results. However, a substantial signal change is not always guaranteed, especially in the early phase of the training when self-regulation ability is not high. Even if a participant could regulate the brain activation, a task-induced fMRI signal change is often smaller than the noise variance; TCNR is usually lower than 1.0 [65,66]. For detecting such a low contrast signal change, noise reduction should be applied as much as possible. Indeed, the simulation for the neurofeedback task data demonstrated that the noise variance ratio was still significant with a selfregulation task. Furthermore, if a noise-contaminated neurofeedback signal was used, the training could depend on the regulation of noise-associated motion or physiological states. Thus, the caveat of the present results could be applied to general neurofeedback training.
While we presented beneficial RTP components for the voxel-wise and dynamic FC signals on average across the brain regions, the optimal combination with minimum regressors will depend on the target regions. For example, the cardiac noise effect concentrated on the areas overlapping with large cerebral blood vessels (figure 7), which was removed by RETROICOR regression. If the target is out of these regions, RETROICOR might not be needed. Although the present results help consider which kinds of noise could affect which brain areas and what RTP helps reduce it, the optimal RTP pipeline needs to be designed for each RTP application with a simulation analysis [36,64].
Regarding the sliding-window dynamic FC, the RTP effect could be different depending on window width. The previous study [32] indicated that the correlation between motion and dynamic FC increased with the window width so that a narrow window was more robust to noise. Also, a narrow window could be preferred in a neurofeedback application to provide a timely feedback signal; a wide-window dynamic FC reflects less current brain activation and depends more on the long history of past brain states. We observed this limitation in the results of the 30-TR sliding-window connectivity (supplementary material, p 17). However, several studies with an offline analysis have shown that a wide window is required to evaluate dynamic FC robustly [64], and that even with a wide window (e.g. 2 min), dynamic FC is challenging to detect with a single session measurement [65]. Thus, in choosing the window width for a dynamic FC neurofeedback, we need to consider the tradeoff between the robust evaluation of brain state and timely and less noisy feedback.
While the optimal choice could be different between individual applications, we consider that a narrow window can be approved in a neurofeedback application. We should also note that a TR-wise estimate of regional brain activation is not robust, similar to a narrow-window dynamic FC, and many studies have proven that even such a fragile estimate helped participants learn brain self-regulation. Thus, even if a real-time estimate cannot meet the requirement of the offline statistical analysis, it still could be useful in training participants to self-regulate their brain state. Indeed, studies have demonstrated the utility of narrow window connectivity for neurofeedback training [32][33][34], and the two-point dynamic FC neurofeedback, which is one extreme window width focusing on feedback timeliness, also demonstrated successful training in self-regulating FC [30,66]. As both a narrow or wide window can be approved for dynamic FC neurofeedback, the choice should depend on the task and target in each individual application. Refer also to the supplementary material (p. 17) for the caveats using a wide-window dynamic FC for neurofeedback.
While using a control ROI or control connectivity was not considered in the present study, it could be an economical noise reduction regressor. However, we should remember that the control selection should be specific to the target region and is not a trivial task. The present simulation demonstrated that the global signal and the white matter/ventricle average signals could not be a substitute for the RETROICOR, suggesting that just taking the average of the regions apart from the target would not be enough to suppress the noise. Also, if the selected control area had a unique functional response, subtracting the control signal from the target could induce an artifactual signal fluctuation. The self-regulation task usually invites brain activations in many regions, not limited to the target [67][68][69][70][71]. Thus, we must be very cautious in choosing an appropriate control signal, which needs an extensive independent evaluation to optimize noise reduction performance.

Conclusions
The present study demonstrated that extensive realtime noise reduction had a minimal cost of computation time with significant benefit in reducing motion and physiological noises for rtfMRI. We believe our comprehensive evaluation will promote a widespread utilization of extensive real-time noise reduction for neurofeedback applications. The high-quality realtime estimate of brain activity would pave the way for robust and replicable neurofeedback training, as well as the urgently needed development of novel and robust interventions for several mental and neurological disorders, as well as a potential application in presurgical mapping.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.