Using carpet plots to analyze transit times of low frequency oscillations in resting state fMRI

A “carpet plot” is a 2-dimensional plot (time vs. voxel) of scaled fMRI voxel intensity values. Low frequency oscillations (LFOs) can be successfully identified from BOLD fMRI and used to study characteristics of neuronal and physiological activity. Here, we evaluate the use of carpet plots paired with a developed slope-detection algorithm as a means to study LFOs in resting state fMRI (rs-fMRI) data with the help of dynamic susceptibility contrast (DSC) MRI data. Carpet plots were constructed by ordering voxels according to signal delay time for each voxel. The slope-detection algorithm was used to identify and calculate propagation times, or “transit times”, of tilted vertical edges across which a sudden signal change was observed. We aim to show that this metric has applications in understanding LFOs in fMRI data, possibly reflecting changes in blood flow speed during the scan, and for evaluating alternative blood-tracking contrast agents such as inhaled CO2. We demonstrate that the propagations of LFOs can be visualized and automatically identified in a carpet plot as tilted lines of sudden intensity change. Resting state carpet plots produce edges with transit times similar to those of DSC carpet plots. Additionally, resting state carpet plots indicate that edge transit times vary at different time points during the scan.


Full SHAG carpet plots for DSC-MRI and rs-fMRI
SHAG carpet plots, when first reordered, tend to result in voxels grouped in the upper or lower portions of the plots which do not contain signal patterns consistent with that seen in the majority of voxels (See Fig. S1 for full SHAG carpet plots). These grouped voxels were removed from each carpet plot in order to perform a linear fit to edges present in the plots (See Fig. S2 for percentage of voxels removed). The distribution of discarded voxels in GM, WM, and CSF are shown in Figure S3. Figure S1: The full SHAG carpet plots for DSC and rs-fMRI for eight subjects. The DSC-MRI SHAG carpet plots are shown in the first and third rows while the corresponding rs-fMRI SHAG carpet plot are shown in the second and fourth rows. Figure S2: The percentage of cropped portion of carpet plots for the analysis throughout the paper. 99.2 ± 0.7 % of all voxels for DSC (blue) and 71.3 ± 5.2 % for rs-fMRI (red) are used. Figure S3: A histogram for each subject presenting the distribution of discarded voxels in GM, WM, CSF, and the selected voxels in total. The means and the standard deviations of the percentage of discarded voxels in GM, WM, and CSF were 17%±6%; 16%±5%; 8%±3%, respectively.

The effects of spatial smoothing and motion artifact on the transit time and delay time calculation
Given the fact that spatial smoothing will mix signals in voxels with their neighborhood voxels, potentially mixing sLFOs with other signals (e.g., Mayer waves), we compared the rs-fMRI SHAG carpet plots and their corresponding transit times with and without spatial smoothing being applied. Signal-to-noise ratio was improved and fewer voxels were discarded when the spatial smoothing was applied (See Fig. S4). There were no significant differences on the transit time in cropped carpet plots in full, GM, WM, and CSF, respectively (ANOVA test). Dynamic change patterns within the subject also remain very similar (See Fig. S5). For these reasons, it was determined that spatial smoothing is still necessary in the preprocessing step.
Even after motion correction, there is an obvious artifact observed in one subject's rs-fMRI SHAG carpet plot (See Fig. 2). A simulation was done to evaluate the effect of a single motion artifact on the delay time calculation (See Fig. S6). In the simulation, the original BOLD was shifted by 3 seconds to form a new BOLD signal (shifted BOLD). The head motion was added to both signals at the same time point (See Fig. S6(A)). The size of the signal spike due to the head motion was simulated as 0 to 5 times (with 0.1 increment) the standard deviation of the given BOLD signal. After the preprocessing steps, signal interpolation and filtering used in the paper, cross-correlation was performed between the sLFOs from both original BOLD and the shifted BOLD. The resulting delay times and the corresponding maximum cross-correlation coefficients (MCCCs) were shown in Figure S6(B). The MCCC values are very high in all cases (0.970-0.975). The corresponding delay times did deviate from the true value (3s), which are 2.9 s (0-4.9 times the SD) and 2.8 s (5 times the SD). Even with the large motion artifact (5 times the SD), the deviation is within 0.2 s. Therefore, we demonstrated that head motion occurring instantaneously across the whole brain had limited effects on the delay time calculations.   shows simulated original BOLD (in blue color) and shifted BOLD (in red color; shifted by 3 sec) with the motion artifact occuring at the same time (the magnification factor of the motion spike is 5 times the standard deviation of the original BOLD without motion artifact). (B) displays the cross-correlation results between the original BOLD and shifted BOLD with the magnification factors from 0 to 5 times the standard deviation of the original BOLD, the delay time in seconds (blue circle) and the corresponding maximum crosscorrelation coefficients (MCCC, red cross).

Spatial comparison of the delay map from rs-fMRI and the TTP map from DSC-MRI
Positive correlations between the sLFO delay time of rs-fMRI and TTP of DSC-MRI were found in all subjects ( Fig. S7(a)). The plot was divided into two portions. In the first portion shown in blue (about -4.7 s to -1.7 s), the fitted slope of the line between the TTP from the DSC-MRI and the sLFO delay time from rs-fMRI is roughly flat. Those voxels largely reside in central gyrus and gray matter (as shown in Fig. S7(b)). In the second portion shown in red (about -1.7 s to 2 s), the fitted slope of the line between the TTP and the sLFO delay time is close to 1, which is expected. Those voxels largely reside in some parts of gray matter, white matter, and large draining veins (as shown in Fig. S7(c)). These results demonstrated that the spatial blood flow patterns derived from the rs-fMRI using the sLFO and the DSC-MRI are positively correlated with each other in a large extent. However, we can also see the discrepancy between the delay map from rs-fMRI and the TTP map from the DSC-MRI (See Fig. S8). The potential reasons can be found in the main paper.

Slope-detection algorithm
The slope-detection algorithm supports the choice of two input parameters: n, the number of edges to detect, and edge, which determines whether edges will be detected on the rising edge (black-to-white; edge = -1) or falling edge (white-to-black; edge = 1). The chosen carpet plot image ( Fig. S3.A) is input as the image matrix I. To increase signal to noise ratio (SNR), I is first smoothed using the 2D blurring filter

Slope and Transit Time Computation
The smoothed carpet plot Ismooth is filtered using the 2D horizontal derivative filter hd2, where The value nvoxels represents the number of voxels included in the (cropped) carpet plot image input into the slope-detection algorithm, slope is the computed edge slope, and TR is the repetition time for the fMRI scan protocol.
We note that there is some degree of error associated with the computation of the transit times due to time resolution limitations. The computed transit time error can be thought of as a function of three variables: (1) the true underlying transit time, (2) the TR of the imaging sequence, and (3) the time point at which sampling of the true signal begins. The potential error for a given edge is inversely proportional to the true transit time duration, as a longer true duration can be better sampled with a given TR. In the extreme case where the true underlying transit time is less that the TR, then the transit time computation can have an error up to the magnitude of the TR. However, as the true transit time becomes greater than the TR, then the potential error shrinks. In summary, the higher the transit time, the more accurate the measured value. Since we primarily expect transit times in the range of 4-6 seconds, we can expect that the potential error in measurements is much lower than 1.5 seconds.
Tests were performed by creating artificial "perfect" images with a temporal resolution of 0.1 s and tilted edges based on varying true transit times. These artificial plots were then sampled every 15 time points to simulate a TR of 1.5 s, and input into the slope detection algorithm to see how the algorithm performed at measuring the transit time. An example comparison of the computed transit time with the true underlying transit time is shown in Figure S4.A. As listed as variable (3) above, we also note that the transit time computation will vary based on when this sampling begins within the original image, since (in this case) there are 15 possible "starting points" for the sampling. This represents how the propagation of a sLFO could begin at any point within the given TR. Figure S4.B shows the effect of different sampling starting points on the error of the computed signal for a true transit time of 4.5 s. This shows that we can expect that the error of the transit time computation ranges between approximately -0.6 to 0.3 s for a true underlying time of 4.5 s. It was also observed that this type of error tends to cause the transit time to be underestimated slightly more often than it is overestimated. This is a limitation of the transit time computation method, however we believe that these errors are not significant enough to change the general conclusions we draw in the paper; in particular, these errors are much smaller than the variations in rs-fMRI transit times observed within a single subject.

Noise Model Evaluations
The slope-detection algorithm was evaluated for robustness by constructing a simple noise model to analyze how consistently the algorithm performed under variation of two parameters: intensity contrast and noise level. Within this model, carpet plot images were assumed to consist of two components: an underlying clean image and an added image of Gaussian noise. The contrast value for a given carpet plot was estimated by computing the difference between the average value over selected regions of highest values (lightest) and that of selected regions of lowest value (darkest). The noise level of a carpet plot was estimated as follows: the underlying clean image was estimated by applying a 2D smoothing filter (identical to that of the slope-detecting algorithm) to the original image. The noise image was then estimated by subtracting the clean image estimate from the original. The standard deviation over all values in the estimated noise image was then computed and formed the quantified noise level of the given carpet plot. Rs-fMRI carpet plots produced contrasts in the range 0.735-1.471 and noise standard deviations in the range 0.737-0.817. DSC-MRI carpet plots produced contrasts in the range 2.899-3.226 and noise standard deviations in the range 0.197-0.330. Due to the processing and scaling of the fMRI data these measurements are simply in units of signal intensity of the processed carpet plot data.
This model was developed in order to evaluate the reliability of the slope-detection algorithm, completed by scaling a chosen DSC carpet plot to varying contrast (ranging from 0.5 -3.5) and noise levels (ranging from 0.1 -1.0) and conducting repeated trials (n = 30) of random noise addition at each pairing of contrast and noise level to determine the consistency of transit time computation. Scaling of the chosen DSC carpet plot to varying contrast and noise levels was completed as follows: a single cropped DSC carpet plot was chosen to serve as the baseline image and was smoothed using the filter h described above. Then, for a given trial, the smoothed image was scaled to a specified contrast level (completed by scaling the image by a factor defined by the ratio between the desired contrast level and the contrast level of the selected DSC carpet plot) and was added to an equally-sized image of random Gaussian noise at a specified standard deviation. 30 trials of each combination of contrast (ranging from 0.5-3.5) and random noise at a given standard deviation (ranging from 0.1-1.0) were completed. The results of this analysis are given in Figure 7 in the main document.
We note that this noise model represents an oversimplification of the noise observed in rs-fMRI data. The noise which corrupts the BOLD signal is not simply Gaussian in nature, and it may not be ultimately correct to assume perfect and linear underlying edges exist beneath the noise in carpet plot images. For instance, Figure S11 displays a DSC carpet plot image which has been scaled to the contrast and noise level of RS carpet plots compared to an unaltered RS carpet plot image. We note that despite being "equal" in contrast and noise level to the RS image, the DSC image still appears to contain a bolder and cleaner edge that the RS image, indicating that our noise model does not fully capture the characteristics which impact the quality of edges seen in carpet plots. These observations limit our ability to yet claim causal physiological differences between edges of different transit time. Given these limitations, an additional test was run to evaluate the consistency of slope detection where the noise model was used to "reconstruct" the input carpet plots. To complete this "reconstruction", a given carpet plot was first smoothed using the filter h described above. Then an image of Gaussian noise at standard deviation defined by that carpet plot's noise level (computed as described previously) was added back into the smoothed image. The goal of these operations was to recreate a given carpet plot with an equivalent contrast and random noise level so that the resulting transit time computations could be analyzed. 30 trials were completed for each DSC and RS carpet plot, where this reconstruction process was applied and the slopedetection algorithm was then used to find a single edge's transit time (that of the most prominent edge as detected by the algorithm). Results are displayed in Figure S12. We note that application of the reconstruction process resulted in a slight increase in transit time for some carpet plots, while it resulted in a slight decrease in transit time for others, indicating no obviously clear pattern.
The reconstruction trials resulted in a wider variance of transit times for RS carpet plots compared to those from DSC imaging, which is to be expected due to higher SNR in DSC data. While the size of this variance varied slightly from subject to subject, the resulting standard deviations were all small relative to the variance in transit times observed across different time points for a given subject, indicating that the variance in transit time for different edges within the same carpet plot was likely due to true underlying differences as opposed to algorithm error. Figure S12. Noise model reconstruction of DSC and RS carpet plots. For each subject's DSC and RS dataset, the carpet plot was filtered to remove noise, and the standard deviation of removed noise was computed. Random white noise at the computed standard deviation was then added back to the filtered image (creating the reconstructed image) and the transit time recomputed. This was completed for 30 trials per carpet plot. Red points indicate the original transit time for the first slope detected in the unaltered carpet plot. Vertical bars for reconstructed trials (blue) are twice as long as the standard deviation of resulting transit times and centered at the mean for that subject.