Modular segregation drives causality of the dynamic oscillatory network responses during threat processing

Abstract Physiological responses to threat and stress stimuli entrain synchronized neural oscillations among cerebral networks. Network architecture and adaptation may play a critical role in achieving optimal physiological responses, while alteration can lead to mental dysfunction. We reconstructed cortical and sub-cortical source time series from high-density electroencephalography, which were then fed into community architecture analysis. Dynamic alterations were evaluated in terms of flexibility, clustering coefficient and global and local efficiency, as parameters of community allegiance. Transcranial magnetic stimulation was applied over the dorsomedial prefrontal cortex during the time window relevant for physiological threat processing and effective connectivity was computed to test the causality of network dynamics. A theta band-driven community re-organization was evident in key anatomical regions conforming the central executive, salience network and default mode networks during instructed threat processing. Increased network flexibility entrained the physiological responses to threat processing. The effective connectivity analysis showed that information flow differed between theta and alpha bands and were modulated by transcranial magnetic stimulation in salience and default mode networks during threat processing. Theta oscillations drive dynamic community network re-organization during threat processing. Nodal community switches may modulate the directionality of information flow and determine physiological responses relevant to mental health.


MRI data acquisition
For all participants magnetic resonance images were acquired at the Neuroimaging Center

Heart rate estimation
In the EEG signals, volume conduction is thought to be linear and instantaneous, and it is expected that the sources of cardiac signals are not generally time locked to the sources of EEG activity, reflecting the activity of cortical neurons. 1 The lCA can accurately identify the time courses of activation and scalp topographies of relatively large and temporally-independent sources from simulated scalp recordings, even in the presence of a large number of low-level and temporally-independent source activities. 2 For heart rate detection analysis, the rows of the input matrix y are the EEG signals recorded at the 256 electrodes, the rows of the output data matrix ϑ = Xy are time courses of activation of the lCA components, and the columns of the inverse matrix, X −1 , give the projection strengths of the respective components onto the scalp sensors.
In general, and unlike PCA, the component time courses of activation will be non-orthogonal.
Corrected EEG signals can then be derived as y ′ = (X) −1 ϑ′, where ϑ′ is the matrix of activation waveforms,ϑ , with rows representing sources of cardiac artifacts which are then extracted for further estimations from each participant. In total for experiment 1 we concatenated the 36 CS+ trials to take a total of 180 seconds and 24 CS-trials to take 120 seconds. For experiment 2 we concatenated the 54 CS+ trials to take 270 seconds and 36 CS-trials to take 180 seconds.

Reliability check of the EEG signals using inter-trial phase coherence (ITPC) analyses
Single trial data were first decomposed into their time-frequency representation by using the multitaper method. 3,4 In this method the spectrum is estimated by multiplying the data with K different windows (i.e. tapers). In this study, K = 7 orthogonal tapers were used with good leakage and spectral properties, the discrete prolate spheroidal sequences (DPSS) are applied. 5 The electrodes were grouped in line with the five lobes to have a global robust measure over a v certain lobe instead of selecting certain electrodes. The ITPC was estimated in theta (4)(5)(6)(7) and alpha (8)(9)(10)(11)(12)(13) Hz) frequency bands separately. This measure is stimulus-locked and independent of amplitude changes. A value of 0 represents absence of synchronization and a value of 1 indicates perfect synchronization. The baseline activity was taken as a reference and was calculated as the average at each frequency band and across conditions, from -250 to 0 ms before visual stimuli. The change with respect to the baseline interval at 6 windows of 250 ms each after the visual stimulus was then extrapolated up to 1500 ms. Finally, the ITPC difference between the CS+ and CS-was estimated. The significance levels of the ITPC are assessed using surrogate data by randomly shuffling 1000 times the single-trial spectral estimates from different latency windows during the baseline period.

Reconstruction of brain activity
The forward problem is the computation of the scalp potentials for a set of neural current sources. An established procedure was used by estimating the lead-field matrix with specified models for the brain; a volume conduction model with a finite-element method (FEM) was used. 6 For the forward modelling the surfaces of the compartments like the skin, skull, CSF, gray matter, and white matter extracted from the individual T1 MRI, and individual electrode locations were used. The forward modeling and the source analysis were done in FieldTrip. 7 The lead-field matrix (LFM) contains information about the geometry and conductivity of the model. The complete description of the solution for the forward problem has been described previously. 8,9 A full description of the beamformer linear constrained minimum variance spatial filter is given elsewhere. 10, 11 The output of the beamformer at a voxel in the brain can be defined as a weighted sum of the output of all EEG channels. The weights determine the spatial filtering characteristics of the beamformer and are selected to increase the sensitivity to signals from a voxel and reduce the contributions of signals from (noise) sources at different locations. The frequency components and their linear interaction are represented as a cross-spectral density (CSD) matrix. In order to visualize power at a given frequency range, a linear transformation was used based on a constrained optimization problem, which acts as a spatial filter. 10 The spatial filter assigned a specific value of power to each voxel. For a given source the beamformer weights for a location of interest are determined by the data covariance matrix and the LFM. A voxel size of 5 mm was used in this study, resulting in 6676 voxels covering the entire brain. The created source model was then interpolated on the brain regions defined according to the Harvard Oxford cortical-subcortical regions of interest (ROIs) defined in the MNI space. For each frequency band (theta and alpha) the activated voxels were selected by a within-subject surrogate analysis to define the significance level, which was then used to identify voxels in the regions as activated voxels. Once the brain region voxels were identified, their activity was extracted from the source space. In a further analysis, all the original source signals for each Harvard-Oxford 12 region with several activated voxels were combined by estimating the second order spectra and employing a weighting scheme depending on the analyzed frequency range to form a pooled source signal estimate for each region as previously described separately for both stimulus (Cs+, CS-). 13-15 Finally, the time series difference between the two conditions (CS+, CS-) was obtained for all the following analyses.

Investigating causal relationships between network nodes
The effective connectivity analysis was performed on all the nodes separately, for each of the three new communities and each frequency separately. Using time-frequency causality, is possible not only to focus on a particular frequency, but also to analyze the dynamics of the causality at that frequency. The time-frequency causality estimation using the TPDC is based on dual extended Kalman filtering (DEKF), 16,17 and allows time-dependent auto regressive (AR) coefficients to be estimated. One EKF estimates the states and feeds this information to the other; the second EKF estimates the model parameters and shares this information with the first. By using two Kalman filters working in parallel with one another, we can estimate both states and model parameters of the system at each time instant. After estimating the timedependent multivariate (MVAR) coefficients, the next step is to use those coefficients for the calculation of causality between the time series. By calculating the time-dependent MVAR coefficients at each time point, we can also calculate partial directed coherence (PDC) at each time point. The frequency bands taken into account were the theta and alpha. After estimating the TPDC values the significance level was calculated from the applied data using a bootstrapping method 18 . In short, we divide the original time series into smaller nonoverlapping windows and randomly shuffle the order of these windows to create a new time series. Data length was taken from -250 milliseconds to 1500 milliseconds i.e. 1750 milliseconds. Number of data points was computed by taking the product of the data length and sampling rate (1.750secs X 250 = 437.5 data points). Data length was at least 438 data points (multiplied with number of signals) for sub-samples. Model order used was 30 and data length was long enough for this model order. The model order was estimated for the first sub-interval and then it was fixed for the subsequent samples. We used Akaike's final perdition error as our model estimation criterion. Based on the literature, the model order can be estimated using Schwarz Bayesian Criterion (SBC) or Akaike's final prediction error (FPE). However, the FPE tends to over-estimate the model order than SBC because SBC favors space solutions and penalizes the number of parameters more strongly than FPE. 19 Moreover, it is also generally preferred to select model order that is slightly larger large than the other way around. 19 Hence, we used Akaike's model. To note, it has been shown that a small variation of model order do not significantly alter the connectivity patterns. 20 Time-varying MVAR coefficients were calculated using sub-intervals so that the dynamics of MVAR coefficients can be captured. The MVAR model is fitted to the shuffled time series and TPDC is estimated. The bootstrapping is performed 1000 times and the average TPDC value is taken as the significance threshold for all our connections. This process is performed separately for each participant. The resulting value is the significance threshold value for all our connections. This process is performed separately for each subject. In this study the open source Matlab package autoregressive fit (ARFIT) 21,22 was used for estimating the autoregressive coefficients from the spatially filtered source signals of the identified nodes in the three new communities. We applied time reversal technique TRT 23 as a second significance test on the connections already identified by TPDC using a data-driven bootstrapping surrogate significance test.

Threat processing related inter-trial phase coherence changes
The frontal theta showed a significant inter-trial phase coherence (ITPC) increase in between the baseline (-250-0 ms) and the four subsequent time windows from (0 to 1000 ms) (experiment 1 p < 0.001; experiment 2 p < 0.001; supplementary Fig. 1). The occipital alpha showed decreased ITPC between the two conditions (CS+ and CS-) but the difference was only significant between the baseline and the subsequent three time windows from (T1 to T3, 0-750 ms) in both experiment 1 (p < 0.001) and experiment 2 (p < 0.001). The ITPC in the theta frequency showed an inverted U-shape like temporal pattern, increase in the first two windows

TMS induced inter-trial phase coherence changes
In experiment 2, the ITPC at the theta band was perturbed by the TMS in the dmPFC at 1000 ms and induced a significant increase (p < 0.01) in the frontal theta but did not affect the occipital alpha (supplementary Fig. 1). We were able to modulate the ITPC using TMS at the dmPFC in the frontal lobe during the time interval of threat processing. The control experiment of TMS at 80 ms did not induce any significant change in the frontal lobe indicating choosing the correct temporal window for perturbation is vital.

Topological characteristics of dynamic community alliance in the alpha band
In central executive network (CEN), the alpha band did not exhibit any significant differences in neither experiment nor analyzed network measures.
For the salience network (SN) (Supp. Fig. 2) no significant decreases with respect to baseline were detected for T1 to T4 in neither experiment. However, flexibility significantly decreased between T5 and T4 in both experiments (p < 0.001). The clustering coefficient also showed no changes in respect to baseline in neither experiment. However, a decrease in clustering coefficient in T5 in comparison to T4 was significant in both experiments (p < 0.001). For local and global efficiency, no significant changes with respect to baseline were detected in neither experiment. The decrease in global efficiency in T5 compared to T4 was significant in both experiments (p < 0.001). The decrease in local efficiency in T5 compared to T4 was significant in both experiments (p < 0.001).

Control experiment 3 with TMS at 80 milliseconds
The three network parameters were estimated for the control experiment as described in the main manuscript. The parameters did not differ significantly (p > 0.05 for all three parameters between experiment 1 (without TMS) and experiment 3 (with TMS at 80 ms) ( Supplementary   Fig. 7). However, the T5 window (1000-1250 ms) differed significantly between the experiment 2 (with TMS) and experiment 3 for all the three parameters (p < 0.01) (refer to Supplementary   Fig. 1). The control experiment results indicated that the TMS at 80 ms did not change the course of the dynamic network adaptation. On the other hand, this experiment validates that the correct temporal window for TMS is 1000 ms to modulate the network response in each of the three newly formed communities.

Inter-trial coherence as a substrate of threat processing
We found an increase of the theta inter-trial phase coherence in the frontal lobe and a simultaneous decrease in the inter-trial phase coherence for the occipital alpha, relative to stimuli presentation. This suggests that threat processing is not a purely autonomous response to stimulus presentation, but rather that it facilitates interactions between regions 24, 25 and for specific temporal conditioned stimuli. 26 Previous studies on memory function have shown that there is a relationship between brain response to external stressors and the phase of the synchronized oscillations, 27, 28 which can be prolonged by exciting a small number of neurons 29 that participate in such oscillatory behavior. Low-frequency oscillations, such as theta (4)(5)(6)(7) and alpha (8)(9)(10)(11)(12) can be recorded in different specific anatomical regions and especially facilitate communication between hippocampus, 30 amygdala, 31 and prefrontal cortex. The specificity of attention in instructed threat studies suggests that these oscillations provide a temporal window for inter-regional communication. 32 Intra-regional functional communication has been found for interactions involving the fronto-occipital circuit during directed attention to visual stimuli. 33,34 It has also previously been shown that frontal theta phase consistency reflects coordination of information transfer between distant brain areas. 35,36 On the contrary, the decrease in the inter-trial phase for the alpha oscillations in the occipital lobe could be related to alterations due to pre-visual threat processing. Earlier studies 37,38 have shown disruption in phase consistency over successive trials in the occipital lobe, which suggest that the inter-trial coherence of these oscillations drives the physiological response during instructed threat processing. Our data support this hypothesis and localize it differentially in both the frontal and occipital lobes. Here we use the dynamics of theta driven alterations for the application of TMS pulses to dmPFC, and apply TMS at a time relevant for threat processing (1000 ms) and a time point before this (80 ms). Using this approach, we were able to achieve different effects at the network behavior level. Specifically, a TMS pulse before active processing leads to a community independent increase of network flexibility and clustering without a preservation of inter-network interactions. A TMS pulse applied during a time point physiologically relevant for processing mirrors the network response and intercommunity information transfer.

Validation using different brain atlas
To evaluate if the results obtained in the study are dependent on the choice of the brain atlas, we performed the community analysis and computed a network measure -Flexibility using a brain atlas with 100 regions of interests 39 mapped to seven functional networks. 40 We showed that the results are replicable (suppemmentary Fig. 8) for three functional networks dorsal attention, salience ventral attention and default mode corresponding to the ones we used in the study -central executive, salience and default mode network. There was a significant increase in flexibility for the three networks respectively compared to baseline all the other following windows (factors condition (F1,18 = 8.48, p = 0.01; F1,18 = 10.15, p < 0.001; F1,18 = 9.28, p < 0.001) and time (F6,108 = 3.68, p = 0.009; F6,108 = 5.09, p = 0.004; F6,108 = 6.42, p = 0.001)), but increased flexibility in experiment 2 (factors condition (F1,25 = 5.45, p = 0.03; F1,25 = 5.47, p = 0.004; F1,25 = 5.20, p = 0.003) and time (F6,150 = 2.86, p = 0.01; F6,150 = 4.89, p = 0.007; F6,150 = 6.29, p = 0.003)). We also found, replicating our findings, a significant increase in the flexibility for the windows at T5 compared to T4 in both experiments (p < 0.001) for the two networks dorsal attention and salience ventral attention. In the default mode network, the flexibility decreased as in our original analyses for the theta frequency band between the windows comparing T5 and T4 for both experiments (p < 0.001). From these analyses, it is clear the network architecture shows similar dynamic processing (Flexibility) irrespective of the used atlas in this study.