Matching Rules for Collective Behaviors on Complex Networks: Optimal Configurations for Vibration Frequencies of Networked Harmonic Oscillators

The structure-dynamics-function has become one of central problems in modern sciences, and it is a great challenge to unveil the organization rules for different dynamical processes on networks. In this work, we study the vibration spectra of the classical mass spring model with different masses on complex networks, and pay our attention to how the mass spatial configuration influences the second-smallest vibrational frequency () and the largest one (). For random networks, we find that becomes maximal and becomes minimal if the node degrees are point-to-point-positively correlated with the masses. In these cases, we call it point-to-point matching. Moreover, becomes minimal under the condition that the heaviest mass is placed on the lowest-degree vertex, and is maximal as long as the lightest mass is placed on the highest-degree vertex, and in both cases all other masses can be arbitrarily settled. Correspondingly, we call it single-point matching. These findings indicate that the matchings between the node dynamics (parameter) and the node position rule the global systems dynamics, and sometimes only one node is enough to control the collective behaviors of the whole system. Therefore, the matching rules might be the common organization rules for collective behaviors on networks.


Introduction
Various dynamical processes on complex networks [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16] have attracted a great deal of interest in many disciplines, ranging from physical, chemical, and biological sciences, to even social sciences and engineering technology, with the central target to decipher and even utilize the emergent features coming out of the interplay between the structure and function via dynamics. Since a function is a dynamic property, the structure-dynamics-function has been one of common problems in modern sciences [1][2][3][4][5]. These dynamical behaviors can be very rich including synchronization in coupled nonlinear systems, consensus or grouping in multi-agent systems, self-organized criticality, epidemic spreading, traffic jam, stochastic resonance, and many others. They pose astonishing challenges for us to unveil the organization rules for all these different dynamical processes on networks.
So far the dynamics on extensive systems composed of many coupled identical units has been extensively studied, for example, the complete synchronization in coupled chaotic systems based on the stability analysis of the synchronous manifold and in the framework of the master stability function [17][18][19][20][21][22][23]. However, the assumption of identical units is not very realistic for many realistic systems, in which the units composing the ensemble always present a disparity in the values of some characteristic parameters. Under these conditions, the system (parameter) diversity becomes very important [24][25][26][27][28][29]. It is, however, still unknown how these nodes with different characteristic parameters settled on networks determine the systems dynamics, especially when some key parameters are properly placed on some special positions of networks. This optimal configuration problem has never been carefully studied, to the best knowledge of the authors, although it certainly could help to better understand the collective behaviors of coupled nonidentical systems, which are determined by not only the interplay between the temporal information (the node dynamics with different parameters) and the spatial information (these nodes' positions on networks), but also their matchings.
For this purpose, we study the optimal configuration problem for the vibration spectrum (or frequencies) of coupled nonidentical harmonic oscillators on complex networks. Each vertex is occupied by a mass point and they are connected by an edge (spring) each other with a uniform spring constant. This massspring (linear) model and the normal mode analysis originally developed in classical mechanics [30,31], have been found very useful in many other disciplines, such as lattice vibrations and associated phonon excitations in solid-state physics [32], and protein dynamics in structural biology [33,34]. As we are interested in the diversity effect, we only consider the case for each node having a different mass; this is the key distinction with all previous studies on coupled harmonic oscillators [35][36][37][38].

Model
The motion equation for the systems can be written as where m j represents the j-th oscillator's mass, k is the coupling strength (or spring coefficient), and N represents the number of mass points. a ij~1 if nodes i and j are connected and a ij~0 otherwise. For simplicity, we take k:1 throughout the paper. The masses m j for all j are different. The system's motion can be viewed as N normal modes with different frequencies: x j~P N i~1 y j e iv i t , where v i and y j denote one independent frequency and the amplitude, respectively. v 1 ƒv 2 ƒ . . . ƒv N .
Clearly v 1~0 , representing the translational motion. For the other modes, usually the second smallest (slowest) vibration frequency v 2 characterizes the most global vibrational motion of the systems, whereas the largest vibration frequency v N reflects the most tightly packed and constrained local motion; both of them are the most representative. Therefore, below we only consider the values of v 2 and v N for different spatial configurations of the masses.
The above equations [Eqs.
(1)] can be written in a compact form based on the normal mode analysis: where the mass matrix M~diag(fm j g), the Laplacian matrix L~(d ij d i {a ij ), where d ij~1 for i = j and d ij~0 otherwise, d i is the degree of node i, and Y is the N-dimensional nonzero column vector with the components y i . Solving Eqs.
(2) is equivalent to finding v satisfying det(v 2 M{L)~0, where ''det'' denotes the determinant. Since M is nonsingular, it is convenient to change the mass-weighted Laplacian matrix to a symmetric one and get a faster numerical result. Therefore, we have where l denotes the eigenvalue of the corresponding matrix. And we obtain v based on solving the eigenvalues of the symmetrical matrix M {1=2 LM {1=2 for a given network. In numerics, we have studied various networks, such as the random Erdos-Renyi (ER) [39], scale-free (SF) [40], and small-world (SW) [41] networks. For each network, an ensemble of masses is produced from a random distribution (1ƒm j ƒ10) and these values are fixed and scattered on networks. For each configuration, its vibration frequencies can be easily calculated from the above eigenvalue analysis [Eq. (3)].
We are interested in how the correlation between degrees and masses influences the normal mode frequencies v 2 and v N . Figure 1(a) illustrates the numerical results of v 2 in an ascending order for completely random configurations (i.e., all masses can be arbitrarily shuffled on complex networks); a random ER network (N = 200 and the average degree SdT~30) is chosen here. The corresponding distribution re-scaled by the peak value is given in Fig. 1(b). Clearly v 2 shows a very wide distribution from v 2 &1:2 to v 2 &1:65. However, we may also choose configurations purposefully, such as point-to-point-positive correlation between the node masses and the node degrees (namely, we may arrange the masses in an ascend order fm 1 ƒm 2 ƒ . . . ƒm N g and the degrees in an ascend order also fd 1 ƒd 2 ƒ . . . ƒd N g, and put these arranged masses on the network in a point-to-point-positive correlation order fd 1 =m 1 ,d 2 =m 2 , . . . ,d N =m N g). We really find its v 2 for this configuration is maximal, as shown the open square in Figs. 1(a) and 1(b). Correspondingly, for the configuration of pointto-point-negative correlation fd 1 =m N ,d 2 =m N{1 , . . . ,d N =m 1 g illustrated by a solid square, it seems that v 2 is very small.

Simulation results
Below let us see if the above strong correlation condition for the point-to-point matching could be loosened. Figure 1(c) shows the results for the single-point-positive correlation fd 1 =m 1 g (upper curve) and fd N =m N g (lower curve), with all other N{1 masses randomly placed. Here we use fd 1 =m 1 g to denote the configuration of putting the lightest mass on the lowest-degree node, whereas fd N =m N g for that of the heaviest on the highest-degree node. It is discernible that the curve for fd 1 =m 1 g is slightly above the counterpart for fd N =m N g, which is quite similar to that for the random configurations in 1(a). Moreover, Figs. 1(d) and 1(e) show the results for the single-point-negative correlation fd N =m 1 g and fd 1 =m N g, respectively, with all other N{1 masses randomly placed again. Clearly the distribution in 1(d) is still similar to that for the random configurations in 1(a). However, very remarkably the values of v 2 for fd 1 =m N g (with one and only one heaviest mass on the lowest degree) have been greatly squeezed into a much smaller regime: 1:204vv 2 v1:210 and they are indeed very small [noticing the different ordinates and their scales used in Figs. 1(a) and 1(e)]. This range is only one-hundredth of the original range for the random distribution or other single-point matching distributions. We also find that v 2 for the point-to-pointnegative correlation in Fig. 1(a) is not minimal, which is re-plotted in Fig. 1(e) with the same solid square. This indicates that the single-point matching fd 1 =m N g has already caught the essential feature of the collective behavior in such a model of harmonic oscillators on random networks, regardless of all other nodes. As a result, for v 2 , the single-point matching fd 1 =m N g should give rise to (v 2 ) min , whereas the point-to-point matching fd 1 =m 1 , d 2 =m 2 , . . . ,d N =m N g should lead to (v 2 ) max , although these observations still need proof in theoretical analysis.  Fig. 2(a), give rise to the minimal and maximal values of v N , respectively. They are just opposite to those for v 2 . Again v N for the configuration of point-to-point-positive correlation is very far from the distribution in 2(b). In Figs. 2(c) and 2(e), the single-point matching distributions for fd 1 =m 1 g, fd N =m N g, and fd 1 =m N g give no significant difference with that for the random distribution. However, for the single-pointnegative correlation fd N =m 1 g in 2(d), we do find v N 's have been squeezed into a much smaller regime; this time it corresponds to the part of the largest v N . All these findings show the key qualitatively unchanged features with those for v 2 .
To show the generality of matching rules for any random networks, we give some more examples.   Figs. 4(a) and 4(b), the two points keep as extreme, locating far away the distribution for random configurations. However, with decreasing random connections for smaller P, as shown in 4(c) and 4(d), these relations become broken. Based on these comparisons, we understand that the rules of point-to-point matching and the single-point matching may only work for sufficiently random networks, such as ER or SF networks, and these rules could gradually be broken with the weakening of the network random connections, such as the SW networks.

Analysis
So far some unusual effects of matching rules have been well revealed. Below let us demonstrate them in a rigid way with the aid of mathematical analysis. As we have known, the normal mode frequencies are solely determined by the mass-weighted Laplacian matrix [Eq. (3)]. Therefore, some approximation results for the weighted Laplacian matrix developed in the matrix algebra and applied in the chaos synchronization study should be very valuable for our problems here [20,21,42]; for the second smallest and largest eigenvalues l 2 and l N , we have where c and c 0 can be approximated by 1 for most complex networks, and S i , the intensity function of a node, is defined as S i~P N j~1 W ij a ij , where W ij is the element of the weighted matrix. Since W ij~1 =m i here, we have S i~d i m i and further When the network is sufficiently random (d min &1 and d min & ffiffiffiffiffiffiffiffi ffi SdT p ), Figure 2. The results of v N for a random ER network. Similar to Fig. 1 for the results of v N with the same ER network considered instead. Again the effect of one is enough appears, but this time v N becomes maximal for the single-point-negative correlation configuration fd N =m 1 g in (d). For more details, see the text. doi:10.1371/journal.pone.0082161.g002 Therefore, for a sufficiently large SdT, As a result, l 2 and l N (and hence v 2 and v N correspondingly) should be completely determined by S min and S max , respectively. These linear relations have been perfectly proved by our simulations, for example, the plot of S min versus v 2 2 in Fig. 1(f) and that of S max versus v 2 N in 2(f). Based on these analyses [Eq. (5)], we immediately know that the configuration of single-point-negative correlation fd 1 =m N g (with the heaviest mass on the lowest-degree node) should always give rise to (S min ) min and (v 2 ) min irrespective of all other masses. Here one node match is dominant.
The result for point-to-point-positive correlation for (v 2 ) max can also be easily understood. Without losing generality, suppose that the masses and the degrees have been perfectly placed, i.e., Here uvv, d u ƒd v and m u ƒm v . If any one pair of nodes (u-th and v-th) are not placed in order, namely, node with mass m u is placed on the node of degree d v while m v on d u instead, we have the new configuration Based on the following inequalities du mv ƒ dv mu , du mv ƒ du mu , du mv ƒ dv mv , ð10Þ giving rise to we can derive that any one permutation from the original, perfect configuration of the point-to-point-positive correlation [Eq. (8)] would not let S min increase. Therefore, the original configuration would be globally optimal, corresponding to (v 2 ) max . On the other hand, we may use the same idea to the analysis of S max for (v N ) min . Finally after mastering these features of coupled harmonic oscillators, we may even control the systems collective behavior by a slight manipulation [43][44][45]. E.g., a single node with the heaviest mass m = 10 (1ƒm j ƒ10) has been added and connected to any nodes of an ER network, as shown in Fig. 5

Discussion
In conclusion, we have studied the parameter diversity effect of coupled harmonic oscillators with different masses on complex networks, and found that the values of v 2 and v N highly depend  on the configurations of these masses on the network spatial structure. Especially, two key matching rules including the pointto-point matching and single-point matching determine their extreme values. These findings might be helpful for explaining some interesting phenomena, such as biological swarming and flocking in nature, where the match of the head's ability with its position in hierarchy is always crucial for the behaviors of the whole group. In addition, the optimal configuration not only exhibits the importance of matching between the node dynamics and the node position and their impact on the collective behaviors of coupled systems, but also provides a potential control method for the manipulation of the systems dynamics. A possible candidate for this might be the point mutation in the protein dynamics study. We expect that the feature of ''one is enough'' may arouse general interest in studying the structure-dynamics-function relation in complex systems, and the matching rules may also become the common organization rules for various dynamical processes on complex networks. Figure 5. Controllability of network. Control of v 2 by adding one node (m = 10) and connecting it to any one node (lower part) and any two nodes (higher part) in an ER network, the same as in Figs. 1 and 2, in contrast to the uncontrolled value of v 2 (dashed line). doi:10.1371/journal.pone.0082161.g005