Event-marked Windowed Communication: Inferring activity propagation from neural time series

Neural communication, Functional connectivity, MEG, connectome 15 16

17 Tracking signal propagation in nervous systems is crucial to our understanding of brain 18 function and information processing. Current methods for inferring neural communication 19 track patterns of sustained co-activation over time, making them unsuitable to detect discrete 20 instances of signal transmission. Here, we propose Event-marked Windowed Communication 21 (EWC), a new analytical framework to infer functional interactions arising from discrete 22 signalling events between neural elements, in otherwise continuous time series data. In 23 contrast to conventional measures of functional connectivity, our method utilises an event-24 based subsampling of neural time series, which allows it to capture the statistical analogue 25 of activity propagation. We test EWC on simulations of neural dynamics and show that it is 26 capable of retrieving ground truth motifs of directional signalling , over a range of model 27 configurations. Critically, we demonstrate that EWC’s subsampling approach affords profound 28 reductions in computation times, compared to established network inference methods such 29 as transfer entropy. Lastly, we showcase the utility of EWC to infer whole -brain functional 30 networks from MEG recordings. Networks computed using EWC and transfer entropy were 31 highly correlated (median r=0.821 across subjects), but EWC inference was approximately 6.5 32 times faster per epoch. In summary, our work presents a new method to infer signalling from 33 time series of neural activity at low computational costs. Our framework is flexible and can be 34 applied to activity time series captured by diverse functional neuroimaging modalities, 35 opening up new avenues for the study of neural communication. 36 37 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 2

38 Communication between neural elements — neurons, neural populations and grey matter 39 regions, plays a central role in the functioning of the brain. Understanding the principles by 40 which signals are dynamically and flexibly transmitted in networked nervous systems remains 41 an open challenge and active area of research 1–3. Efforts in this direction have produced a 42 vast number of models of brain communication 1,4,5, ranging from information transfer via 43 neural oscillations6–10 to network measures of connectome communication 11–13. Despite the 44 abundance in models, efforts in model validation have been lacking, and it remains unclear 45 which methods faithfully describe empirical patterns of neural dynamics and signalling. A key 46 barrier to progress in model validation is our current inability to infer events of signal 47 transmission from recordings of neural activity. 48 49 Functional connectivity (FC) is commonly used as a proxy for neural communication. Broadly 50 defined, FC quantifies statistical dependencies in time series of neural activity using measures 51 such as the Pearson correlation or Mutual Information , capturing the extent to which the 52 dynamics of two brain regions are synchronised over a period of time. FC estimation has been 53 extensively used in model formulation and validation as well as in experimental and clinical 54 studies3,14–24. While useful, these FC measures are limited to capturing stable, sustained 55 statistical associations over time, which can dilute or mask discrete events of directional signal 56 transmission from one neural element to another. Information -theoretical measures such as 57 Transfer Entropy stand to address this issue, but are typically data - and computation-heavy 58 for whole-brain network inference25. 59 60 Explicitly tracking how activity (endogenous/exogenous perturbations) or signals propagate 61 across the underlying anatomical substrate can provide important insights into the principles 62 of inter-areal communication and help validate existing computational models . Current 63

in neuroimaging technology have seen activity propagation tracing often applied 64 to smaller controlled stimulation experiments in microscale neuronal networks 26–28. Its 65 application at the w hole-brain scale is still limited 29–32, although approaches combining 66 existing data to yield robust insights seem to show promise33. In this work, our aim was to 67 develop an efficient method of inferring activity propagation patterns by combining the 68 practicality and efficiency of undirected FC estimation with the dynamical resolution afforded 69 by activity propagation tracing. We term this modified FC estimation protocol Event-marked 70 Windowed Communication (EWC ). Specifically, EWC involves identif ying salient events in 71 neural recordings (“stimuli”) and the subsequent estimation of FC within short temporally 72 ordered windows/subsamples of the signal. In this manner, EWC captures an indirect 73 statistical analogue of activity propagation as opposed to spatiotemporally localised 74 stimulation-response measurements. We first validate and demonstrate the utility of the EWC 75 implementation using a simple in-silico network motif with ground truth signalling embedded 76 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 3 in its dynamics. We also show that the EWC approach allows us to use undirected FC 77 measures to capture directional interactions. We then study the computational tractability of 78 the approach by comparing network inference runtimes as a function of number of nodes. 79 Finally, we demonstrate a real-world application of the method on source -localised resting-80 state magnetoencephalography (MEG) recordings. 81 82

83 84 Figure 1 – A pipeline to es1mate communica1on pa5erns from regional 1me series. (A) Regional ac1vi1es are divided into 10 85 second epochs, and each epoch is processed independently. (B) For each epoch, the regional ac1vi1es are z-scored to iden1fy 86 instances of devia1on from mean behaviour, which are termed “supra-threshold events”. (C) Each region is successively chosen 87 as a “source”. For each event of a source, a communica1on window that encompasses the ac1vity over a dura1on of 1 second, 88 is defined star1ng at the 1me point of the event. Similar communica1on windows are defined over the ac1vi1es of all other 89 regions (termed “targets”), star1ng at 1mepoints delayed with respect to the source event, in propor1on to the Euclidean 90 distance between the source and target. (D) For each event for a chosen source, the Event-marked Windowed Communica1on 91 (EWC) is es1mated between the source and all possible targets by compu1ng the condi1onal mutual informa1on between 92 the ac1vi1es contained in the communica1on windows of the source and target. The ac1vity of the target over 1 second in 93 the past is used as the condi1oning variable . The pairwise EWC values over all events of a source (within an epoch) are 94 averaged, to populate a row of the EWC matrix corresponding to the source index. Steps B-D are repeated over all epochs to 95 give epoch-level EWC matrices that capture the dynamic communica1on pa5erns over the en1re scan dura1on. 96 In this study, we develop a new method to quantify functional interactions between neural 97 elements based on the detection of directional signalling events in their activity time series 98 (Fig.1). The Event-marked Windowed Communication (EWC) protocol is composed of 3 main 99 elements: 1) Event -identification 2) Temporal ordering based on conduction delays, and 3) 100 Windowing or subsampling. Specifically, our approach involves first identifying significant 101 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 4 deviations from the mean of activity time series, which we term “supra-threshold events” (For 102 brevity, we will refer to them as “events” for the remainder of this paper). We then trace the 103 downstream effects of these events on the activities of other neural elements , while 104 considering the possible delays in signal conduction due to inter -element distances. 105 Importantly, we restrict the estimation of communication to a short “communication window”, 106 that starts at the timepoint of an event, for the source, and at a future timepoint for the target 107 (proportional to the source-target distance). Once a communication window is specified, the 108 magnitude of the signalling event is quantified by the statistical association between the 109 within-window activity time series of the source and target . EWC is flexible and this 110 association can be computed using various measures of FC, such as the partial correlation or 111 Mutual Information. Lastly, we derive an aggregate communication measure by averaging the 112 FC across communication windows. 113 114 The analyses of this paper can be divided into three main sections. First, we tested our 115

in simulated time series produced by a simple connectivity motif in which node 116 activity was governed by Linear Stochastic Model (LSM) dynamics (Fig 2A). The simple 117 dynamical landscape of the LSM allowed us to impose ground truth signalling between 118 nodes. We assessed our method’s ability to identify these communication patterns in relation 119 to previously proposed measures of FC. Second, we characterise gains in computational time 120 afforded by estimating FC using EWC. Third, we apply our method on source localised MEG 121 recordings, and gauge the agreement between the FC inferred v ia EWC and bivariate 122 Transfer Entropy. 123 In the remainder of the paper, we will use the suffix “-Full” after the FC measure to indicate 124 a conventional estimation technique using the entire signal, and “ -EWC” to indicate our 125 implementation. 126 127 Asymmetric signalling over a network mo tif 128 We tested the EWC protocol in a four-node motif with Linear Stochastic Model (LSM) 129 dynamics (Fig.3). Three of the nodes in this network were connected in a linear topology 130 (nodes 1, 2, and 3), whereas node 4 was isolated (i.e. not influenced by the dynamics of other 131 nodes). In addition to the LSM dynamics, the nodes 1, 3, and 4 also pulsated at random times 132 according to a Poisson process. These pulses emulated punctual events of directional 133 signalling, with known sources at nodes 1, 3, and 4. 134 135 EWC was implemented using two measures of within -window FC, partial correlation (PC -136 EWC) and conditional mutual information (cMI -EWC). Note that, traditionally, these are 137 symmetric measures that cannot resolve directionality of functional interactions. We 138 benchmarked EWC against Transfer Entropy (TE -Full), which was conventionally estimated 139 based on the entire time series. We compared the ability of these measures to retrieve the 140 ground truth signalling motif in simulations with increasing noise levels. To enable comparison 141 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 5 of the different measures, we estimated the standardised contrast between the FC in true 142 (existing) and false (absent) connections. 143 144 Figure 2 - Communica1on over a network mo1f. (A) (inset) We test our method of es1ma1ng communica1on pa5erns on a 145 simple 3-node network mo1f with a linear topology and an isolated node (4 -cyan). The ac1vi1es of the individual nodes in 146 this network are described by a Linear Stochas1c Model. The dynamics of the red, blue and cyan nodes have an addi1onal 147 Poisson process, causing it to spike at an average rate of 0.2 Hz, emula1ng communica1on events. 𝛿xy is the delay between 148 nodes x and y, in ms. The noise amplitude of the LSM, s , is varied rela1ve to the fixed Poisson pulse amplitude of the sources. 149 (A) Transfer entropy (in nats) between nodes, es1mated over the full signal. (B) Condi1onal Mutual Informa1on (in nats) 150 between nodes, es1mated as per the EWC protocol (C) Par1al Correla1on (Pearson-R) between nodes, es1mated as per EWC. 151 We show only the absolute Pearson correla1on. (D -F) Standardised contrast (difference between EWC of True and False 152 connec1ons, normalised by their pooled standard devia1on) plots associated with the FC es1mates. A non-zero contrast value 153 indicates discriminability. Shading in the plots represent ±SEM (20 trials). All plots maximally smoothed (i.e. using all available 154 points) to clearly reveal trends. 155 156 We considered the performance of our method as a function of the ratio between the noise 157 level in the system and the amplitude at which the source regions pulsate. The delay between 158 the sources (nodes 1, 3, and 4 ) and the target ( node 2) was set at 15ms, and the sources 159 pulsated according to a Poisson process with a mean frequency of 0.2Hz (Fig.2A - inset). The 160 ground truth comprises connections from node 1 to node 2 and 3 to 2. We observed that 161 estimating the TE-Full resulted in a good representation of the ground truth, with strong 162 interactions from nodes 1 and 3 to 2 (Fig.2A). Importantly, the corresponding contrast 163 distribution indicated that the interaction between node 4 (isolated) and 2 could be clearly 164 discriminated from the true interactions , providing evidence of good specificity (Fig.2D). 165 When the FC was instead measured using cMI-EWC, we observe d an interaction between 166 regions 4 and 2 , although this false positive interaction was considerably weaker than the 167 ground truth interactions from node 1 to 2, and 3 to 2, particularly at low noise levels (Fig.2B). 168 Interestingly, PC-EWC performed well when compared to cMI -EWC. Like TE-Full, PC-EWC 169 also resulted in a good representation of the ground truth (Fig.2C). The contrast between 170 estimated interactions show ed that true connections could be distinguished from absent 171 Partial Correlation (EWC) 1 3 24 42 𝜎 𝜈 = 0.2Hz = 15ms42== A B C D E F .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 6 connections even at large levels of noise (Fig.2F), although not as well as TE-Full. Importantly, 172 in the case of cMI -EWC and PC-EWC, asymmetric interactions were estimated, despite the 173 FC measures being undirected/symmetric. This asymmetry is a result of the temporal ordering 174 within the EWC implementation. 175 176 Additionally, we also tested whether a common source equidistant from two targets would 177

in a spurious functional link between them (a closed triangle problem). For this we set 178 region 2 as the pulsating source and estimated the PC-EWC. We observed that our method 179 correctly identifie d communication from node 2 to both 1 and 3, with negligible 180 communication in the opposite direction, or between the targets (Fig.S1). We also tested the 181 EWC implementation over a range of conduction delays, and mean firing rates of the sources 182 (Fig.S2). 183 184 To summarise, our comparative analyses demonstrate that EWC can accurately capture 185 directional interactions over a range of noise levels, delays, and firing rates. We emphasize 186 that PC and cMI, being undirected/symmetric measures, were only able to capture directional 187 signalling due to the temporal ordering implemented by the EWC protocol. Despite these 188 positive results, based on our observations of the standardised contrast, we find that TE-Full 189 still provides the strongest discriminability between the absence and presence of signall ing 190 patterns. For this reason, we continue using TE-Full as a representative benchmark of directed 191 FC in the subsequent analyses. 192 193 Computational tractability of the EWC protocol 194 In this section, we compared the computational tractability of the EWC protocol by 195 comparing full versions of TE and MI, and EWC versions of TE, cMI and the PC. We used 196 these methods to estimate functional networks comprising increasing number of nodes. 197 Nodal time series were sampled from source -localised MEG recordings (see Methods) . 198 199 As expected, the computation time increased with network size in all cases. The conventional 200 implementation of the TE was the most computationally intensive approach (Fig.3A), whereas 201 PC-EWC estimation required ≈75% less time relative to TE -Full. However, the EWC 202 implementation of TE was considerably quicker than the conventional approach and required 203 just as long as the cMI-EWC. Interestingly, MI estimated from the full signal requires less time 204 to compute than the closely related cMI -EWC. This is however due to an increase in the 205 number of computations in the EWC approach, as we will discuss later. 206 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 7 207 Figure 3 –Time taken for network inference as a func1on of network size. Signals spanning 200s (Twenty 10s epochs) sampled 208 from source-localised MEG recordings of 3 subjects. Sampling repeated 10 1mes per subject for each network size (Each box 209 plot spanning 30 data points). Circles represent outliers. (Inset) Percentage of 1me taken to compute PC-EWC, rela1ve to the 210 1me taken to compute TE-Full. For a range of network sizes, PC -EWC takes ≈75% less 1me to compute for a 200s scan, or 211 equivalently, is ≈4 1mes faster. 212 213 In short, we find that both methods of estimation (Full or EWC) and the within-window EWC 214 FC measure affect computation time. Importantly, using the EWC protocol generally leads to 215 a decrease in computation time. In addition, the gain in computational time grows as a 216 function of network size, indicating that EWC may confer pronounced benefits for the 217 estimation fine-grained functional networks (e.g. >1000 nodes). Combined with our results 218 from the previous section, we find that PC -EWC can be computed ≈4 times faster than TE-219 Full, while also resolving asymmetries in signalling. 220 221 Inferring whole -brain interaction patterns from MEG recordings 222 Having compared both the accuracy and computational efficiency of the EWC to the 223 conventional approach in-silico, in our final set of analyses, we extended the comparison onto 224 empirical neuroimaging data . Specifically, we compute and compare PC -EWC and TE-Full 225 from resting-state source-localised MEG recordings of 30 subjects , for the left hemisphere 226 (right-hemisphere results Fig.S3). The temporal ordering was based on delays proportional 227 to the inter-regional Euclidean distance (see Methods). 228 229 For each subject, we computed the correlation between the PC -EWC and TE-Full matrices 230 and found that the two measures led to highly correlated estimates of functional connectivity 231 0 200 400 600 800 1000 1200 1400 1600Network inference time for a 200s recording (s) 3 10 20 50 100 N TE-Full cMI-EWC TE-EWC MI-Full PC-EWC -100 -80 -60 -40 -20 Relative time taken to compute PC-EWC (% of time to compute TE-Full) 10 20 N 503 100 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 8 (Fig 4B). The median correlation coefficient across subjects was r≈0.82 (p<0.0001) (Fig 4C 232 top). Critically, although PC-EWC and TE-Full were very similar, the EWC protocol offered a 233 speedup per epoch by a factor of ≈6.5 (Fig.4D bottom). In addition to the high subject-level 234 agreement, the EWC implementation also allows us to study how signalling directions and 235 strengths change across the scan, by estimating the degree of asymmetry and variance of 236 estimates at different scales (event-, epoch-, and subject-levels) (Fig.S4). 237 238 To summarise, in this section we demonstrate a real-world application of the EWC, using PC-239 EWC and TE-Full. We show that the FC estimated as per the EWC shows good agreement 240 with our conventional reference using the TE, and does so at a fraction of the computational 241 cost. This lends support to the EWC as a viable technique for practical applications , 242 particularly involving multiple functional components. 243 244 Figure 4 – Applica1on of EWC on source-localised MEG recordings. (A) Res1ng-state MEG scans of 30 subjects from the Human Connectome Project (HCP) were pre-processed, source localised, orthogonalized, and epoched into 10 second segments. (B) (Top) The subject-level FC for the leb hemisphere was es1mated as PC-EWC and TE-Full. (Bo5om) Edgewise correla1on distribu1on between the subject-level FC matrices. The red line marks the median correla1on. (C) (Top) Sca5erplot of edge weights for a representa1ve subject (closest to the median correla1on). Correla1on aber excluding outliers (> 4 standard devia1on) - r≈0.82, p<0.0001. (Bo5om) inference 1me per epoch across subjects, for each of the methods. .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 9

245 Understanding how the brain routes information is an open problem in neuroscience with 246 crucial implications to our knowledge of perception and cognition. In this work, we focus on 247 the first step to wards addressing this problem – reliably inferring communication from 248 neuroimaging data. We introduce a time series analytical technique, termed Event -marked 249 Windowed Communication (EWC), to capture dynamic statistical relationships driven by inter-250 regional communication. This method quantifies inter-regional relationships using functional 251 connectivity (FC) measures estimated on select subsamples of the signal. We demonstrated 252 the merits and limitations of this method in silico, using different FC measures, and compared 253 it to established methods of inferring directional relationships (Fig.2). We then applied it to 254 source-localised MEG data, where high correlation was evident between subject-level FC 255 patterns derived using EWC and conventional approaches (Fig.4). Importantly, we also 256 showed that the EWC protocol is computationally efficient (≈6.5 times faster per 10s-epoch 257 when using the PC ) (Fig.3,4C-bottom), and additionally capable of capturing changes in 258 asymmetric relationships over the scan (Fig.S4). 259 260 Studying how activity or information flows between neural elements is crucial to our 261 understanding of the mechanisms of communication. In addition to EWC, multiple previous 262 works have explored this important topic from diverse standpoints, such as neuronal 263 avalanches26,27,31, direct stimulation -response measurements 33, functional hierarchies 21, FC 264 dynamics32,34, and information -theoretic methods 25,35,36, to name a few. EWC differs from 265 previous work primarily in its handling of time series . As opposed to measures/features 266 computed using the entirety of recordings or parts of them (as in sliding window approaches), 267 EWC restricts FC estimation to select subsamples of the signal, to infer activity propagation 268 between neural elements. By targeting the FC estimation to select subsamples of the scan , 269 EWC effectively computes a resting-state “event-related potential ”. Importantly, being a 270 general and flexible framework, the EWC implementation can be easily incorporated into 271 previously developed methodologies involving neural time series e.g. using feature vectors 272 instead of raw time series35. 273 274 Estimating EWC comprises three steps: 1) event identification 2) temporal ordering 3) 275 windowing/subsampling. Each of these elements have associated advantages and 276 limitations. We discuss these points below. 277 278 Typically, FC is estimated for the entire scan, or in the case of dynamical-FC, for consecutive 279 windows of the scan24,37,38. A limitation of such an approach to capture activity propagation 280 is that inter -regional communication may occur only at select time epochs and this 281 information may be obscured by internal dynamics. The purpose of the event-identification 282 step in the EWC is to identify salient features in regional dynamics, that then serve as a 283 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 10

point from which downstream effects are gauged. In task-based paradigms, these 284 segments are generally marked by controlled stimulus triggers. On the other hand, in task -285 free recordings, change-point detection or point-process analysis-based tools can be used to 286 capture salient changes in regional dynamics39–42. EWC reuses a simple method employed in 287 previous works that traces activity propagation to parts of the signal that deviate from mean 288 behaviour, through z-scoring31. Limiting the signal that is analysed to the proximity of the 289 events additionally makes the EWC computationally tractable when compared to 290 conventional approaches (Fig.3). 291 292 The temporal ordering aspect, which is designed to account for delays in signal conduction 293 over the network, incorporates directionality into the protocol, irrespective of the FC measure 294 used. This allows us to use undirected/symmetric measures like PC or cMI, to resolve transient 295 asymmetric FC relationships (Fig.2B,C,E,F). Any symmetric FC measure (coherence, phase 296 locking value, phase lag index etc.43–45) can similarly be rendered directional if implemented 297 in this manner. Care must however be taken to ensure that the considered delays are within 298 a reasonable range, to avoid capturing effects that might not likely be caused by the observed 299 event. 300 301 Limiting the FC estimation to a window/subsample of the signal proximal to the events offers 302 multiple advantages. It is well suited for scenarios where there is limited information regarding 303 the spatiotemporal span of the effects, allowing us to capture extended effects as opposed 304 to a local one (at a single time point) . Furthermore, averaging the results over multiple 305 windows (events) increases the robustness of the estimate. It also reduces the size of the 306 signal for which the FC must be estimated, decreasing the computational time. 307 308 In contrast to conventional FC implementations however, EWC typically involves a higher 309 number of number of computations – instead of a single computation per epoch, there are 310 now computations associated with each event within the epoch (albeit involving shorter 311 segments of data). This effect is apparent in Fig.3, where MI-Full is computed considerably 312 faster than its EWC counterpart – cMI-EWC. 313 314 Although we find that there is a strong agreement between the subject-level PC-EWC and 315 TE-Full, we noted that agree ment in terms of their asymmetry or directionality was not 316 statistically significant. This could be due to the PC and TE captur ing different statistical 317 dependencies in the data. This difference may have been subtle in the simple network-motif 318 (resulting in both the PC and TE capturing the true signalling asymmetries), but more 319 pronounced at larger scales. 320 321 Methodological considerations 322 323 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 11 We designed EWC to capture communication-driven statistical relationships or as a statistical 324 analogue of activity propagation . Its development involved the use of several simplifying 325 assumptions to ensure computational and analytical tractability. These points need to be 326 considered to fully appreciate the scope of this work and limitations. For example, a z-value-327 based event identification as we used in this study, will be sensitive to changes in the mean 328 activity – this is evident in Fig.S2F, where an increase in the source firing rate results in firing 329 events not being identified, resulting in diminished estimates. Similarly, errors in delay 330 estimation, which affect the temporal ordering, can result in the effect of the event being 331 missed entirely. The window length is also a crucial component and must be chosen so that 332 it captures the immediate effects of the event, while also being long enough to ensure that 333 there are enough datapoints for an accurate FC estimate. The window and epoch lengths 334 also determine the resolution at which the dynamics of asymmetric relationships can be 335 observed (event-/epoch-scale vs. subject-scale). 336 337 It is also important to remember that EWC is not a direct measure of communication but is 338 instead a measure of the statistical consequence of communication. It can only be used to 339 infer communication. Care must be taken while choosing the FC measure, since different 340 measures capture different features of signal similarity (correlation, phase synchrony etc.), 341 while also being valid only on certain derivatives of the original signal e.g. Mutual information 342 can be computed between any signal (with different interpretations), whereas the phase lag 343 index is computed between the instantaneous phase time series. The effect of FC measure 344 choice is evident in Fig.2B,C,E,F, where the cMI indicates strong interaction from node 4 to 345 2, despite the lack of any interaction by construction. The PC on the other hand , does not 346 capture this functional link. In the case of non-negative FC measures like the cMI, estimates 347 can sometimes be inflated due a positive bias (despite filtering based on significance testing). 348 This can be corrected by subtracting a null-derived estimate from the observed (significant) 349 estimate. 350 351 Additionally, it must be noted that although our protocol restricts the estimation of the EWC 352 to supra-threshold events, we do not claim that communication does not take place during 353 the sub-threshold segments. The significant deviations are used to systematically reduce the 354 analysis space, since it is more difficult to establish whether EWC estimates are due to 355 communication as opposed to some other regional process or noise , in the sub -threshold 356 segments. 357 358

359 In conclusion, our work presents a new method to infer neural communication patterns 360 through a restricted FC estimation approach based on activity propagation events. Our 361

yields FC estimates comparable to that of well -established techniques at a fraction 362 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 12 of the computational time, opening up new avenues for the investigation of signalling in large 363 networks of neural elements. 364 365 Code availability 366 All the code used to generate activity propagation maps using the EWC is available at 367 https://github.com/vmadanmohan/EWC . 368 369

370 Dataset 371 Resting state MEG scans of 30 subjects (22 -35 years, 17F), along with associated MEG 372 anatomical data, 3T structural MRI data, and empty-room recordings, were obtained from the 373 Human Connectome Project (Van Essen et al., 2013), through the ConnectomeDB platform. 374 The MEG scans varied in duration from 5-6 minutes, at a sampling rate ( 𝑆𝑅) of 2035Hz, and 375 anti-aliasing low pass filtered at 400Hz. 376 377 Processing 378 The MEG scans were processed entirely using the MATLAB 46 -based Brainstorm software47, 379 and in accordance with pipeline described in Brainstorm's HCP -MEG tutorial 48. MEG 380 recordings were first coregistered to the subject's structural MRI using the MEG anatomical 381 data. A notch filter (60, 120, 180, 240 and 300Hz), followed by high pass filter (0.3Hz) were 382 applied to resting state and empty -room recordings, to filter out power -supply and slow -383 wave/DC-offset artifacts respectively. Each subject's recording was then visually inspected, 384 along with the channel power spectral density, to weed out bad channels and bad time 385 segments. The ECG and EOG recordings were then used to identify heartbeats and eye -386 blinks, after which associated artifacts were removed using their signal space projections 387 (SSPs)49. 388 389 The source-level activities at 8000 points, defined by the fsLR4k mesh, were then estimated 390 from the sensor -level recordings. This involved first computing the head model using 391 overlapping spheres and constrained dipoles normal to the cortical surface, and estimating 392 the noise covariance from the empty -room recordings. Source -level activities were then 393 estimated using the dSPM method 50 available in Brainstorm. The sources were then 394 parcellated using the Schaefer-Yeo 7-network 100 atlas51, with the parcel activity computed 395 as the mean of the constituent source activities. 396 397 The centroid coordinates of the Regions of Interest (ROIs) were used to define a Euclidean 398 distance (ED) matrix between all possible pairs of regions. Assuming a conduction velocity of 399 10m/s, inter -regional distances were proportionally converted to time delays ( 𝛿 ) as 400 𝛿!" = 𝐸𝐷!" 𝑣⁄ 401 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 13 Where 𝑣 is the conduction velocity (CV) of neural signals, set to 10 𝑚 𝑠⁄ . 402 Parcellated source-localised resting-state MEG recordings were corrected for source leakage 403 effects by removing zero-lag correlations as per 52, using the OHBA Software Library (OSL) 404 package in Python53. 405 406 Estimating inter-areal communication 407 “-EWC” implementation 408 The three main elements of EWC are – 1) event-identification, and 2) Temporal ordering 409 based on conduction delays, and 3) estimation of the FC between event -marked 410 windows/subsamples. 411 Source-level parcellated time series were first epoched into 10𝑠 segments. Epochs that 412 contained bad segments were removed, and each epoch was then processed independently. 413 The data within an epoch was first z-scored: This helped identify the spread of regional 414 activities around their respective mean behaviours. For all regions, the timepoints at which 415 the |𝑧| > 3 were labelled as “events” and marked strong significant deviations from mean 416 behaviour. We argued that these significant events should cause measurable effects on the 417 activity of the rest of the brain. Additionally, finite signal conduction times would imply that 418 these downstream effects be delayed. We also focused the estimation of inter -areal 419 communication to a segment of the time series starting at the event timepoint – termed a 420 “communication window”, spanning a duration of 1𝑠. 421 422 The following steps were executed iteratively, over all ROIs, for each event: 423 424 1. A region was chosen as a source, and a window of 1𝑠 was defined starting at the 425 event. Any additional events that fell within this window were removed and not 426 considered in the subsequent iterations, to avoid excessive computation. The window 427 length was chosen to be short enough to capture the immediate effects of the 428 communication events and minimise contributions from internal dynamics. 429 2. A window of similar length was placed at the time series of all other ROIs (targets) 430 within the same hemisphere , at timepoints delayed in proportion to the Euclidean 431 distance between the source and the target (see Processing). The 1𝑠 length of the 432 window can also accommodate the effects of source signals travelling at velocities 433 much slower than 10𝑚/𝑠 (used in the calculation of inter-regional delays). 434 3. The FC measure (PC, cMI) was computed between all pairs of the source and target 435 signals within the respective windows, conditional upon the target's immediate past 436 activity over the duration of one window, up to the event. This conditioning was done 437 to remove biases due to internal dynamics (not associated with communication). The 438 FC estimates were checked for significance, with a Bonferroni-corrected threshold of 439 𝛼 = !.!# $!"!#$% , and set to zero for p -values above this threshold. The final measure is 440 termed the Event-marked Windowed Communication (EWC). 441 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 14 4. The event -level EWC values from the source to all possible targets were then 442 averaged to yield the source’s epoch-level EWC. 443 444 Once the above steps were completed for all possible sources, we were left with an epoch -445 level 𝑁 × 𝑁 EWC matrix, that captures the communication between all source -target pairs, 446 averaged over all the events in the epoch. For a given subject whose scans were divided into 447 𝑀 epochs, the communication protocol yields an 𝑁 × 𝑁 × 𝑀 EWC matrix over the entire scan 448 duration. This matrix is then averaged across epochs to give the subject-level EWC matrix. 449 450 “-Full” implementation 451 Source-level parcellated time series were first epoched into 10𝑠 segments. Epochs that 452 contained bad segments were removed, and each epoch was then processed independently. 453 The FC was computed in a pairwise manner, between all regions, between the entire epoch 454 time series. Each FC estimate was significance tested with a threshold of 𝛼 = 0.01. Like the 455 EWC implementation, the pairwise FC was computed for all epochs, yielding an 𝑁 × 𝑁 × 𝑀 456 matrix (for 𝑀 epochs) over the entire scan duration. 457 458 Partial correlation 459 The partial correlation (PC) quantifies the linear association between two random variables, 460 while discounting the effects of a third variable. In the EWC implementation, we set the third 461 variable as the target’s past over a single window duration, up to the timepoint of the event. 462 We compute the PC using MATLAB’s “partialcorr” function. The significance level was set to 463 𝛼 = !.!# $!"!#$% . The partial correlation is a continuous function varying from -1 to 1 , indicating 464 perfectly negative to perfectly positive linear correlation respectively. 465 466 Conditional Mutual Information 467 The conditional Mutual Information (cMI) is an information theoretic measure that quantifies 468 the shared information between random variables (representing regional activities) in the 469 context of the activity of some exogenous additional regions (represented as a conditioning 470 variable). It is a variant of the widely used Mutual Information (MI)25,54–56. 471 Given random variables X and Y, containing the activities of two regions, and a third random 472 variable Z, containing the activity of a third region (or a set of multiple regions), the cMI is 473 measured in terms of entropies as: 474 𝐼(𝑋; 𝑌|𝑍) = 𝐻(𝑋|𝑍) + 𝐻(𝑌|𝑍) − 𝐻(𝑋, 𝑌|𝑍) 475 476 Where 𝐻(𝑋|𝑍) and 𝐻(𝑌|𝑍) are the conditional entropies associated with variables X and Y, 477 and 𝐻(𝑋, 𝑌|𝑍) is the conditional joint entropy. 478 The cMI is a symmetric measure, that varies from 0 when the two random variables are 479 independent of each other, to ∞ when they are identical. 480 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 15 The target’s past activity was used as the conditional variable, to remove any biases to the 481 measured communication from internal dynamics. 482 483 Transfer Entropy 484 The Transfer Entropy is a widely used information-theoretic measure of directional influences 485 2,3,25,36,57,58. Specifically, for two random variables X and Y, it quantifies how much uncertainty 486 in the future of Y (target) is reduced through the knowledge of the past of X (source), 487 conditional on the past of Y. It is technically a special case of the cMI, where the shared 488 information is estimated between relatively time-shifted random variables: 489 𝑇%→'(𝑘, 𝑙) = 𝐼 @𝑋( (*); 𝑌(,#A𝑌( (-)B , 490 Where 𝑋( (*) = {X./0/#, … , X./#, X.} . There are multivariate variants of the TE, which also 491 condition it on other possibly mediating variables. In this work, we use its simplest form – the 492 bivariate TE. Multivariate forms of the TE are generally better than other techniques at 493 capturing ground truth asymmetries but are considerably more data - and computation -494 intensive25. 495 496 We used the Gaussian estimator in the Java Information Dynamics Toolkit (JIDT)59 to compute 497 all information theoretic quantities. Using the Gaussian estimator allowed each estimate to 498 be significance -tested analytically from the χ1 -distribution. The estimators were run with 499 default properties, except for the source -destination delay in the TE estimator, which was 500 changed from the default value of 1, to the Euclidean distance -based inter-regional delay 501 (see Processing). All information theoretic quantities (including the EWC) are measured in 502 nats. 503 504 Demonstration on a network model 505 In order to identify the merits and limitations of our approach prior to its application to 506 neuroimaging data, we tested it on the activities of a four-node system with Linear Stochastic 507 Model60 dynamics. The dynamics of a node 𝑖 is described as: 508 𝑥2̇ = −𝑥3 + 𝐾 L 𝐶34 𝑥4N𝑡 − 𝛿34Q 4∈$ + 𝜎𝒩(0,1) 509 Where 𝐶34 is the connectivity strength between nodes 𝑖 and 𝑗 , 𝐾 is a global coupling 510 parameter ( 𝐾 = 1), δ34 is the time delay between 𝑖 and 𝑗 (proportional to the Euclidean 511 distance), 𝒩(0,1) denotes random standard normal noise , and σ is the noise amplitude. 512 Three of the nodes were connected in a linear chain-like topology, and the fourth node was 513 isolated from the rest of the system i.e. 𝐶#1 = 𝐶1# = 𝐶16 = 𝐶61 = 1; 𝐶#6 = 𝐶6# = 0; 𝐶73 = 𝐶37 =514 0 ∀ 𝑖 ∈ {1,2,3}. 515 In addition to the LSM dynamics, a Poisson process caused both the terminal nodes to 516 “pulse” randomly, with an average frequency ν, and with a pulse amplitude of 0.1. This 517 resulted in the dynamics of the central node to be a mix of its internal dynamics, and random 518 inputs from the two terminal nodes. While Poisson spiking is a convenient means of exercising 519 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 16 control over ground-truth signalling, it can result in an inflation of nodes’ mean activities from 520 the “baseline” at high ν, due to rapid fluctuations. 521 The communication protocol was tested across systematic variations of σ, 𝛿16, and ν. 522 523 For each set of parameter values, we simulated the dynamics for 205𝑠, at a sampling rate of 524 2035𝐻𝑧, using an Euler integrator. Epoch duration and communication window durations 525 were set to 10𝑠 and 1𝑠 respectively. 20 trials of each simulation were carried out. 526 527 Communication was inferred using the TE (Full implementation), cMI and PC (EWC 528 implementation) (see Protocol), on the activity time series of all four nodes. 529 530 Computation time benchmarking 531 To compare the computation times for different protocols as a function of network size, 532 activity propagation was inferred as TE -Full, MI-Full, TE-EWC, cMI-EWC, and PC -EWC for 533 time-series data of 200s. To ensure that the estimated computation times were realistic, the 534 recordings were randomly sampled source localised MEG recordings (see Processing). The 535 following steps were repeated for a network size 𝑁 (for 10 trials): 536 1. 𝑁 regions are randomly chosen, and the associated regional time series were trimmed 537 to 200s (sampling frequency = 2035Hz). 538 2. Inter-regional delays were set to 0.015ms. 539 3. Activity propagation was estimated between all pairs of regions as TE -Full, MI-Full, 540 TE-EWC, cMI-EWC, and PC-EWC. 541 4. Repeat with a new random selection of regions. 542 543 Acknowledgment s 544 V.M.M. is supported by the Melbour ne Research Scholarship , University of Melbourne . 545 R.F.H.C. is funded by a NHMRC Emerging Leadership Investigator Grant ( Grant number: 546 2017527). C.S. acknowledges support from the Australian Research Council ( Grant number: 547 DP170101815). A.Z. is supported by an ARC Future Fellowship and the Rebecca L. Cooper 548 Foundation. Data were provided [in part] by the Human Connectome Project, WU -Minn 549 Consortium (Principal Investigators: David Van Essen and Kamil Ugurbil; 1U54MH091657) 550 funded by the 16 NIH Institutes and Centers that support the NIH Blueprint for Neuroscience 551 Research; and by the McDonnell Center for Systems Neuroscience at Washington University. 552 This research was supported by The University of Melbourne’s Research Computing Services 553 and the Petascale Campus Initiative. 554 555 556 557 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 17

558 1. Avena-Koenigsberger, A., Misic, B. & Sporns, O. Communication dynamics in complex 559 brain networks. Nature Reviews Neuroscience vol. 19 17–33 Preprint at 560 https://doi.org/10.1038/nrn.2017.149 (2018). 561 2. Kirst, C., Timme, M. & Battaglia, D. Dynamic information routing in complex networks. 562 Nat Commun 7, (2016). 563 3. Palmigiano, A., Geisel, T., Wolf, F. & Battaglia, D. Flexible information routing by transient 564 synchrony. Nat Neurosci 20, 1014–1022 (2017). 565 4. Hahn, G., Ponce-Alvarez, A., Deco, G., Aertsen, A. & Kumar, A. Portraits of 566 communication in neuronal networks. Nat Rev Neurosci 20, 117–127 (2019). 567 5. Seguin, C., Sporns, O. & Zalesky, A. Brain network communication: concepts, models and 568 applications. Nat Rev Neurosci (2023) doi:10.1038/s41583-023-00718-5. 569 6. Fries, P . A mechanism for cognitive dynamics: Neuronal communication through neuronal 570 coherence. Trends Cogn Sci 9, 474–480 (2005). 571 7. Fries, P . Rhythms for Cognition: Communication through Coherence. Neuron 88, 220–572 235 (2015). 573 8. Jensen, O. & Mazaheri, A. Shaping functional architecture by oscillatory alpha activity: 574 Gating by inhibition. Front Hum Neurosci 4, (2010). 575 9. Bonnefond, M., Kastner, S. & Jensen, O. Communication between brain areas based on 576 nested oscillations. eNeuro 4, (2017). 577 10. Hahn, G., Bujan, A. F., Frégnac, Y., Aertsen, A. & Kumar, A. Communication through 578 Resonance in Spiking Neuronal Networks. PLoS Comput Biol 10, (2014). 579 11. Goñi, J. et al. Exploring the Morphospace of Communication Efficiency in Complex 580 Networks. PLoS One 8, (2013). 581 12. Seguin, C., Van Den Heuvel, M. P . & Zalesky, A. Navigation of brain networks. Proc Natl 582 Acad Sci U S A 115, 6297–6302 (2018). 583 13. Avena-Koenigsberger, A. et al. A spectrum of routing strategies for brain networks. PLoS 584 Comput Biol 15, (2019). 585 14. Bazinet, V., Vos de Wael, R., Hagmann, P ., Bernhardt, B. C. & Misic, B. Multiscale 586 communication in cortico-cortical networks. Neuroimage 243, (2021). 587 15. Schipul, S. E., Keller, T. A. & Just, M. A. Inter-regional brain communication and its 588 disturbance in autism. Frontiers in Systems Neuroscience Preprint at 589 https://doi.org/10.3389/fnsys.2011.00010 (2011). 590 16. Chapeton, J. I., Haque, R., Wittig, J. H., Inati, S. K. & Zaghloul, K. A. Large-Scale 591 Communication in the Human Brain Is Rhythmically Modulated through Alpha 592 Coherence. Current Biology 29, 2801-2811.e5 (2019). 593 17. Papadopoulos, L., Lynn, C. W., Battaglia, D. & Bassett, D. S. Relations between large-594 scale brain connectivity and effects of regional stimulation depend on collective 595 dynamical state. PLoS Comput Biol 16, (2020). 596 18. Schnitzler, A. & Gross, J. Normal and pathological oscillatory communication in the brain. 597 Nat Rev Neurosci 6, 285–296 (2005). 598 19. Shafiei, G., Baillet, S. & Misic, B. Human electromagnetic and haemodynamic networks 599 systematically converge in unimodal cortex and diverge in transmodal cortex. PLoS Biol 600 20, (2022). 601 20. Suárez, L. E., Markello, R. D., Betzel, R. F. & Misic, B. Linking Structure and Function in 602 Macroscale Brain Networks. Trends Cogn Sci 24, 302–315 (2020). 603 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 18 21. Vázquez-Rodríguez, B., Liu, Z. Q., Hagmann, P . & Misic, B. Signal propagation via cortical 604 hierarchies. Network Neuroscience 4, 1072–1090 (2020). 605 22. Lynall, M. E. et al. Functional connectivity and brain networks in schizophrenia. Journal of 606 Neuroscience 30, 9477–9487 (2010). 607 23. Han, S. et al. Resting state functional network switching rate is differently altered in 608 bipolar disorder and major depressive disorder. Hum Brain Mapp 41, 3295–3304 (2020). 609 24. Hutchison, R. M. et al. Dynamic functional connectivity: Promise, issues, and 610 interpretations. Neuroimage 80, 360–378 (2013). 611 25. Novelli, L. & Lizier, J. T. Inferring network properties from time series using transfer 612 entropy and mutual information: Validation of multivariate versus bivariate approaches. 613 Network Neuroscience 5, 373–404 (2021). 614 26. Beggs, J. M. & Plenz, D. Behavioral/Systems/Cognitive Neuronal Avalanches in 615 Neocortical Circuits. (2003). 616 27. Beggs, J. M. & Plenz, D. Neuronal avalanches are diverse and precise activity patterns 617 that are stable for many hours in cortical slice cultures. Journal of Neuroscience 24, 5216–618 5229 (2004). 619 28. Lerner, T. N., Ye, L. & Deisseroth, K. Communication in Neural Circuits: Tools, 620 Opportunities, and Challenges. Cell 164, 1136–1150 (2016). 621 29. Shriki, O. et al. Neuronal avalanches in the resting MEG of the human brain. Journal of 622 Neuroscience 33, 7079–7090 (2013). 623 30. Mitra, A. & Raichle, M. E. How networks communicate: Propagation patterns in 624 spontaneous brain activity. Philosophical Transactions of the Royal Society B: Biological 625 Sciences 371, (2016). 626 31. Sorrentino, P . et al. The structural connectome constrains fast brain dynamics. Elife 10, 627 (2021). 628 32. Griffa, A. et al. Transient networks of spatio-temporal connectivity map communication 629 pathways in brain functional systems. Neuroimage 155, 490–502 (2017). 630 33. Seguin, C. et al. Communication dynamics in the human connectome shape the cortex-631 wide propagation of direct electrical stimulation. Neuron (2023) 632 doi:10.1016/j.neuron.2023.01.027. 633 34. Faskowitz, J., Esfahlani, F. Z., Jo, Y., Sporns, O. & Betzel, R. F. Edge-centric functional 634 network representations of human cerebral cortex reveal overlapping system-level 635 architecture. Nat Neurosci 23, 1644–1654 (2020). 636 35. Nguyen, A., McMullin, O., Lizier, J. T. & Fulcher, B. D. A feature-based information-637 theoretic approach for detecting interpretable, long-timescale pairwise interactions from 638 time series. (2024). 639 36. Novelli, L., Wollstadt, P ., Mediano, P ., Wibral, M. & Lizier, J. T. Large-scale directed 640 network inference with multivariate transfer entropy and hierarchical statistical testing. 641 Network Neuroscience 3, 827–847 (2019). 642 37. Schulz, D. & Huston, J. P . The sliding window correlation procedure for detecting hidden 643 correlations: existence of behavioral subgroups illustrated with aged rats. J Neurosci 644

121, 129–137 (2002). 645 38. Vergara, V. M., Abrol, A. & Calhoun, V. D. An average sliding window correlation method 646 for dynamic functional connectivity. Hum Brain Mapp 40, 2089–2103 (2019). 647 39. Aminikhanghahi, S. & Cook, D. J. A survey of methods for time series change point 648 detection. Knowl Inf Syst 51, 339–367 (2017). 649 40. Chen, G., Lu, G., Shang, W. & Xie, Z. Automated Change-Point Detection of EEG Signals 650 Based on Structural Time-Series Analysis. IEEE Access 7, 180168–180180 (2019). 651 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Madan Mohan et al., 2024 19 41. Tagliazucchi, E., Balenzuela, P., Fraiman, D. & Chialvo, D. R. Criticality in large-scale brain 652 fmri dynamics unveiled by a novel point process analysis. Front Physiol 3 FEB, (2012). 653 42. Tagliazucchi, E., Balenzuela, P., Fraiman, D., Montoya, P. & Chialvo, D. R. Spontaneous 654 BOLD event triggered averages for estimating functional connectivity at resting state. 655 Neurosci Lett 488, 158–163 (2011). 656 43. Colclough, G. L. et al. How reliable are MEG resting-state connectivity metrics? 657 Neuroimage 138, 284–293 (2016). 658 44. Lachaux, J. P ., Rodriguez, E., Martinerie, J. & Varela, F. J. Measuring phase synchrony in 659 brain signals. Hum Brain Mapp 8, 194–208 (1999). 660 45. Stam, C. J., Nolte, G. & Daffertshofer, A. Phase lag index: Assessment of functional 661 connectivity from multi channel EEG and MEG with diminished bias from common 662 sources. Hum Brain Mapp 28, 1178–1193 (2007). 663 46. The MathWorks Inc. MATLAB version: 9.13.0 (R2022b). Preprint at (2022). 664 47. Tadel, F., Baillet, S., Mosher, J. C., Pantazis, D. & Leahy, R. M. Brainstorm: A user-friendly 665 application for MEG/EEG analysis. Comput Intell Neurosci 2011, (2011). 666 48. Niso, G. et al. Brainstorm pipeline analysis of resting-state data from the open MEG 667 archive. Front Neurosci 13, (2019). 668 49. Uusitalo, M. A. & Ilmoniemi, R. J. Communication Signal-Space Projection Method for 669 Separating MEG or EEG into Components. Biol. Eng. & Cornput vol. 35 (1997). 670 50. Dale, A. M. et al. Dynamic Statistical Parametric Mapping: Combining fMRI and MEG for 671 High-Resolution Imaging of Cortical Activity. Neuron 26, 55–67 (2000). 672 51. Schaefer, A. et al. Local-Global Parcellation of the Human Cerebral Cortex from Intrinsic 673 Functional Connectivity MRI. Cerebral Cortex 28, 3095–3114 (2018). 674 52. Colclough, G. L., Brookes, M. J., Smith, S. M. & Woolrich, M. W. A symmetric multivariate 675 leakage correction for MEG connectomes. Neuroimage 117, 439–448 (2015). 676 53. Quinn, A. J., van Es, M. W. J., Gohil, C. & Woolrich, M. W. OHBA Software Library in 677 Python (OSL). Preprint at (2023). 678 54. Hulata, E., Segev, R. & Ben-Jacob, E. A method for spike sorting and detection based on 679 wavelet packets and Shannon’s mutual information. J Neurosci Methods 117, 1–12 680 (2002). 681 55. Wilmer, A., de Lussanet, M. & Lappe, M. Time-Delayed Mutual Information of the Phase 682 as a Measure of Functional Connectivity. PLoS One 7, (2012). 683 56. Palus, M. Detecting Phase Synchronization in Noisy Systems. Physics Letters A vol. 235 684 (1997). 685 57. Schreiber, T. Measuring Information Transfer. Phys Rev Lett 85, 461–464 (2000). 686 58. Shorten, D. P ., Spinney, R. E. & Lizier, J. T. Estimating transfer entropy in continuous time 687 between neural spike trains or other event-based data. PLoS Comput Biol 17, (2021). 688 59. Lizier, J. T. JIDT: an information-theoretic toolkit for studying the dynamics of complex 689 systems. Front Robot AI 1, (2014). 690 60. Galán, R. F. On how network architecture determines the dominant patterns of 691 spontaneous neural activity. PLoS One 3, (2008). 692 693 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Event 2 Event 3Event 1 Source Target Target 1 2 3 𝛿12 Communication window 𝛿13 Average over all events Epoching Event localization (Activity propagation pattern between all pairs of regions, in an epoch) EWC 1 1 2 2 3 3 Targets Sources Repeat over epochs Epoch 1 Epoch 2 Epoch M ... A B C D .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Partial Correlation (EWC) 1 3 24 42 𝜎 𝜈 = 0.2Hz = 15ms42= = A B C D E F .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint 0 200 400 600 800 1000 1200 1400 1600Network inference time for a 200s recording (s) 3 10 20 50 100 N TE-Full cMI-EWC TE-EWC MI-Full PC-EWC -100 -80 -60 -40 -20 Relative time taken to compute PC-EWC (% of time to compute TE-Full) 10 20 N 503 100 .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint Corr. Biv. Transfer Entropy (Full) Sub. 1 Sub. 2 Sub. N Partial correlation (EWC) Sub. 1 Sub. 2 Sub. N 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 |Pearson-R| 0 2 4 6 8 10 12Probability Density Median correlation = 0.82095 Human Connectome Project (HCP) MEG subjects (n=30) . . . MEG scans (Sensor space) Regional activities (Source space) Source localisation Parcellation Delay estimation Orthogonalisation Source Target Communication window Information flow A B C Partial corr. (EWC) TE (Full) 0 10 20 30 40 50 60Inference time per epoch (s) .CC-BY-NC 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2024. ; https://doi.org/10.1101/2024.07.30.605466doi: bioRxiv preprint