Inferring the physical connectivity of complex networks from their functional dynamics
© Ta et al. 2010
Received: 5 September 2009
Accepted: 26 May 2010
Published: 26 May 2010
Skip to main content
© Ta et al. 2010
Received: 5 September 2009
Accepted: 26 May 2010
Published: 26 May 2010
Biological networks, such as protein-protein interactions, metabolic, signalling, transcription-regulatory networks and neural synapses, are representations of large-scale dynamic systems. The relationship between the network structure and functions remains one of the central problems in current multidisciplinary research. Significant progress has been made toward understanding the implication of topological features for the network dynamics and functions, especially in biological networks. Given observations of a network system's behaviours or measurements of its functional dynamics, what can we conclude of the details of physical connectivity of the underlying structure?
We modelled the network system by employing a scale-free network of coupled phase oscillators. Pairwise phase coherence (PPC) was calculated for all the pairs of oscillators to present functional dynamics induced by the system. At the regime of global incoherence, we observed a Significant pairwise synchronization only between two nodes that are physically connected. Right after the onset of global synchronization, disconnected nodes begin to oscillate in a correlated fashion and the PPC of two nodes, either connected or disconnected, depends on their degrees.
Based on the observation of PPCs, we built a weighted network of synchronization (WNS), an all-to-all functionally connected network where each link is weighted by the PPC of two oscillators at the ends of the link. In the regime of strong coupling, we observed a Significant similarity in the organization of WNSs induced by systems sharing the same substrate network but different configurations of initial phases and intrinsic frequencies of oscillators.
We reconstruct physical network from the WNS by choosing the links whose weights are higher than a given threshold. We observed an optimal reconstruction just before the onset of global synchronization.
Finally, we correlated the topology of the background network to the observed change of the functional activities in the system.
The results presented in this study indicate a strong relationship between the structure and dynamics of complex network systems. As coupling strength increases, synchronization emerges among hub nodes and recruits small-degree nodes. The results show that the onset of global synchronization in the system hinders the reconstruction of an underlying complex structure. Our analysis helps to clarify how the synchronization is achieved in systems of different network topologies.
Biological networks, such as protein-protein interactions, metabolic, signalling, transcription-regulatory networks and neural synapses, are representations of large-scale dynamic systems. Network maps are graphical representations of dynamic systems in life. A network with a static connectivity is dynamic in the sense that the nodes implement so-called functional activities that evolve in time. In biological context, these activities may represent the concentration of a molecule, phosphorylation state of enzyme, expression level of a gene, depolarization of a neuron or circadian rhythm.
Despite the remarkable diversity of networks in nature, their architecture is governed by a few simple principles that are common to most networks of major scientific and technological interest [1, 2]. For decades network systems have been modeled either as chains, grids, lattices and fully-connected graphs which are completely regular or as random Erdős-Rényi network whose node degrees follow a Poisson distribution . However, a number of recent findings indicate that real networks including large communication systems [4–6], biological systems [7–9], and a variety of social interaction structures  are characterized by a power-law degree distribution, P(k) ∝ k -γ , where degree k is the number of neighbours of a given node. In these so-called scale-free (SF) networks, most of the nodes have only a few links, whereas a few nodes have a very large number of links, which are often called hubs .
Synchronization is a common nonlinear phenomenon in a broad class of systems ranging from physics and chemistry to biology and social sciences [12–14]. This is a dynamical process by which the activities of two or more individuals change coherently, almost in unison. For example, thousands of cells synchronize their activities to make our heart beat rhythmically. An excellent example of a realistic complex network system is the brain, where thousands of neurons fire synchronously to respond to external stimuli. A sudden and unexpected synchronization among a large population of neuronal units may cause some diseases , such as epileptic seizures. Therefore, understanding the path to synchronization can be of aid to diagnosing human disease.
Investigating how dynamical activities arise from complex network topology is of fundamental importance to understanding the functions of real-world systems [1, 16, 17]. Given a model of dynamical elements and the wiring diagram among them, what can we conclude of the dynamics such as synchronization arising from the model? Conversely, how can we infer the presence or absence of individual connections of a network system from observations of its behaviours or measurements of its functional dynamics? For example, how does an electroencephalography (EEG) pattern (which is the recording of electrical activities in different positions within the brain) reflect the details of the axons among cortical neurons?
Most previous studies on structure-dynamics relations analyzed the ability to obtain idealized complete synchronization or collective coherent oscillations, using regular networks with local coupling or in networks with global coupling . It has been shown that the response dynamics of a regular network system when applying external driving signals depends on both the driving and the network connectivity, yielding the ability to reconstruct the network connectivity . Recently, many works have focused on synchronization induced by more realistic complex network systems with the small-world property, the scale-free degree distribution, and the modular hierarchical structure [11, 19, 20]. The impact of these topological features of a network on the dynamical processes has been intensively studied. For example, it has been shown that the global synchronization is enhanced by the small-world property [21, 22] but hindered by the modularity of the network [23–25]. The evolution of global synchronization inside a network system depends much on the heterogeneity of the underlying network [26, 27].
In this paper, our aim is to elucidate how the network structure drives the functional dynamics, and consequently, analyze the ability to reconstruct the physical connectivity of an underlying network. We adopt a network of Kuramoto phase oscillators, a useful model to display a large variety of synchronization patterns while being sufficiently exible to be adapted to many realistic systems . The system was investigated comprehensively at both weak and strong regimes of coupling with the background network interpolating between regular and scale-free topologies. The collective behaviour of oscillators in this model is conventionally represented via the global order parameter which is the phase coherence of the population of oscillators  (see Methods). Evolution of the global order parameter as a function of coupling strength reflects the path from an incoherent to coherent state of the system. The global order parameter, however, fails to describe where the synchronization emerges and how it propagates inside the system. We calculated the pairwise phase correlation (PPC) for all pairs of oscillators and built a weighted network of synchronization (WNS), an all-to-all functionally connected network where each link is weighted by the PPC of two oscillators at the ends of the link (see Methods). PPC of two phase oscillators shows how dependent their motions are (see Methods). Therefore, the plot of PPC versus the coupling strength K and the product of node degrees helps to find which pairs of nodes synchronize first. The changes in the organization of WNS, as the coupling strength increases, help to explore how synchronization is achieved. To be able to observe if the dynamical process is purely driven by the background structure, we compare the WNSs obtained by the systems sharing a common substrate network but different initial configurations.
In spite of dealing with the reconstruction of physical connectivity as the work in  does, our work focuses on the reconstruction of physical connectivity from functional dynamics induced by a complex network system without applying any external inputs to the system. Namely, the links are predicted from an averaged WNS based on their weights (see Methods). Whereas the reconstruction of physical connectivity in  is implemented only in the stationary state of dynamics with regular background networks, our reconstruction is tested at both regimes of global incoherence and coherence with a family of networks having the same number of nodes and links but different topologies ranging from regular to scale-free. Interestingly, we observed a bimodal distribution of the links weights right before the onset of the global synchronization, irrespective of topologies of the background network. Moreover, the value by which two modes of the distribution are well separated can distinguish very successfully physical links from non-physical links, yielding an optimal reconstruction at this regime of weak coupling. Our findings are useful for practical applications because the weak-coupling regime is biologically more realistic.
The PPC between two nodes depends on the product of their degrees. As the coupling strength increases, the development of the structure of the skeleton of the WNS indicates that synchronization emerges among hub nodes and propagates into peripheral nodes in the system. This phenomenon can be explained as follows. Node i is considered to interact with the network via its effective coupling strength, (Figure 1d). When K is strong, r i approaches one (Figure 1c), is mostly proportional to the degree k i and, thus, hubs have stronger effective coupling. Therefore the hubs are synchronized more easily.
Figure 7b shows the evolution of the global order parameter as a function of coupling strength K for several network topologies ranging from regular to scale-free. The onset of global synchronization first occurs for the SF network. The more heterogeneous the network is, the smaller value of K is needed for the onset of the global synchronization. Conversely, the path to the complete synchronization is faster for networks with a homogeneous degree distribution. Whereas the nodes in homogenous network system suddenly become globally synchronized as the coupling strength increases, the path to synchronization of heterogeneous network system is nontrivial: synchronization starts among hubs and recruits low-degree nodes.
Figure 7d shows the behaviour of PPC for physically connected nodes (upper panel) and disconnected nodes (lower panel) in different types of networks. In the regime of global incoherence (K <K c ), for all network topologies, the mean PPC of physically connected nodes is Significant whereas that of two disconnected nodes is approximately zero, indicating that the dynamical process in the systems at a weak coupling is driven by physical connectivity. To clarify this point, we show the distribution of all links in terms of their weights in Figure 6. In this regime of weak coupling, we observed a bimodal distribution where the right-side (left-side) mode is well-consistent with the distribution of physical (non-physical) links, irrespective of network topologies. Namely, the value that separates two modes of the distribution can distinguish very successfully physical links from non-physical links. Therefore, reconstruction of an underlying network from the functional activities is highly possible for all network topologies, at this weak coupling regime (Figure 7c).
Right at the onset of global synchronization (K ≈ 0.02 for σ(κ) ≈ 3.0, 0.03 for σ(κ) ≈ 1.4 and 0.04 for σ(κ) = 0.0), mean PPC for disconnected nodes deviates from zero for all network topologies (lower panel, Figure 7d), indicating that pairwise correlation between disconnected nodes starts at the emergence of global synchronization. Moreover, this happens earlier in more heterogeneous networks. In Figure 6, as σ(κ) increases, the overlap between two PPC distributions for connected and disconnected nodes at K = 0.05 becomes more and more significant. The distribution of all links shows a unimodal distribution, causing a difficulty to determine the threshold while implementing the reconstruction of physical connectivity. This is due to the fact that in heterogeneous networks, such as SF network, disconnected hub nodes are strongly pairwise synchronized whereas some connected small-degree nodes remain weakly correlated. This dependence of the PPC of two nodes on their degree is shown in Figure 1a. The existence of pairwise synchronization between the nodes that are not physically connected hinders the inference of a background network from the functional activities around K = 0.05 for all network topologies (Figure 7c).
At the strong coupling, K = 0.140, the behaviour of PPC changes when the network topology varies between homogeneous and heterogeneous (Figure 6). In homogeneous networks (σ(κ) = 0.0 and σ(κ) ≈ 0.4), both of the distributions of PPC for connected and disconnected nodes are very sharp compared to those in heterogeneous networks. This indicates that, at a strong coupling strength, PPC of two nodes in homogeneous networks purely depends on whether they are physically connected or not. Only for homogeneous networks does the distribution of all links show a bimodal distribution, facilitating the inference of the underlying network for homogeneous topology at the strong coupling (Figure 7c). In heterogeneous networks, the more neighbours a node has, the higher effective coupling strength via which it interacts with the network. This results in the dependence of the PPC of two nodes on their degrees. At a strong coupling strength, hub nodes are strongly pairwise synchronized whether they are connected or disconnected. In heterogeneous networks, the distributions of PPC both for connected and disconnected nodes are broader, and there is a Significant overlap between the two populations. This results in a lower capacity for reconstructing the background network for heterogeneous networks in the regime of strong coupling.
There have been many works aiming to discover the relationship between the dynamical process in the system and the underlying physical connectivity [21–27]. Many studies performed so far have investigated the impact of background network topologies on the collective behaviour of the system as represented by the global order parameter [23–27, 35]. In this paper, we have comprehensively investigated a dynamical network system at both regimes of weak and strong coupling with a family background networks interpolating between regular and scale-free topologies.
Our study on the collective behaviour of the system proves that heterogeneous network systems are easier to be synchronized whereas homogeneous network systems go to the complete synchronization faster. By calculating pairwise phase correlation for every pair of either connected or disconnected nodes, we showed how the synchronization emerges and propagates inside the system. The homogeneous distribution of node degree results in the fact that the pairwise correlation of two nodes in a homogeneous network is purely dependent on whether they are connected or not. In heterogeneous network systems, the path to synchronization of heterogeneous network system is nontrivial. Synchronization starts among hubs and propagates toward low-degree nodes.
The strong relationship between the physical connectivity and the functional dynamics suggests an ability to solve one of inverse problems: reconstruct physical connectivity from observed functional dynamics. In spite of dealing a similar reconstruction of physical network, we use the synchronization pattern induced by a complex network system without external inputs whereas Timme M. in  use the method of analyzing the response dynamics of the network system upon perturbation. More comprehensively, our reconstruction of physical connectivity was tested for the systems with network topologies ranging from regular to complex at both weak and strong regime of coupling. We pointed out that the regimes of weak coupling, right before the onset of global synchronization, facilitate a successful reconstruction of physical connectivity, irrespective of network topologies. Clearly, our method to reconstruct physical connectivity is applicable to a wide range of networks from regular to scale-free topologies.
Our study of functional dynamics in complex networks has important implications for medicine. The sudden emergence of synchronization among the elements of a system might cause serious damages. For example, excessive synchronization of neuronal activity in basal ganglia cortical loops is the hallmark of activity in Parkinson's disease . An epileptic seizure is assumed to be associated with abnormal excessive or synchronous neuronal activity in the brain . If we know, which locations inside the system become synchronized first, we can anticipate the systems' damages prior to their occurrence.
where i = 1, ..., N, θ i and ω i are the phase initiated randomly within [-π, π] and the intrinsic frequency distributed uniformly within [-0.1, 0.1] of node i, K is the coupling strength. To make the model more realistic, we add Gaussian white noise, ξ i (t) with intensity D satisfying ⟨ξ i (t)⟩ = 0 and ⟨ξ i (t)ξ j (t')⟩ = 2Dδ ij δ(t - t'). In this study, we use D = 0.01. A ij is an element of the adjacency matrix, which is 1 if nodes i and j are connected and 0 otherwise. The coupling term, therefore, is summed over the neighbours of node i. We solved the eqs. (1) by using the improved Euler method with a step size of 0.01.
where the brackets ⟨...⟩ t signify the time averaging. The global order parameter measures the extent of global synchronization of the population of N oscillators.
We calculate the global order parameter at different regimes of coupling strength to study how the collective behaviour of all the oscillators changes from full-desynchronized and full-synchronized states.
Plotting RN β/νversus K for various network size, one can find the value of β/ν that gives a well-defined crossing point at K c . In this work, N = 128, 256, 512 and 1024. After obtaining K c and β/ν, one can use eq. (4) to determine the value of (1- β)/ν, which combined with the known value of β/ν, yields the values of β and ν.
C ij , which shows how dependent the motions of two oscillators at nodes i and j are, is equal to 0 or 1 corresponding the full incoherence or coherence between nodes i and j, respectively. This measure of phase coherence was used in .
We monitored the changes in PPC of all pairs of nodes as K increases to find which pairs of nodes synchronized first. In Figure 1a, we plot log (1/(1 - C ij )) versus the coupling strength K and the product of k i and k j to see how the PPC depends on the neighborhood connectivity. The range of color from blue to brown corresponds to the value of C ij changing from 0 to nearly 1.
We define a weighted network of synchronization as an all-to-all functionally connected network where each link is weighted by PPC of two oscillators at the ends of the link. In a WNS, physical and non-physical links are between the nodes that are connected and disconnected in the background network, respectively. WNS, therefore, keeps the information of the organization of synchronization of the system.
In (7), r i is the local order parameter which measures the coherence of k i neighbors of node i. The oscillator at node i interacts with the network via an effective coupling strength = Kk i r i (Figure 1d).
The reconstruction of physical connectivity from functional dynamics is depicted in Figure 4. The eqs. (1) are solved with different configurations of initial phases and intrinsic frequencies of oscillators. A so-called averaged WNS is built where each link is weighted by the PPC averaged over several initial configurations. In the averaged WNS, links whose weights are higher than a given threshold are predicted as physical links. On the contrary, predicted non-physical links have weights that are lower than the threshold.
The threshold can be determined by presenting the distribution of all links in terms of their weights. At the regimes of weak coupling where the reconstruction performs the best, the weight distribution of all links shows a bimodal distribution. The value by which two modes of the distribution are well separated can be used as the threshold. For example, at K = 0.010, a bimodal distribution of all links is observed for all the network topologies and the threshold can be around 0.130 (Figure 6).
The performance of reconstruction is analyzed by a Receiver Operating Characteristic (ROC) curve, which depicts relative trade-offs between true positive (benefits) and false positive (costs). We compute true positives (TP), false positives (FP), true negatives (TN) and false negatives (FN). TP denotes true predicted links. FP denotes non-physical links that are predicted as physical links by our method. TN denotes true predicted non-physical links. FN denotes the physical links that are predicted as non-physical links by our method. In the plot of the ROC curve, the x-axis represents false positive rate (FPR), that is FP/(TN+FP), and the y-axis represents true positive rate (TPR), TP/(TP+FN) as the threshold gradually varies.
In this study, the best reconstruction would yield a point in the upper left corner or coordinate (0,1) of the ROC space, representing 100% true predicted physical links without false positives. The (0,1) point corresponds to an optimal reconstruction. It the physical network was reconstructed by randomly choosing links from the WNS, this would give a point along a diagonal line from the left bottom to the top right corners.
The reconstruction performance is quantified by the area under the ROC curve (AUC), which is usually used for model comparison in machine learning . This measure can be interpreted as the probability that when we randomly pick one positive and one negative example, the classifier will assign a higher score to the positive example than to the negative. An AUC of 0.5 reflects a reconstruction of physical network by randomly choosing links from the WNS, while AUC = 1 implies a perfect reconstruction.
The authors are grateful to Kim WS, Kim S, Kook H, Khang BN, Jeong H, Oymak H, and Kim PJ for invaluable comments and discussions. This work has been supported by the Korea Science and Engineering Foundation grant (R15-2004-033-007001-0) funded by the Ministry of Science of Technology, Korea. Hung X Ta acknowledges the support from Marie Curie grant (MRTN-CT-2005-019475). Seung Kee Han acknowledges the support from the BK21 program at Chungbuk National University
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.