{'fraction_non_alphanumeric': 0.045202476141346405, 'fraction_numerical': 0.028157510102312783, 'mean_word_length': 4.513747333491349, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 2, 'https://': 1, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Time-delay structure predicts clinical scores for patients with disorders of consciousness using resting-state fMRI

Highlights • TDE and PFE can characterize the temporal organization of iBA in DOC patients.• DOC patients showed abnormal TDP and PFP mainly in DMN and SMN.• DOC patients had a shortened range of TD value relative to the healthy controls.• The temporal propagation structure of iBA can be used to predict clinical scores for DOC patients.


Introduction
Disorders of consciousness (DOC), including vegetative state/unresponsive wakefulness syndrome (VS/UWS) and minimally conscious state (MCS), refers a continuum of disruption in the wakefulness and awareness systems resulted from severe brain injury (Giacino et al., 2002;Laureys et al., 2010;Lee et al., 2020). As two components of consciousness, wakefulness refers to a physiological state of being awake and awareness is the perception of self and environment (Laureys, 2005). The VS/UWS patients could keep wakefulness but show no signs of awareness (Laureys, 2005). The MCS patients are characterized by inconsistent but reproducible signs of awareness (Giacino et al., 2002). However, it is difficult to determine the state of consciousness for DOC patients due to the lack of overt behavior response from patients (Schnakers et al., 2009). Recently, the functional magnetic resonance imaging (fMRI) technique is used to explore intrinsic brain activity (iBA) Abbreviations: BOLD, blood oxygen level dependent; CCF, cross-covariance function; CPM, connectome-based predictive modeling; CRS-R, Coma Recovery Scale-Revised; DMN, default mode network; DOC, disorders of consciousness; FC, functional connectivity; iBA, intrinsic brain activity; IFG, inferior frontal gyrus; MCC, middle cingulate cortex; MCS, minimally conscious state; PF, probabilistic flow; PF P , lag projection of probabilistic flow; PoCG, postcentral gyrus; PreCG, precentral gyrus; PFE, probabilistic flow estimation; ROI, region of interest; TD, time delay; TD P , lag projection of time delay; TDE, time delay estimation; VS/UWS, vegetative state/unresponsive wakefulness syndrome.
in DOC patients and the altered iBA is considered to assist the clinical assessment for DOC patients (Cao et al., 2019;Song et al., 2018). The iBA refers to the neural states that are produced spontaneously by the brain and not as responses to external stimulus (Boly et al., 2008). The iBA is believed to play a central role in human brain function and can be analyzed using resting-state fMRI (rs-fMRI) data (Raichle, 2015).
Previous studies analyzed inter-regional functional connectivity (FC) based on the rs-fMRI data and found abnormal iBA in the thalamocortical network (Schiff, 2008), frontoparietal network , and default mode network (DMN) (Vanhaudenhuyse et al., 2010) in DOC patients. FC analyses assumed that the iBA among brain regions is synchronized or zero-time lag correlation (Preti et al., 2017). However, recent studies suggested that iBA exhibits dynamic fluctuation (Allen et al., 2014;Vidaurre et al., 2017). The abnormality of fluctuation of iBA has been identified in DOC patients (Cao et al., 2019;Demertzi et al., 2019;Huang et al., 2020;Luppi et al., 2019). For example, the reduced dynamic FC was found in the posterior DMN, sensory network and somatomotor network in DOC patients compared with healthy controls (Cao et al., 2019;Luppi et al., 2019).
In addition to the findings on fluctuation of iBA, accumulating evidence (Ferezou et al., 2007;Mitra et al., 2017b;Mohajerani et al., 2013;Mohajerani et al., 2010) showed an intrinsic brain propagation activity in the whole brain. Using electrodes, several animal studies found that the iBA is asynchronous in the rat brain (Ferezou et al., 2007;Mohajerani et al., 2013;Mohajerani et al., 2010) and inferred a time delay (TD) to describe the iBA propagation from one region to other regions. Based on the concept of TD, Mitra et al. (2014) introduced TD estimation (TDE) to assess propagation structure of iBA between a pair of regions in human brain. Different from dynamic FC analysis, TDE emphasizes the temporal ordering of iBA rather than the temporal fluctuation of iBA. The studies of TDE showed that iBA propagated across whole-brain brain networks on a time scale of ~1s in healthy controls (Mitra et al., 2015a;Mitra et al., 2014). Subsequently, based on TDE and correlation analysis, Mitra et al. (2020) introduced probabilistic flow estimation (PFE) to describe propagation probability of iBA by calculating the probabilistic flow (PF) between the timeseries of each voxel in whole brain. Totally, TDE and PFE offer the possibility to describe the features of intra-and inter-regional temporal organization of iBA. Previous TDE studies have reported the altered propagation structure in autism (King et al., 2018;Mitra et al., 2017a), post-traumatic stress disorder (PTSD) (Weng et al., 2018), and epilepsy patients Shah et al., 2019). As for DOC patients, the impact of severe brain injury on the temporal organization of iBA is remain unclear.
Machine learning has been widely applied to neuroimaging data to aid the clinical assessment for various brain disorders (Edlow et al., 2021;Rehme et al., 2015;Shen et al., 2017). Several studies (Campbell et al., 2020;Cao et al., 2019;Song et al., 2018) used machine learning to explore the imaging indicator that can promote the diagnosis and/or prognosis of DOC patients. Song et al. (2018) proposed a prognostic model by taking the FC as features in a classification algorithm, and discriminated the DOC patients who would later recover consciousness from these not recover with an accuracy of 88%. However, we are still lacking the accurate and efficient neuroimaging indicator for accessing the clinical scores of DOC patients. In this study, we asked the question whether the TD or PF could be used as a neuroimaging indicator to predict clinical scores of DOC patients. Thus, we applied connectomebased predictive modeling (CPM), a data-driven machine-learning method (Shen et al., 2017), to assess the prediction performance of TD and PF for the clinical scores of patients. The CPM was selected because it is a validated approach for extracting and pooling the most relevant features from neuroimaging data to predict the clinical score (Finn et al., 2015;Shen et al., 2017).
Our aims were to reveal the altered propagation structure and propagation probability of iBA and then to evaluate the predictive performance of the temporal organization for clinical scores in DOC patients. We hypothesized that the severe brain injury may affect the propagation structure or propagation probability of iBA in several critical brain regions of DMN and FPN in DOC patients, and the temporal metrics (TD and PF) of iBA may be related to the clinical scores of DOC patients.

Subjects
A total of 45 DOC patients were recruited from the Liuhuaqiao Hospital in Guangzhou City, China. The patients were assessed by two medical doctors (QX and RY) according to the Coma Recovery Scale-Revised (CRS-R) scale, a 23-item scale for measuring the behavioral responses of DOC patients. The scale includes 6 subscales that measure the auditory, visual, motor, oromotor/verbal, communication and arousal levels of the patient, respectively. We also recruited 17 healthy volunteers who had no history of neurological or psychiatric illnesses as the control group. The study protocol was approved by the Ethics Committee of the Liuhuaqiao Hospital. Written informed consent was obtained from each patient's guardian and healthy volunteers.

Data acquisition
The MRI data were acquired on a GE 3 T MR scanner with an eightchannel phased-array head coil. The rs-fMRI data were obtained using a single-shot multi-slices gradient-echo EPI sequence with the following parameters: repetition time (TR) = 2,000 ms, echo time (TE) = 26 ms, flip angle (FA) = 90 • , field of view (FoV) = 240 × 240 mm 2 , data matrix = 64 × 64, 36 axial interleaved slices covering the whole brain, and 240 volumes obtained in about 8 min. Subjects were requested to keep their eyes open when subjects were in the rs-fMRI scanning. In practice, we had checked to make sure the patient eyes open before and after the scan. If a patient cannot keep eyes-opening, we stopped the scan procedure and resumed it later. The high-resolution brain structural images were also acquired using a T1-weighted 3D Fast Spoiled Gradient Recalled (FSPGR) sequence (TR = 8.86 ms, TE = 3.52 ms, FoV = 240 × 240 mm 2 , data matrix = 256 × 256, FA = 90 • , voxel size = 0.94 × 0.94 × 1 mm 3 , and 176 sagittal slices). Additionally, the 2D T1-weighted and 2D T2-FLAIR images were also acquired for clinical diagnosis. Fig. 1a shows the visual quality checking steps. All brain images under analysis were visually inspected and those showing an insufficient level of quality are excluded. The exclusion criteria for the patients were 1) focal brain damage, 2) artifacts, and 3) excessive head movement. In the calculation, for each subject, we have checked the head motion of the fMRI datasets and excluded the whole datasets from the further analysis if the mean framewise displacement mean (FD) > 0.2 (Power et al., 2012). We also excluded the patients who were diagnosed with locked-in syndrome (LIS), or emerged minimally conscious state (EMCS), or had disease course more than one year. According to the exclusion criteria, we excluded 28 patients' data from the recruited 45 patients, and only retained 17 patients. We included 17 healthy subjects for the further analysis. The demographic and clinical profiles of the patients are described in Table 1.

Data preprocessing
The functional images were preprocessed using SPM12 (htt p://www.fil.ion.ucl.ac.uk/spm) and DPABI (Yan et al., 2016). For each individual's rs-fMRI data, we removed the first ten volumes, performed slice-timing correction, corrected head movement. We coregistered functional images to the individual structural images and normalized them in the MNI standard space. Finally, the functional images were smoothed with a 6 mm full-width-at-half-maximum (FWHM) Gaussian kernel and filtered to 0.01-0.1 Hz. In the calculations, the nuisance variables, including 24 head movement parameters (Friston et al., 1996) and the mean signals of brain white matter (WM), cerebrospinal fluid, and whole brain, were regressed out. To reduce the computational load, we resampled the functional data to (6 mm) 3 isotropic voxels and acquired n = 8,434 voxels in the MNI space for each individual brain.

Relationship between FC, TD, and PF
FC and TD are independent metrics to describe the spatial and temporal relationship between a pair of blood oxygen level dependent (BOLD) signals, respectively. Both these metrics could be used to describe systems-level organization of iBA in large-scale network level (Mitra et al., 2020). The conventional FC, by measuring the zero-lag Pearson's correlation coefficient between pairs of timeseries, represents the synchronization between spatially distinct neurophysiological events (Friston et al., 1993). However, the FC may be unsuitable to reveal the directionality of inter-regional FC in the whole brain. In fact, there are several approaches can be taken to infer the directionality of BOLD signal propagation, including dynamic causal modelling (DCM) and Granger causality analysis (GCA). However, DCM is limited in the number of ROIs and specific hypotheses, and thus is unsuited to an Fig. 1. The workflow to illustrate the exclusion criteria for the patients with disorders of consciousness (DOC) and the estimation of time delay (TD) and probabilistic flow (PF). (a) Exclusion criteria for the fMRI datasets of DOC patients. A total of 28 patients' fMRI datasets were excluded because of the following reasons: 1) excessive brain lesion or atrophy, 2) heavy artifact, 3) overlarge head movement, 4) diagnosed as patients with locked-in syndrome (LIS) or emerged minimally conscious state (EMCS), and 5) other types of imaging quality problems. This study retained the resting-state fMRI datasets of 17 patients, including 11 patients with vegetative state/unresponsive wakefulness syndrome (VS/UWS) and 6 patients in a minimally conscious state (MCS), for the further analysis. (b) Illustration of time delay estimation (TDE). We computed the lagged crosscovariate function (CCF) between pairs of timeseries in two different voxels (i and j) and determined the maximum of CCF by using a parabolic interpolation.
The τ i,j was defined as TD value when CCF reached maximum. (c) Illustration of probabilistic flow estimation (PFE). The values of correlation coefficient (r) and the time delay (TD) were normalized to [-1, 1]. The PF value was computed by combining the normalized r and the normalized TD. (d) Calculation of the lag projection of TD and PF (TD P and PF P ). We obtained the TD P (PF P ) by averaging each column of the TD (PF) matrix for each subject. (e) Illustration of connectome-based predictive modeling (CPM). After acquiring the ROI-wise TD/PF matrix based on Shen-268 atlas, we first calculated Pearson's correlation coefficient between the edges (TD/PF value) and the clinical scores, and then summed the significant edges as the feature for each subject. Using the LOOCV, we split the all patients into the training set to build the linear regression model and the test set to predict the clinical score of the remaining patient. At last, the significance of this model was evaluated by 1,000 nonparametric permutation tests. exploratory analysis at the whole-brain level (Stephan et al., 2010). GCA could be used in whole-brain analysis, but the theoretical assumption of GCA may not fit the properties of fMRI data (a low sampling rate and a long delay between the BOLD signal and neural activity) (Friston et al., 2013). Thus, we took TDE and PFE to estimate the direction of interregional interactions.

Time delay estimation
We performed TDE by calculating the lagged CCF between BOLD signals at two voxels (Mitra et al., 2014). The lagged CCF represents the correlation between two signals with all possible lags (Mitra et al., 2014) and is given by where τ represents the lag between the two BOLD signals x i and x j at two spatial indexes i and j, and the τ value is defined as TD when C xixj (τ)reaches the maximum. The zero-lag (τ = 0) CCF is the conventional inter-regional FC. As shown in Fig. 1b, the maximum of C xixj (τ) was obtained using a parabolic interpolation. The space index n represents the total voxels of the whole brain. In this study, the window size is the total scan time t = 230 TRs (TR = 2 s) of the functional data. The τ value was set not to exceed 6 s since BOLD response is typically modeled as the convolution of a stimulus and the hemodynamic response function (HRF) which typically peaks at 5-6 s after human brain respond to the stimulus delivery (Aguirre et al., 1998;Saad et al., 2003). Due to the aperiodicity of BOLD timeseries , the TD value is single and well-defined. The TD between any pairs of timeseries was estimated across all voxels (n = 8,434) in whole brain according to Eq.1, and the resulting set of τ which represents TD value was assembled into a matrix as follows: where TD matrix is a square anti-symmetric matrix. The off-diagonal elements τ i,j = -τ i,j , and the diagonal elements are zeros.
For each subject, we computed the lag projection of time delay (TD P ) by the following formula: (3) where the voxel-wise TD P is calculated through averaging the elements in each column of the TD matrix, which indicates whether the BOLD signal in a given voxel is either earlier or later than the rest of the brain. Fig. 1d illustrates the steps for calculating the voxel-wise TD P for each subject. In each group, the group-level TD P map was obtained by averaging voxel-wise TD P across all DOC patients or healthy controls, indicating the propagation structure of iBA at the group-level. In this way, we obtained the voxel-wise TD P for each subject and obtained the grouplevel TD P map.

Probabilistic flow estimation
The probabilistic flow (PF) was calculated according to Mitra et al. (2020). For a pair of voxel-wise timeseries, we first estimated the TD (τ) value by using lagged CCF and a parabolic interpolation. Meanwhile, the correlation value (r) at τ was acquired using lagged CCF and was termed the peak correlation ( Fig. 1c). This calculation was repeated for all pairs of voxels. To diminish the influence of noise signals on PF, we set the threshold |r| = 0.10 in calculating the peak correlation (Raut et al., 2019b) and set PF = 0 when |r| < 0.10. Next, both r and TD values were separately normalized to [-1, 1] to make sure their magnitudes on the same scale. The PF value, propagation probability for a pair of timeseries at voxels, was given by (Mitra et al., 2020): where r represents the peak correlation coefficient, and τ represents the TD value.
For each subject, we first acquired the asymmetric PF matrix from all pairs of timeseries at voxels, which is given by Mitra et al. (2020) Then we calculated the voxel-wise lag projection of probabilistic flow (PF P ) according to the following equation: where PF i p = i/n n j=1 PF i,j shown in Fig. 1d, n = 8,434 is the number of voxels. The group-level PF P map was obtained by averaging the voxel- wise PF P values across all the DOC patients or the healthy controls. To build probabilistic distributions for the PF, we normalized the cumulative sum of all positive and negative prob to one, separately. The detailed theory and MATLAB code of TDE and PFE can be found in (Mitra et al., 2020;Raut et al., 2019b).

Statistical analysis for TD P and PF P between groups
A nonparametric permutation t-test with 5,000 iterations was applied based on FSL Randomize tool (Chen et al., 2018b;Winkler et al., 2016) to test between-group differences in whole-brain TD P and PF P . The significance threshold p-values were set to 0.05 and was corrected for multiple comparisons using threshold-free cluster enhancement (TFCE). In the calculations, the mean FD, gender, and age were taken as nuisance covariates. Fig. 1e shows the procedures of the CPM analysis. Consistent with previous studies on CPM (Finn et al., 2015;Rosenberg et al., 2016), we selected Shen-268 atlas (Shen et al., 2013) to extract the timeseries of each ROI by averaging all voxels' signals. Next, we estimated the values of TD and PF for all pairs of ROIs by using TDE and PFE, resulting in a 268 × 268 TD matrix and a 268 × 268 PF matrix for each subject. Since the TD and PF matrices are anti-symmetric, we only selected the upper triangle elements as the edges (i.e., TD or PF) for the following analysis.

Connectome-based predictive modeling for the CRS-R scores
We used the leave-one-out cross validation (LOOCV) approach to construct linear regression model in the training set and to predict the clinical score of the remaining patient in the test set. A total of 16 of 17 DOC patients were taken as the training set, and the remaining one patient was selected as the test set. We calculated Pearson's correlation coefficient between CRS-R score (the sum of subscale score) and each edge in the TD (PF) matrix for each patient in the training set. Then, edges that significantly positively and negatively correlated with the clinical score (threshold p edge < 0.001) were retained. These edges will make up the positive (negative) network (Fig. 5c). Next, we summed the strength of edges within the positive (negative) network for each subject as the positive (negative) features.
In the training set, we took the positive and negative features and the CRS-R scores of 16 subjects to build a linear regression model by: where Y is the CRS-R scores, X pos is the positive feature, X neg is the negative feature, w 1 and w 2 are the weight for X pos and X neg . Then the predicted score was calculated using the linear regression model based on the features of the remaining patient in test set. The procedure was repeated 17 times (because of 17 patients) and the predicted CRS-R score of each patient was acquired at last. In the calculation, LOOCV can control the influence of false positives in feature selection because when there was a high proportion of false positives at the feature selection step in the training set, the model cannot generalize well in the new data in the test set (Finn et al., 2015;Rosenberg et al., 2016). We acquired the model performance (r perf ) by calculating the Pearson's correlation coefficient between the predicted and observed CRS-R scores. The p-value was calculated by using 1,000 nonparametric permutation tests. The CPM scripts are available from (https://www.nitrc. org/projects/bioimagesuite/) and visualization tool for the results presentation is available from (http://bisweb.yale.edu/connviewer/).
To test the effect of different thresholds on the accuracy and generalizability of the CPM result (Shen et al., 2017), we set the edges with additional thresholds to select features for each subject. In total, we set five thresholds (threshold p edge = 0.01, 0.005, 0.001, 0.0005, and 0.0001), which include the default threshold p edge = 0.001, to determine the optimal threshold. The predictive procedures were repeated for each threshold.

Demographical and clinical information
The rs-fMRI data of 17 DOC (13 M/4F, aged 40.6 ± 16.4 years old, 11VS/6MCS) and 17 health controls (9 M/8F, aged 32.9 ± 8.6 years old) was analyzed in this study. No significant between-group differences were found in age (t = 1.43, p = 0.09), gender (χ 2 = 2.69, p = 0.15), and mean FD (t = 0.05, p = 0.27). Fig. 2a and 2c show TD P map for the DOC patients and the healthy controls, respectively. The voxel color-coded in cool (warm) indicate that the BOLD signal in the voxel is earlier (later) than that of the rest voxels of the whole brain. Fig. 2e shows the t-map in TD P between the DOC patients and the healthy controls. Fig. 2b and 2d show PF P map for the DOC patients and the healthy controls, respectively. The cool hue represents that a voxel receives BOLD signals from the rest voxels while the warm hue represented that a voxel sends out BOLD signals to the rest voxels. The deeper the hue is, the more likely it is that the voxel will act as the sender/receiver. The sum of all the cool (or warm) hues equals one. Fig. 2f shows the t-map in PF P between the DOC patients and the healthy controls. Fig. 3 shows each cluster with significant TD P difference between the two groups. Relative to the TD P in the healthy controls, the DOC patients showed significant differences in 1) the right precentral gyrus (PreCG), 2) right precuneus, 3) right middle cingulate cortex (MCC), 4) right inferior frontal gyrus (IFG), 5) left PreCG, and 6) right postcentral gyrus (PoCG). Specifically, the right PreCG, right precuneus, right MCC, right IFG, left PreCG show early iBA in the healthy controls but not in the DOC patients. In addition, the right PoCG shows late iBA in the healthy controls but not in the DOC patients. Fig. 4shows the clusters with significant between-group differences in PF P . Compared with the healthy controls, the DOC patients showed significantly different PF P in five clusters, which were located in 1) the right precuneus, 2) right PreCG gyrus, 3) left precuneus, 4) right MCC, and 5) right PoCG. Specifically, the right precuneus, right PreCG gyrus, left precuneus, and right MCC are high probability sender of iBA in the healthy controls but not in the DOC patients, whereas the right PoCG is a high probability receiver in the healthy controls but not in the DOC patients.

Voxel-wise time delay estimation and probabilistic flow estimation
We also calculated the top and bottom quartile of TD P to measure the range of TD value for the DOC patients and the healthy controls, separately. In the calculations, we took the mean FD, age, and gender as covariates. We found that the top quartile of TD P was smaller (t = -3.33, p = 0.002) in the DOC patients (0.147 ± 0.025 s) than that in the healthy controls (0.184 ± 0.026 s), and the bottom quartile of TD P was larger (t = 3.92, p < 0.001) in the DOC patients (− 0.147 ± 0.024 s) than that in the healthy controls (− 0.187 ± 0.027 s). Meanwhile, we also repeated the TDE calculation without regressing out the global signal and found the range of TD in DOC patients was still shorter than that in the healthy controls. Specifically, the top quartile of TD P was significantly smaller (t = -1.97, p = 0.028) in the DOC patients (0.202 ± 0.050 s) than that in the healthy controls (0.242 ± 0.049 s), and the bottom quartile of TD P was significantly larger (t = 2.45, p = 0.009) in the patients (-0.206 ± 0.070 s) than that in the controls (-0.258 ± 0.068 s). Thus, the shortened range of TD in DOC patients was not associated with the global signal. Table 2 lists the detailed model performance (r perf value) among five thresholds based on whole-brain ROI-wise TD matrix and PF matrix, respectively. Specifically, for default threshold p edge = 0.001, we found that the features based on TD matrix (r perf = 0.662, p = 0.006) could be used for significantly predicting the clinical scores of DOC patients. However, the model performance based on TD matrix reduces when the p edge is too large (p edge = 0.005 and 0.01) or too small (p edge = 0.0005 and 0.0001, see Fig. 5a), which indicated that additional/insufficient edges may do affect the models' predictive performance (Ren et al., 2021). In addition, the model performance based on PF matrix is poor (r perf from − 0.639 to − 0.053) compared with based on TD matrix. Because the model performance based on the TD matrix is better than that based on the PF matrix and the best performance based on the TD matrix is acquired when threshold p edge = 0.001, we selected the result based on TD matrix in default threshold p edge = 0.001 to perform the further analysis. For a better visualization, we converted the CRS-R scores by z score = x− μ δ (Beaty et al., 2018), where × is the raw observed (predicted) clinical score of each patient, μ is the mean of the observed (predicted) clinical scores of all patients, δ is the standard deviation of the observed (predicted) clinical scores of all patients. Fig. 5b shows the scatterplots of standardized observed and predicted clinical scores. In addition, to test the robustness of CPM result based on TD matrix, we also repeated the calculations by taking the Brainnetome atlas (BN-274) (Fan et al., 2016) to construct TD/PF matrices at the optimal threshold (p edge = 0.001, Fig. 5a) that determined by Shen-268 atlas. We found that the TD matrix (r perf = 0.505, p = 0.031) rather than the PF matrix (r perf = 0.05) can accurately predict the clinical scores of DOC patients. This result is consistent with that derived from Shen-268 atlas. When using p edge = 0.001, the CPM results shows that a total of 42 edges in TD matrix are significantly correlated with clinical scores. Fig. 5c shows vertical and left view of the glass brain on the distribution Fig. 2. TD P , PF P , and t-map on brain surface for the disorders of consciousness (DOC) patients and the healthy controls. (a) TD P map of the DOC patients (averaged TD P across all DOC patients). The clock-like color bar signifies that the propagated intrinsic brain activity (iBA) is along the clockwise direction. The hour hand and minute hand represent the bottom and top quartile in the range of propagated iBA in the DOC patients or healthy controls. The cool (warm) hues represent that the lag structure of this region is early (late) relative to the rest of the brain. (b) PF P map of the DOC patients (averaged PF P across all DOC patients). The cool hues represent region that send out signals to the other region and the warm hues represent region that receive signals in the whole brain. The sum of all the cool (or warm) hues equals one. (c) TD P map of the healthy controls (averaged TD P across all healthy controls). (d) PF P map of the healthy controls (averaged PF P across all healthy controls). (e) Between-group comparison in TD P . (f) Between-group comparison in PF P . of 42 edges. To characterize the functional roles of predictive networks which consist of the 42 significant edges, we distributed the all 268 ROIs into ten macroscale brain resting state networks (RSNs) according to previous studies (Beaty et al., 2018;Shen et al., 2017) and checked the number of significant edges between RSNs. Fig. 5d shows the distribution of 42 edges in predictive networks. Due to not all brain regions equally involved in the prediction, we only display the nodes connected with more than one edge, including 11 nodes with total 24 edges. Table 3 lists detailed information on the 11 nodes. Fig. 5e shows the propagation direction of these 11 nodes with other nodes in whole brain. Among these 11 nodes in the predictive networks, 5 nodes in the prefrontal cortex (PFC) joins 8 significant edges; 3 nodes in the cerebellum joins 7 significant edges; 1 node in the limbic region joins 3 significant edges; 1 node in the subcortical region joins 2 significant edges; 1 node in the parietal region joins 2 significant edges. Specifically, we found 8 of 11 nodes are distributed in the cerebellum and PFC, and dense edges are the connections between the cerebellum and PFC.

Discussion
In the current study, we first analyzed the altered propagation structure of intrinsic brain activity (iBA) by using the time delay estimation (TDE), and then detected abnormal propagation probability of iBA by using the probabilistic flow estimation (PFE) for the patients with disorders of consciousness (DOC). With the connectome-based predictive modeling (CPM), we also predicted the clinical scores for the patients according to the ROI-wise time delay (TD) and probabilistic flow (PF) matrix of iBA.

Clusters with altered TD P and PF P in DOC patients
After calculating TDE and PFE, we compared the spatial distribution of the clusters with significant between-group differences in TD P and PF P (Fig. 2e and 2f). Specifically, we found four shared clusters with abnormal TD P and abnormal PF P , including the precentral gyrus, precuneus, middle cingulate cortex, and postcentral gyrus. However, the Fig. 4. Clusters with significant between-group differences in PF P . Each cluster is indicated with a circle in the t-map. Bar plot indicates the group-level PF P of the DOC patients and the healthy controls for the corresponding cluster. The PF P > 0 (or < 0) represents that this voxel is a receiver (or sender) relative to the rest voxels. The higher the value of |PF P |, the higher probability the voxel as a sender or receiver. Abbreviations: PoCG, postcentral gyrus; PreCG, precentral gyrus; MCC, middle cingulate cortex; R (L), right (left) hemisphere. range of spatial distribution for each shared cluster with difference in PF P is relatively smaller than that of the cluster with difference in TD P .
We observed abnormal TD P and PF P in the precuneus (or PCC/precuneus) between groups (Figs. 3 and 4), which showed the altered temporal ordering of iBA in precuneus in the DOC patients. The precuneus, a hub region of the DMN (Buckner et al., 2008), plays a vital role in maintaining self-referential processing and introspection (Cavanna and Trimble, 2006). Previous studies reported the abnormal functional interactions of precuneus in DOC patients by using either static FC (Qin et al., 2015;Stender et al., 2014;Wu et al., 2019) or dynamic FC (Cao et al., 2019;Luppi et al., 2019). Several recent dynamic causal modelling (DCM) studies (Chen et al., 2018a;Crone et al., 2017;Crone et al., 2015) also found the disrupted effective connectivity between PCC/precuneus and other regions (e.g., medial frontal cortex, inferior parietal lobule, and globus pallidus) in DOC patients, reflecting the abnormal information flow related to precuneus at ROI level.
Moreover, our studies found that the precuneus possibly acted as a source in propagation structure of iBA in healthy controls, which was in line with previous voxel-wise TDE studies (Mitra et al., 2015a;Mitra et al., 2014). However, the role of precuneus in propagation structure was lost in the DOC patients (Fig. 3), which may reflect that a steady propagation activity of iBA in precuneus is associated with the  Table 2 The model performance (r perf ) and its significant level (p) in CPM predicting CRS-R scores for DOC patients under different thresholds of the edge selection. The permutation test will be not calculated when the model performance r perf < 0.2. *, p < 0.01 (Bonferroni correction, p < 0.05/5 = 0.01).

Table 3
The 11 nodes connected with more than one edge in CPM according to the TD matrix. maintenance of consciousness. We inferred that the alteration in precuneus may be related to the disrupted brain WM connection in DOC patients (Fernández-Espejo et al., 2012;Lant et al., 2016;Zheng et al., 2017). Using the diffusion tensor imaging (DTI) technique, previous studies revealed the disrupted WM tracts between the precuneus and other regions (e.g., thalamus, hippocampus, and temporoparietal junctions) in DOC patients (Fernández-Espejo et al., 2012;Tan et al., 2019). The DOC patients showed the altered TD P and PF P in both the precentral gyrus (PreCG) and postcentral gyrus (PoCG) compared with the controls (Figs. 3 and 4). In the healthy controls, we found that the timing of PreCG was earlier than other regions, indicating that it possibly acted as a source for the information transmitting, while the timing of PoCG was later than other regions, indicating that it possibly acted as a destination for the information reception and integration (Mitra, Kraft, et al., 2018;. Their roles in the brain information flow may reflect the brain physiological function. Specifically, both the PreCG and PoCG belong to the sensorimotor network (SMN), the PreCG is responsible for the body movement (Lemon, 2008), which may need to transmit the information to other regions. The PoCG is responsible for the processing of peripheral sensory information (Kropf et al., 2019), which may correspond to the information reception and integration in propagation structure. The altered temporal ordering in PreCG and PoCG may be related to the impaired sensorimotor function in clinical performance for the DOC patients. Previous fMRI studies also reported abnormal brain activity in the PreCG and PoCG in DOC patients (Demertzi et al., 2015;Ovadia-Caro et al., 2012;Yao et al., 2015). Yao

Table 4
Summary of human brain resting-state fMRI studies based on the time delay estimation conducted with lagged cross-covariance function (Jan. 2014-Dec. 2020). ±0.5 s 1) Significant altered TD P was found in PCC/precuneus, middle prefrontal cortex, right angular, and left pre-and postcentral cortex in PTSD patients.
2) There was a correlation between clinical scores and TD P in PTSD. (Weng et al., 2018) 9 Specific propagation structure of iBA between the DMN and the DAN in HC HC: 1,376 TR = 3 s 3 T ±1 s 1) The iBA in DMN propagated from the retrosplenial cortex to prefrontal cortex.
2) The iBA in DAN propagates from the frontal eye fields to the parietal cortex.  10 Relationship between propagation structure of iBA in PE patients and laterality PE: 26 TR = 2.07 s 3 T

±1 s
The TD of PE patients was correlated with seizure foci lateralization. (Shah et al., 2018) 11 Propagation structure of iBA between the cerebellum and cerebral cortex in HC HC: 10 TR = 2.2 s 3 T ±0.8 s The iBA in cerebellum lagged the brain cortex about 0.1-0.4 s in HC. (Marek et al., 2018) 12 Potential propagation structure of iBA in ASD patients ASD: 23 HC: 23 TR = 2.2 s 3 T ±0.5 s 1) Significant changed TD P was found in frontopolar cortex, occipital cortex, and putamen in ASD patients.
2) There was a correlation between clinical scores and TD P in ASD. (Mitra et al., 2017a) 13 Similarity in propagation structure of iBA between sleeping infants and awake/sleeping adults Infant: 45 Adult: 63 TR = 2/2.5 s 3 T

±0.75 s
The propagation structure of iBA for infants was more similar to sleeping adults rather than awake adults. (Mitra et al., 2017b) 14 Differences in propagation structure of iBA between the hippocampus and brain cortex under awake and SWS HC: 38 TR = 2.08 s 3 T

±1.5 s
The iBA propagated in opposite directions between the hippocampus and cerebral cortex under awake and SWS.  15 Propagation structure of iBA in HC under wake and SWS HC: 63 TR = 2.08 s 3 T ±0.75 s Significant differences in propagation structure between cerebral and subcortical cortex were found under wake and SWS. (Mitra et al., 2015b) 16 Describing lag processes of iBA in HC HC: 1,376 TR = 3 s 3 T ±1 s Multiple, highly reproducible, temporal sequences of iBA were reported in HC, which termed lag thread. (Mitra et al., 2015a) 17 Temporal features of propagation structure of iBA in HC HC: 692 TR = 3 s 3 T ±0.5 s 1) The iBA propagated across the whole brain on a timescale of about 1 s. 2) Propagation activity of iBA reflected neuronal processes rather than hemodynamic delay. (Mitra et al., 2014) Abbreviations: ASD, Autism spectrum disorder; BOLD, blood oxygen level-dependent; CCF, cross-correlation function; DAN, dorsal attention network; DMN, default mode network; DOC, disorders of consciousness; HC, healthy controls; iBA, intrinsic brain activity; PE, pediatric epilepsy; PTSD, Posttraumatic stress disorder; TD P , lag projection of time delay; RE, Rolandic epilepsy; SWS, slow wave sleep; TLE, Temporal lobe epilepsy; TR, repetition time; PCC, posterior cingulate cortex. et al. (2015) found a decreased amplitude of low-frequency fluctuations (ALFF) in the bilateral PreCG in DOC patients. Ovadia-Caro et al. (2012) found a decreased FC between the bilateral PoCG in DOC patients compared with healthy controls.

Decreased range of TD in DOC patients
The current study found that the range of TD value for the patients was shorter than that of the healthy controls by comparing the quartile between groups (Fig. 2c). We observed that the TD P was in a range of 1 s across the controls (Fig. 2c), consisting with most previous TDE studies listed in Table 4. A possible explanation was that the impaired propagation structure of iBA in critical regions lead to the shortened TD P range in the patients. The precuneus, PreCG, and PoCG were a source or destination of the propagation activity in the controls but not in the DOC patients (Fig. 3). Another reason for the shortened TD P range may be the low level of brain activity in DOC patients (Demertzi et al., 2019;Fridman et al., 2014;Sinitsyn et al., 2018). Previous studies found a global decreased metabolic regulation in DOC patients compared with healthy controls (Di Perri et al., 2016;Fridman et al., 2014).

CPM predicting clinical scores based on TD and PF matrices
Using CPM, we found that the ROI-wise TD matrix can be used to predict clinical scores for DOC patients (Fig. 5a). The predictive results expanded the application of CPM in predicting the clinical scores through using the inter-regional temporal information (e.g., TD) in addition to the inter-regional spatial information (e.g., FC). However, the ROI-wise PF matrix showed a worse predictive performance compared with the TD matrix in predicting the clinical scores for DOC patients. This may be associated with the distinct properties of TD and PF. The PF was built on the TD and the peak correlation. To reduce the influence of low correlations from noise signals on the PF, we set the PF between pairs of regions as zero when the peak correlation value r < 0.1 by following Mitra et al., (2020). These zero values in PF matrix may distort Pearson's correlation coefficient between the PF value and clinical scores at feature selection, which may further hinder linear model construction in training set and the model prediction in the test set.
For the results derived from TD matrix, we found that the nodes playing a predictive role were mainly located in the cerebellum and PFC (Fig. 5e). In fact, the cerebellum is believed involving in many cognitive functions in addition to maintaining movement, such as working memory, task switching, and language processing (Kansal et al., 2017). Cao et al. (2019) found that the dynamic FC between the cerebellum and somatomotor network could be used to predict clinical scores of DOC patients. In addition, the PFC is believed subserving for decision making, planning complex behavior, and moderating social behavior (Carlen, 2017;Miller, 2000). Previous studies indicated that WM tracts between the PFC and thalamus was critical in sustaining consciousness (Schiff, 2010) and the WM tracts have been found to be correlated with clinical scores (Lant et al., 2016).
In addition, previous studies showed that the cerebellum and PFC play key roles in propagation structure of iBA (Marek et al., 2018;Mitra et al., 2015a). For example, Marek et al. (2018) showed that the BOLD signals of cerebellum were lagged behind the entire cerebral cortex at both individual-level and group-level in healthy controls based on the Midnight Scan Club (MSC) data. Mitra et al. (2014) analyzed 692 healthy subjects' fMRI data and reported that the cerebellum and PFC were located at the late stage in the propagation structure. Based on the above description, our conjecture is that the cerebellum and PFC are the receiver of the iBA, and the impeded degree of propagation structure in these regions may be helpful for accessing the clinical scores of DOC patients. Further studies are needed to uncover the mechanism behind the relationship between the CRS-R scores and the identified predictive networks before the findings can be translated into clinical practice.

Limitations
Several limitations should be addressed. First, the small sample size may affect the generality of the results. This study included only 17 DOC patients as the trade-off between the sample size and imaging quality, although we obtained the rs-fMRI and high-resolution brain structural images from 45 DOC patients. Further studies could test the results of TD P , PF P , and CPM in a larger DOC patients sample.
Second, this study revealed that the precuneus acts as a source of propagation structure in the healthy controls, which was in line with previous TDE (Mitra et al., 2015a;Mitra et al., 2020) and GCA studies (Deshpande et al., 2011;Garg et al., 2011). However, we also noticed that the precuneus was a receiver of information in several studies Deco et al., 2021;Seguin et al., 2019). This discrepancy may be caused by the different analysis level including voxel level or ROI level. Our result derived from a voxel-wise TDE, while the other results were mainly estimated on the ROI-wise analysis. Crone et al. (2015) estimate the interaction relationship between 4 ROIs, including the medial frontal cortex (MFC), PCC, left inferior parietal lobule (IPL), and right IPL by using DCM and found that the PCC acts as the main driven hub in healthy controls but not in DOC patients. Based on the voxel-wise TDE, however, Mitra et al. (2014) analyzed the directionality of BOLD signals in whole-brain voxels and showed the precuneus was a source of the propagation structure in healthy controls. Thus, we also estimated propagation structure of iBA for either the DOC patients or healthy controls at a ROI level to analyze the temporal ordering of BOLD signal between the PCC/precuneus and the other regions. The results showed that the iBA in PCC/precuneus was not earlier than all the rest ROIs in the healthy controls (see Supplemental Materials).
Third, this study only applied rs-fMRI to analyze the BOLD signals of DOC patients, which lacked the verification of multimodal imaging techniques. In the future, the more reliable results about propagation structure of iBA could be obtained by combining multimodal noninvasive techniques (Peterson et al., 2020), such as electroencephalogram (EEG) (Claassen et al., 2019;Gui et al., 2020;Pan et al., 2020) or magnetoencephalography (MEG) (Boto et al., 2018) or transcranial magnetic stimulation (TMS) (Rosanova et al., 2012) or transcranial direct-current stimulation (tDCS) (Bourzac, 2016).
TDE has been used in the studies of healthy subjects  and patients (Mitra et al., 2017a;Shah et al., 2019) for characterizing temporal organization of iBA (Mitra et al., 2017a;Raut et al., 2019a). To better understand the application of TDE, we summarized the last 7 years (Jan. 2014-Dec. 2020) healthy controls and patient's fMRI studies on TD in Table 4. The PFE could be used to reflect the probability of a given region in receiving or sending out information with other regions. In addition to constructing a PF matrix by following Mitra et al., (2020), we suggested that the TD and FC matrix may be used to build a directional FC network for reflecting the inter-regional interaction directionality across the brain for patients with DOC or other brain diseases in the future work.

Conclusion
In the current study, we estimated time delay (TD) and probabilistic flow (PF) of BOLD signals to infer the propagation structure and propagation probability of intrinsic brain activity in DOC patients. We found that the DOC patients showed not only the altered TD and PF in DMN and SMN but also a shortened range of TD compared with the healthy controls. The findings of the aberrant temporal properties of intrinsic brain activity in unconsciousness may contribute to understanding the pathophysiology of DOC patients. With CPM, we predicted CRS-R scores of DOC patients by using the whole-brain ROI-wise TD and PF matrices. The CPM analysis showed that the TD rather than the PF had a superior predictive performance to the CRS-R scores, in which the cerebellum and prefrontal cortex were the main predictive regions. The results suggested that the propagation structure of iBA in DOC patients may provide a potential neuroimaging indicator to aid clinical applications.

Funding
The study was supported by grants from the National Natural Science Foundation of China (Grant numbers: 82171914, 81871338, 81371535, 82171174, and 81974154) the National Key R&D Program of China (Grant number: 2018YFC1705006), and Key Realm R&D Program of Guangzhou (Grant number: 202007030005). The funding organizations played no further role in study design, data collection, analysis and interpretation, or paper writing.