Simulated evolution of protein-protein interaction networks with realistic topology

One mechanism by which PPI networks change is gene duplication (DU) Ohno_1970 ; Zhang_2003 ; Xiao_2008 . In DU, an existing gene is copied, creating a new, identical gene. In our model, duplications occur at a rate $d$, which is assumed to be constant for each organism. All genes are accessible to duplication, with equal likelihood. For simplicity, we assume that one gene codes for one protein. One of the copies continues to perform the same biological function and remains under the same selective pressures as before. The other copy is superfluous, since it is no longer essential for the functioning of the cell Wagner_2003 .

The superfluous copy of a protein/gene is under less selective pressure; it is free to lose its previous function and to develop some other function within the cell. Due to this reduced selective pressure, further mutations to the superfluous protein are more readily accepted, including those that would otherwise have been harmful to the organism Koch_1972 ; Taylor_2004 . Hence, a superfluous protein diverges rapidly after its DU event Lynch_2000 ; Maslov_2004 . This well-known process is referred to as the post-duplication divergence. Following Vazquez_2003 , we assume that the link of each such superfluous protein/gene to its former neighbors is deleted with probability $\phi$. The post-duplication divergence tends to be fast; for simplicity, we assume the divergence occurs within the same time step as the DU. The divergence is asymmetric Kellis_2004 ; Gu_2005 : one of the proteins diversifies rapidly, while the other protein retains its prior activity. We delete links from the original or the duplicate with equal probability because the proteins are identical. As discussed in the supporting information (SI), this is closely related to the idea of subfunctionalization, where divergence freely occurs until redundancy is eliminated. In our model, $\phi$ is an adjustable parameter.

In many cases, the post-duplication divergence results in a protein which has lost all its links. These ‘orphan’ proteins correspond to silenced or deleted genes in our model. As discussed below, our model predicts that the gene loss rate should be slightly higher than the duplication rate in yeast, and slightly lower in flies and humans.

We simulate a gene duplication event at time $t$ as follows:

1a) Duplicate a randomly-chosen gene with probability $d\Delta t$.

2a) Choose either the original (50%) or duplicate (50%), and delete each of its links with probability $\phi$.

3a) Move on to the next time interval, time $t+\Delta t$.

Our model also takes into account that DNA can be changed by random mutations. Most such mutations do not lead to changes in the PPI network structure. However, some protein mutations lead to new interactions with some other protein (which we call the target protein). The formation of a novel interaction is called a neofunctionalization (NE) event. NE refers to the creation of new interactions, not to the disappearance of old ones. Functional deletions tend to be deleterious to organisms Lynch_1998 . We do not account for loss-of-function mutations (link deletions) except during post-duplication divergence because damaged alleles will, in general, be eliminated by purifying selection. In our model, NE mutations occur at a rate $\mu$, which is assumed to be constant. All proteins are equally likely to be mutated.

How does the mutated protein choose a target protein to which it links? We define a probability $q$ that any protein in the network is selected for receiving the new link from the mutant protein. To account for the possibility of homodimerization, the mutated protein may also link to itself Wagner_2003 ; Ispolatov_2005 . Random choice dictates that $q=1/N$ (see SI).

Many PPI’s are driven by a simple geometric compatibility between the surfaces of the proteins Jones_1996 . The simplest example is the case of PPI’s between flat, hydrophobic surfaces Tovchigrechko_2001 , a type of interaction which is very common Wu_2007 . These PPI’s have a simple planar interface, and the binding sites on the individual proteins are geometrically quite similar to one another. One consequence of these similar-surface interactions is that if protein A can bind to proteins B and C, then there is a greater-than-random chance that B and C will interact with each other. We refer to this property as transitivity: if A binds B, and A binds C, then B binds C. The number of triangles in the PPI network should correlate roughly with transitivity. As discussed below, the number of triangles (as quantified by the global clustering coefficient) is about 45 times higher in real PPI networks than in an equally-dense random graph. This suggests that transitivity is quite common in PPI networks. Another source of transitivity is gene duplication. If A binds B, then A is copied to create a duplicate protein A’, then A’ will (initially) also bind B. If A interacts with A’, then a triangle exists. However, duplication is unlikely to be the primary source of transitivity; recent evidence shows that, due to the post-duplication divergence, duplicates tend to participate in fewer triangles than other proteins Vinogradov_2009 .

A concrete example of transitivity is provided by the evolution of the retinoic acid receptor (RAR), an example of neofunctionalization which has been characterized in detail Escriva_2006 . Three paralogs of RAR exist in vertebrates (RAR$\alpha$, $\beta$, and $\gamma$), as a result of an ancient duplication. The interaction profiles of these proteins are quite different. Previous work indicates that RAR$\beta$ retained the role of the ancestral RAR Escriva_2006 , while RAR$\alpha$ and $\gamma$ evolved new functionality. RAR$\alpha$ has several interactions not found in RAR$\beta$. RAR$\alpha$ has novel interactions with a histone deacetylase (HDAC3) as well as seven of HDAC3’s nearest-neighbors (HDAC4, MBD1, Q15959, NRIP1, Q59FP9, NR2E3, GATA2). None of these interactions are found in RAR$\beta$. The probability that all of these novel interactions were created independently is very low. RAR$\alpha$ has 65 known PPI’s and HDAC3 has 83, and the present-day size of the human PPI network is a little over 3000 proteins. Therefore, the chance of RAR$\alpha$ randomly evolving novel interactions with 7 of HDAC3’s neighbors is less than 1 in a billion. This strongly suggests that when a protein evolves an interaction to a target, it has a greater-than-random chance of also linking to other, neighboring proteins.

How do similar-surface interactions affect the evolution of PPI networks? First, consider how an interaction triangle would form. Suppose proteins A and B bind due to physically similar binding sites. Protein X mutates and evolves the capacity to bind A. There is a reasonable chance that X has a surface which is similar to both A and B. If so, protein X is likely to also bind to B, forming a triangle. Denote the probability that two proteins interact due to a simple binding site similarity by $a$. The probability that A binds B (and X binds A) in this manner is $a$. Assuming these probabilities are identical and independent, the probability that X binds B is $a^{2}$.

So far, we have discussed transitivity as it affects the PPI’s in which protein A is directly involved (A’s first-neighbors). We now introduce a third protein to the above example, resulting in a chain of interactions: protein A binds B, B binds C, but C does not bind A. Protein X mutates and gains an interaction with A (with probability $a^{2}$). What is the probability that X will also bind C? The probability that B binds C due to surface similarity is $a$. Thus, X will bind C (A’s second-neighbor) with probability $a^{3}$. In general, the probability that X will bind one of A’s $j^{\text{th}}$ neighbors is $a^{j+1}$. We refer to this process as assimilation, and the ‘assimilation parameter’ $a$ is a constant which varies between species. As discussed in SI, it is primarily mutliple-partner proteins which bind to their partners at different times and/or locations which are affected by this process; consequently, at most one link is created by assimilation at the first-neighbor level, second-neighbor level, etc. Assimilation is assumed to act on a much shorter time scale than duplication and neofunctionalization; in our model, it is instantaneous.

Our hypothesized assimilation mechanism makes several predictions that could be tested experimentally: (1) the probability of a protein assimilating into a new pathway should be $a^{2}$ (at the first-neighbor level), $a^{3}$ (at the second-neighbor level), and so on, where $a$ is a constant which varies between species; (2) weak, nonspecific binding and planar interfaces should be overrepresented in interaction triangles (and longer cycles) between non-duplicate proteins; (3) competitive inhibitors should be overrepresented in interaction triangles; and (4) domain shuffling should be associated with assimilation. (See SI for discussion of (3) and (4).)

We simulate a neofunctionalization event at time $t$ as follows:

1b) Mutate a randomly-chosen gene with probability $\mu\Delta t$.

2b) Link to a randomly-chosen target protein.

3b) Add a second link to one of the target’s first-neighbor proteins, chosen randomly, with probability $a^{2}$.

4b) Add a link to one of the target’s second-neighbor proteins, with probability $a^{3}$, etc.

5b) Move on to the next time interval, time $t+\Delta t$.

A flowchart of how PPI networks evolve in our model is shown in Figure 1. To simulate the network’s evolution, one of the two mechanisms above is used at each time step, using Gillespie_1977 . We call each possible time series a trajectory. We begin each trajectory starting from two proteins sharing a link (the simplest configuration that is still technically a network). Each simulated trajectory ends when the model network has grown to have the same total number of links, $K$, as found in the experimental data, $K_{\text{data}}$. Here, we perform sets of simulations for three different organisms: yeast (Saccharomyces cerevisiae), fruit flies (Drosophila melanogaster), and humans (Homo sapiens). Because evolution is stochastic, there are different possible trajectories, even for identical starting conditions and parameters. We simulated 50 trajectories for each organism. Our figures show the median values of each feature as a heavy line, and individual trajectories as light lines.

For a given data set, the number of links ($K_{\text{data}}$) is known. We estimate the duplication rate $d$ from literature values. There have been several empirical estimates of duplication rates, mostly falling within an order of magnitude of each other Lynch_2000 ; Gu_2002 ; Gao_2004 ; Lynch_2003 ; Osada_2008 ; Lynch_2003 ; Cotton_2005 ; Pan_2007 . We averaged together the literature values to estimate $d$ for each species (Table 1).

The quantity $\mu$ is not as well known. Its value relative to $d$ has been the topic of considerable debate Wagner_2003 ; Beltrao_2007 ; Lynch_2008 ; Gibson_2009 . Although, in principle, $\mu$ is a measurable quantity, it has proven difficult to obtain an accurate value, in part because the fixation rate of neofunctionalized alleles varies with population size Kimura_1957 ; Walsh_1995 . In the absence of a consensus order-of-magnitude estimate, in our model, we treat $\mu$ as a fitting parameter. Consistent with the findings of He_2005 and Beltrao_2007 , our best-fit values of $\mu$ are within an order of magnitude of each other for yeast, fruit fly, and human networks. Best-fit parameter values are given in Table 1.

One test of an evolutionary model is its predictions for present-day PPI network topologies. Current large-scale PPI data sets have a high level of noise, resulting in significant problems with false positives and negatives Deane_2002 ; Deeds_2006 . To mitigate this, we compare only to ‘high-confidence’ experimental PPI network data gathered in small-scale experiments (see Methods). We computed 10 topological features, quantifying various static and dynamic aspects of the networks’ global and local structures: degree, closeness, eigenvalues, betweenness, modularity, diameter, error tolerance, largest component size, clustering coefficients, and assortativity. 8 of these properties are described below (see SI for others).

The degree $k$ of a node is the number of links connected to it. For protein networks, a protein’s degree is the number of proteins with which it has direct interactions. Some proteins interact with few other proteins, while other proteins (called ‘hubs’) interact with many other proteins. Previous work indicates that hubs have structural and functional characteristics that distinguish them from non-hubs, such as increased proportion of disordered surface residues and repetitive domain structures Patil_2010 . The high degree of a protein hub could indicate that protein has unusual biological significance Jeong_2001 . The network’s overall link density is described by its mean degree, $\langle k\rangle$ (Table 2). The degree distribution $p(k)$ is the probability that a protein will have $k$ links. PPI networks have a few hub proteins and many relatively isolated proteins. The heavy tail of the degree distribution shows that PPI networks have significantly more hubs than random networks have. Simulated and experimental degree distributions are compared in Figure 2. (For quantitative comparisons, see SI.)

Component refers to a set of reachable proteins. If any protein is reachable from any other protein (by hopping from neighbor to neighbor), then the network only has one component. If there is no path leading from protein A to B, then A and B are in different components. The fraction of nodes in the largest component ($f_{1}$) is a measure of network fragmentation (Table 2 and Figure S2). Note that, although silent genes (proteins with no links) exist in real systems, these genes do not appear in data sets consisting only of PPI’s. Therefore, calculations of $f_{1}$ for all models exclude orphan proteins (proteins with $k=0$).

Gene loss, the silencing or deletion of genes, is known to play an important role in evolution. The loss of a functioning gene will damage an organism, making the gene loss unlikely to be passed on. The exception is if the gene is redundant. Consistent with this reasoning, evidence suggests that many gene loss events are losses of one copy of a duplicated gene Ku_2000 ; Kellis_2004 . Although empirical estimates of the gene loss rate varied considerably, a consistent finding across several studies is that the rates of gene duplication and loss are of the same order-of-magnitude Lynch_2000 ; Cotton_2005 ; Gao_2004 . This broad picture is in good agreement with our model. In our model, a gene is considered lost when it has degree zero. Our model predicts that the ratio of orphan to non-orphan proteins is $1.6\pm 0.4$ in yeast, $0.58\pm 0.06$ in flies, and $0.67\pm 0.09$ in humans. The gene loss rate has been previously estimated to be about half the duplication rate in both flies and humans Lynch_2000 ; Cotton_2005 , consistent with our model’s prediction.

The distance between nodes $i$ and $j$ is defined as the number of node-to-node steps that it takes along the shortest path to get from node $i$ to $j$. The closeness centrality of a node $i$, $\ell_{i}$, is the inverse of the average distance from node $i$ to all other nodes in the same component. The diameter, $D$, of a network is the longest distance in the network. Simulated closeness distributions are compared to experiments in Figure 3. Interestingly, proteins have about ‘six degrees of separation’, similar to social networks Travers_1969 ; Leskovec_2008 . The closeness distributions $p(\ell)$ have peaks around $1/\ell\approx 5-7$.

Another property of a network is its modularity Yook_2004 . Networks are modular if they have high densities of links (defining regions called modules), connected by lower densities of links (between modules). One way to quantify the extent of modular organization in a network is to compute the modularity index, $Q$ Newman_2004 ; Newman_2006 :

where $k_{i}$ and $k_{j}$ are the degrees of nodes $i$ and $j$, $u_{i}$ and $u_{j}$ denote the modules to which nodes $i$ and $j$ belong, $\delta(u_{i},u_{j})=1$ if $u_{i}=u_{j}$ and $\delta(u_{i},u_{j})=0$ otherwise, and $A_{ij}=1$ if nodes $i$ and $j$ share a link, and $A_{ij}=0$ otherwise. $Q$ quantifies the difference between the actual within-module link density to the expected link density in a randomly connected network. $Q$ ranges between $-1$ and $1$; positive values of $Q$ indicate that the number of links within modules is greater than random. The numerical value of $Q$ required for a network to be considered ‘modular’ depends on the number of nodes and links and method of computation. To calibrate baseline $Q$ values given our particular network data, we used the null model described in Maslov_2002 . Our non-modular baseline values are $Q=0.603$ for the human PPI net, $Q=0.590$ for yeast, and $Q=0.722$ for flies (see SI). As shown in Table 2, PPI networks are highly modular, and our simulated $Q$ values are in good agreement with those of experimental data.

The clustering coefficient, $C_{i}$, for a protein $i$, is a measure of mutual connectivity of the neighbors of protein $i$. $C_{i}$ is defined as the ratio of the actual number of links between neighbors of protein $i$ to the maximum possible number of links between them,

In a PPI network, clustering is thought to reflect the high likelihood that proteins of similar function are mutually connected vonMering_2002 . The average (or global) clustering coefficient, $\langle C\rangle$, quantifies the extent of clustering in the network as a whole. As shown in Table 2, PPI networks have large global clustering coefficient values; the yeast PPI network, for example, has a value of $\langle C\rangle$ which is 45 times higher than that of a random graph of equivalent link density. In flies and humans, our simulated networks have $\langle C\rangle$ values in excellent agreement with the data; in yeast, our predicted value is slightly low.

A network is said to be ‘hierarchically clustered’ if the clustering coefficient and degree obey a power-law relation, $C\sim k^{-\xi}$ Ravasz_2002 (Figure 4), indicating that nodes are organized into small-scale modules, and the small-scale modules are in turn organized into larger-scale modules following the same pattern Barabasi_2004 . By plotting each node’s clustering coefficient against its degree, we observed a trend consistent with hierarchical clustering, although data in the tail is very limited.

The betweenness of a node measures the extent to which it ‘bridges’ between different modules. Betweenness centrality, $b$, is defined as:

Betweenness has been proposed as a uniquely functionally-relevant metric for PPI networks because it relates local and global topology. It has been argued that knocking out a protein that has high betweenness may be more lethal to an organism than knocking out a protein of high degree Joy_2005 . Betweenness distributions are shown in Figure 5.

If a network’s well-connected nodes are mostly attached to poorly-connected nodes, the network is called disassortative. A simple way to quantify disassortativity is by determining the median degree of a protein’s neighbors ($n$) as a function of its degree ($k$). Previous work has found that yeast networks are disassortative Maslov_2002 . It has been argued that disassortativity is an essential feature of PPI network evolution, and recent modeling efforts have heavily emphasized this feature Zhao_2007 ; Wan_2010 . However, it was noted by Aloy_2002 that disassortativity may simply be an artifact of the yeast two-hybrid technique, and Hakes_2008 pointed out that this trend is quite different among different yeast datasets, and in some cases is completely reversed, resulting in assortative mixing, where high degree proteins prefer to link to other high-degree proteins. As shown in Figure 7 and Table S1, the empirical data shows no evidence of disassortativity in flies or humans, and even the trend in yeast is quite weak. This conclusion is based solely on analysis of the empirical data, and casts further doubt on the role of disassortative mixing in PPI network evolution.

Some higher-order features of the network are simply a result of its degree sequence, and other features might be important in their own right. As discussed in Maslov_2002 , it is possible to isolate the effects of the degree sequence by ‘rewiring’ (detaching then reattaching links) the network at random, subject to the restriction that the degree sequence must be preserved. If a property contains extra information about the network’s structure, then it should be different in the rewired network. On the other hand, if the network is rewired many times, and the property is always the same, then it is likely to just be a result of the degree sequence. We randomly rewired the empirical network $4K_{\text{data}}$ times.¹¹1Network rewiring was carried out using a script downloaded from http://www.cmth.bnl.gov/$\sim$maslov/matlab.htm. As expected, modularity is decreased by random rewiring. Upon rewiring, we find $Q=0.603\pm 0.002$ in humans, $Q=0.590\pm 0.003$ in yeast, and $Q=0.722\pm 0.007$ in flies (median $\pm$ standard deviation from 50 repeats of the rewiring algorithm). Rewiring also shrinks the diameters of PPI networks to $D=13\pm 0.9$ in humans, $D=12\pm 1.0$ in yeast, and $D=15\pm 1.1$ in flies. These results suggest that these features contain important structural information about the network, and are not merely consequences of the degree sequence.

One reason we are interested in calculating $Q$ and $D$ is simply to check that the values are comparable between the simulated and experimental networks. However, on a more qualitative level, we would also like to have some idea of what the threshold is for a network to be considered ‘modular’ or ‘small-diameter’. The rewired $Q$ and $D$ values are useful because these features are dependent on the size of the network (number of nodes $N$ and links $K$); given an identical network construction method, $Q$ and $D$ will generally be different in sparse versus dense networks. We use these $Q$ and $D$ values as baseline values with which the experimental and simulated networks can be compared; we considered $Q$ and $D$ values differing from the rewired values by more than a standard deviation to be significantly different.

Comparisons of simulated and experimental eigenvalue spectra and error tolerance curves are shown in SI (Figures S1 and S4). As discussed in SI, the various per-node network properties we have analyzed are largely uncorrelated (Figure S6).

We now consider the question of how PPI networks evolve in time. The present-day networks show a rich-get-richer structure: PPI networks tend to have both more well-connected nodes and more poorly connected nodes than random networks have. In our model, the rich-get-richer property has two bases: duplication and assimilation. The equal duplication chance per protein means the probability for a protein with $k$ links to acquire a new link via duplication of one of its interaction partners is proportional to $k$. Likewise, the probability of a protein to receive a link from the first-neighbor assimilation probability $a$ is proportional to its degree $k$. ‘Rich’ proteins get richer because the probability of acquiring new links rises with the number of existing links.

First, we discuss two dynamical quantities for which experimental evidence exists: the rate of gene loss, and the relation between a protein’s age and its centrality. Gene losses in our model correspond to ‘orphan’ proteins which have no interactions with other proteins. As shown in Figure S2, the fraction of orphan proteins grows quickly at first, then levels off. This is consistent with the findings of Cotton_2005 : in humans, while the overall duplication rate is higher than the loss rate, when only data from the past 200 Myr are considered, the loss rate is slightly higher than the duplication rate. In our model, after the initial rapid expansion, the rate of gene loss stabilizes relative to the duplication rate.

We define the ‘age’ of a protein in our simulation according to the order in which proteins were added to the network. Our model shows that a protein’s age correlates with certain network properties. Consistent with earlier work Woese_1987 ; Krapivsky_2001 ; Qin_2003 ; Zhu_2011 , we find that older proteins tend to be more highly connected. We plotted the ‘age index’ of a protein (the time step at which the protein was introduced) versus its centrality scores. As shown in Figure 10, the age index negatively correlates with degree, betweenness, and closeness centralities: older proteins tend to be more central than younger proteins. Figure 10 shows our model’s prediction that a protein’s age correlates with degree, betweenness, and closeness centrality. We confirmed this prediction by following the evolutionary trajectories of individual proteins (Figure S3). These results are consistent with the eigenvalue-based aging method described in Zhu_2011 (Figure S5). Phylogenetic protein age estimates indicate that older proteins tend to have a higher degree Woese_1987 ; Zhu_2011 , which our model correctly predicts. Interestingly, the eigenvalue-based scores are only modestly correlated with other centrality scores (0.36 degree, 0.47 betweenness, and 0.10 closeness correlations). Using the eigenvalue method in tandem with our centrality-based method could provide stronger age-discriminating power for PPI networks than either method alone.

The correlation between centrality and age suggests that static properties of present-day networks may be used to estimate relative protein ages. Suppose each normalized centrality score ($k^{\prime}\equiv k/\max(k)$, $\ell^{\prime}\equiv\ell/\max(\ell)$, $b^{\prime}\equiv b/\max(b)$) represents a coordinate in a 3-D ‘centrality space’. We can then define a composite centrality score ($S$) as $S^{2}\equiv(k^{\prime})^{2}+(\ell^{\prime})^{2}+(b^{\prime})^{2}$.

Do older proteins typically have different functions than newer proteins? We classified S. cerevisiae proteins using the GO-slim gene ontology system in the Saccharomyces Genome Database. As shown in Figure 9, GO-slim enrichment profiles were somewhat different between the oldest and youngest proteins (as measured by their $S$ values). Several categories which were more enriched for the oldest proteins were the cell cycle, stress response, cytoskeletal and cell membrane organization, whereas younger proteins were overrepresented in several metabolic processes. Overall, the differences were not dramatic, suggesting that cellular processes generally require both central and non-central proteins to function. Consistent with this, ancient proteins tend to be centrally located with modules, as their betweenness values gradually decline over time (Figure S3). The roughly linear relation between degree and betweenness also suggests that ancient proteins do not occupy structurally ‘special’ positions within the network, such as stitching together separate modules (Table S1 and Figure 6). This may indicate that modules tend to accumulate around the most ancient proteins, which act as a sort of nucleus. Thus, ancient proteins are involved in all kinds of pathways, because they have each nucleated their own pathway.

In contrast to the two dynamical quantities discussed so far, most structural properties of PPI networks have only been measured for the present-day network. Although our model accurately reproduces the present-day values of these quantities, there is no direct evidence that the simulated trajectories are correct; rather, these are predictions of our model. Figure 8 shows that both modularity $Q$ and diameter $D$ increase with time. These are not predictions that can be tested yet for biological systems, since there is no time-resolved data yet available for PPI evolution. Time-resolved data is only currently available for various social networks (links to websites, co-authorship networks, etc.). Interestingly, the diameters of social networks are found to shrink over time Leskovec_2005_2 . Our model predicts that PPI networks differ from these social networks in that their diameters grow over time. In addition to $Q$ and $D$, we tracked the evolutionary trajectories of several other quantities: the evolution of the global clustering coefficient, the rate of signal propagation, the size of the largest connected component (Figure S2), as well as betweenness and degree values for individual nodes (Figure S3). See SI for details.

The number of ‘extra’ links created by each assimilation event should independent of $N$. Proteins with multiple interaction partners may interact with their partners simultaneously (so-called ‘party hubs’) or at different times/locations (‘date hubs’) Han_2004 . Since they bind to several partners simultaneously, party hubs typically have multiple binding sites, each specific to one binding partner Kim_2006 . This specificity suggests that a protein which evolved the capability to bind to a party hub would be unlikely to undergo assimilation. By contrast, the binding sites of date hubs are often disordered regions which are able to form transient interactions with multiple partners Ekman_2006 ; Singh_2007 . If a protein evolves the capability to bind to a date hub, it is likely to share the physical characteristics of the hub’s neighbors, leading to assimilation. However, to avoid competition for the same binding site, the interaction partners of date hubs tend not to be coexpressed Kim_2006 . One consequence of this is that assimilating proteins will likely only bind one of the target protein’s neighbors – whichever neighbor happens to be present at that time and place. Although the capability to bind to the hub protein’s other neighbors may initially be present, these will presumably remain unused in the cell. Our expectation is that the assimilating protein will therefore be unlikely to retain this capability, as it evolves. Similarly, only a single extra link should be generated at the second-neighbor level, third-neighbor level, etc. Consistent with the evidence discussed above, the number of links created by assimilation is approximately independent of the total network size. Party hubs typically are centrally-located within modules, while date hubs often function to stitch together large-scale modules in the cell. It may be that duplication-only models are unrealistically fragmented (Table 2) because their modules are not properly attached with date hubs; instead, the modules are disconnected components.

The DUNE model has four parameters. One parameter, the DU rate $d$, is estimated from empirical data. The other three are adjustable parameters: the NE rate $\mu$, the divergence probability $\phi$, and the assimilation probability $a$. To gain a better understanding of how these parameters affect the final network structure, starting with each organism’s set of parameters, we systematically adjusted all 4 parameters. Results are shown in Figure 11.

As expected, the divergence parameter $\phi$ was positively correlated with both the modularity $Q$ and the diameter. When $\phi\approx 0$, the network quickly reaches a fully-connected state, where all proteins are linked to all other proteins. Consequently, the network is not organized into modules, and all distances in the network are equal to 1. The ‘thinning out’ of duplicate links is therefore essential to generate non-trivial network features. The opposite occurs when the NE rate $\mu$ is too low: $f_{1}\approx 0$, and the network evolves towards a completely disconnected state.

Why does gene duplication lead to modularity? Consider an initially uniform network. When a node is duplicated at random, this causes the original node, the copy, and their immediate neighbors to share more links internally than they do with the rest of the network. Subsequent duplications amplify this effect: if a node that has 10 links within a module but only 2 external links is duplicated, there are now 20 internal and 4 external links (prior to post-duplication divergence).

Interestingly, $Q$ has a weak negative correlation with the assimilation parameter $a$. When $a$ is large, the probability to link to distant neighbors of the target protein is relatively high, and mutated proteins have a non-negligible chance to generate links to proteins outside of their target’s pathway. This causes modules to blur at the edges; their member proteins will share a higher number of links to other modules than for a low $a$ network. Although the modularity is reduced for a high $a$ network, it does not disappear entirely. Similarly, there is a sharp decrease in $Q$ as the NE rate surpasses the DU rate, indicating the important role of the DU mechanism in modular organization.

Diameter is also negatively correlated with $a$. When a single NE event has a significant chance to generate links to the target protein’s neighbors, this tends to reduce the overall separation of proteins in the network.

Genome-wide PPI screens have a high level of noise Deane_2002 , and specific interactions correlate poorly between data sets Deeds_2006 . We found that several large-scale features differed substantially between types of high-throughput experiments (see SI). Due to concerns about the accuracy and precision of data obtained through high-throughput screens, we chose to work with ‘high-confidence’ data sets consisting only of pairwise interactions confirmed in small-scale experiments, which we downloaded from the public HitPredict database Patil_2011 . We found sufficient high-confidence data in yeast (S. cerevisiae), fruit flies (D. melanogaster), and humans (H. sapiens).

All simulations and network feature calculations were carried out in Matlab. Our scripts are freely available for download at http://interacto.me. We computed betweenness centralities, clustering coefficients, shortest paths, and component sizes using the MatlabBGL package. Modularity values were calculated with the algorithm of Blondel_2008 . All comparisons (except the degree distribution) are between the largest connected components of the simulated and experimental data.

Due to the human network’s somewhat larger size, most dynamical features were calculated once per 50 time steps for the human network, but were updated at every time step in the yeast and fly networks. For dynamical plots, the $y$ coordinates of the trend line are medians-of-medians. The amount of time elapsed per time step (the $x$ coordinate) varies between simulations. We binned the time coordinates to the nearest 10 million years for yeast and fly, and 25 million years for human. When multiple values from the same simulation fell within the same bin, we used the median value. We then calculated the median value between simulations. Scatter plot trend lines are calculated in a similar way. The trend line represents the median response variable ($C$, $b$, or $\ell$) value over all nodes within a single simulation with degree $k$. The $y$ coordinate of the trend line is therefore the median (across 50 simulations) of these median response variables. This median-of-medians includes all simulations that have nodes of a given degree.

We downloaded large-scale data sets from BioGRID Stark_2006 , and used the Wilcoxon rank-sum test to compare aggregate statistical features across various experimental types in yeast (S. cerevisiae) and humans (H. sapiens) Han_2004 . As expected, we found that data obtained by affinity capture was significantly different than pairwise experimental data (primarily yeast two-hybrid and in vitro complexation), as the affinity capture interactions represent entire complexes, which is somewhat different information than the pairwise interactions we are attempting to capture using our model. However, more surprisingly, the only feature to show significant agreement between pair-wise techniques was the eigenvalue distribution of the walk matrix ($P>0.05$). Further sub-dividing the individual techniques into smaller data sets containing only results obtained in single experiments, we discovered that, again, only the spectra agreed between different screens.

Note that, due to the small size of the fly network, there may be too many missing links to obtain an accurate description the network’s large-scale topology. Although, by appropriate parameter tuning, our model is able to accurately reproduce the fly network, it is possible that different parameters will be required to match the fly network once it becomes more fully characterized experimentally. The data sets considered here do not include interactions which are enabled through post-translational modifications. Although these data sets are far from complete, and may be susceptible to false-positive detections, these appear to be the most accurate data available at the present time.