The Backbone Network of Dynamic Functional Connectivity

Asadi, N., Olson, I. R. & Obradovic, Z. (2021). The Backbone Network of Dynamic Functional Connectivity. Network Neuroscience. Advance publication. https://doi.org/10.1162/netn_a_00209 Journal: NETWORK NEUROSCIENCE 1 RESEARCH 2 The Backbone Network of Dynamic Functional Connectivity: Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 3 Nima Asadi1 , Ingrid R. Olson2,3 , Zoran Obradovic1 1 4 Department of Computer and Information Sciences, College of Science and Technology, Temple University, Philadelphia, PA 19122 2 5 Department of Psychology, College of Liberal Arts, Temple University, Philadelphia, PA 19122 3 6 Decision Neuroscience, College of Liberal Arts, Temple University, Philadelphia, PA 19122 7 Keywords: Dynamic functional connectivity, Backbone network, Null model, Optimization, Autism Spectrum Disorder ABSTRACT 8 Temporal networks have become increasingly pervasive in many real-world applications, including the 9 functional connectivity analysis of spatially separated regions of the brain. A major challenge in analysis 10 of such networks is the identification of noise confounds, which introduce temporal ties that are 11 non-essential, or links that are formed by chance due to local properties of the nodes. Several approaches 12 have been suggested in the past for static networks or temporal networks with binary weights for 13 extracting significant ties whose likelihood cannot be reduced to the local properties of the nodes. In this 14 work, we propose a data-driven procedure to reveal the irreducible ties in dynamic functional 15 connectivity of resting state fRMI data with continuous weights. This framework includes a null model 16 that estimates the latent characteristics of the distributions of temporal links through optimization, 17 followed by a statistical test to filter the links whose formation can be reduced to the activities and local 18 properties of their interacting nodes. We demonstrate the benefits of this approach by applying it to a 19 resting state fMRI dataset, and provide further discussion on various aspects and advantages of it. AUTHOR SUMMARY  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 20 In this work we propose an optimization-based null model to infer the significant ties, meaning the links 21 that cannot be reduced to the local strengths and properties of the nodes, from the dynamic functional 22 connectivity network. We asses multiple aspects of this approach and demonstrate that it is adaptable to 23 most temporal segmentation methods. We demonstrate that this approach provides several advantages 24 such as taking into account the global information of the network. We also compare the proposed model 25 with several commonly applied null models empirically and theoretically. Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 INTRODUCTION 26 Dynamic functional connectivity (dFC) has been widely used to analyze temporal associations among 27 separate regions of the brain as well as the correlation between functional patterns of connectivity and 28 cognitive abilities(Allen et al., 2014; Jones et al., 2012; Van Dijk et al., 2010; von der Malsburg, Phillps, 29 & Singer, 2010). In order to identify co-activation patterns of dFC over the period of experiment, a 30 temporal segmentation (such as sliding window) is commonly applied on the time courses of BOLD 31 activation of brain regions to divide them into consecutive temporal windows(Allen et al., 2014; 32 Hutchison et al., 2013; Smith et al., 2012). Then, the connectivity between separate regions is measured 33 to generate one graph adjacency matrix per each temporal window(Damaraju et al., 2014). Building on 34 this core framework, several enhancements have been proposed in the past years, such as different 35 temporal segmentation approaches, to increase the power and precision of dFC analysis(Chang & Glover, 36 2010; Heitmann & Breakspear, 2018; Hindriks et al., 2016; Kiviniemi et al., 2011). 37 However, a major challenge in analysis of dynamic functional connectivity is to distinguish and 38 address the existing noise confounds in the data, which influence the brain connectivity measures and the 39 structure of the dFC network(Birn, Smith, Jones, & Bandettini, 2008; Chang & Glover, 2009; Kalthoff, 40 Seehafer, Po, Wiedermann, & Hoehn, 2011; Shmueli et al., 2007). This issue especially intensifies with 41 the increase in spatial resolution of the analysis as well as in resting state fMRI data(Birn, 2012; 42 Hallquist, Hwang, & Luna, 2013; Kalthoff et al., 2011). There are several possible sources of noise in 43 resting state fMRI data, including displacements, even as small as a millimeter or less, which could add 44 random noise to the generated time series, and therefore decrease the statistical power in resting state 45 functional connectivity (rsFC) analysis(Van Dijk, Sabuncu, & Buckner, 2012). Even more challenging, it –2–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 46 can result in false positive or negative activation if the displacements are correlated with the 47 stimuli(Lydon-Staley, Ciric, Satterthwaite, & Bassett, 2019; Patanaik et al., 2018; Savva, Kassinopoulos, 48 Smyrnis, Matsopoulos, & Mitsis, 2020). Cardiovascular and respiratory signals are also widely identified 49 as a source of noise, causing synchronized fluctuations in MRI signal(Glover & Lee, 1995). 50 Due to these challenges, neuroscientists often face the concern of analytical models being Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 51 noise-induced(Choe et al., 2017; Gorgolewski, Storkey, Bastin, Whittle, & Pernet, 2013; Murphy, Birn, 52 & Bandettini, 2013). A number of correction techniques have been suggested in the past to reduce the 53 influence of these confounds, including modelling fMRI signal variations using independent measures of 54 the cardiac and respiratory signal variations(Behzadi, Restom, Liau, & Liu, 2007; Bollmann et al., 2017; 55 Bright & Murphy, 2017; Murphy et al., 2013). However, the effect of various sources of noise on the 56 dynamic connectivity of fMRI data is yet to be addressed through a data-driven and systematic 57 framework(Beall & Lowe, 2007; Kundu, Inati, Evans, Luh, & Bandettini, 2012). 58 Moreover, temporal ties that can be reduced to node properties can exist between nodes due to the 59 nature of the data itself. Highly active regions could in principle form a larger number of trivial ties with 60 other regions, and reciprocally, the information of ties that regions with lower activity form can be lost in 61 common analytical procedures(Gemmetto, Cardillo, & Garlaschelli, 2017; Kobayashi, Takaguchi, & 62 Barrat, 2019). In general, if the network representation of a real-world system can be inferred based on 63 local properties of the nodes, such as their activity level or degree, the true interaction and functional 64 homologies between the nodes can not be detected(Gemmetto et al., 2017). 65 Therefore, the objective of this work is to put forward a data-driven approach to distinguish the 66 significant ties that construct the functional connectivity of the brain from ties that are the result of 67 random observational errors or chance. The latter group of temporal links are known as reducible ties, 68 whereby they can be fully attributed to intrinsic node-specific features such as degree or strength of their 69 link weights. On the other hand, the temporal ties that are incompatible with the null hypothesis of links 70 being produced at random are known as irreducible or significant ties, and the network of such significant 71 ties is known as the backbone network. Therefore, the goal of this study is to develop a data-driven 72 framework to infer the two-dimensional backbone network from the multilayer network of dynamic 73 functional connectivity. –3–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 74 Multiple approaches have been proposed to extract the significant ties in a network through statistical 75 means, most of which target static networks(Alvarez-Hamelin, Dall’Asta, Barrat, & Vespignani, 2005; 76 Casiraghi, Nanumyan, Scholtes, & Schweitzer, 2017; Gemmetto et al., 2017; Kobayashi et al., 2019; Ma, 77 Ma, Zhang, & Wang, 2016; Nadini, Bongiorno, Rizzo, & Porfiri, 2020; Serrano, Boguná, & Vespignani, 78 2009; Tumminello, Micciche, Lillo, Piilo, & Mantegna, 2011; Yan, Jeub, Flammini, Radicchi, & 79 Fortunato, 2018). Across these approaches, a key step towards inferring the backbone network is the Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 80 formulation of a reliable null model to characterize the reducible fraction of the temporal interactions, 81 and to steer the process of filtering that fraction of network links. Several null models have been 82 suggested in the literature whose focus is on static networks, spanning from basic weight thresholding of 83 multilayer networks to more advanced techniques(Cimini et al., 2019; Kobayashi et al., 2019; Li et al., 84 2014; Tumminello et al., 2011). 85 One of the main disadvantages with weight thresholding approaches is that they commonly fail to 86 control for the difference in intrinsic attributes of the nodes, thus favor highly active nodes or nodes with 87 other strong local properties, which can potentially have a large number of reducible links. A number of 88 null models have been used to evaluate the statistical significance of dFC based on generating null data 89 using randomization frameworks. Two main approaches of this type include autoregressive 90 randomization (ARR) and phase randomization(PR)(Allen et al., 2014; Chang & Glover, 2010; 91 Handwerker, Roopchansingh, Gonzalez-Castillo, & Bandettini, 2012; Zalesky, Fornito, Cocchi, Gollo, & 92 Breakspear, 2014). In this category of time series-based approaches, null hypothesis testing is then 93 applied by comparing statistics from the original data against those from the generated null data. A 94 backbone approach named significant tie filtering (ST filter) for dynamic networks was proposed by 95 Kobayashi et. al, based on a network modeling concept named activity-driven network (ADN) where 96 individual propensity of generating connections over time is determined by a latent nodal parameter 97 commonly known as activity, and the probability of creating a link at a specific time instant between two 98 nodes is the product of the individual latent activities of interacting nodes(Perra, Gonçalves, 99 Pastor-Satorras, & Vespignani, 2012; Starnini & Pastor-Satorras, 2014; Zino, Rizzo, & Porfiri, 2017). 100 Due to their analytical flexibility and interpretability, activity-driven network models have gained 101 popularity in explaining features of real networks in various areas of research(Liu, Perra, Karsai, & 102 Vespignani, 2014; Rizzo, Frasca, & Porfiri, 2014; Zino et al., 2017). However, in mentioned studies, a –4–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 103 binomial or Poisson distribution is considered for the temporal connections over time, which limits the 104 approach to unweighted networks(Kobayashi et al., 2019; Nadini et al., 2020), whereas many relational 105 networks based on real data, including various types of fMRI-based networks have continuous weights 106 containing significant information regarding the interactions between the nodes as well as the local and 107 global properties of the network. Therefore, inspired by the work of Kobayashi et. al, we propose an 108 approach for extracting the significant ties for temporal networks with continuous weights that meet the Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 109 characteristics of normality and independence of temporal ties, which are discussed in the methodology 110 section. We demonstrate that this methodology controls for intrinsic local node attributes, with a null 111 model that not only takes into account the global structure of the network, but also the temporal 112 variations of the dynamic connectivity links. In the next section we explain the proposed approach in 113 detail, followed by the experimental results on a real dataset of resting state fMRI. We then present an 114 analysis of the results and discuss the advantages and shortcomings of our approach. METHODOLOGY 115 In this section, we outline the methodological framework for identifying the significant links from the 116 networks of dynamic resting state functional connectivity. A key step towards extracting the backbone 117 network is the formulation of a valid and robust null model. For the sake of simplicity, we name our 118 proposed approach the weighted backbone network (WBN). A null model assumes that all connections 119 are formed randomly, meaning that the probability of an interaction between two nodes at a specific time 120 window and the weight of interactions between them could be explained by chance(Gemmetto et al., 121 2017; Kobayashi et al., 2019; Nadini et al., 2020). The objective of inferring the backbone network is 122 thus to detect links that are not compatible with the null hypothesis, meaning that their formation or 123 strength is not driven by chance. 124 The null model that we present can be interpreted as a temporal fitness model, which is characterized 125 by latent parameters that shape its distribution. In this vein, the first step is to estimate these parameters 126 which are not directly observed from the data. For this purpose, we use a maximum likelihood estimation 127 approach that exploits the global and temporal information of the network of dynamic connectivity. We 128 discuss the details of this methodology in the next section. –5–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 129 Estimation of latent distribution variables 130 We consider a dynamic network of N nodes with links evolving over τ observation windows of size ∆ 131 such that t = 1, ..., τ . At each time step t, a weighted undirected network is formed whose adjacency 132 matrix At stochastically varies in time, and the weights of temporal links (links that are formed at time 133 step t) between each pair of nodes i and j form a Gaussian distribution over the τ time steps. Normality of the distribution of weights between each pair of nodes over time τ is concluded based on Central Limit Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 134 135 Theorem and the assumption that the distribution of temporal weights has a finite variance(Dudley, 1978, 136 2014; Haller & Bartsch, 2009; Smith, 2012). Moreover, an empirical assessment of normality of the 137 distribution of temporal weights on a real dataset of resting state fMRI is provided in the result section. 138 We define a temporal null model in which each node i is assigned two intrinsic variables ai , bi ∈ (0, 1], 139 that rule the probability of mean µ and standard deviation σ of the temporal distribution of its interactions 140 with other nodes over τ time steps, such that: µi,j = ai × aj (1) σi,j = bi × bj 141 Therefore, a each parameter of the distributions of temporal ties between each pair of node i and j is 142 the realization of a Bernoulli variable. The null model thus lays out a baseline for the expected mean and 143 standard deviation of the distribution of interactions between two nodes over τ time given their intrinsic 144 variables, if interacting nodes are selected at random at each time step. 146 To uncover significant links with regards to the null model described above, we proceed in two steps. 147 First, given a set of weighted undirected temporal networks with N nodes, we estimate the intrinsic 148 variables a∗ = (a∗1 , ..., a∗N ) and b∗ = (b∗1 , ..., b∗N ) by calculating the maximum likelihood estimation of the 149 set of parameters for each node. For this purpose, we consider the joint probability function over τ time t 150 intervals and edge weights w ∈ (wi,j ; i, j ∈ 1, ..., N ; t ∈ 1, ..., τ ) for the entire temporal connections of 151 the network. Therefore we have: –6–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 145 Figure 1. A schema of the Backbone network inference procedure –7–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: τ 1 2 √ e−(wi,j −µi,j ) /2σi,j t 2 Y Y t f (wi,j |µi,j , σi,j ) = (2) i,j,i6=j t=1 σi,j 2π 152 Where µi,j and σi,j denote the mean and standard deviation of the distribution of temporal edges 153 between nodes i and j observed over τ time steps in the null model. t 154 The log-likelihood function for the empirical data wi,j (weight of the link between i and j at time Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 155 interval t) with replacing the values of µi,j = ai .aj and σi,j = bi .bj will lead to: X t log(f (wi,j |µi,j , σi,j )) = [−n log (bi bj ) i,j,i6=j τ (3) n X 1 − log 2π − (wt − ai aj )2 ] 2 t=1 2(bi bj )2 i,j 156 By differentiating the log-likelihood function with respect to the first parameter, ai , and setting it to 157 zero we have: τ o X X wi,j X [a∗i a∗j − ]= [a∗i a∗j − wi,j o ] = 0, ∀i = 1, ..., N (4) j,j6=i t=1 τ j,j6=i 158 Similarly, by differentiating the log-likelihood function with respect to bi and setting it to zero we have: τ X X o (wi,j − a∗i a∗j )2 [−(b∗i b∗j )2 + ] = 0, , ∀i = 1, ..., N (5) j,j6=i t=1 τ 159 In which the the maximum likelihood estimation of a∗i for every node i was calculated from equation 160 4. Therefore, for a temporal network with N nodes, the pair of latent variables ai , bi for each node i can 161 be estimated by solving the system of N nonlinear equations 4 and 5. The system of nonlinear equations 162 can be solved through a standard numerical algorithm such as the Newton method. The initial values of 163 ai and bi are calculated by dividing the temporal degree of node i averaged over τ time steps by the 164 doubled number of total temporal edges as follows: τ XX s X ai = wi,j /τ 2 ∗ wi,j /τ (6) j,j6=i t=1 i<j –8–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 165 The general schema of the proposed methodology is provided in figure 1. Note that the proposed 166 maximum likelihood approach incorporates the global information of the network as the weights of all 167 temporal links are considered in the system of equations. 168 Selection of significant ties 169 After estimating the latent distribution variables for each node, we then compute for each pair of nodes i Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 170 and j the probability distribution of their interaction over the τ time steps in the null model: 1 2 e−(wij −µi,j ) /2σi,j t ∗ ∗ 2 t g(wij |µ∗i,j , σi,j ∗ )= ∗ √ (7) σi,j 2π 171 ∗ where µ∗i,j and σi,j are the mean and standard deviation based on the estimated latent variables 172 ai , bi , aj , bj through maximum likelihood estimation and equation 1. In order to determine the t 173 reducibility of a temporal link wi,j between nodes i and j at time t, it is compared against the c-th c 174 percentile (0 ≤ c ≤ 100) weight wij of the maximum likelihood estimated distribution of temporal links t c 175 between i and j. If the empirical value wi,j is larger than the wij , then it cannot be explained by the null c t 176 model at significance level α = 1 − . Therefore, the link wi,j is determined to be a significant tie. The 100 177 significance level α is given as an input to the model, providing a systematic adjustable filtering 178 mechanism. The significance threshold α can also be assigned with Bonferroni correction, in which it 179 can be adjusted by dividing by the sum of weights of edges to control for false positives. The P-value of 180 the test is thus given by: wc ij X p=1− g (8) wij =0 181 In order to determine whether a significant tie wij exists between nodes i and j, we simply count the 182 number of times that the temporal link between them was found to be significant according to the 183 threshold value α. If the count of significant links between i and j falls above half of the τ time intervals, 184 i.e. count > τ2 , the link wij is retained in the backbone network. Note that the final backbone network is 185 a binary network, meaning that the weight of links is 1 if the link between two nodes is determined to be 186 significant, and 0 otherwise. However, a weighted network of significant ties can be easily created –9–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 196 Figure 2. Figure a: an example of the effect of estimated a (distribution mean parameter), on admissibility of an empirical temporal link w where the 197 threshold α is 0.1 (90th percentile). A link with a low weight can be admitted to the backbone network as long as its estimated mean is sufficiently low 198 (the blue distribution), therefore, controlling for the effect of high intrinsic distribution weights on acceptance to the backbone network. Figure b: correlation 199 between estimated latent distribution mean variable for each node i (ai ) and the aggregated dFC weights corresponding to the node over τ time steps for the 200 left hippocampus. The weighted degree-estimated a pair values are averaged across all subjects within the study dataset. 187 through various error measures such as averaging the difference between the weights of the temporal t c 188 links wij and the c-th percentile weight wij of the distribution. 189 An important property of the proposed null model is that the tie between two nodes at time t can be t 190 significant even if the weight of temporal link wi,j is small, with the condition that their individual latent 191 variables a, b, and in turn the mean and standard deviation of their temporal distribution, are sufficiently 192 low. On the contrary, ties with large weights might not be deemed significant by WBN if their estimated 193 a, b are large. This property is illustrated in figure 2, where large estimated µ = ai .aj shifts the c − th 194 percentile threshold to the right side of the distribution such that it becomes increasingly difficult for 195 temporal links to meet the threshold. 203 Moreover, strong correlation exists between the MLE based estimated values of distribution means 204 (ai .aj ) and the degree of the nodes, calculated as the sum of weights of the edges over τ time intervals. 205 An example of such correlation is shown in figure 2 for left hippocampus (283 voxels) (further empirical 206 results are provided in the supplementary information) where the aggregated node degree-estimated –10–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 201 Figure 3. Correlation between share of significant ties connected to each node and the MLE estimated latent variable a for right and left hippocamous 202 regions. 207 latent variable a were averaged across all subjects of the study data. Also, as figure3 shows, there exists a 208 weak negative correlation between the share of significant ties that is connected to each node i, and the 209 MLE estimated variables ai corresponding to it. The share of significant ties are calculated as the number 210 of ties connected to node i that are admitted to the final backbone network divided by the total edges 211 connected to it (N − 1). These results establish the property that, based on WBN model, the admissibility 212 of an edge wi,j to the irreducible network is not attributed merely to its degree, therefore controlling for 213 the effect of local strengths of nodes. 214 In the next section, we assess and compare the backbone networks detected based on WBN with 215 autoregressive randomization (ARR) as well as phase randomization(PR). ARR and PR are different 216 from ST filter and WBN in the sense that they are applied to the fMRI time series of each temporal 217 window before drawing the connectivity maps of the brain regions, and are used to explain the fluctuation 218 in generated FC links. ARR assumes that the fMRI data at time t is a linear combination of the fMRI data 219 from the previous p time points: p X xt = Al xt−l + t (9) l=1 –11–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 220 where p ≥ 1, xt is the N × 1 vector of fMRI data at time t, and  corresponds to zero-mean Gaussian 221 noise, and Al is an N × N matrix of model parameters which contains the linear dependencies between 222 each time t and its previous time point. ARR first estimates the model parameters for each time point 223 (A1 , ...Ap ) from the fMRI data. Each null fMRI time series is generated by randomly selecting p 224 successive time points from the original data, and then applying the AR model to generate T p new time 225 points until time series of length T are generated. Naturally, significant deviation of the original data Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 226 from ARR null data means null hypothesis being rejected. 227 The PR procedure initiates the null time series generation by performing Discrete Fourier Transform 228 (DFT) of each time course, and then adds a uniformly distributed random phase for each frequency, and 229 same random phase is added across all variables. Finally, an inverse DFT is performed to obtain the null 230 time series. PR generates data with linear, weak-sense stationarity (WSS), and Gaussian properties 231 whose auto-covariance sequence R0 , ..., RT −1 is similar to those of the original time series. A rejection of 232 the null hypothesis based on the two mentioned null models could be due to the fMRI time series not 233 possessing either one of the three properties of the null data or a combination of them. The experimental 234 results for WBN as well as the mentioned baseline approaches are provided in more detail in the next 235 section. EXPERIMENTAL RESULTS 236 In order to assess the proposed methodological framework, we apply it to a resting state fMRI data set of 237 300 subjects from the Autism Brain Imaging Data Exchange (ABIDE) database, including 150 subjects 238 diagnosed with Autism Spectrum Disorder (ASD)(Di Martino et al., 2014). This dataset was selected 239 from the C-PAC preprocessing pipeline and was slice time and motion corrected, and the voxel intensity 240 was normalized using global signal regression. The automated anatomical labeling atlas (AAL) was then 241 used for parcellation of regions of interest(Tzourio-Mazoyer et al., 2002). Then, the temporal links 242 between each pair of nodes were extracted based on the Pearson correlation between their BOLD 243 activation time series within each temporal window t, and were then rescaled based on min-max feature 244 scaling to have continuous values within the range [0, 1]. The implementation code for the methodology 245 in this work is available in 246 https://github.com/ThisIsNima/Weighted-Backbone-Network. –12–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 247 After extracting the backbone networks, we probed several aspects and measures of them which will 248 be discussed in this section. In particular, we provide a closer assessment of backbone networks on four 249 brain regions, namely the left and right hippocampus and the left and right amygdalas. We also provide 250 part of the experimental results for the cerebellar regions in the main manuscript and the rest in the 251 supplementary document. The reason for choosing these regions is the extensive focus of prior literature 252 related to diagnosis and pattern discovery in functional connectivity among ASD patients on them and Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 253 the fact that several types of abnormality have been discovered related to these regions among this group 254 of patients(Cooper et al., 2017; Guo et al., 2016; Ramos, Balardin, Sato, & Fujita, 2019; Rausch et al., 255 2016; Shen et al., 2016). 256 As the first step of our analysis, we examined the normality of the distribution of temporal links 257 between each pair of nodes i, j across the experiment time τ . For this purpose, we used the 258 Kolmogorov-Smirnov test on temporal ties between each pair of nodes for four different window sizes 259 ∆ ∈ {5, 10, 15, 20}. Table 1 demonstrates the average p values of the normality tests for the distribution 260 of temporal ties between every pair of voxels across 300 subjects for four separate regions. These results 261 demonstrate that the p values are below the 0.005 common threshold for rejecting the null hypothesis. 262 Furthermore, the p value tends to increase as the size of temporal windows decreases, which can be 263 attributed to the increase in total number of temporal windows τ . Beyond the theoretical basis of Central 264 Limit Theorem, these results further highlight that the assumption of normality for the distribution of 265 temporal edges in our resting state fRMI data is reasonable. 266 The backbone networks of the right hippocampus (region 38 per AAL atlas) for one control subject 267 based on four different threhsolds (c = 1 − α) are provided as heatmaps in figure 4, where each cell 268 represents a voxel, and white cells represents the significant ties. Note that self links are removed from 269 these networks, thus the value of the diagonals of the heatmaps are set to zero. For this analysis, time 270 courses were segmented into 20 temporal windows through the sliding window approach, with an overlap 271 of 5 time points between consecutive windows (this is the default setup for the other parts of the 272 experiments. Otherwise, we denote the temporal window size set up). The visualizations in figure 4 273 indicate that the number of admitted links decreases by increasing the threshold c. Moreover, the links 274 between voxels in the vicinity of the diagonal line tend to endure the increase in threshold c, which 275 highlights the strength of links between spatially close voxels. –13–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: Brain region ∆=5 ∆ = 10 ∆ = 15 ∆ = 20 L hippocampus 2.1004e−11 2.4801e−9 1.1422e−22 2.5488e−32 R hippocampus 8.1161e−13 2.5102e−10 1.1084e−9 1.1100e−9 L Amygdala 6.3835e−14 1.3560e−9 1.1609e−9 1.1102e−9 R Amygdala 3.3875e−10 5.0045e−7 1.4108e−7 3.0545e−7 Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 276 Table 1. p values of Kolmogorov-Smirnov test for normality of the distribution of temporal links. The p values presented in the table are averaged across all 277 links ((N (N − 1))/2 edges for N nodes) of the network. 278 The relation between the threshold c and number of significant ties is further inspected in figure 8-c in 279 SI, where the threshold increases from 0.5 to 1 with a fixed step resolution of 10−2 . The number of 280 significant ties for the network within each time window t = 1, ..., τ is also provided in figure 8-a and 281 8-b, where the red bars show the number of edges admitted to the final backbone network. As noted 282 earlier, only the ties that meet the significance threshold in over 50% of the time steps τ qualify to be 283 included in the final backbone network (red bar), thus the number of admitted links is usually smaller 284 than the significant ties within various temporal windows. However, as figure 8 in SI demonstrates, the 285 number of significant ties does not demonstrate a large variation across different temporal windows. 286 For the next step of the analysis, we examined and compared the backbone networks of the two cohorts 287 (control and ASD) within our experimental dataset with similar temporal segmentation as the previous 288 step. For this analysis, the value of α was set to 0.2, i.e. c-th percentile = 0.80. Figures 1 and 2 in the 289 supplementary information present the networks of significant ties extracted from the dynamic 290 connectivity of the left and right hippocampus from four subjects, including two subjects diagnosed with 291 ASD, and figures 3 and 4 in SI show the extracted significant ties from the left and right Amygdalas for 292 eight subjects, four of whom were diagnosed with ASD. Moreover, in order to provide a more 293 comprehensive perspective of the irreducible networks of the mentioned regions, the averaged backbone 294 networks of the two cohorts (Control and ASD) across the entire dataset are presented in figures 5 and 6 295 in the supplementary information. 296 As mentioned in the introduction, several null models have been applied to fMRI connectivity data in 297 the past based on null time series generation. Among these models, autoregressive randomization (ARR), –14–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 298 and phase randomization (PR) have been two of the most widely focused approaches. Therefore, we 299 compare the backbone networks based on those two methods with WBN (Handwerker et al., 2012; 300 Liegeois, Laumann, Snyder, Zhou, & Yeo, 2017; Liegeois, Yeo, & Van De Ville, 2021). A comparison of 301 the backbone network extracted through the WBN approach with ARR and PR null models is provided in 302 figure 7 in supplementary information where the averaged backbone networks of the two cohorts for the 303 right hippocampus are provided based on each null model. We can see that compared to the backbone Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 304 networks in figure 5, despite the fact that the backbone networks based on ARR and PR demonstrate a 305 higher density of weights around the diagonals, their averaged values are dispersed across the regions 306 with lower average values. This means that ARR and PR demonstrate a lower consistency of null 307 hypothesis rejection across the subjects in this study compared to WBN. Moreover, WBN demonstrates 308 higher accuracy in detection of randomly injected edges, which will be discussed in the next sections. An 309 explanation for these results can be the fact that in WBN the global and spatial information of the 310 network is considered in latent parameters of each node due to their dependency on the parameters of 311 every other node in the network, which can result in a more stable null model. Another reason can be the 312 fact that the resting state time courses of different regions can demonstrate variations in statistical 313 properties (Guan et al., 2020; Gultepe & He, 2013). Moreover, stationary linear Gaussian (SLG) models 314 might lack the ability to explain more complex aspects of fMRI dynamics. These issue can particularly 315 intensify in case studies with higher spatial resolution such as voxel-level analysis. 316 Furthermore, we assessed the effect of the length of temporal windows on the extracted significant ties. 317 For this purpose, we measured the difference between backbone networks of dFC based on four different 318 window sizes: ∆ ∈ {5, 10, 15, 20}, where the overlap between consecutive windows was 2 time points 319 for the smallest window (∆ = 5), and 5 time points for the other three window sizes. As the 320 measurement of dissimilarity, we used the mean percentage error (MPE) of the voxel-wise difference 321 (between the values of corresponding matrix cells) between the backbone networks averaged across 300 322 subjects. The results of this analysis are provided in Figure 5 for two threshold values of 0.5 and 0.9. As 323 this analysis demonstrates, the dissimilarity between the extracted backbone networks calculated as MPE 324 is negligibly small for both temporal resolutions, which indicates the consistency of the backbone 325 network against variations of the temporal window size. –15–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 326 In order to evaluate the effect of the choice of temporal segmentation criteria on the backbone 327 networks, we compare the networks based on sliding-window criteria as well as a change point detection 328 (DCR) approach for single-subject data (Cribben, Wager, & Lindquist, 2013). The DCR approach 329 proposed by Cribben et al. detects the data partitions with the smallest combined Bayesian Information 330 Criterion (BIC) score to obtain the candidate change points (Cribben et al., 2013). For this analysis, we 331 assigned the value of ∆ (the minimum possible number of time points between adjacent change points) Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 332 to be 10 time points. By comparing figures 5 (based on sliding window) and 9 (based on DCR) in the SI, 333 we can note an overall similar backbone structure with between the networks based on the two 334 segmentation approaches. 335 As the last part of the voxel-level experiments, we examined the correlation between the empirical 336 weight of the links and degree of the nodes in dynamic functional connectivity network with the 337 backbone link wights and estimated latent variables a, b. In figures 6 a and b, the average backbone 338 network of the right amygdala of 300 subjects as well as their average dFC over τ windows are presented. 339 Additionally, the correlations between node degrees of the dFC network, calculated as the sum of the 340 weights of temporal links for each node, and their estimated a, b as well as the correlation between the 341 average backbone link weights of 300 subjects and the average weight of their corresponding dFC links 342 over τ windows is provided in figure 6-c. Results for additional regions are provided in supplementary 343 information. As these results demonstrate, there is a weak correlation between the weight of the dFC 344 links and the average weight of backbone links (note that averaging binary backbone links results in 345 continuous weights). Additionally, there is a relatively small negative correlation between node degree of 346 dynamic functional connectivity and estimated distribution latent variables a, b. In line with the argument 347 provided in the methodology, these empirical results further illustrate that WBN considers global and 348 temporal information of the network beyond the local node degree and the weight of the links in the dFC. 349 Full brain Analysis 350 Just like the voxel-level analysis, a full brain backbone network of rsFC can be extracted where each 351 node is a region of interest (ROI). For this purpose, the time courses within each region based on the 352 AAL atlas was averaged. The averaged full-brain networks of significant ties for 150 control and 150 353 ASD subjects are demonstrated in figures 7, 8, 9, and 10, where the seed regions are the left and right –16–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 354 amygdalas and hippocampus, and the average backbone links with weights below 0.05 where filtered out 355 to facilitate easier presentation. In these figures, the link weights (illustrated by thickness in the figure) 356 correspond to the count of their corresponding links appearing in the binary backbone networks across 357 each cohort. We can observe from the figures 9, and 10 that the amygdalas and hippocmapus develop a 358 larger number of significant ties with other regions among the control group compared to the ASD cohort 359 as the network of the latter cohort is more sparse. The width of links in figures 7 and 8 also represent the Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 360 strength of the average correlations. We can therefore also observe a relatively high averaged backbone 361 link between the right and left hippocampus among both the control and ASD group. However, certain 362 differences can be detected between the two cohorts, including a stronger average backbone tie between 363 the hippocampus and the cerbellums as well as right and left olfactories among control subjects 364 compared to the ASD cohort. Moreover, a higher average backbone tie can be noticed between the left 365 and right amygdalas and the superior and the middle temporal gyrus among control subjects. Further 366 related experimental results are provided in the supplementary information, which include the average 367 backbone connectivity with several cerebellar regions being the seed area. These results demonstrate the 368 benefit of the weighted temporal backbone network in revealing the differences in irreducible ties 369 between different regions of interest. 370 Also, figure 11 depicts the averaged backbone connectivity of the cerebellums (18 regions per AAL) 371 and the vermis (8 regions per AAL) of the two cohorts in this study, which indicates higher connectivity 372 level among the control group compared to the ASD group. These results, along with the experimental 373 results provided in supplementary information (figures 5-10), can indicate that the increased 374 cerebro-cerebellar functional connectivity detected in some studies can be driven by a large number of 375 links that fail to be incompatible with the null hypothesis of links being produced at random. In other 376 words, despite the lower connectivity detected in cerebro-cerebellar subnetwork among the control group 377 in terms of number of links or their weights, the number of meaningful and irreducible links in that 378 subnetwork among the healthy cohort tend to be larger compared to the ASD cohort(Khan et al., 2015; 379 Mostofsky et al., 2009; Ramos et al., 2019). 380 Detection of random links –17–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 381 To compare the performance of the proposed approach with other proposed null models, we considered 382 three different models, namely binary ST filter, ARR, and PR. 383 we created a simulated dataset by injecting random weights to a subset of edges of the real rsFC 384 networks of our dataset. For this purpose, 100 random weights were injected into 100 links of rsFC of the 385 left hippocampus, and the precision of WBN as well as the ST filter approach, proposed by Kobayashi et Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 386 al., in excluding them from the final network were calculated(Kobayashi et al., 2019). Due to the fact that 387 the ST filter operates on temporal binary weight networks, in order to evaluate it we converted the rsFC 388 link weights as well as the randomly injected weights into binary links by drawing a temporal link 389 between each pair of node whose weight in the original rsFC network was above the entire 390 network'average. The result of this experiment is provided in figure 12, where WBN demonstrates an 391 advantage over the ST filter in random link detection precision. Similar experiment with other regions of 392 interest was conducted, which is provided in the supplementary information (figure 20 in SI). The 393 evaluation measure for this analysis were calculated by comparing the detection of injected random link 394 weights with the ground truth. Part of the superior performance of WBN can be attributed to the fact that 395 the process of conversion to binary network for the ST filter setup results in loss of information and 396 precision, which is an inherent disadvantage of backbone network detection approaches that are designed 397 for binary networks. DISCUSSION 410 In this work, we proposed a new approach for detecting the significant ties between nodes on voxel and 411 ROI level networks of resting state dFC. The proposed framework entails two computational steps; first, a 412 maximum likelihood optimization is performed to calculates the latent variables that characterize the 413 optimal Gaussian distribution of the temporal links between each pair of nodes across τ time steps. Then, 414 the empirical link weights between each pair of node within each temporal window are compared to the 415 c − th percentile of the Gaussian distribution to detect the significant links that form the backbone 416 network. This process is performed for every pair of nodes in the temporal network of dFC. Aside from 417 providing a systematic filtering framework for weighted temporal networks such as resting state dFC, this 418 approach has several analytical advantages over other prior filtering approaches that we discuss in this –18–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 398 Figure 4. Derived backbone networks of the right hippocampus from one control subject given four threshold values. 399 Figure 5. A comparison based on different window sizes using mean percentage error (MPE) of the voxel-wise difference between the backbone networks 400 of dFC averaged across 300 subjects –19–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 401 Figure 6. Figure a: average backbone network of the right Amygdala for 300 subjects. Figure b: average dFC network of the same region for the same Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 402 subjects across τ = 20 time intervals. Figure c: correlations between average node degree and estimated a, b backbone link weights as well as average dFC 403 link weight over τ = 20 intervals for 300 subjects. –20–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 404 Figure 7. The averaged full brain backbone networks of 150 control subjects based on four seed regions –21–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 405 Figure 8. The averaged full brain backbone networks of 150 ASD subjects based on four seed regions –22–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 406 Figure 9. Averaged full brain backbone networks of 150 ASD and 150 control subjects with the left and right Amygdalas as the seed regions –23–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 407 Figure 10. Averaged full brain backbone networks of 150 ASD and 150 control subjects with the left and right hippocampus as the seed regions –24–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 408 Figure 11. Averaged full brain backbone networks of cerebro-cerebllar and vermis regions across 150 ASD and 150 control subjects 409 Figure 12. The AUC of detection of injected random weights based on the ST filtering, ARR, PR and WBN in the left hippocampus. –25–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 419 section. We also discuss the limitations of the proposed methodology along with possible suggestions for 420 improvement and future plans. 421 As mentioned previously, inclusion of a temporal link in the backbone network is determined by 422 testing the hypothesis that the link can be explained by the null model that links are created uniformly at 423 random. This comparison is applied to every link in the dFC individually, i.e. between every pair of Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 424 nodes and within every temporal window. Therefore, temporal properties and variations of the network 425 structure over time are taken into account in backbone network inference. This property is an advantage 426 of the proposed methodology over some of the prior approaches that consider a constant intrinsic activity 427 value for the nodes over time. It also offers the power of determining a cut-off percentage of ties having a 428 larger weight over the c − th percentile, which was decided to be %50 in this study. 429 Another advantage of the suggested approach is the fact that it considers the interplay of global and 430 local information of the network in estimating the latent variables a and b. In other words, the 431 significance of temporal ties cannot be attributed merely to node properties such as degree or centrality 432 measures, because each equation in the system of N equations of equations 4 and 5 takes into account the 433 combination of weights over time for each link for node i as well as their combination with other links 434 between i and every other node in the network. This property has been discussed in more detail in 435 methodology section and evaluated in the results section. 436 The refinement of parameters of the distribution through maximum likelihood optimization requires 437 solving the system of N equations for N nodes (one set of N equations for each of the two parameters), 438 which can be solved through several optimization approaches such as gradient-based optimization, search 439 methods, or the Newton method. Solving these equations does not require any hyper parameter tuning as 440 the only parameters that need to be selected as input is the threshold value α and the percentage of times 441 that the weight of the link meets the c = 1 − α percentile of the distribution, which offer the flexibility for 442 a comprehensive assessment of the temporal ties in the dynamic connectivity network. 443 Unlike some of the null models suggested in the past for binary networks based on binomial or Poisson 444 distributions, the methodology put forward in this work does not assume a strictly positive weight 445 between interacting nodes. This property provides the flexibility for ties that are generated through –26–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 446 various approaches such as correlation measures to be considered in the null model, as negative 447 correlation is a possibility between interacting nodes. 448 Another advantage of the proposed approach is the fact that the backbone networks are learned for 449 each subject individually. As explained in the methodology section, the input for WBN is the weighted 450 dynamic connectivity network of a subject, and its output is the network of irreducible ties corresponding to the subject. This property has the benefit of taking into account the individual differences when Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 451 452 inferring the backbone network in an isolated fashion. 453 The suggested methodological framework can be used in studies with various scales and resolution of 454 dFC networks, meaning that instead of voxel-level analysis, dFC networks consisting of regions of 455 various scales as nodes can benefit from this approach as well. Moreover, this approach is independent of 456 temporal segmentation step, as long as the statistical properties of independence and normality are met. 457 Limitations 458 Despite the mentioned advantages, the proposed approach bears certain limitations which we highlight in 459 this section. 460 As discussed in the methodology, the first step of the suggested framework entails estimation of latent 461 variables a and b, which rule the propensities to generate a distribution of links with a certain average and 462 standard deviation. However, these variables are estimated and compared across the experiment time τ , 463 i.e. the length of the fMRI signals. In other words, the mean and standard deviation of the distributions, 464 and in turn the backbone network calculations, can vary depending on the length of the experiment. 465 The structural characteristics of dFC can be influenced by temporal fluctuations in the data throughout 466 the course of the experiment. In other words, reducible links might not have the same statistical features 467 at any time during the observation as node properties might not be constant over time. Therefore, more 468 improvements need to be applied to WBN to take such variations into account. 469 Despite its adaptability with regards to different temporal segmentation approaches such as sliding 470 window or DCR, WBN requires equal number of temporal windows across the entire region of interest 471 for calculating the latent distribution variables for each node due to the number of optimization equations 472 that it solves. –27–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 473 Another limitation of the suggested approach is the assumption of normality for larger temporal 474 window sizes. As the empirical tests demonstrated, an increase in size of the temporal windows could in 475 principle weaken the normality assumption of the distribution of the temporal links. Despite the evidence 476 of normality for reasonable and common window sizes in the literature, this assumption needs to be 477 further explored for various different datasets. The MLE optimization for estimating the intrinsic variables a, b plays the largest role in the Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 478 479 computational complexity of the methodology presented in this work. The computation time depends on 480 the number of nodes, i.e. spatial resolution, and the number of time intervals that the signal is segmented 481 to. By definition of the approach, the spatial resolution plays a more significant role in the computational 482 complexity (refer to equations 4 and 5). In this study, the system of N equations were solved through the 483 trust-region-dogleg method, whose computation time for regions below 1000 voxels was 10 minutes for 484 8 GB of RAM memory. However, more efficient approaches can be employed for this purpose. 485 Alleviating the mentioned limitations requires further methodological explorations and analytical 486 studies on various datasets. As future work, our objective will include assessment of the backbone 487 network of resting state dFC of other cohorts and data from various neurological conditions and to study 488 different group differences. Furthermore, assessment of significant temporal structures and graph 489 communities and motifs as well as exploring the effect of different preprocessing pipelines and temporal 490 sample size on the outcome of the proposed approach can be fruitful paths for further experiments in the 491 area of dynamic functional connectivity. ACKNOWLEDGMENTS 492 This work was supported by National Institute of Health grants to I. R. Olson [R01HD099165; RO1 493 MH091113; R21 HD098509; and 2R56MH091113-11]. SUPPORTING INFORMATION 494 Supporting information (also referred to as SI in this manuscript) document contains the figures from 495 several steps and aspects of experiments which are referred to in this manuscript. Moreover, extra 496 experiments on different regions of the brain are provided in SI. –28–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: COMPETING INTERESTS 497 The authors declare no competing interests 498 499 REFERENCES 500 Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 501 Allen, E. A., Damaraju, E., Plis, S. M., Erhardt, E. B., Eichele, T., & Calhoun, V. D. (2014). Tracking whole-brain 502 connectivity dynamics in the resting state. Cerebral cortex, 24(3), 663–676. 503 Alvarez-Hamelin, J. I., Dall’Asta, L., Barrat, A., & Vespignani, A. (2005). K-core decomposition of internet graphs: 504 hierarchies, self-similarity and measurement biases. arXiv preprint cs/0511007. 505 Beall, E. B., & Lowe, M. J. (2007). Isolating physiologic noise sources with independently determined spatial measures. 506 Neuroimage, 37(4), 1286–1300. 507 Behzadi, Y., Restom, K., Liau, J., & Liu, T. T. (2007). A component based noise correction method (compcor) for bold and 508 perfusion based fmri. Neuroimage, 37(1), 90–101. 509 Birn, R. M. (2012). The role of physiological noise in resting-state functional connectivity. Neuroimage, 62(2), 864–870. 510 Birn, R. M., Smith, M. A., Jones, T. B., & Bandettini, P. A. (2008). The respiration response function: the temporal 511 dynamics of fmri signal fluctuations related to changes in respiration. Neuroimage, 40(2), 644–654. 512 Bollmann, S., Kasper, L., Vannesjo, S. J., Diaconescu, A. O., Dietrich, B. E., Gross, S., . . . Pruessmann, K. P. (2017). 513 Analysis and correction of field fluctuations in fmri data using field monitoring. Neuroimage, 154, 92–105. 514 Bright, M., & Murphy, K. (2017). Cleaning up the fmri time series: Mitigating noise with advanced acquisition and 515 correction strategies. NeuroImage, 154, 1–3. 516 Casiraghi, G., Nanumyan, V., Scholtes, I., & Schweitzer, F. (2017). From relational data to graphs: Inferring significant 517 links using generalized hypergeometric ensembles. In International conference on social informatics (pp. 111–120). 518 Chang, C., & Glover, G. H. (2009). Relationship between respiration, end-tidal co2, and bold signals in resting-state fmri. 519 Neuroimage, 47(4), 1381–1393. –29–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 520 Chang, C., & Glover, G. H. (2010). Time–frequency dynamics of resting-state brain connectivity measured with fmri. 521 Neuroimage, 50(1), 81–98. 522 Choe, A. S., Nebel, M. B., Barber, A. D., Cohen, J. R., Xu, Y., Pekar, J. J., . . . Lindquist, M. A. (2017). Comparing 523 test-retest reliability of dynamic functional connectivity methods. Neuroimage, 158, 155–175. 524 Cimini, G., Squartini, T., Saracco, F., Garlaschelli, D., Gabrielli, A., & Caldarelli, G. (2019). The statistical physics of Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 525 real-world networks. Nature Reviews Physics, 1(1), 58–71. 526 Cooper, R. A., Richter, F. R., Bays, P. M., Plaisted-Grant, K. C., Baron-Cohen, S., & Simons, J. S. (2017). Reduced 527 hippocampal functional connectivity during episodic memory retrieval in autism. Cerebral Cortex, 27(2), 888–902. 528 Cribben, I., Wager, T., & Lindquist, M. (2013). Detecting functional connectivity change points for single-subject fmri data. 529 Frontiers in computational neuroscience, 7, 143. 530 Damaraju, E., Allen, E. A., Belger, A., Ford, J. M., McEwen, S., Mathalon, D., . . . others (2014). Dynamic functional 531 connectivity analysis reveals transient states of dysconnectivity in schizophrenia. NeuroImage: Clinical, 5, 298–308. 532 Di Martino, A., Yan, C.-G., Li, Q., Denio, E., Castellanos, F. X., Alaerts, K., . . . others (2014). The autism brain imaging 533 data exchange: towards a large-scale evaluation of the intrinsic brain architecture in autism. Molecular psychiatry, 19(6), 534 659–667. 535 Dudley, R. M. (1978). Central limit theorems for empirical measures. The Annals of Probability, 899–929. 536 Dudley, R. M. (2014). Uniform central limit theorems (Vol. 142). Cambridge university press. 537 Gemmetto, V., Cardillo, A., & Garlaschelli, D. (2017). Irreducible network backbones: unbiased graph filtering via 538 maximum entropy. arXiv preprint arXiv:1706.00230. 539 Glover, G. H., & Lee, A. T. (1995). Motion artifacts in fmri: comparison of 2dft with pr and spiral scan methods. Magnetic 540 resonance in medicine, 33(5), 624–635. 541 Gorgolewski, K. J., Storkey, A. J., Bastin, M. E., Whittle, I., & Pernet, C. (2013). Single subject fmri test–retest reliability 542 metrics and confounding factors. Neuroimage, 69, 231–243. 543 Guan, S., Jiang, R., Bian, H., Yuan, J., Xu, P., Meng, C., & Biswal, B. (2020). The profiles of non-stationarity and 544 non-linearity in the time series of resting-state brain networks. Frontiers in Neuroscience, 14, 493. –30–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 545 Gultepe, E., & He, B. (2013). A linear/nonlinear characterization of resting state brain networks in fmri time series. Brain 546 topography, 26(1), 39–49. 547 Guo, X., Duan, X., Long, Z., Chen, H., Wang, Y., Zheng, J., . . . Chen, H. (2016). Decreased amygdala functional 548 connectivity in adolescents with autism: a resting-state fmri study. Psychiatry Research: Neuroimaging, 257, 47–56. 549 Haller, S., & Bartsch, A. J. (2009). Pitfalls in fmri. European radiology, 19(11), 2689–2706. Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 550 Hallquist, M. N., Hwang, K., & Luna, B. (2013). The nuisance of nuisance regression: spectral misspecification in a 551 common approach to resting-state fmri preprocessing reintroduces noise and obscures functional connectivity. 552 Neuroimage, 82, 208–225. 553 Handwerker, D. A., Roopchansingh, V., Gonzalez-Castillo, J., & Bandettini, P. A. (2012). Periodic changes in fmri 554 connectivity. Neuroimage, 63(3), 1712–1719. 555 Heitmann, S., & Breakspear, M. (2018). Putting the “dynamic” back into dynamic functional connectivity. Network 556 Neuroscience, 2(02), 150–174. 557 Hindriks, R., Adhikari, M. H., Murayama, Y., Ganzetti, M., Mantini, D., Logothetis, N. K., & Deco, G. (2016). Can 558 sliding-window correlations reveal dynamic functional connectivity in resting-state fmri? Neuroimage, 127, 242–256. 559 Hutchison, R. M., Womelsdorf, T., Allen, E. A., Bandettini, P. A., Calhoun, V. D., Corbetta, M., . . . others (2013). 560 Dynamic functional connectivity: promise, issues, and interpretations. Neuroimage, 80, 360–378. 561 Jones, D. T., Vemuri, P., Murphy, M. C., Gunter, J. L., Senjem, M. L., Machulda, M. M., . . . others (2012). Non-stationarity 562 in the “resting brain’s” modular architecture. PloS one, 7(6), e39731. 563 Kalthoff, D., Seehafer, J. U., Po, C., Wiedermann, D., & Hoehn, M. (2011). Functional connectivity in the rat at 11.7 t: 564 Impact of physiological noise in resting state fmri. Neuroimage, 54(4), 2828–2839. 565 Khan, A. J., Nair, A., Keown, C. L., Datko, M. C., Lincoln, A. J., & Müller, R.-A. (2015). Cerebro-cerebellar resting-state 566 functional connectivity in children and adolescents with autism spectrum disorder. Biological psychiatry, 78(9), 625–634. 567 Kiviniemi, V., Vire, T., Remes, J., Elseoud, A. A., Starck, T., Tervonen, O., & Nikkinen, J. (2011). A sliding time-window 568 ica reveals spatial variability of the default mode network in time. Brain connectivity, 1(4), 339–347. –31–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 569 Kobayashi, T., Takaguchi, T., & Barrat, A. (2019). The structured backbone of temporal social ties. Nature 570 communications, 10(1), 1–11. 571 Kundu, P., Inati, S. J., Evans, J. W., Luh, W.-M., & Bandettini, P. A. (2012). Differentiating bold and non-bold signals in 572 fmri time series using multi-echo epi. Neuroimage, 60(3), 1759–1770. 573 Li, M.-X., Palchykov, V., Jiang, Z.-Q., Kaski, K., Kertész, J., Miccichè, S., . . . Mantegna, R. N. (2014). Statistically Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 574 validated mobile communication networks: the evolution of motifs in european and chinese data. New Journal of Physics, 575 16(8), 083038. 576 Liegeois, R., Laumann, T. O., Snyder, A. Z., Zhou, J., & Yeo, B. T. (2017). Interpreting temporal fluctuations in 577 resting-state functional connectivity mri. Neuroimage, 163, 437–455. 578 Liegeois, R., Yeo, T. B., & Van De Ville, D. (2021). Interpreting null models of resting-state functional mri. bioRxiv. 579 Liu, S., Perra, N., Karsai, M., & Vespignani, A. (2014). Controlling contagion processes in activity driven networks. 580 Physical review letters, 112(11), 118702. 581 Lydon-Staley, D. M., Ciric, R., Satterthwaite, T. D., & Bassett, D. S. (2019). Evaluation of confound regression strategies 582 for the mitigation of micromovement artifact in studies of dynamic resting-state functional connectivity and multilayer 583 network modularity. Network Neuroscience, 3(2), 427–454. 584 Ma, L.-l., Ma, C., Zhang, H.-F., & Wang, B.-H. (2016). Identifying influential spreaders in complex networks based on 585 gravity formula. Physica A: Statistical Mechanics and its Applications, 451, 205–212. 586 Mostofsky, S. H., Powell, S. K., Simmonds, D. J., Goldberg, M. C., Caffo, B., & Pekar, J. J. (2009). Decreased connectivity 587 and cerebellar activity in autism during motor task performance. Brain, 132(9), 2413–2425. 588 Murphy, K., Birn, R. M., & Bandettini, P. A. (2013). Resting-state fmri confounds and cleanup. Neuroimage, 80, 349–359. 589 Nadini, M., Bongiorno, C., Rizzo, A., & Porfiri, M. (2020). Detecting network backbones against time variations in node 590 properties. Nonlinear Dynamics, 99(1), 855–878. 591 Patanaik, A., Tandi, J., Ong, J. L., Wang, C., Zhou, J., & Chee, M. W. (2018). Dynamic functional connectivity and its 592 behavioral correlates beyond vigilance. NeuroImage, 177, 1–10. –32–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 593 Perra, N., Gonçalves, B., Pastor-Satorras, R., & Vespignani, A. (2012). Activity driven modeling of time varying networks. 594 Scientific reports, 2, 469. 595 Ramos, T. C., Balardin, J. B., Sato, J. R., & Fujita, A. (2019). Abnormal cortico-cerebellar functional connectivity in autism 596 spectrum disorder. Frontiers in systems neuroscience, 12, 74. 597 Rausch, A., Zhang, W., Haak, K. V., Mennes, M., Hermans, E. J., van Oort, E., . . . Groen, W. B. (2016). Altered functional Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 598 connectivity of the amygdaloid input nuclei in adolescents and young adults with autism spectrum disorder: a resting state 599 fmri study. Molecular autism, 7(1), 1–13. 600 Rizzo, A., Frasca, M., & Porfiri, M. (2014). Effect of individual behavior on epidemic spreading in activity-driven 601 networks. Physical Review E, 90(4), 042801. 602 Savva, A. D., Kassinopoulos, M., Smyrnis, N., Matsopoulos, G. K., & Mitsis, G. D. (2020). Effects of motion related 603 outliers in dynamic functional connectivity using the sliding window method. Journal of Neuroscience Methods, 330, 604 108519. 605 Serrano, M. Á., Boguná, M., & Vespignani, A. (2009). Extracting the multiscale backbone of complex weighted networks. 606 Proceedings of the national academy of sciences, 106(16), 6483–6488. 607 Shen, M. D., Li, D. D., Keown, C. L., Lee, A., Johnson, R. T., Angkustsiri, K., . . . Nordahl, C. W. (2016). Functional 608 connectivity of the amygdala is disrupted in preschool-aged children with autism spectrum disorder. Journal of the 609 American Academy of Child & Adolescent Psychiatry, 55(9), 817–824. 610 Shmueli, K., van Gelderen, P., de Zwart, J. A., Horovitz, S. G., Fukunaga, M., Jansma, J. M., & Duyn, J. H. (2007). 611 Low-frequency fluctuations in the cardiac rate as a source of variance in the resting-state fmri bold signal. Neuroimage, 612 38(2), 306–320. 613 Smith, S. M. (2012). The future of fmri connectivity. Neuroimage, 62(2), 1257–1266. 614 Smith, S. M., Miller, K. L., Moeller, S., Xu, J., Auerbach, E. J., Woolrich, M. W., . . . others (2012). Temporally-independent 615 functional modes of spontaneous brain activity. Proceedings of the National Academy of Sciences, 109(8), 3131–3136. 616 Starnini, M., & Pastor-Satorras, R. (2014). Temporal percolation in activity-driven networks. Physical Review E, 89(3), 617 032807. –33–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 618 Tumminello, M., Micciche, S., Lillo, F., Piilo, J., & Mantegna, R. N. (2011). Statistically validated networks in bipartite 619 complex systems. PloS one, 6(3), e17994. 620 Tzourio-Mazoyer, N., Landeau, B., Papathanassiou, D., Crivello, F., Etard, O., Delcroix, N., . . . Joliot, M. (2002). 621 Automated anatomical labeling of activations in spm using a macroscopic anatomical parcellation of the mni mri 622 single-subject brain. Neuroimage, 15(1), 273–289. Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 623 Van Dijk, K. R., Hedden, T., Venkataraman, A., Evans, K. C., Lazar, S. W., & Buckner, R. L. (2010). Intrinsic functional 624 connectivity as a tool for human connectomics: theory, properties, and optimization. Journal of neurophysiology, 103(1), 625 297–321. 626 Van Dijk, K. R., Sabuncu, M. R., & Buckner, R. L. (2012). The influence of head motion on intrinsic functional 627 connectivity mri. Neuroimage, 59(1), 431–438. 628 von der Malsburg, C. E., Phillps, W. A., & Singer, W. E. (2010). Dynamic coordination in the brain: From neurons to mind. 629 MIT Press. 630 Yan, X., Jeub, L. G., Flammini, A., Radicchi, F., & Fortunato, S. (2018). Weight thresholding on complex networks. 631 Physical Review E, 98(4), 042304. 632 Zalesky, A., Fornito, A., Cocchi, L., Gollo, L. L., & Breakspear, M. (2014). Time-resolved resting-state brain networks. 633 Proceedings of the National Academy of Sciences, 111(28), 10341–10346. 634 Zino, L., Rizzo, A., & Porfiri, M. (2017). An analytical framework for the study of epidemic models on activity driven 635 networks. Journal of Complex Networks, 5(6), 924–952. TECHNICAL TERMS 636 Temporal segmentation The process of slicing the fMRI time courses into consecutive temporal 637 windows within which connectivity matrices are formed based on the correlation between the fMRI time 638 series. A common temporal segmentation is the sliding window approach that computes a succession of 639 pairwise correlation matrices using the time courses from a given parcellation of brain regions. 640 Backbone network The network that is formed by the significant ties between nodes by filtering out 641 randomly generated edges or noise-induced links. –34–  == D R A F T == Journal: NETWORK NEUROSCIENCE / Title: The Backbone Network of Dynamic Functional Connectivity Authors: 642 Null model Null models are formulated as a baseline for comparison with the system to verify whether 643 the system displays properties that would not be expected on a random basis or as a consequence of 644 certain constraints. 645 Newton method for solving system of nonlinear equations An approach for solving a system of 646 nonlinear equations by finding the roots of differentiable functions. Downloaded from http://direct.mit.edu/netn/article-pdf/doi/10.1162/netn_a_00209/1962697/netn_a_00209.pdf by guest on 15 September 2021 647 Temporal ties The link between a pair of nodes whose weight might vary across the time of 648 experiment. 649 Latent model variables Model variables that are not directly observed or assigned, but are inferred 650 from other measurements from the data or variables that are observed. 651 weak-sense stationary A random process whose mean function and its autocovariance function do not 652 fluctuate by variations in time. –35–