Dynamical states, possibilities and propagation of stress signal

The stress driven dynamics of Notch-Wnt-p53 cross-talk is subjected to a few possible dynamical states governed by simple fractal rules, and allowed to decide its own fate by choosing one of these states which are contributed from long range correlation with varied fluctuations due to active molecular interaction. The topological properties of the networks corresponding to these dynamical states have hierarchical features with assortive structure. The stress signal driven by nutlin and modulated by mediator GSK3 acts as anti-apoptotic signal in this system, whereas, the stress signal driven by Axin and modulated by GSK3 behaves as anti-apoptotic for a certain range of Axin and GSK3 interaction, and beyond which the signal acts as favor-apoptotic signal. However, this stress system prefers to stay in an active dynamical state whose counterpart complex network is closest to hierarchical topology with exhibited roles of few interacting hubs. During the propagation of stress signal, the system allows the propagator pathway to inherit all possible properties of the state to the receiver pathway/pathways with slight modifications, indicating efficient information processing and democratic sharing of responsibilities in the system via cross-talk. The increase in the number of cross-talk pathways in the system favors to establish self-organization.


Results
Regulatory biochemical network model (Fig. 1), which allows cross-talk between p53 and Wntch (Wnt and Notch) signaling pathways, are proposed based on some experimental reports (see Methods) to investigate complex signaling processes among them and behaviors of possible dynamical states those could correspond to certain cellular states.
Dynamical states driven by stress signal. The dynamical state of a system, we define here, is the dynamics of the system governed by a time series, T(t) which follow Mandelbrot's Multifractal behavior 38 throughout the time series, where, the fractal function F i (κ i ) follows simple power law F i (κ i ) ~ κ D with fractal dimension D. Variation in concentration of nutlin trigger DNA damage [20][21][22] (larger the nutlin concentration in the system, stronger the DNA damage), as a consequence stress signal is imparted to the system via p53 which is very sensitive to the stress signal. p53 dynamic become stabilized for low concentration of nutlin (small value of k 35  , which may correspond to nearly normal dynamical state of the system ( Fig. 2A upper panel). Increase in k 35 (0.006 < k 35 < 0.31) drives the p53 dynamics in the mixed dynamical state, damped for a certain range of time and then stabilized, which can be expressed as, , where, F s (κ) is for certain  range of time (time period and wavelength are constant, but amplitude decays). This may indicate that the system is activated first due to induced stress, and upon expended the stress the system comes back to normal. As k 35 increases the range of T s (κ) also increases. For significant strong stress (0.32 < k 35 < 0.52), the p53 dynamics become sustain oscillation dynamical state (range of T s becomes infinite), which can be represented by, , κ is an integer, and T(t) becomes time period in this case. This dynamical state may correspond to activated state of the system where active interaction of the molecular species is involved. Further increase in k 35 (0.58 < k 35 < 5.1) again brings back to the mixed dynamical state. Excess stress due to nutilin (k 35 > 5.3) drives the dynamics to stabilized dynamical state, which may correspond to apoptosis. Thus we could able to find five such dynamical states as listed in Table 1. The system switch to these various dynamical states in response to the stress signal.
Another way to impart stress signal to the system could be by changing Axin concentration (changing k 22 ) in the system which is reflected in the dynamics of p53 (Fig. 3A). Similar dynamical states, as obtained in the case of nutlin, can also be seen driven by Axin.
Complexity in the dynamical states. The complexity of the dynamical states (Figs 2A,H,L and 3A,H) driven by stress inducer nutlin and Axin can be measured by calculating permutation entropy (S i , i → p53, Notch, Axin) of these states (see Methods). The calculated S i magnitudes are in the order: S i (first stabilized state) < S i (first mixed state) < S i (sustain oscillation state) < S i (second mixed state) < S i (second stabilized state) (Figs 2B,I,M and 3B,I). The first stabilized dynamical state may correspond to normal, which is an ordered state, giving the minimum value of S i . Whereas, the second stabilized dynamical state corresponds to apoptotic state, which is a disorder state, providing maximum value of S i .
Multifractal due to long range correlation. The stress driven p53 time series are multifractal due to short and long correlations in the time series and they might follow some probability distributions. The fluctuation variation function F s (Fig. 2C upper panel) follows power law with fluctuation parameter s (log-log plots show approximately straight line) for all five dynamical states ( Fig. 2C) with varied slopes which correspond to Hurst exponents H q . The H q magnitudes with respect q for different dynamical states are as follows: H q (sustain oscillation) > H q (second stabilized state) > H q (first and second mixed state)H q (first stabilized state). Since large H q value corresponds to large fluctuations introduced in the system due to active molecular interaction driven by stress signal or internal molecular mechanisms constituted, active (sustain oscillation) dynamical state associates maximum fluctuations on an average, whereas stabilized dynamical states (both first and second) have comparatively much lower fluctuations. Since natural systems are far from equilibrium 3 , fluctuations driven systems might incorporate the fluctuations in positive way 2 to establish their own adaptable state for better existence 5 . Further, H q variation is prominent mainly in negative q regime as compared to the positive q regime, indicating multifractal nature in the time series is contributed mainly by long range correlation. This behavior is also reflected in the singularity function F α as the function of α.
The multifractal properties of the p53 time series for various dynamical states induced by Axin show similar behavior as obtained in the case of nutlin driven p53 dynamical states (Fig. 3C). The only slight change in the nature of singularity function F α which show similar behavior for all dynamical states.
Complex stress management. The amplitude death scenario in p53 dynamics (A p53 ) as a function of k 35 ( Fig. 2D) could be signatures of normal (first stabilized dynamical state) and apoptotic (second stabilized dynamical state) cellular states. Monotonic is increasing and decreasing of A p53 as a function of k 35 correspond to first and second mixed states, where removing the stress may come back to either normal or apoptotic state. Slow change in A p53 corresponds to sustain oscillation dynamical state, where the system is strongly activated due to active molecular interaction induced by strong stress in the system. GSK3, the main mediator of the Notch-Wnt-p53 cross-talk, could probably regulate the transition of the system at these dynamical states by managing the stress signal. Permitting the increase in GSK3 concentration through biochemical mechanism in the system allows the system to stay in active state (sustain oscillation dynamical state) for larger range of k 35 (Fig. 2D and E), trying to save the system from apoptotic phase. Therefore, GSK3 favors the system to stay in normal state for small stress signal, whereas for significantly larger stress signal, it  forces the system or the system itself prefers to stay in the active state to prevent from apoptotic state (Fig. 2E). Therefore, nutlin induced stress signal along with GSK3 acts as anti-apoptotic signal in the system.
On the other hand, this stress management by GSK3 is quite different for this system when stress in it is introduced by Axin via k 22 (Fig. 3D and E). Even though the dynamical states, which we found in the case of nutlin driven stress p53 dynamics, also obtained in this case, the pattern is quite opposite. In this case, increase in the concentration of mediator GSK3 (k 39 ) the range of k 22 within which sustain oscillation dynamical state can be observed becomes decreased (Fig. 3D). This increase in k 39 force the available, accessible area of active state to diminish (Fig. 3E), favoring the system to go to apoptotic state. Hence, the stress signal modulated by Axin along with GSK3 acts as favored-apoptotic signal.
Complicated self-organization. The time series of the dynamical states of stress p53 induced by nutlin ( Fig. 2A) can be transformed into their respective complex networks (Fig. 2G) using visibility graph algorithm (see Methods), where, properties of the system, reflected in p53 dynamical states, can be studied using topological properties of the constructed networks. The probability of degree distribution P(k), clustering co-efficient C(k) and neighborhood connectivity C n (k) follow power law behaviour as a function of the degree k, The power law behaviour in these topological properties are verified and confirmed by using the fitting technique to test power law distribution proposed by Clauset et al. 39 , where all statistical p-values for all data, calculated against 2500 random samplings, are found to be larger than 0.1 which is the critical limit, and goodness of fits are found to be less than and equal to 0.3. The power law behaviours of these three topological parameters are the signatures of Hierarchical features 40 in these networks exhibiting multifractal nature in their structures. The lowest value of γ = 2.3 in sustain oscillation dynamical state exhibits importance of, not only system level organization of modules, but also few hubs in the network regulation. Networks corresponding to other dynamical states also follow these hierarchical properties, but the network constructed from second stabilized state (apoptotic state) lacks most of the network properties (γ = 4.3, α = 0.03, β = 1.5) due random distribution in the dynamics.
Centrality measurements, namely, betweenness C B , closeness C C , and eigen-vector C E centralities (see Method), which provide the nature of information processing and to identify most influencing nodes in the network, again follow power low behaviors as a function of k (Fig. 2G), We again verify the qualification of the respective power laws to the corresponding data with the statistical technique applied above 39 . The increase in these centrality values as a function of degree (k) indicate that hubs are most influencing nodes, and take significant roles in information processing in the network. Even though the absence of these hubs do not cause network breakdown, they and their interaction among them have better responsibilities in regulating the network in each dynamical state.
The constructed networks from the dynamical states of stress p53 triggered by Axin ( Fig. 3G) also show similar hierarchical properties as obtained in nutlin induced stress p53 networks. The calculated values of exponents of power law topological parameters are, The multifractal properties of dynamical states of Notch due to stress signal propagated from p53 either induced by nutlin ( Fig. 2J and N), or Axin ( Fig. 3J) show very similar nature as compared to that of p53 dynamical states (Fig. 2C). Further, topological parameters of the networks corresponding to dynamical states of stress receivers Notch and Wnt (due to stress inducer nutlin) ( Fig. 2K and O); and Notch (due to stress inducer Axin) (Fig. 3K) exhibit similar properties as compared to the properties of the networks corresponding to dynamical states of stress propagator p53. This closely similar properties of the stress signal propagators and receivers reveals the inheritance of fundamental properties of the stress signal from propagators to receivers. Hence, the responsibilities, beneficials and injuries due to stress signal are approximately equally distributed among the important candidate proteins and/or genes in the interacting pathways. However, depending on the nature of stress signal and magnitude, and stress receivers, the inherited properties are modified according to their needs and comfortability.
Since the response of the stress propagated among the cross-talking pathways are received excellently via various mediators, the responsibilities of protecting the system as whole from any change due to external perturbations or failures in internal mechanisms might be distributed approximately equally among these pathways.

Signature of assortivity.
Hierarchically organized networks generally qualify most of the features of self-organization, namely, fractal behaviours of topological parameters, system level organization of modules and absence of central control mechanism (removal of hubs do not cause network breakdown) 41,42 . However, neighborhood connectivity in these networks constructed from the dynamical states of stress p53 driven by either nutlin or Axin follows, C n (k) ~ k β , which is power law of positive exponent (Figs 2G and 3G). The fractal behaviour is the signature of importance of few hubs and their interaction exhibiting assortivity nature of the network topology 43,44 . Even though the few hubs do not have the capability of full control of the network, their significantly strong cross-talks can have the possibilities of regulating the network up to some extent. In these time series of dynamical states, there could be few hubs corresponding to few time states which are responsible for signal processing and management.
The propagation of this stress signal from the propagator p53 pathway to the receivers Notch and Wnt pathways associates the assortivity characteristics. This indicates that in the time series corresponding to dynamical states of the system, interaction of few hubs (formation of the rich-club) is probably essential for efficient signal propagation and stress management.
Active state is preferred state. The active dynamical state (sustain oscillation state) generally has maximum H q values (Figs 2C,J,N and 3C,J middle pannels) due to large fluctuations than the other dynamical states.
Since the system is perturbed with stress signal, the fluctuations coming from active molecular interaction due to stress signal becomes optimal, which could probably be utilized by the system in a constructive way 2 , to establish simple fractal law where the system can stay comfortably with this possibility. The network, corresponding to this dynamical state, exhibit topological properties closest to hierarchical network properties 40,43 , However, the noise associated with this dynamical state is far smaller than the noise at second mixed and stabilized (apoptotic) states, as evident from permutation entropy S p53 calculations.
Cross-talks favor to establish self-organization. Stress p53 dynamics driven by nutlin exhibits optimized active dynamical state as number of cross-talks of pathways are increased (Fig. 4A-E). S i for p53 in the p53 pathway (direct stress is introduced via nutlin) is maximum, then the S i is pulled down drastically as it just cross-talks to Wnt pathway. When Notch-Wnt-p53 cross-talk is established S i is slightly increased, may be trying to optimize the fluctuations imparted by active molecular interaction due to stress in the system. Calculated H q value is minimized as number of cross-talks of pathways is increased (Fig. 4C), but amplitude A p53 is optimized (Fig. 4D).
The constructed networks from these p53 dynamics corresponding to different combinations of cross-talks of the pathways (Fig. 4E) show that the topological properties of the network corresponding to Notch-Wnt-p53 interaction exhibits better organized and closer to hierarchical features, The assortivity property, indicating the importance of hubs in the network, is still one of the important properties of the network. The important responsibilities of these hubs in information processing in the network are reflected and become clear in the fractal nature of the centrality measurements.

Discussion and Conclusion
Natural systems always try to be self-organized, subjected to a number of possibilities to change and corresponding states which are far from equilibrium, try to adapt to the preferred changes beneficial to their own organization for better and comfortable survival. During this struggle, the system tries to optimize any change (such as fluctuations due to stress in the system) in a constructive manner, utilize and distribute it, and look for better stability. If the system cannot able to face the change, it will break down.
Stress induced in Notch-Wnt-p53 cross-talk by stress inducing molecules (nutlin and Axin in our case), which are exhibited in stress p53 dynamics, drives the system at various dynamical states defined by different fractal laws, and the system switch to these dynamical states depending on the amount of stress induced. This stress system prefers to stay in an active dynamical state which has simple fractal rule subjected to the optimal fluctuations available due to active molecular interaction driven by stress. However, the system still associates a group of few hubs (assortive topology), but not in dependent manner (absence of these hubs do not cause system's breakdown), for better signal processing and system regulation. Then this stress signal is propagated throughout the pathways, and found to inherit all the properties of the propagator pathway to the receiver pathways may be with slight modifications in them. This excellent co-ordination in cross-talk helps the system to save it from one directional apoptosis (once the system falls in this phase, it can never come back to normal situation) by regulating available active molecular interaction. This regulating mechanism could be different depending on the type of stress induced in the system (nutlin and Axin in our case). The increase in the number of cross-talk pathways allows the system to optimize the imparted fluctuations, and therefore favours to establish self-organization in the whole system. This cross-talk regulates various biological functions out of which is regulating apoptosis. In the case of nutlin induced stress system, nutlin along with cross-talk mediator GSK3 acts as anti-apoptotic signal. Whereas, Axin induced stress along with GSK3 acts in two fold, first as an anti-apoptotic within a certain range of their interaction range, second favour-apoptotic signal beyond that regime. Since this Notch-Wnt-p53 cross-talk is involved in many important cellular processes in normal and cancerous cells, one needs to investigate these cellular activities with a number of open questions to understand how they work at a fundamental level.

Methods
Coupled Wntch (i.e. Wnt/Notch) Signaling in development. The notion Wntch signaling was introduced by Martinez Arias in 2008, for the elements of integrated system controlling cell fate i.e. Wnt and Notch 45 . The most important part of multicellular organisms is development which requires the proliferation of cell that decides its fate 46 . This decision for fate determination converge on the enhancement of specific genes that create a combination of transcription factors that determine the state and behaviour of the cell 47 . According to this analysis obtained from the study of Dorsophila genetics it can be proposed that Wnt is involved in regulating the probability of cell fate adoption 48 . The current understanding of Wnt signalling gives two striking features as reviewed by Martinez Arias and Hayward 2006, first is that many elements of this pathway take part in many other signaling processes, and second, is the ability of Wnt signaling to co-operate with other transcription factors and modulate their effect in other signaling pathways 48 . Thus, it acts as a noise filter, i.e making fluctuation in gene expression less strong 49 . However, there is one pathway that is consistent with Wnt, which is generally known as Notch signaling pathway 48 . Since Notch is involved in lateral inhibition, it together with Wnt plays crucial role in patterning the Notch inhibitory process 47 . Since, the development of wings in Dorsophila was first established by Couso and Martinez Arias in 1994, it also gave us a strong proof of interaction of notch and Wnt with strong sensitivity 50 . These evidences indicate that the elements of Notch and Wnt form a closed pathway for information processing, where Notch is available for NICD and Wnt regulate the Notch traffic using the β-catenin and activity of GSK3 in specific cells. Thus, Wntch signaling gives a clear picture of transition of states in cell during development 51 .
Notch-Wnt-p53 cross-talk biochemical network. Notch-Wnt cross-talk is a well-studied model [16][17][18]25,52,53 which regulates variety of biological functions in normal and cancer cells, from lower to higher level organisms. On the other hand, p53 is very sensitive in responding stress in the system, may be via DNA damage or due to other molecular mechanisms, and drives the cell at various cellular states and lets the cells to decide their fate of survival 19,54 . The cross-talk between Notch-Wnt and p53 pathways are particularly important to understand, (i) how does p53 interfere Notch-Wnt interaction, and vice versa, (ii) nature of stress signal propagation in interacting pathways, and (iii) regulation of apoptosis. We present here Notch-Wnt-p53 regulatory biochemical network (Fig. 1) by incorporating simplified versions of the above models and other experimental findings.
Model Description. Wnt signaling pathway allows Wnt to interact with Dsh (Desheilved protein) to inhibit the formation of ubicutinaceous complex between Gsk3β and Axin2, which de-phosphorylate β-catenin (k 6 ) to enter the nucleus (with a rate k 4 ) and initiates the expression of Axin2-mRNA (k 16 ). Then Axin2-mRNA is transported into the cytoplasm as part of the negative feedback loop between Wnt and Axin2. Free Gsk3β available in the system will allow to communicate Wnt with p53 by competing with the binding of p53 and Mdm2, with the formation of complex Gsk3β with p53 (k 39 ). The process then activates the transcription of Mdm2-mRNA (k 41 ), and it is translated (k 27 ) and translocated to the cytoplasm where it binds to unbound p53 to degrade (k 33 ). On the other hand, β-catenin present in the nucleus binds to Lef1 (k 49 ) to increase the expression of mRNA of Delta gene (k 54 ), followed by translation of Delta (k 51 ) and translocated to the cytoplasm (k 58 ). Now, this Delta protein interacts with Notch-NICD binary complex present on the cell membrane 55 to initiate the release of NICD to form ternary complex (k 67 ) and facilitating its transport to the nucleus (k 73 ). The release of NICD in the cytoplasm is due to the ligand Delta that performs two activities cis-inhibition (intracellular i.e degradation of Notch) and trans-activation (intercellular) of Notch 53,55 . In our case cis-inhibition and trans-activation goes simultaneously within the cell, where cis-inhibition degrades Notch with rate constant k 65 and trans-activation of Notch releases NICD in Cytoplasm (see Supplementary Fig. 1) followed by the nuclear transportation 55 . Further, NICD in the nucleus triggers the expression of Lfringe (mRNA) (k 76 ) 56 followed by its translation (k 79 ). After translation, Lfringe goes to bind to its site on Notch followed by the inhibition of formation Notch delta complex to inhibit its own (Lfringe) transcription (k 76 ) exhibiting negative feedback loop between Notch and Lfringe 57 . Lfringe does not prevent the binding of DSL ligands to Notch, while it potentiates both via Notch 58 . On the other hand, Lfringe decreases the binding of Notch and Delta by modifying (O-linked fucose glycosylation) the extracellular region of Notch2 59 . However, Wnt directly inhibits the formation of Lfringe by direct transfer of information via interaction of Dsh and NICD present in cytoplasm (k 86 ). The periodic expression of Lfringe is essential for somite formation 56 . Notch pathway interacts with Wnt via a mediator molecule Dsh, that bind to the cytoplasmic form of NICD (NICDc) (k 85 ). This Dsh inhibits the formation of the destruction complex (formed by the association of GSK3β, Axin and β-catenin). GSK3β from this complex binds to p53 in the presence of Mdm2 forming binary complex p53-GSK3β (k 39 ). This binary complex activates the transcription of Mdm2 with a rate constant of 0.042 min −1 (k 41 ), increasing the concentration of Mdm2 inside the cell. This hike in the Mdm2 decreases the p53 concentration (k 33 ). The bound p53 is restored by Nutilin (k 36 ), a small molecule that compete to the binding site of p53 on Mdm2. On the other hand activated β-catenin is translocated into the nucleus (k 4 ) and bind with the Lef1 (k 49 ) to activate the transcription of Delta in the nucleus (Delta mRNA). This Delta mRNA is translated in cytoplasm to the protein that go to Notch to promote the cleavage of NICD (k 51 ). The interaction of Notch-Wnt-p53 pathways allows to regulate apoptosis as well as stress propagation.
Mathematical framework of the network. The Notch-Wnt-p53 regulatory network model (Fig. 1) is defined by N = 28 molecular species (Table 2) corresponding to the reaction network description provided in Table 3. The state of the system at any instant of time 't' is given by the state vector,  x (t) = (x 1 , x 2 , … , x N ) T , where, N = 28 and 'T' is the transpose of the vector. By considering feedback mechanism of in p53, Wnt and Notch oscillators and coupling reaction channels of the two oscillators, we could able to reach the following coupled ordinary differential equations (ODE), where, the functions in the equation (8) {F i (x 1 , x 2 , … , x N )}, i = 1, 2, … , 28 are given by,                        We used standard Runge-Kutta method (order 4) of numerical integration to simulate equation (8) to find the solution of the variables listed in Table 2 for the parameter values given in Table 3. We then analyzed the constructed mathematical model to get possible approximate analytical solutions of the variables (slow variables) using quasi-steady state approximation.
Multifractal DFA approach. Fractal properties in non-stationary time series, and associated important correlations can studied using Multifractal detrended fluctuation analysis (MF-DFA) 60 . Important fractal parameters which characterize the time series, namely, Hurst exponent (H), generalized dimension (D) etc can be calculated numerically using a method adopted by Kantelhardt et al. 11 as summarized below. First, the time series signal {x j } of length N is taken as random walk, and can be represented by the profile, where, 〈 x〉 is the mean value of the signal, and i = 1, 2, … , N. Second, the profile Y(i) is now divided into = ( ) N int s N s equal non-overlapping equal segments of size s. To take into account all data points, 2N s segments are considered by counting starting from both ends of the data. Third, the following variance is determined, where, ν = N s + 1, … , 2N s , and y ν (i) is the fitting polynomial in segment ν. Fourth, the q th order fluctuation function is estimated by averaging over all segments, where, H q is the generalized Hurst exponent, which represents the measure of self-similarity and correlation properties of the signal. Then H q is related to the classical scaling exponent τ(q), q and from the definition of Holder exponent, α = τ d dq , the singularity function f(α) 60 is given by, Then, generalized fractal dimension of the signal is measured by, Now, D 0 , for q = 0, is the fractal or Hausdorff dimension, D 1 is information dimension and D 2 represents correlation dimension 60  Visibility graph of time series. This technique maps a time series to a network 61 , where each observation in time series is translated to a node and an edge between any two nodes is introduced when the following visibility condition is satisfied i.e. two nodes corresponding to observations x(t a ) and These networks are undirected due to symmetry in visibility condition. Since the properties of the time series are inherited to the corresponding network, the studies of this network provide useful information which can't be observed in traditional time series data.
Topological properties of networks. The following topological properties of the networks are studied to study the important behavior of the networks.
Degree distribution. The degree k of a node indicates the number of links the node connects with other nodes in the network. Consider a network defined by a graph G = (N, E), where N and E are number of nodes and edges respectively. The probability of degree distribution (P(k)) of the network is the probability that any chosen node will have a degree k, which is given by, k where, n k is the number of nodes having degree k. P(k) in random and small-world networks follow Poisson distribution, whereas, it obeys power law P(k) ~ k −γ in scale-free and hierarchical networks depending on the value of γ which indicates the importance of hubs or modules in the network 62 .
Clustering co-efficient. Clustering co-efficient of a network characterize how strongly a node's neighborhood nodes are interconnected. For an undirected network, clustering co-efficient (C(k i )) of ith node is the ratio of number of its nearest neighborhood edges to the total possible number of edges of degree k i , and can be calculated by, where, e i is the number of connected pairs of nearest-neighbor of ith node, and k i is its degree. C(k) in scale free networks is independent of k, whereas in hierarchical network it follows a power law, C(k) ~ k −α , with α ~ 1.
Neighborhood connectivity. Neighborhood connectivity of a node of a network is a measure of the average connectivities of the nearest neighbors of the node in the network 63 , and is given by, where, P(q|k) is conditional probability that a link belonging to a node with connectivity k points to a node with connectivity q. The power law nature of C N (k), C N (k) ~ k −β is a signature of hierarchical topology in the network 43 . However, the positive power dependence of C N (k) could be an indicator of assortivity in the network topology 44 .
Betweenness centrality. Betweenness centrality of a node characterize the ability to (i) extract benefits from information flows in the network 64 , and (ii) extent to which the node has control over the other nodes in the network through communication 65,66 . If d ij (v) indicates the number of geodesic paths from node i to node j passing through node v, and d ij represents number of geodesic paths from node i to j, then betweenness centrality (C B (v)) of a node v can be measured by, If M denotes the number of node pairs excluding v, then normalized betweenness centrality is given by, Closeness centrality. Closeness centrality (C C ) estimates how fast information is spread from a node to other nodes reachable from it in the network 67 . C C of a node i is the harmonic mean of geodesic distance between the node and all other nodes connected to it in the network, where, d ij is geodesic path length between nodes i and j, and n is the total number of nodes in the network connected to node i.
Eigenvector centrality. Eigenvector centrality of a node i (C E (i)) in a network is proportional to the sum of i′ s neighbor centralities 68 , and it is measured by, where, nn(i) indicates nearest neighbors of node i in the network. λ is the eigenvalue of the eigenvector v i given by, Av i = λv i , where, A is the adjacency matrix of the network. The principal eigenvector of A, which corresponds to maximum eigenvalue λ max , is taken to have positive eigenvector centrality score 69 . Since node's eigenvector centrality function smoothly varies over the network and depends on its neighbors, node with high eigenvector centrality is embedded in the locality of nodes of high eigenvector centralities, and chance of having isolated nodes in and around the locality is very low 68 . Hence, eigenvector centrality can be used as an indicator of node's spreading power in the network.
Algorithm for calculating permutation entropy. Important information contained in a time series can be measured by calculating permutation entropy of the time series 70,71 . Permutation entropy H of a time series of a dynamical variable x(t) of a system can be calculated as follows. The time series x(t) can be mapped onto a symbolic sequence of length N: x(t) = {x 1 , x 2 , … , x N }. This sequence is then partitioned into M number of short sequences of equal length U i.e. x(t) = {w 1 , w 2 , … , w M }, where ith window is given by w i = {x i+1 , x i+2 , … , x i+U }. This window is allowed to slide along x(t) with maximum overlap. Permutation entropy of a window w i can be calculated by defining a short sequence of embedded dimension r, S i = {x i+1 , x i+2 , … , x i+r } in r-dimensional space, finding all possible inequalities of dimension r and mapping the inequalities along the ascendingly arranged elements of w i to find the probabilities of occurrence of each inequality in w i . Since q out of r! permutations are distinct, one can define a normalized permutation entropy as