Probing neural networks for dynamic switches of communication pathways

Dynamic communication and routing play important roles in the human brain in order to facilitate flexibility in task solving and thought processes. Here, we present a network perturbation methodology that allows investigating dynamic switching between different network pathways based on phase offsets between two external oscillatory drivers. We apply this method in a computational model of the human connectome with delay-coupled neural masses. To analyze dynamic switching of pathways, we define four new metrics that measure dynamic network response properties for pairs of stimulated nodes. Evaluating these metrics for all network pathways, we found a broad spectrum of pathways with distinct dynamic properties and switching behaviors. We show that network pathways can have characteristic timescales and thus specific preferences for the phase lag between the regions they connect. Specifically, we identified pairs of network nodes whose connecting paths can either be (1) insensitive to the phase relationship between the node pair, (2) turned on and off via changes in the phase relationship between the node pair, or (3) switched between via changes in the phase relationship between the node pair. Regarding the latter, we found that 33% of node pairs can switch their communication from one pathway to another depending on their phase offsets. This reveals a potential mechanistic role that phase offsets and coupling delays might play for the dynamic information routing via communication pathways in the brain.


Introduction
Over the past decades it has been shown that the brain, facing a specific task or not, exhibits well-structured functional connectivity [1][2][3][4][5]. This has been specifically investigated for resting-state networks [6][7][8][9], but also for other networks when the brain is performing different tasks [8,10]. These findings lead to the idea that networks of dynamically coupled neural populations reflect an inherent functional organization of the brain which is optimized to perform a wide range of tasks it encounters frequently [11,12]. If faced with a task that requires synchronization between brain areas not typically coupled at rest, this organization has to be altered temporarily in order to perform that task efficiently [13][14][15].
One mechanism that has been proposed for such dynamic re-organization is top-down suppression of distractor areas (or task-irrelevant areas) by a slow, oscillatory rhythm compared to the neural rhythm at which the signal is being processed [16][17][18]. While this suppression mechanism acts on individual nodes in the network, dynamic re-organization by changes in the phase relationships could act directly on the edges of the network by modulating the effective connectivity between nodes [19][20][21]. Here, we focus on the latter and consider the dynamic re-organization as a self-organizing property of the network that flexibly adapts functional coupling to changes in extrinsic stimulation. According to the communicationthrough-coherence theory [19,20], the phase relationship between two neural oscillators affects their success in communication. A neural population which tries to establish communication with other target populations at a certain frequency is more successful if the signals arrive at the excitable phase of the target populations. In other words, successful communication between two neural oscillators is reflected by their amount of entrainment. The functional coupling pattern that arises from this process could be temporarily altered by advancing or slowing down the phase of each of the neural oscillators. This way, changes in phase relationships between brain areas could serve as a mechanism for the re-organization of effective connectivity.
Within a complex network like the human brain, multiple structural pathways exist between most pairs of nodes given a sufficiently high spatial resolution. Since functional coupling patterns are established via synchronization along such pathways, we set out to identify general principles of how these pathways interact with each other during synchronization processes. Several studies have placed emphasis on the importance of information transmission delay for synchronization processes as well as its role in the formation of functional clusters in the brain [5,20,[22][23][24][25][26][27][28][29][30][31][32]. In their fundamental work on the dynamic implications of discrete delays in networks of rate-based neurons, Roxin and colleagues were able to demonstrate that the incorporation of delays enriches the dynamic repertoire of the system, allowing for various periodic and aperiodic behaviors [33,34]. Using connectome models composed of delay-coupled Kuramoto oscillators, Petkoski and colleagues showed that delay distributions in heterogeneously coupled networks can influence whether nodes tend to synchronize in phase or anti-phase [32]. Furthermore, Gosh and colleagues showed in a connectome model of coupled Fitz-Hugh Nagumo systems that the incorporation of realistic network delays allowed the noise-driven model to show a dynamic repertoire similar to experimentally reported dynamics of resting-state functional connectivities [35]. Here, we investigate the interplay between coupling delays and synchronization phase offsets in a model of the human connectome. In the case of two delay-coupled neural oscillators, a tight relationship between their optimal phase relationship and coupling delay can be expected [20]. Given a complex network with multiple communication pathways between node pairs, this raises the question of whether the set of active communication pathways depends on the phase offset between the nodes. Finding a difference in this set of utilized communication pathways over varying phase offset strongly implies that communication pathways express preferences for the phase relationship between the nodes they connect. These preferences will be influenced by the time delays of the communication pathways, which are determined by axonal signal transmission delays as well as rise and decay properties of the post-synaptic response. Here, we focus on the former and assume the latter to be relatively homogeneous across the cortico-cortical connections considered in our model, since they all resemble glutamatergic synaptic connections. In particular, we predict that two brain regions trying to communicate at a certain frequency with a given phase offset will use only a fraction of their available communication paths. Furthermore, we predict the selection of communication paths to be influenced by their interaction time delays and the phase offset at which the communication is initially attempted.
To test these hypotheses, we investigate cortico-cortical communication processes via a computational whole-brain network model of delay-coupled, oscillating nodes. We introduce a set-up with two oscillatory drivers to probe the network for changes in interactions between the stimulated pair of nodes. Importantly, our stimulation approach relies on the entrainment of the given pair of nodes to oscillate at the same frequency, but with a certain phase lag relative to each other. While Fig 1A illustrates the overall extrinsic stimulation setup, Fig 1B visualizes why it is necessary to stimulate at different phase lags between the stimulation signals. In particular, Fig 1B shows that even in the absence of any interaction through the network, there is the externally induced trivial coherence between the two stimulated nodes. Thus, the coherence is measured for many different stimulation phase offsets and the measurement with the lowest coherence is chosen as the baseline. Any deviation in the coherence from this baseline can be attributed to induced changes in the coupling between the two stimulated nodes through the network. As indicated in Fig 1C and 1D, this may happen due to a switching in the communication pathways. Comparing the coherence along different pathways over different stimulation phase offsets then reveals the phase preferences for different routes. We The simultaneous stimulation at the two nodes could induce a trivial coherence (green) between the two nodes even in the absence of any interaction between the two nodes through the network. (C) In this example a stimulation with no phase offset induces a network interaction between the two stimulated nodes through path 1 (blue). (D) In contrast, when stimulating with a phase offset of π another path is activated (red). speculate that these differences in phase preferences at different pathways could act as a switching and gating mechanism used by the brain to establish communication between remote brain areas when needed. Our method allows the investigation of these mechanisms by probing the network for these dynamic switches in communication pathways.

Computational model
Our computational model (Fig 2) is based on the widely used Jansen-Rit neural mass model [36] which employs a mean-field approach to model the interaction between a pyramidal cell population (green), an excitatory interneuron population (blue), and an inhibitory interneuron population (red), as illustrated with the relevant equations in the Figure. The function σ(V) transforms the average membrane potentials to average firing rates via a parameterized sigmoid (depicted in cyan). The standard parametrization originally proposed by Jansen and Rit reflects cortical oscillatory activity in the alpha frequency band. These parameters were chosen based on experimental findings in the neuroscience literature and are reported in the Materials and Methods section. Since the aim of this article is the investigation of the effect of pathway time delays on neural synchronization processes and not the effect of information processing delays at the node level (arising from time constants of cell membranes, synapses etc.), we decided to use this standard parametrization for each node in our network [36]. For the purpose of investigating networks of extrinsically perturbed Jansen-Rit models, the following two extensions were added: First, we coupled multiple Jansen-Rit nodes via delayed, weighted connections between their pyramidal cell populations (depicted in yellow). In our model this pyramidal cell population is interpreted as a lumped representation of all projection cells within a modeled brain area, whereas the excitatory and inhibitory interneuron populations describe local feedback loops of that pyramidal cell population. We avoid the original layer-specific (infragranular, granular, and supragranular) interpretation offered by Jansen and Rit [36] because inter-areal connectivities would be subject to layer-specific profiles [37] which would be hard to estimate reliably for a whole-brain network model. Note that this approach is in accordance with previous simulation studies of networks of multiple coupled Jansen-Rit models that employed pyramidal to pyramidal cell coupling as well [38][39][40][41][42]. Secondly, weak external drivers were applied at two stimulation sites influencing the average membrane potential of the pyramidal cells with phase offset Δφ between the two drivers (depicted in purple). By applying the stimulation to the pyramidal cells, it enters at the same population which also receives the network input. It can, thus, either be conceived as input of a neural origin not explicitly incorporated in our model or as an extrinsic input that mainly affects the excitability of the pyramidal cells.

Results
We will first present our results of the simulations in simple networks of only 2 or 3 nodes. Specifically, we show how the coherence between a pair of stimulated nodes depends on the transmission delay of their connections and on the relative phase offsets between the two stimulation signals. Subsequently, we move on to a model of the cortico-cortical synchronization processes within a single hemisphere of the human brain. Cortico-cortical coupling strengths and delays were informed by an approximation of the human connectome obtained from diffusion tensor imaging (DTI) data. Additionally, we fit the model to match functional connectivities calculated from electroencephalography (EEG) recordings. In this model, we show how the coherence between stimulated nodes changes over phase lags and how this effect relates to the connectedness and distance of the nodes. In a final step, we identify the pathways responsible for the interaction between the stimulated nodes, analyze their phase lag preferences and identify cases of phase-related switching between pathways. For all simulations, we use the computational model defined in the previous section.

Simple models with 2 or 3 nodes
The idea behind the extrinsic stimulation approach can be well explained using a simple toymodel of 2 directly coupled Jansen-Rit nodes, where each node is stimulated with a f ext = 11 Hz sinusoidal signal with strength c ext = 0.1 mV. Fig 3 shows the coherence between the driven nodes for systematic changes in the phase offset between the stimuli and the distance between the coupled nodes. While uni-directionally coupled nodes can have preferences for any stimulation phase offset, as shown in Fig 3A, bi-directionally coupled nodes are more susceptible for stimulation at in-or anti-phase (see Fig 3B). The latter is in line with other studies that showed neural synchronization of remote, delay-connected neurons or neural populations to typically appear at in-or anti-phase [31,32,43]. This shows that the communication between coupled pairs of nodes can be modulated by stimulation and that communication pathways can have characteristic stimulation phase offset preferences, depending on their length [44].
To quantify the modulation of communication, we define the pathway-synchronizationfacilitation (PSF), measuring for a given pair of weakly stimulated network nodes k i and k j how their interaction is dependent on specific phase offsets: where coh(k i , k j , Δφ) is the coherence between network nodes k i and k j for stimulation phase offset Δφ. The PSF is high for node pairs if their coherence is high for one stimulation phase offset and low for another, i.e., the relative phase of the stimulation at the two sites matters strongly. The PSF curves in Fig 3A and 3B show that in both cases there is a PSF effect (PSF > 0) and in the case of bi-directionally coupled nodes the strength of this effect depends on the distance between the nodes. To extend this idea to communication via indirect pathways, we investigated synchronization between 2 nodes connected only indirectly via a third intermediate node. We used bidirectional couplings for both connections and both end nodes were stimulated as described previously. As can be seen in Fig 3C and 3D, the interaction between the two weakly stimulated nodes not only depended on the length of the connection, but mostly on the relative position of the third node on the indirect path. Thus, phase lag preferences of pathways connecting two target nodes arise primarily from the distance between intermittent nodes along the pathway. In a heterogeneous network like the connectome, we hence expect to find node pairs connected by pathways with distinct phase lag preferences. Depending on the stimulation phase offset at the target nodes, we expect to be able to induce a switch in the pathway the nodes employ to entrain each other.

Connectome model without atimulation
As shown in the previous section, the coherence in a network of only three nodes can already exhibit very complex dependencies on the stimulation phase offset. To extend our analysis to network communication patterns in the case of a complex network with multiple competing pathways, we used a model of 33 delay-connected nodes, representing one hemisphere of the human connectome [30]. The structural connectivity matrix was obtained from DTI-based tractography data as described in the Materials and Methods section. Fig 4A shows the sparse connectivity matrix C mj used to connect the 33 regions and Fig 4B the corresponding average fiber lengths L mj between region pairs. These two matrices are used for connecting the different populations of excitatory pyramidal cells (green in Fig 2) in a delay-coupled network of a single hemisphere.
For the same 33 regions, EEG resting-state recordings from the same subjects were used to calculate pairwise coherences in the 10 Hz range as shown in Fig 4C. To resemble this empirical resting-state functional connectivity, we used a grid-search over the two remaining free parameters that were not defined in [36]: the global scaling of the connection weights c net and the global velocity v that is used to scale the distance data resulting in pair-wise signal transmission delays. We simulated the 33 connected neural-mass-models and processed the timeseries of the pyramidal cells in the same way as the EEG data, i.e. we band-pass filtered the Probing neural networks for dynamic switches of communication pathways time series at 10 Hz, applied the Hilbert transform, and calculated the pair-wise coherence. This yielded a 33 x 33 functional connectivity matrix which we compared to the empirical functional connectivity by calculating the Pearson correlation coefficient.
The selection of parameters was based on the rationale to match the functional connectivity observed in the network model as well as possible to empirical EEG-based functional connectivity from human subjects. By fitting the velocity, we ensured that our pathway delays reflect realistic, empirically observed timescales of cortico-cortical interactions ( Fig 4D). We selected the parameters c net = 14 and v = 2.6 m/s for subsequent analyses which show a high correlation to the empirical functional connectivity (r = 0.64, p < .0001). Notice that this correlation is comparable with values of other bottom-up models reported in the literature [30,45], which is remarkable considering that we set a substantial amount of structural connections to 0. The simulated functional connectivity using these selected parameters is shown in Fig 4E. A comparison of these different connectivity matrices (Fig 4A, 4B, 4C and 4E) revealed that the empirical functional connectivity correlated more strongly with the inter-regional fiber distances (r CB = −0.70) than with the structural connectivity (r CA = 0.55), while the simulated functional connectivity correlated more strongly with structural connectivity (r EA = 0.78) than with the inter-regional fiber distances (r EB = −0.71).
As a last analysis of the unperturbed network, we evaluated the power spectral density of each network node ( Fig 4F). All 33 nodes showed a frequency peak around 11 Hz.

Connectome model with stimulation
Based on this model of cortical activity, we used the stimulation approach to investigate how pathways facilitate synchronization between network nodes at certain phase lags between the nodes. Specifically, we weakly stimulated different pairs of cortical regions with varying phase offsets between the two stimulation signals while measuring the coherence between the stimulated nodes at each phase offset. As argued above, finding differences in the coherence over stimulation phase offsets would indicate phase-specific communication modulation between the stimulated nodes. Before analyzing PSF effects in the connectome model, it was necessary to determine the optimal stimulation frequency and strength for this model, because complex synchronization patterns can arise in networks of delay-coupled non-linear elements [33], especially when perturbed extrinsically [41]. This was performed in two steps. First, we determined the optimal frequency by extrinsic stimulation of a single network node. Second, we determined the optimal stimulation strength by stimulation of two network nodes at that frequency.
Since our main analysis will focus on coupling effects through different network paths between two stimulated nodes, we aimed to find the stimulation frequency at which the extrinsic stimulation has a high penetration depth into the network (as measured by the mean coherence to all network nodes). Notice that the stimulation frequency with the optimal network penetration depth may not necessarily coincide with the optimal stimulation frequency to entrain the stimulated network node. To determine the optimal stimulation frequency, we stimulated a single region in our network with a stimulus of varying frequency (4-30 Hz) and strength (0.03-4 mV) while evaluating the mean coherence between the stimulus and all nodes in the network. As Fig 5A shows, this average coherence to the full network shows a strong peak at 11 Hz, which is the intrinsic frequency of the unperturbed network (see Fig 4F). The mean coherence of the stimulation signal to the stimulated node and its direct neighbours can be observed in Fig 5B, while the coherence to only the stimulated node can be observed in Fig  5C (always averaged over simulations in which the extrinsic stimulation was applied to each of the 33 regions). The latter reflects the well-studied relationship between a driver and an oscillator described by the so-called Arnold tongue [46], except that this coherence to the stimulated node was higher at driving frequencies that are different to the intrinsic network frequency of 11 Hz (see Fig 5C). This reflects the competition between entrainment via the extrinsic stimulation and via the network. At the intrinsic frequency of the network (11 Hz) and its harmonics, the driver has to compete with the network in entraining the stimulated node at a certain phase relationship. This competition decays with the deviation of the driving frequency from the networks' intrinsic frequency and its harmonics. Thus, the Arnold tongue of a single driven node embedded in the connectome model shows coherence peaks shifted away from its intrinsic frequency. However, as can be seen in Fig 5B, entrainment of the direct neighbours of the stimulated node is still strongest at the networks' intrinsic frequency, which explains the mismatch between the optimal stimulation frequency for penetrating the network (Fig 5A) vs. entraining a single node (Fig 5C). To ensure that the network is not pushed into a fully synchronized state by the extrinsic perturbation, we also evaluated the Kuramoto order parameter of the network (see Fig 5D). This metric reflects the level of synchronization in the network, with a value of 1 implying full synchronization and a value of 0 implying a fully asynchronous state. As can be seen in Fig 5D, the Kuramoto order parameter varies between 0.1 and 0.25 across extrinsic stimulation parametrizations and expresses a relatively low value in this range for an 11 Hz driver. Based on these results, we set the frequency of our stimulus to 11 Hz, at which the network (and not only the directly stimulated node) was most susceptible for entrainment by an external stimulation, but was not driven into a fully synchronized state.
In a second step, we stimulated pairs of nodes with 11 Hz stimuli. We varied the stimulus strength (0.25-1 mV) and the relative phase offset between the stimuli (0 − 2π) while evaluating the coherence between the stimulated nodes. All other parameters were chosen to be the same as for the previous simulation. The variability in the coherence between stimulated region pairs that we observed over stimulation phase offsets (as depicted for 2 example region pairs in Fig 6A and 6B) shows that the stimulated region pairs interacted with each other and that the interactions expressed a characteristic profile of phase offset preferences. As can be seen in Fig 6A and 6B, the variance of the coherence over phase offsets depended on the stimulation strength. Based on visual inspection of the coherence patterns of 20 different region pairs, we chose our stimulus scaling to be c ext = 0.5 mV, leaving the variance of the post-synaptic potential of the neural masses in a biologically plausible range and such that the external driver is relatively weak in comparison to the internal network dynamics (the membrane potential fluctuations of a single Jansen-Rit node are in the order of 1-10 mV). This gave us the final set of global model parameters which were used throughout all subsequent simulations and are reported in Table 1.
To statistically confirm the variance in the coherence between stimulated region pairs over phase offsets, we ran simulations with subsequent stimulations of each possible node pair in  the network. Again, we varied the phase offset between the two stimuli (16 equally spaced phase offsets between 0 and 2π) and evaluated the coherence between the stimulated nodes for each phase offset. Subsequently, those coherences were used to calculate the PSF for each region pair as defined in Eq (1). Using a one-sample t-test, we found the PSF effect to be significantly larger than zero (mean = 0.10021, CI = [0.092875, 0.10755], t = 26.8268, p < .0001). Hence, we were able to show with our extrinsic stimulation approach that pathways facilitated synchronization between cortical nodes and that the facilitatory strength depended on the phase lag between the region's average PSPs.
With the PSF effect established, we continued by investigating its dependence on certain features of the underlying structural connectivity graph. For this purpose, we searched for all possible pathways between each pair of stimulated nodes based on the structural connectivity matrix reported in Fig 4A. Since every stimulated pair of nodes was connected by at least one path via at most 6 edges, we restricted the search to pathways including 6 edges at maximum. With these pathways at hand, we started out by evaluating how the PSF effect changed with increasing network distance. An analysis of variance showed that the effect of shortest path length (minimum number of edges seperating a pair) on log(PSF) was significant, F(5,521) = 75.22, p < .0001. As can be seen in Fig 6C, we observed the trend that the PSF effect decreases with the number of nodes separating the stimulated nodes. Furthermore, as depicted in Fig  6D, this trend was supported by a significant correlation between the PSF effect and the length of the shortest pathway between the stimulated nodes (r = -0.57, p < .0001), a measurement that is strongly related to both interregional distance and minimal number of separating edges. Thus, we conclude that there is a tendency for a decrease in the interaction of stimulated node pairs with increasing network distance, which can be measured either as the number of interposed edges or as the summed up length of the edges of the shortest pathway connecting the nodes.
Next, we investigated the dependence of the PSF effect on the connectedness between the stimulated nodes. In this regard, we characterize the connectedness either by the connection strength along the shortest path (evaluated in Fig 6E) or by the number of different pathways connecting the two nodes (evaluated in Fig 6F). In the first case, we found a high correlation between the Dijkstra distance, calculated from the inverted connection strengths 1/C mj , and the PSF (r = -0.73, p < .0001) which is depicted in Fig 6E. This confirms previous findings of the dependence of the (directed) interactions between delay-coupled oscillators on their connection strength [29,32].
In the second case (Fig 6F), we measure the connectedness by the number of different pathways connecting the two nodes. It was necessary to limit this analysis to include only pathways with up to 6 edges, because of the combinatorial growth of the number of possible pathways with more edges. In other words, considering pathways with more edges could lead to scenarios in which node pairs show a large number of connected paths, even though those pathways all have to pass through many intermediate nodes and might not contribute to the interaction of these nodes through the network at all. Thus, to render this measure of connectedness less noisy, we decided to exclude paths with more than 6 edges. An analysis of variance showed that the effect of the number of connecting paths (only counting paths with 5 edges or less, and all nodes with more than 5 connecting paths were pooled into one level) on log(PSF) was significant, F(5,501) = 7.42, p < .0001. Together, these results (Figs 6 and 3) reveal how the variance of the interaction of network nodes over different phase lags (as measured by the PFS) depends both on the strength and delay of their connecting pathways.

Evaluation of pathway activation
To identify, how different pathways between stimulated node pairs contribute to the PSF effect, we next aimed to quantify how much a particular pathway was involved in the node interaction at a given phase offset. For this analysis, we define the pathway activation (PA) for a pathway through n nodes k i with i = 1‥n at a phase offset Δφ as the product of the pairwise coherences between neighboring pathway nodes: In other words, if communication fails at any point along a pathway, leading to a reduced coherence between the involved nodes, this is considered to be a bottleneck for the information flowing through that pathway. This is similar to the definition of efficiency of neural communication via communication pathways provided in [21]. Thus, for the remainder of this article, we will use changes in information flow and changes in pathway activation interchangeably. Furthermore, this metric ensures that short pathways are preferred over longer pathways, since the coherence is bound to the interval between 0 and 1. We evaluated the pathway activation (PA) for all pathways of up to n = 5 nodes connecting a given pair of stimulated nodes for all stimulation phase offsets. Doing this for each stimulated node pair, we found different classes of pathway interactions: Some pairs show only a very small selectivity for the stimulation phase offset (Fig 7A), while other node pairs were connected by paths with PA values with a strong dependence on the phase offset (Fig 7B and 7C). Moreover, some of these node pairs switched their interaction between different pathways depending on the stimulation phase offset. This is shown in Fig 7C, while the pathways between which the switching occurs are visualized in Fig 7D and 7E.
There are several alternative ways to define the pathway activation. A similar analysis like in pathway activation did not change the overall results in a qualitative way, meaning that using any of these alternative definitions there are always all three types of node pairs in the network: 1) pairs with no phase selectivity in their strongest pathway (like in Fig 7A), 2) pairs with phase selectivity without switching (like in Fig 7B), 3) pairs with pathway switching (like in Fig 7C).
To further analyze how the communication via specific pathways depends on the stimulation phase offset, we define the pathway phase selectivity (PPS) of a pathway P 1 similar to the PSF as Pathways with relatively constant PA values for all stimulation phase offsets have a low PPS (example in Fig 7A), while pathways with a high variation in the PA values have a high PPS (example in Fig 7B and 7C). In general, the metric is bound to the interval between 0 and 1. The evaluation of PPS values for all node pairs results in the distribution shown in Fig 7F. The activation of pathways with a PPS close to 0 is very hard to influence with phase offsets between the nodes at their ends. However, the further the PPS grows towards 1, the more sensitive the PA is towards changes in these phase offsets. As indicated by the long tail of the PPS distribution, a significant amount of pathways in the connectome model show such phase preferences.
In the next step, we analyzed the relationship of pathway-specific phase preferences (as shown in Fig 7A-7C) to the phase preferences of the stimulated nodes (as shown in Fig 5C  and 5D). We chose the most active pathway per node pair (averaged over all stimulation phase offsets) and calculated the phase difference between the stimulation phase offset with the highest coherence between the stimulated nodes and the stimulation phase offset with the highest PA. The histogram of these phase differences is significantly different from uniform, χ 2 (15, N = 514) = 155.31, p < .001, and has a peak at 0 (Fig 7G). A similar analysis for the second strongest pathway (excluding all pathways with overlapping sections with the strongest path), results in a histogram that differs only slightly from a uniform distribution, χ 2 (15, N = 514) = 26.82, p = 0.03 (Fig 7H). Hence, we found that the pathway with the strongest PA shows a similar phase preference as the coherence between the two stimulated nodes. We take this as Probing neural networks for dynamic switches of communication pathways evidence that the interaction between the nodes through the network is strongly modulated by this pathway.
Finally, we quantified the process of switching between the strongest and second strongest pathway per node pair. To this end, we define the pathway switching index (PSI) between pathways P 1 and P 2 as The PSI is positive if the two pathways switch their activation depending on the stimulation phase offset, meaning that at one phase offset the first path is more active and at another phase offset the second path is more active. We found that 33% (170 of 514) of node pairs have a positive PSI between their non-overlapping strongest and second strongest pathways (Fig 7I). These results suggest that in this network of 33 nodes of the human connectome many node pairs have the capacity to switch their communication between at least two different pathways with a PA characteristic similar to the example shown in Fig 7C.

Discussion
We have carried out a computational study of cortico-cortical synchronization processes that strongly emphasizes the role of phase relationships for dynamic switches in communication pathways. In this process, we introduced a novel method to detect network interactions between pairs of cortical regions via an extrinsic stimulation scheme. Using our method, we were able to quantify the influence of different pathways on cortico-cortical coupling between all pairs of 33 brain regions. We could further identify the pathways those region pairs use to interact with each other. These pathways represent communication channels with distinct preferred interaction time lags. One of the main findings of our work is that the pathway activation depends on the phase lag between the nodes they connect. This finding is in line with the communication-throughcoherence theory that predicts neural communication to depend on oscillatory phase differences and has received support from various experimental results [20]. In a transcranial magnetic stimulation study, Elswijk and colleagues demonstrated the effect of the stimulation on primary motor cortex to depend on the oscillatory phase of the latter [47]. Similarly, Helfrich and colleagues found the performance in a visual detection task to be modulated by the phase of a 10 Hz transcranial alternating current stimulation applied to the parieto-occipital cortex [48]. Together, these experimental results support the idea of the communication-throughcoherence theory that the effectiveness of neural communication between a source and a target population is modulated by the oscillatory phase of the latter. In a study using multielectrode recordings in frontal eye field and area V4 of monkeys, Gregoriou and colleagues were able to show that neural populations at both recording sites synchronized at characteristic time lags when processing a stimulus in their joint receptive field [49]. Based on the communicationthrough-coherence theory this may be explained by the axonal and synaptic delays of the communication pathways between these two areas, which render synchronization at this time lag most effective. The work presented in this article contributes to this notion via its systematic investigation of the dependence between pathway delays and neural synchronization. Our observation that the synchronization of two brain regions via simultaneous extrinsic stimulation has a preference for particular phase offsets between the two stimulation signals demonstrates the mechanistic role of phase lags for neural communication within a computational model. We were able to extend this finding to individual communication pathways, i.e., we showed that these pathways can express more or less strong preferences for the phase lag between the brain regions they connect. Together with the results of our toy-model simulations, which showed a direct dependence of neural synchronization on propagation delays in various connectivity motifs, this emphasizes the importance of connection delays for pathway switching. Moreover, it suggests that the effectiveness of information exchange via certain pathways depends on the phase relationship between the communicating sites.
The functional significance of these phase lag preferences is underlined by our finding that different pathways may be employed at different stimulation phase offsets between the communicating regions. Exploiting this mechanism, brain regions could dynamically switch between communication pathways by alterations in their phase lags. This flexibility in the choice of pathways could reflect a major self-organization mechanism of the brain. Specifically, the processing of transient sensory events may trigger strong phase shifts at certain sites in the brain that render certain pathways inaccessible due to their inherent delays. Our results suggest that in such a situation, the brain could employ a different set of routes to still enable communication between populations that previously relied on the now inaccessible pathway. As suggested in [12,50], flexible switching between different processing pathways could also allow for the dynamic binding of remote neural representations into different combined representations. Both external stimuli and network intrinsic signals could act as phase-resetting mechanisms at the communicating sites and thus switch from one pathway to another within a few oscillatory cycles [51].
In our study, we considered the case of simultaneous, sinusoidal stimulation of two network nodes with a relative phase offset between the stimulation signals. In the following, we discuss how such a phase offset could arise in a biological network without external stimulation. On the one hand, a network internal stimulation could originate from a common neural driver that projects to several target sites with slightly different delays. Our results suggest that such delays could favour some communication pathways between these target sites over others. Thus, the common neural driver would not only be able to project its information to its target sites, but also choose certain pathways via which information about the input is processed. Given different neural drivers that project to the same target sites with distinct delays, different routes for integrating the information could be associated with the driving delays. On the other hand, the stimulation could also originate from an external sensory stimulus. In this scenario, different features of the stimulus are processed along different neural pathways that will converge along the sensory hierarchy. Our results suggest that the convergence of these pathways is sensitive to differences in the processing delays across pathways. Differences in processing delays could for example arise from variations in the prominence of certain stimulus features for which a certain pathway is sensitive. Depending on the prominence of the feature, the neural responses would be more or less synchronized which directly translates to differences in the processing delay further up the processing hierarchy [21]. Alternatively, the feature prominence could affect the processing delay by advancing the phase of an oscillation, by increasing the oscillation frequency, or by reducing the time to entrain other oscillators [20,21]. This has for example been demonstrated for neural signal transmission between primary and secondary visual cortex in macaque monkeys [52]. Along these lines, Voloh and Womelsdorf provide a convincing overview with respect to different experimental scenarios in which stimulus-triggered phase resets lead to re-organizations of functional networks [53].
Different scenarios could be imagined in which input-dependent changes in the phase of a specific node affect the pattern of active output pathways of that node. One example would be the scenario described in [19,20], where two bottom up neural drivers encoding different stimuli try to entrain a target population at different phase lags. While our results suggest that entrainment at different phase lags could affect the pathways along which the winning stimulus is processed further, this scenario of competing input stimuli is not explicitly considered in this study. Thus, future work could investigate pathway switching due to transient inputs delivered at different phases of the ongoing oscillations. To this end, a major challenge will be the required online read-out of the phase of the source node to correctly time the input. This scenario would allow to further test our proposal of dynamic switching of communication pathways in both computational models and experimental setups. A potential candidate for the latter would be brain stimulation studies, employing for example transcranial stimulation or optogenetics in combination with online recordings of neural activity [54,55]. For a review on how such phase shifts can be induced in models of single cells and networks of interconnected cells, see [54].
We would like to propose two further future directions along which our line of research could be extended and at the same time, address two important simplifications in our model. First, we only considered discrete delays in our model, resembling average lengths of DTI-tractography based cortico-cortical tracts. However, from experimental studies it is known that long-range interactions in cortex express less sharp delay profiles [56]. Indeed, computational studies have found that the variance of delay distributions affects synchronization properties of delay-coupled systems in various ways [57,58]. Thus, a systematic study of how stimulus phase-dependent effects on cortical routing (as reported in this study) may change with the variance of distributed delays could shed light on the validity of our results for biologically more realistic delay-coupling profiles.
Another simplification of our study concerns the Jansen-Rit model we used to represent the cortical dynamics at each of our network nodes. This model was developed to describe dynamics of average membrane potential fluctuations in three interconnected cell populations. However, the model equations have not been derived from a single cell representation of such a network [36]. As such, the model provides limited information about the underlying single cell states, such as the within-population synchronization degrees. This means that the model does not allow investigating the dependence of its synchronization with other network nodes on the degree of neural synchronization within its underlying neural populations. However, how strongly the excitability of a population varies within its oscillatory cycle clearly depends on this degree of within-population synchronization [20]. Recently developed mathematical descriptions of the macroscopic dynamics of networks of spiking neurons have overcome these problems and allow deriving the synchronization of the neural population directly from its mean-field description [59]. Unfortunately, models of cortical microcircuits have yet to be developed for these novel mean-field descriptions. With such models at hand, it would be feasible to investigate interactions between within-and betweenpopulation synchronization processes at a macroscopic scale and to investigate whether communication pathways in delay-connected networks express similar phase preference profiles as we found in this study.
Taken together, our results suggest a potential mechanism the human brain might have developed to use the physiological constraints imposed by coupling delays to its computational advantage. We have laid out several future directions of research that could help to advance and confirm these findings. From an experimental point of view, the stimulation method applied in our computational model could guide the development of brain stimulation protocols to probe dynamic switching between communication pathways in the brain. Additionally, our stimulation method and graph metrics are applicable in future theoretical studies characterizing the dynamic properties of network graphs.

Materials and methods
In the following, we present the detailed parameters of the neural mass model, the pre-processing of the structural connectivity, and the functional connectivity.

Neural mass model
In the Jansen-Rit model, signal transmission between cell populations is realized via a convolution of average pre-synaptic firing rates with a post-synaptic response kernel. These convolution operations are mathematically described by the coupled ordinary differential equations in Fig 2 (each line describes a synaptic convolution operation) and turn pre-synaptic firing rates into post-synaptic potentials. The simple exponential form of the response kernel is in line with empirically measured post-synaptic responses [60] and provides a sufficient approximation of synaptic response time scales for our purposes. To translate post-synaptic responses back into firing rates, an instantaneous sigmoidal transform is used as shown in Fig 2 which introduces non-linearity to the model behavior.
The parameters of the neural mass model are shown in Table 2.
They reflect the standard parametrization proposed by Jansen and Rit for modeling waxing-and-waning cortical alpha oscillations at around 10 Hz [36]. While some of those parameters (such as the connection strength ratios) were informed by experimental findings, others can be considered free parameters in bio-physiologically plausible ranges. The dynamic dependence of the Jansen-Rit model on its parameters has been investigated in various theoretical studies [36,[61][62][63]. For our study, it is important that each network node is kept in an oscillating dynamic regime in which the oscillation amplitude is sensitive to the network input. Since the model has several control parameters that can lead to aperiodic network behavior (see [36,62,63]), we decided to stick with the standard parametrization proposed by Jansen and Rit (Table 2) for each node in our network.
Thus, the only model parameter that is varied in our simulations is the average input to the model arising from the background noise, network interactions and external stimulation. As shown by Spiegler and colleagues, the Jansen-Rit model undergoes various bifurcations due to alterations of the input strength that can drive the model into bi-stable periodic or aperiodic regimes [63]. To investigate whether such undesired states could be entered under the conditions imposed by our simulations, we calculated the minimum and maximum of the average input a single node may receive in our simulations. Besides the extrinsic stimulation signal, which is centered around 0 mV, there are two rate-based inputs at each node. (1) The sub-cortical noise is modeled by a Gaussian process with a mean of m ¼ 220:0 1 s . (2) The network input is ranging between 0 1 s and e 0 c ¼ 100 1 s . Feeding (1) and (2) into a single excitatory synapse with efficacy H = 3.25mV and τ = 0.01s, we find that it produces a synaptic input ranging between 7.5mV and 10mV. As shown by [63] the Jansen-Rit model expresses a mono-stable dynamic regime in this input range with a limit cycle being the only stable equilibrium. Variation of the input strength within this range mainly results in a modulation of the oscillation amplitude. Thus, the models considered in this study can conceptually be viewed as networks of (sparsely) coupled non-linear oscillators with discrete time delays.
If not reported otherwise, simulation results reported for those models refer to 16 minutes of simulated network behavior. We solved the model equations using the Runge-Kutta algorithm of 4th order (RK4) with an integration step-size of 1 ms. An analysis of the influence of the simulation step size on the accuracy of the simulation can be found in supporting material S1 Fig.

Structural connectivity and distance estimates
In the first step to building a bottom-up model of cortical activity, we needed to approximate the structural connections between different brain regions. As mentioned in the introduction, this can be done via DTI recordings. However, there are several technical limitations of the extent to which human SC can be approximated based on DTI, one of them being the systematic underestimation of inter-hemispheric connections [64,65]. Thus we decided to restrict our analysis to the cerebral cortex of a single hemisphere. To this end, we used the same structural imaging data, pre-processing and probabilistic tracking pipeline as reported by Finger et al. [30], but restricted subsequent processing to the 33 regions of interest (ROIs) of the left hemisphere. This data set included diffusion-and T1-weighted images acquired from 17 healthy subjects (7 female, age mean ± s.d. = 65.6y ± 10.9y) with a 3 Tesla Siemens Skyra MRI scanner (Siemens, Erlangen, Germany) and a 32-channel head coil. The 33 ROIs were registered individually for each subject based on the 'Desikan-Killiany' cortical atlas available in the Freesurfer toolbox (surfer.nmr.mgh.harvard.edu) [66]. The incoming connections to each region were normalized such that they summed up to 1. Since we were only interested in synchronization along indirect pathways, we needed some connections in our model to be strictly 0. Otherwise, it would be difficult to exclude potential synchronization along very weak direct connections. Hence, we chose to set all connections below a strength of 0.1 to zero. Afterwards, we re-normalized the input to each region such that they summed up to 1. The resulting SC matrix as well as the pair-wise distances calculated from the average tract lengths are visualized in Fig 4A and 4B.

Empirical resting-state functional connectivity
Based on those SC and distance information we aimed to build a model of cortical activity able to reflect empirically observed synchronization behavior. Thus we needed empirical observations of cortical activity to evaluate our model. For this purpose, we acquired EEG data from the same 17 subjects as described above. This was done with 63 cephalic active surface electrodes arranged according to the 10/10 system (actiCAP R Brain Products GmbH, Gilching, Germany) for eight minutes of eyes-open resting-state. Again, data acquisition and pre-processing followed the same procedure as reported by Finger et al. [30]. EEG time-series from the surface electrodes were projected onto the centers of the ROIs via a linear constraint minimum variance spatial beam former [67]. The resulting source-space signals were band-pass filtered at 10 Hz and turned into analytic signals using the Hilbert transform. Subsequently, functional connectivity was evaluated as the coherence between all pairs of ROIs [68]. This resulted in the 33 x 33 functional connectivity matrix that can be observed in Fig 4C in the main paper and served as optimization target for our model.