Measurement of dynamic task related functional networks using MEG

The characterisation of dynamic electrophysiological brain networks, which form and dissolve in order to support ongoing cognitive function, is one of the most important goals in neuroscience. Here, we introduce a method for measuring such networks in the human brain using magnetoencephalography (MEG). Previous network analyses look for brain regions that share a common temporal profile of activity. Here distinctly, we exploit the high spatio-temporal resolution of MEG to measure the temporal evolution of connectivity between pairs of parcellated brain regions. We then use an ICA based procedure to identify networks of connections whose temporal dynamics covary. We validate our method using MEG data recorded during a finger movement task, identifying a transient network of connections linking somatosensory and primary motor regions, which modulates during the task. Next, we use our method to image the networks which support cognition during a Sternberg working memory task. We generate a novel neuroscientific picture of cognitive processing, showing the formation and dissolution of multiple networks which relate to semantic processing, pattern recognition and language as well as vision and movement. Our method tracks the dynamics of functional connectivity in the brain on a timescale commensurate to the task they are undertaking.

2 separate tensors, one representing simulated networks , and the second representing interference .
To generate , four spatially distinct networks were constructed based on a previous MEG study . The spatial patterns of connectivity were each represented by an ( ( × ( ) adjacency matrix, , , (where = 1 − 4 and ( is the number of AAL regions). These matrices, alongside a 3D visualisation of the network, are shown in columns 1-4 of Figure S1 and reflect visual, sensorimotor, superior frontal and frontoparietal networks. In addition to the spatial signature, we simulated the time evolution of dynamic connectivity in each network. (i.e. we allow each of these networks to form and dissolve on a timescale commensurate to brain activity). In order to do this, Each network was assigned a simulated timecourse 3 ( ) which was given by 3 ( ) = 93 ( ) + <3 ( ), [S1] In order to generate the network contribution, , to the simulated adjacency , The 4 separate network matrices were combined in a weighted sum, such that ( ), the single adjacency matrix at single time point , is given by: Concatenation over time generated , which represents transient formation of our 4 simulated networks over a timecourse spanning 60 minutes.
In order to add realistic interference to our simulated adjacency tensor we added a second adjacency tensor, , generated using empty room MEG data. 60 minutes of empty room data were recorded at 600 Hz using a 275-channel MEG system (MISL, Coquitlam, BC, Canada) in synthetic third order gradiometer configuration.
Empty room data were processed in the same way as the data in the main paper (save for correcting for leakage) to form . Note here that because we do not perform leakage correction, common mode signals recorded at the MEG sensors in the empty room recording will correlate over AAL regions. This manifests as networks of connectivity which will necessarily change in time (according to fluctuating interference). The interference tensor, , thus effectively represents networks of no interest. [Note that in reality networks of no interest could be formed through environmental noise, biological interference, or the brain itself. Although where that SNR represents the effective signal to noise of the data. A derivation showing the origins of this combination of correlation coefficients is given in our appendix below.

Testing the effectiveness of ICA
We applied the ICA decomposition described in the main paper to in order to test whether, in the presence of interference, we would accurately reconstruct the spatial signatures and timecourses of the stimulated networks. The similarities between simulated and reconstructed data were measured by Pearson correlation between 1) simulated and reconstructed timecourses and 2) simulated and reconstructed adjacency matrices. To summarise the performance of ICA in a single run, for a single network, we developed can range between 0 and 1, where 1 would correspond to a single independent component representing a network with no other component contributing. Importantly, the figure of merit penalises both poor representation and degeneracy, that is to say when all components can equally explain a single simulated network.

Testing SNR and temporal overlap
We extracted 7 components ( OX = 7), and SNR was allowed to range between 0 and 3 in steps of 0.02. For these initial simulations, the onset of the Hanning window was as described above (i.e. on the 3rd, 18th, 33rd and 48th second of every minute for the visual, sensorimotor, frontal and fronto-parietal networks respectively).
We also tested the effect of changing the overlap between temporal representation of networks. In the initial simulation, the time lag between networks peak connectivity was Δ=15 s. However, we wanted to investigate the effect of reducing Δ to the point when two networks timecourses would overlap (i.e. networks were no longer independent). To do this, the simulation was simplified to use only two networks (1 and 2) and the interference removed. Δ was changed from 0 s to 15 s in 1 s steps; ICA was asked to return only 2 components.

Results
: Figure   A) The trial averaged independent component timecourses. B) 3D representations of the independent component (i.e. corresponding columns of the mixing matrix) thresholded to 50% of the maximum value. C) Bar plots showing temporal (left) and spatial (right) correlation between the reconstructed independent components and the original simulated networks. The columns are colour coded to match the colours of the timecourses, and show there is one clear match, with little residual mixing. Figure S3 shows our temporal figure of merit plotted against SNR, for all four simulated networks. At high SNR, all 4 simulated networks are reconstructed faithfully using ICA. However, decreasing SNR sees a sharp transition meaning that, below some minimum 'threshold' value, a specific network is no longer reconstructed. Importantly, these thresholds differ across networks, showing, for example, that the frontal or fronto-parietal networks can be reconstructed using data at a lower SNR compared to the visual or sensorimotor networks. This interesting dissociation likely results from the number of connections in the network; It can be seen clearly from the matrices in Figure S1 that the visual and sensorimotor networks have less connections than the frontal and fronto-parietal networks, and it is these networks that require higher SNR. This relationship is further shown by the graph in Figure S3B. Here, number of connections is represented by the Frobenius norm of the network adjacency matrix. Note that Frobenius norm (given by where O3 represents the element in the i th row and j th column of ) usually represents both number and strength of connections, however since maximum strength of any connection (maximum value of O3 ) in each simulated matrix was made equal, Frobenius norm, in this case, is solely representative of spatial extent. Figure S3B shows clearly that SNR is a monotonically decreasing function of norm, meaning that networks with large extent can be reconstructed using lower SNR data. Figure S3: The effect of the signal to noise ratio (SNR). A) The mean temporal figure of merit for each network plotted as a function of SNR; note that for high SNR all networks are well characterised. However, a critical SNR exists below which our ability to successfully reconstruct any one network drops. Panel B shows the relationship between the Frobenius Norm of the simulated adjacency matrix (representing spatial extent of the network) and critical SNR (defined as the SNR required for a temporal figure of merit of 0.5). Note that adjacency matrices with a high Frobenius norm can be reconstructed using lower SNR data.
Finally, Figure S4 shows the effect of increasing the temporal overlap between simulated network timecourses (i.e. decreasing Δ). It is well known that ICA attempts to pull out temporally independent timecourses meaning that, when timecourses of separate networks begin to overlap, the spatial representations of those networks will become mixed. This is a fundamental assumption of ICA, and thus the way in which networks are characterised using our method. Figure S4A shows the temporal Figure of merit (for a 2 network simulation, see above) plotted as a function of Δ. As expected when timecourses become mixed, the networks are no longer reconstructed separately. Panels B, C and D show both temporal and spatial representations of the independent components for Δ = 6 s, Δ = 3 s and Δ = 0 s respectively. As expected, with no temporal overlap (Δ = 6 s) both the spatial and temporal signature of the two extracted components represent a single simulated network. However, with increasing overlap those spatial and temporal properties combine, meaning that the two networks are no longer separable.

Discussion
We have shown in simulation that ICA can successfully separate functional networks based on their connectivity timecourses. However, importantly, we show that success of this algorithm is a function of both signal to noise ratio and the number of connections in a network. As would be expected, the ability of ICA to uniquely extract networks is diminished at low signal to noise ratio. [A18] So this is how the signal and noise adjacency matrices should be combined. Note that for an SNR of 1, you simply average the 2 correlation coefficients, i.e = 0.5( f + O ), and the value is always bounded by -1 and 1. For extremely high SNR, the term