{'fraction_non_alphanumeric': 0.031900291000933396, 'fraction_numerical': 0.008126063800582001, 'mean_word_length': 4.526092233009709, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Master equation for the degree distribution of a Duplication and Divergence network

Network growth as described by the Duplication-Divergence model proposes a simple general idea for the evolution dynamics of natural networks. In particular it is an alternative to the well known Barab\'asi-Albert model when applied to protein-protein interaction networks. In this work we derive a master equation for the node degree distribution of networks growing via Duplication and Divergence and we obtain an expression for the total number of links and for the degree distribution as a function of the number of nodes. Using algebra tools we investigate the degree distribution asymptotic behavior. Analytic results show that the network nodes average degree converges if the total mutation rate is greater than 0.5 and diverges otherwise. Treating original and duplicated node mutation rates as independent parameters has no effect on this result. However, difference in these parameters results in a slower rate of convergence and in different degree distributions. The more different these parameters are, the denser the tail of the distribution. We compare the solutions obtained with simulated networks. These results are in good agreement with the expected values from the derived expressions. The method developed is a robust tool to investigate other models for network growing dynamics.


I. INTRODUCTION
Physicists are required to build models that extract the essence of observable phenomena seen in nature in order to understand and describe them. As physics endeavors in studying complex phenomena in distinct fields such as biology or social sciences, a commonly applied paradigm is the use of graph theory. In this approach the system under study is described as a network consistent of a set of nodes and a set of links among them.
In the present work we are interested in Network dynamics, i.e. the study of the dynamics behind a growing graph. The motivation lies in the dynamics of biological networks, given that many such systems find a natural description within those models. Understanding the evolutionary selection rules resulting in networks with similar topological characteristics as the real observed ones, may give insights about the underlying biological processes (natural selection) behind these structures. For example, the protein association networks can help understand the evolution of species' genomes [10][11][12][13].
In the present work, the network dynamics is described as a Markovian process. Within this approach, the network state at a given time depends only on its configuration on the previous moment. Given a set of rules that describe how the network changes in each time step, we construct the corresponding master equation representing the evolution of the system's configurations. A similar methodology was used in the work by Ferreira et al. [14]. With this approach, they presented analytical results and simulations of networks growing according to the Barábasi-Albert rule. Here we focus on a complementary approach to model protein-protein network dynamics [15] and explore the evolution of an adapted version of the Duplication Divergence model [16]. The importance in describing the average behavior of stochastic processes in this manner is to know the network behavior for different values of the parameters without the need of long, time consuming, numerical simulations to obtain statistically relevant information.
This article is organized as follows. In the next section we focus in explaining the Duplication and Divergence model. Following the model explanation, we derive an expression for the total number of links as a function of the number of nodes, which gives us a straight forward way to obtain the mean degree of the graph. Then we study the graph growth as a Markovian chain, in which the next degree distribution of the network is a function of the current degree distribution, pondered by the probabilities of all possible occurrences in each time step. Finally, we study the asymptotic limit of the degree distribution.

II. THE MODEL
Given an initial small network (three nodes connected to each other forming a triangle) we study the Markov process where, in each time step a node of the network is randomly chosen to be copied i.e. a new node is created with exactly the same neighbors as the chosen one. In what follows we refer to the copied node as original and its copy as duplicated. After duplication, original and duplicated nodes may diverge, meaning that each link of the original node is lost with probability m o and each link in duplicated node is lost with probability m d . Also, a link between original and the duplicated nodes is always added.
This model is an adaptation of the Duplication and Divergence model, originally developed by Vázquez et. al. [16,17]. In the model proposed by Vázquez a new node is also added by copying an existing node and all its links. New node and ancestor are linked with a probability p. Also, either the link between new node and a third neighbor or the link between ancestor and this neighbor is lost with probability q [17]. Our model treats independently the loss of a link by original and duplicated nodes and sets p = 1.
In Vázquez's model the network growth is based on local rules, that is, rules that require only information on one node instead of rules that require information over all the network. This model of network dynamics applied to protein-protein interaction networks allows all proteins to evolve from a common ancestor through gene copies (represented by duplications) and mutations (divergence). Therefore, it would mimic the entire history of a genome evolution [17].

III. MEAN DEGREE AND TOTAL NUMBER OF LINKS
First, let's evaluate the behavior of the network average node degreek as a function of the number of nodes in the network, t. Given that a node with degree k is chosen to be duplicated, the number of links in the next step changes. The mean change in the number of links is given by [18]: Terms in the right-hand side of Eq. (1) represent, from left to right: existing links, links added due to the duplication of a k-degree node, mean number of lost links of the original node, mean number of lost links for the duplicated node and creation of the original-duplicated link.
Let's define the total mutation parameter M as the sum of the independent mutation parameters m o and m d : M = m o + m d . Now, given the relation between node's degree distribution and number of links in a network, known as the handshake lemma [19], and, considering the time as the number of nodes (since we add a node in each time step), the following map for the network average degree is written This equation can be rearranged to explicit thek t+1 term:k For long times (big values of t), this can be approximated as a continuous differential equation: For 0 ≤ M ≤ 2, the solution of this ODE is: where the constant c is related to network initial conditions. For instance, consider the case of a process without divergence (M = 0). Starting with a triangle (k(3) = 2), the duplication will make every node connected with every other, i.e., a complete graph. In this particular case Eq. (4) results in In the limit t → ∞, the possible asymptotic behaviors for Eq. (4) are: For M > 0.5:: The mean degree converges to For M = 0.5:: The mean degree diverges logarithmically.
We represent these possible situations in Fig. 1, where each point refers to a network evolved through the model dynamics with a different value for the total mutation parameter M . In the vertical axis one has the correspondent mean degree calculated by Eq. (4) (in the limit t → ∞) and in the horizontal axis it is presented the results for the mean degree evaluated from simulated networks growing according to the model rules. The correspondence of the points colors with the parameter M can be read in the color scale. The figure clearly shows that, for high values of the total mutation M (values equal or bigger than 0.8), a network with 20 thousand nodes has already reached its stationary behavior (the term t 1−2M can be neglected). The closer the total mutation gets to M = 0.5, the slower one observes the convergence to the stationary limit. For values of total mutation under 0.5 it is expected that the mean degree diverges and so are the points in this figure departing from the main diagonal.
Using again the handshake lemma in equation (4), we obtain the total number of links as a function of time In Fig. 2 we show results comparing the mean value for the number of links in 5 thousand networks as a function of time with the expected value found through from Eq. (5).

IV. DEGREE DISTRIBUTION FOR FINITE TIMES
Considering the model dynamics, it is possible to establish relations among the number of nodes with degree k at time t, N (k, t), and at time t + 1: the terms in this equation represent, from top to bottom: the number of existing nodes of degree k; the probability of a k-degree node to be chosen to duplicate, leaving this degree; the first sum is the total probability that a node of degree different from k is duplicated and after the divergence process becomes a k-degree node; the second sum is the same total probability for the duplicated node. Finally, the last three terms are the probabilities that the neighbors of the node chosen to be duplicated arrive, from a different degree (either k+1 or k−1) to the degree k and the probability that a k-degree neighbor goes to k + 1 or k − 1.
Since each link, during the divergence process, is lost independently with probability m o or m d , the probabilities inside the sums are binomial distributions representing the probabilities that either the original or duplicated node goes from degree i to degree k − 1 after the divergence, and finally received the original-duplicated link: It is important to note the domain of the above functions (7a) and (7b), i ≥ k−1, that is, the node after passing through duplication and divergence can increase its degree by one unit (which means keeping all its links and adding the copied-duplicated link), keep the same degree (losing one neighbor and adding the copied-duplicated link) or else it will have its degree decreased due to the loss of more than one neighbor.
One can verify the addition of a single node in each time step by summing the master equation for N (k, t), Eq. (6), over all possible degrees ∞ k=0 (N (k, t + 1) − N (k, t)) = 1. As shown in Fig. 3, the numerical solution of the map in equation (6) produces results in excellent agreement with the mean number of nodes with degree k at the time t evaluated from thousands of simulated networks growing through the stochastic process.
As an analytic example, using m o = 0 e m d = 0, no links are lost by neither node, and equation (6) becomes, Eq. (9) can be written as Note that, in the case where the initial network is totally connected, that is, k = t − 1, there is only the flux of all nodes having its degree increased by one in a totally connected network, N (k, t) = tδ k−1,t .

V. ASYMPTOTIC DEGREE DISTRIBUTION
Equation (6) can be conveniently written in its matrix form, defining the column vector N (t) whose components are N k (t) = N (k, t): where the matrix A shows the following elements, constant with respect to time: The coefficients of equation (12) explicit the couplings in equation (11). It is possible to solve equation (11) through a matrix decomposition of the matrix A n×n in its eigenvalues and eigenvectors: Writing A = X −1 DX, putting this decomposition in equation (11) and multiplying from the left by the matrix X, one decouples the equations and is able to solve it for each component of the vector N in the space where the matrix A is diagonal. Returning to the original space one has: Where the λ i are the eigenvalues of matrix A, X i the respective eigenvectors (the columns of the matrix X, properly normalized, are the vectors X i ) and the constants c i depend on the initial conditions. Normalizing the vector N (t) one works with the fraction of nodes with degree k: note that P k (t) is the fraction of nodes of degree k at time t. The solution (13) can be written in terms of the new vector as follows: In the stationary state, condition (8) implies that λ = 1 is an eigenvalue of matrix A (whose left-eigenvector is (1, 1, 1, 1..., 1)) [20], and, to conserve probability, this must be the greatest eigenvalue, and the normalized right-eigenvector which corresponds to the unitary eigenvalue (called main eigenvetor) will be the asymptotic solution of the degree distribution of a network whose growth is governed by the Duplication and Divergence model with mutation rates m o and m d . Fig. 4 compares the main eigenvector numerically obtained with the simulation of networks up to 80 thousand nodes.
Note that the matrix A n×n has to be truncated for taking the limit t → ∞ because its actual size is t × t in Eq. (11), which is the greatest degree possible for (t + 1)node networks. For M > 0.5, this truncation still results in a matrix with an eigenvalue whose value is 1, which is related, as mentioned, to the stationary distribution. When M < 0.5, any truncation results in a matrix A which has all eigenvalues less than 1. So when the mean degree diverges, there is obviously no stationary degree distribution.

VI. PARAMETERS mo AND m d ARE INDEPENDENT
As seen in Eq.(4), the mean degree (stationary or timedependent) is a function that depends only on the total mutation rate, M = m o + m d . However, Eq.(6) cannot be written without considering only the sum of the individual mutations. Therefore, different partitions of the same total mutation lead to different distributions with the same mean. The parameters are interchangeable in both equations which reflects the fact that, after the duplication, the original and copied nodes are indistinguishable.
Finding the stationary distribution numerically for M = 0.80 in two possible scenarios, shown in Fig.5, gives us insights about the network behavior in both cases.
When the mutation parameters are equal (m o = m d = 0.4), one finds a distribution denser near the distribution's average. When the mutation parameters are dissimilar (m o = 0.79; m d = 0.01), the resulting distribution has a denser tail, indicating more nodes with lower and higher degrees than the average value. This is a very reasonable result considering the divergence process. In  6, we show this difference. It is possible to see the time evolution of the degree distribution P (k, t) for degrees in these three cases: lower, near and greater than the stationary average degree (that isk = 20/3, for M = 0.65) obtained numerically by integrating Eq. (6) and simulations.
One can also observe the real part of the eigenvalues of matrix A in both cases, shown in Fig. 7. The imaginary part is symmetric, which only reflects the fact that the elements of the matrix are real, and therefore the complex eigenvalues come in pairs as complex conjugates. From the real part of the eigenvalues one can infer the convergence rate. Therefore, one can conclude that for different values of m o and m d the convergence will take longer, since the eigenvalues are greater. In the symmetrical case, there are smaller eigenvalues than in the asymmetrical case. The rate of convergence to the stationary solution can also be inferred by the second highest eigenvalue of matrix A. In Fig. 8 dots represent the real part of the highest eigenvalues of A. It is possible to see that the second highest eigenvalue of matrix A for total mutation rate higher than 0.55 will be given by 2 − 2M , and therefore P (t) ∝ t 1−2M , which is the same exponent of the convergence of the mean degree of Eq. ( 4).

VII. CONCLUSIONS
In the present work, we derived the master equation for the degree distribution of networks evolving through the duplication and divergence model considering the mutation rates of original and duplicated nodes as independent parameters. Numerical integration of the resulting maps agree well with the average values obtained from simulated networks.
The maps obtained for the network nodes average degree and number of links can be approximated as ordi-  Figure 6. [color on-line] Degree distribution P (k, t) as a function of time, for networks with the same total mutation, but different mutation rates mo and m d . The top graphs are for low degree nodes, the middle ones are near the stationary average and the bottom ones are high degree nodes. nary differential equations and solved analytically. The asymptotic solution for these ODEs, agrees well with simulated data for values of total mutation greater than 0.5, limit for which the network converges to a stationary degree distribution.
Though the average node degree of a network evolving through this process only depends on the total mutation  rate parameter, its degree distribution and rate of convergence will depend on the mutation parameter difference between original and duplicated nodes. The behavior of this distribution can be studied from the spectrum of the matrix A, which contains the transition probabilities between nodes of different degrees in each time step of the Markovian process. The more similar the two parameters are, the faster the distribution converges because, in this case, one has less eigenvalues close to 1, and the less dense the distribution tail will be.
To use a Markovian process in order to describe network evolution, allows one to obtain iterative maps for the time evolution of the graph's properties in a wide variety of models. In the asymptotic limit, the maps can be approximated as ODEs, and therefore solved such that one is able to study qualitatively the dynamics dependence on the model parameters without the need to run computationally intensive numerical simulations. Finally, different network growing models may be treated using the present method. In particular, a work considering a hybrid Barabasi-Duplication/Divergence model is under development.

VIII. ACKNOWLEDGMENTS
We thank the contributions of collaborators from MTC-IF-UFRGS. We also thank CNPq for the financial support.