Gene duplication and hierarchical modularity in intracellular interaction networks

BioSystems 74(2004)51–62Gene duplication and hierarchical modularity in intracellularinteraction networksJennifer Hallinan ∗ARC Centre for Bioinformatics,Institute for Molecular Biosciences,The University of Queensland,Brisbane 4072,Qld,AustraliaReceived 13June 2003;received in revised form 29January 2004;accepted 2February 2004AbstractNetworks of interactions evolve in many different domains.They tend to have topological characteristics in common,possibly due to common factors in the way the networks grow and develop.It has been recently suggested that one such common characteristic is the presence of a hierarchically modular organization.In this paper,we describe a new algorithm for the detection and quantification of hierarchical modularity,and demonstrate that the yeast protein–protein interaction network does have a hierarchically modular organization.We further show that such organization is evident in artificial networks produced by computational evolution using a gene duplication operator,but not in those developing via preferential attachment of new nodes to highly connected existing nodes.©2004Elsevier Ireland Ltd.All rights reserved.Keywords:Modularity;Hierarchical;Network;Intracellular;Yeast;Gene duplication1.IntroductionNetworks of interactions between agents arise in a wide variety of contexts,including social net-works (Newman,2001),the Internet (Albert et al.,1999;Huberman and Adamic,1999)and the world wide web (Kleinberg and Lawrence,2002;Flake et al.,2002),ecological networks (Williams and Martinez,2000),and intracellular interaction networks (Bhalla and Iyengar,1999;Uetz et al.,2000;Sole and Pastor-Santorros,2002;Jeong et al.,2000).Analysis reveals that these diverse networks frequently have topological and dynamic features in common,and it has been suggested that these commonalities arise∗Corresponding author.Tel.:+61-7-3365-7131;fax:+61-7-3365-4388.E-mail address:j.hallinan@.au (J.Hallinan).from similar processes operating during the evolution and development of the networks.Topological characteristics which are common to many naturally occurring networks include a scale-free pattern of connectivity.A scale-free net-work has no characteristic number of connections per node as does a randomly constructed network;the probability P (k )of finding a node with k connections follows a power law:P(k)∝k −γ,(1)where the scaling exponent,γ,varies with the degree distribution of the network.When the degree,k ,of the nodes of a scale-free network is plotted against the probability of occurrence of that degree,P (k ),on a log–log scale the data forms a straight line,the slope of which is γ.Interaction networks may also exhibit “small-world”properties (Watts,1999).A small-world network has a0303-2647/$–see front matter ©2004Elsevier Ireland Ltd.All rights reserved.doi:10.1016/j.biosystems.2004.02.00452J.Hallinan/BioSystems74(2004)51–62 small but significant number of short-cut connectionsbetween otherwise widely separated nodes.This or-ganization leads to characteristic topological features,including a small diameter,where diameter is definedas the longest of the shortest paths between every pairof nodes in the network.Small-world networks alsohave a large average cluster coefficient,C,comparedwith randomly connected networks with the samenumber of nodes and links.Cluster coefficient is ameasure of the extent to which the neighbors of anode are linked to each other:C=1nni=1C iN i(N i−1)/2,(2)where,n is the number of nodes in the network,C i is the number of connections between neighbors of node I,and N i is the number of neighbors of node i(Watts, 1999).Many networks further appear to be organized into a number of modules.A module is generally defined as a subnetwork of a graph,the nodes of which have more connections to other nodes within the module than to external nodes(see,for example,Ancel and Fontana, 2002;Calabretta et al.,1998;Csete and Doyle,2002; Rives and Galitski,2003).The identification of modules within a network is an NP-complete problem(Flake et al.,2002).In prac-tice,a number of algorithms have been used for the identification of modules in networks.One approach involves the analysis offlux modes(the smallest sub-networks enabling the metabolic system to operate in steady state)within the network.Many polynomial time algorithms exist forfinding the maximumflow that can be routed from a source node,to a sink node, while obeying all capacity constraints(Flake et al., 2002;Stelling et al.,2002).This approach requires a “seed node”with which to initialize the algorithm. Another approach to module detection relies upon the identification of nodes or links which lie be-tween modules.Snel et al.(2002)define such link-ers as“orthologous groups with mutually exclusive associations”,and split the network at the linkers to produce modules which appear to be biologically plausible.Similarly,Schuster et al.(2002)split the network at nodes which have more than a threshold number of links,on the contention that such highly connected hubs must be external to the modules.They used a threshold number of links of four,but point out that other values may be useful,depending upon the size of the subnets produced.Module identification can be approached as a form of cluster analysis.Hierarchical clustering algorithms are widely used,even by researchers who are not interested in the cluster tree itself.The cluster tree is simply thresholded at an arbitrary depth in order to determine thefinal clusters(Ravasz and Barabasi, 2003).Girvan and Newman(2002)produced a clus-ter tree by identifying links with high“betweenness”(Freeman,1977)and iteratively removing the links with the highest betweenness to produce a cluster tree.An interesting approach was taken by Holme et al.(2002),who combined the node-removal ap-proach with the betweenness measure to develop an algorithm in which nodes of high betweenness are iteratively removed to deconstruct the network.It has recently been suggested that in addition to a modular organization,biological networks tend to have a hierarchical structure,in which nodes are orga-nized into small modules which are,in turn,organized into larger modules,and so on(Rives and Galitski, 2003).These authors propose a method for the iden-tification of hierarchical modularity which does not require the identification of individual modules.They derive a scaling law for the connectivity of nodes in a hierarchically modular networkC(k)≈k−1where,C(k)is the cluster coefficient defined in Eq.(2).Networks whose C(k)distributionsfit this curve are held to be hierarchically modular.Ravasz et al.(2002)have identified hierarchical modularity in the metabolic networks of43different organisms. While the scaling law provides a simple means of identifying hierarchical modularity in a network,it offers no insights into the form of that modularity and,hence,does not contribute to a detailed analysis of network structure.All of the algorithms discussed above rely upon user judgement either to chose the threshold at which the network is fragmented or to validate the biologi-cal plausibility of the modules.Since the algorithms are topology-based,an objective,topology-based measure of the“goodness”of a module would be a valuable addition to the module detection algorithms. In this paper we describe a new algorithm for theJ.Hallinan/BioSystems74(2004)51–6253detection of modularity,in conjunction with an ob-jective,topology-based measure of the coherence of the modules detected.These tools are combined to produce a coherence profile which can be used to visualize the extent of hierarchical modularity of a network,compare the modular structure of networks, and identify the threshold at which a network has maximum modular coherence.The major evolutionary operators which have been implicated in the evolution of scale-free networks are the preferential attachment of new nodes to highly connected existing nodes(Albert and Barabasi,2000) and the noisy duplication of existing nodes(“gene duplication”;Pastor-Satorras et al.,2002).The rela-tive importance of these operators to the development of real networks is unclear,and probably differs from network to network.Although both of these opera-tors have been demonstrated to produce a scale-free pattern of connectivity,their effect upon the mod-ularity of the network topology has not previously been investigated.We use the coherence profile al-gorithm to examine networks evolved according to several different published algorithms and compare the modularity of the resulting networks with that of the best-characterized biological network,the protein–protein interaction network of the yeast Sac-charomyces cerevisiae.work generation2.1.The yeast protein–protein interaction network Probably the best-characterized subcellular interac-tion network is the protein–protein interaction network of the bakers’yeast,S.cerevisiae.High-throughput methods for the collection of yeast protein–protein interaction data have been developed over the last 5years(Fields and Song,1989),and large interac-tion databases exist on the Web.The data in these databases is known to be both noisy and incomplete (von Meering et al.,2002).Both false negatives(in-teractions which exist in vivo,but have not been picked up by the screens)and false positives(inter-actions which occur under the particular conditions of a yeast two-hybrid screen,but not otherwise)will occur,to an unknown extent.Further,the network which can be constructed from the yeast two-hybrid data is a static snapshot of interactions,with none of the dynamic,temporal qualities of the network in the living cell.These problems mean that considerable care must be taken to choose the most reliable data with which to work,and care must be taken not to over-interpret the results of work done using yeast two-hybrid interaction data.In an effort to use only the most reliable data available,the dataset used for these experiments is the“core”set of S.cerevisiae protein–protein in-teractions identified by Deane et al.(2002)from the Database of Interacting Proteins(DIP database; /).This data is a sub-set of the entire DIP database consisting of those in-teractions which the authors verified using two forms of computational assessment,and is,therefore,less likely to contain false positive relationships than is the DIP database as a whole,although false negatives (missing interactions)undoubtedly occur.The core dataset comprises3003interactions be-tween1788proteins(an average connectivity of1.7). It does not form a single connected component,how-ever;there are139components,of which the largest has1471proteins and2770interactions(average con-nectivity1.9).This largest connected component was used for all investigations(Fig.1).work modelsSince we are interested in the evolution of biolog-ical interaction networks,the yeast protein–protein interaction network was used as the“gold standard”network for these experiments.In order to compare the effects of different evolutionary operators,we used different operators to generate networks with the same general characteristics as the yeast network. The yeast network,and probably many other bio-logical interaction networks,have three major charac-teristics:1.A power law connectivity with a well-defined cut-off.The distribution of connectivity within the network follows a power law.A truly scale-free network obeys this distribution over a wide range of connectivities.Naturally occurring networks, however,tend to deviate from the power law at the extremes of the distribution,probably because of physical factors affecting nodes:people form54J.Hallinan /BioSystems 74(2004)51–62Fig.1.The largest connected component of the curated yeast dataset.In this diagram the circles represent proteins and the lines represent interactions between proteins.relatively fewer new relationships as they age;proteins have physical limitations to the number of binding sites they can support (Amaral et al.,2000),and so on.The yeast network displays such a “cutoff”at the tail of the distribution.2.Sparse average connectivity .Although the range of connectivities is wide,most naturally occurring networks have an average connectivity of around 2.0–2.5.The average connectivity of the core yeast network is 1.9.3.Small-world characteristics .“Small-world”net-works are characterized by a small diameter rel-ative to the number of nodes in the network,and a large cluster coefficient in comparison with that of a randomly connected network of the same size and average connectivity.We generated networks with characteristics as close as possible to the size and average connectiv-ity of the yeast protein–protein interaction network,using two published algorithms which have been demonstrated to produce scale-free networks:gene duplication (Pastor-Satorras et al.,2002);and prefer-ential attachment (Ravasz et al.,2002).In addition,we generated randomly connected networks with ap-proximately the same size and average connectivity as the yeast network.Five networks were generatedusing each algorithm.The size,average connectivity,average diameter,and average cluster coefficient of these networks are described in Table 1.2.2.1.The preferential attachment modelScale-free networks were generated using the algo-rithm described by Albert and Barabasi (2000).In this algorithm,a network grows by the addition of new nodes to an existing node k i with probability,Π,pro-portional to the connectivity k i of node i :Π(k i )=k i +1j (k j +1)(3)Albert and Barabasi’s model produces scale-freenetworks only for a subset of possible values of the parameters,p and q (see Albert and Barabasi,2000for a full analysis of the behavior of the algorithm).The network analysis program Pajek (Batagelj and Mrvar,1998)incorporates an implementation of Albert and Barabasi’s algorithm,with default parameter values which will produce a scale-free network.Starting from these defaults (m 0(initial number of nodes)=3,m (nodes added at each time step)=2,p (probability of adding a link)=0.33333,q =(probability of rewiring a link)0.333335),we iteratively modified the parame-ter values until we obtained scale-free networks which also had an average connectivity as close as possibleJ.Hallinan /BioSystems 74(2004)51–6255Table 1Characteristics of the networks used in the project NetworkNodes Edges Connectivity Diameter Cluster coefficient MeanS.D.Mean S.D.Mean S.D.Mean S.D.Mean S.D.Yeast 1471.00N/A 2770.00N/A 1.88N/A 15.00N/A 0.207800N/A Random1429.0014.382768.20 1.17 1.940.0212.200.980.0022000.001400Preferential attachment 1471.000.002892.2035.95 1.970.02 4.400.49 1.0296000.018700Gene duplication945.4057.002426.40154.162.570.0733.005.400.0000360.000045All values are averaged over five networks,except for the yeast network,which is the largest connected component of the core yeast protein–protein interaction network.to that of the yeast network.The final parameters used were m 0=2,m =1,p =0.333,q =0.334.Fig.2shows a typical example of a network grown using the preferential attachment algorithm.2.2.2.Gene duplicationGene duplication has been an important factor in the evolution of many organisms (Lynch,2002).We used a network generation algorithm based upon gene duplication,in which a gene is interpreted as a node in the network (Pastor-Satorras et al.,2002).At each time step a node is selected at random and duplicated,together with all of its links to other nodes.Links as-sociated with the new node are then added with prob-ability α,or deleted with probability δ.Fig.2.The scale-free network generated using the preferential attachment algorithm.The gene duplication model tends to generate net-works with a large number of single,unconnected nodes.The average connectivity of the largest con-nected component of the network is,therefore,con-siderably higher than the value for the network as a whole.In addition,the highly stochastic nature of the algorithm means that the average connectivity of the network varies considerably from run to run of the algorithm,particularly when generating a relatively small network.Table 2summarizes the results of 100runs of the algorithm with the parameters described in Pastor-Satorras et al.(2002).The average connectivity of the largest connected component in a gene duplication network is depen-dant upon the link deletion parameter δ.It proved56J.Hallinan /BioSystems 74(2004)51–62Table 2Node and link statistics for whole network and largest connected component of the same network,averaged over 100runs of the gene duplication algorithm with δ=0.562,N =2000and k =2.5Nodes Links Average connectivity MeanS.D.Mean S.D.Mean S.D.Whole net 200005075.82020.63 2.5379 1.0103Largest CC490.0162.24797.2162.239.51661.1643impossible to find a value of δwhich would produce a giant component corresponding to that of the yeast network.A very high value of δyields a highly frag-mented network with no single large connected com-ponent,while lower δproduced larger components with connectivity somewhat higher than that of the yeast network.A value of δ(0.75)which would result in a giant component with average connectivity of around 2.5(within the range of average connectivities reported for real networks)was selected empirically.In order to generate a largest connected component of a useful size,networks of 10,000nodes were gener-ated and the largest connected component extracted.Fig.3shows the largest connected component of a typical network generated using the gene duplication algorithm.2.2.3.Random networksControl networks were generated with the appropri-ate numbers of nodes and links,connected atrandomFig.3.A network generated using the gene duplication algorithm.There are 943nodes and 2500links.so that the resulting network has approximately the same average connectivity as the yeast network.These random networks do not have scale-free connectivity,and would not be expected to display any signifi-cant modularity.An example of a random network is shown in Fig.4.3.Quantifying hierarchical modularity 3.1.Iterative vector diffusionThe iterative vector diffusion algorithm operates in the context of a graph G consisting of a vertex set V(G)={v 1...v n }and an edge set E(G)={e 1...e m }where each edge consists of two vertices.The algorithm is initialized by assigning to each vertex a binary vector of length n ,initialized tov i,j =0,i =j1,i =j where,i is an index into the vector and j is the unique number assigned to a given node.This generates an initial set of n orthogonal vectors.The algorithm proceeds iteratively.At each itera-tion an edge from the network is selected at random and the vectors associated with each of its nodes are moved towards each other by adding a small amount,δ,to each element of the vector.This vector diffusion process is iterated until a stopping criterion is met.We chose to compute a maximum number of iterations as the stopping criterion.This number,n ,is dependant upon both the number of connections in the network,c ,and the size of δ,such thatn =c × αδ,where,αis the average amount by which a vector is changed in the course of the run.A value for αofJ.Hallinan/BioSystems74(2004)51–6257Fig.4.A random network with1432nodes and2770edges.0.1was selected empirically in trials on artificially generated networks.At the end of the vector diffusion process the vec-tors,initially mutually orthogonal,are clustered in n-dimensional space.To reduce the dimensionality of the data set,the vectors are then subjected to hi-erarchical clustering using the hierarchical clustering algorithm implemented by Eisen et al.(1998).This algorithm uses the Pearson correlation coefficient as a distance metric.It calculates the distance matrix for all members of the input set of vectors and then uses an agglomerative hierarchical algorithm to create a hierarchical cluster tree,in which the two closest items in the set are joined by a node of the tree,and the two items replaced by a single item representing the new node.The process iterates until only one item remains.3.2.Cluster thresholdingThe output of the cluster algorithm is a binary tree, with a single root node giving rise to two offspring nodes,each of which give rise to two child nodes of their own,and so on.The tree can,therefore,be thresholded at various levels(two parents,four par-ents,eight parents,etc.;see Fig.5)and the modularity of the network at each level can be examined.3.3.Modular coherenceThe problem of identifying modules in a network is essentially an unsupervised cluster analysis task. Nodes are identified as belonging to a given mod-ule,or cluster,on the basis of their closeness to other nodes as assessed using an appropriate metric.There are many cluster analysis algorithms,several of which have been applied to the module detection task,as dis-cussed above.Most clustering algorithms,however, will identify“clusters”in any dataset,whether or not these have any correspondence to real groupings inthe Fig.5.Thresholding a cluster tree.(a)Tree thresholded at parent level2produces two clusters,(b)the same tree thresholded at parent level3has four clusters.58J.Hallinan/BioSystems74(2004)51–62 dataset.In order to validate the output of a clusteringalgorithm,practitioners often examine measures suchas inter-and intra-cluster variance.Such measures arenot easily applied to nodes on a graph.We propose ameasure of modular coherence,which measures therelative proportions of inter-and intra-module linksand assigns a value in the range−1(no coherence)to+1(a fully connected,stand-alone subgraph).The coherence,χ,of a previously identified modulecan be defined asχ=2k in(n−1)−1nnj=1k jik jo+k ji(4)where,k i is the total number of edges between nodesin the module,n is the number of nodes in the net-work,k ji is the number of edges between node j andother nodes within the module,and k jo is the numberof edges between node j and other nodes outside themodule.Thefirst term in this equation is simply the propor-tion of possible links between the nodes comprisingthe module which actually exist;a measure of the con-nectivity within the module.The second term is theaverage proportion of edges per node which are inter-nal to the module.A highly connected node with fewexternal edges will,therefore,have a lower value of χthan a highly connected node with many external edges.χwill have a value in the range(−1,+1).The concept of modularity in a network leads nat-urally to the question of scale.At what scale shouldmodularity be sought?It is important that any charac-teristic scale for modularity in the network arise fromthe data,rather than being imposed by the investiga-tor,since the appropriate scale cannot be determined apriori.This consideration has led to the concept of hi-erarchical modularity:the idea that network modulescan occur at a range of scales,with modules higherup the hierarchy divided into smaller modules,and soon.Holme et al.(2002)consider a fundamental ques-tion in biological network analysis to be:“what thehierarchical organization of subnetworks looks like”.Rather than making an a priori decision about thescale at which network modularity should be analyzed,we propose an approach which provides an overviewof the degree of modularity present in a given networkat every possible scale of modularity.The resultinggraph facilitates visual inspection of the modularity of the network over all possible scales,and permits the selection of a specific“characteristic scale”of the network for further analysis,if required.We call this approach a“coherence profile”.At each level in the hierarchy the number of mod-ules and the average modular coherence of the network was computed.Average coherence was then plotted against threshold level to produce the coherence pro-file summarizing the hierarchical modularity of the network.4.ResultsThe coherence profile for the yeast network is shown in Fig.6.It is immediately apparent that the yeast network has significant positive modular coherence over most of the range of thresholds.At low threshold values, corresponding to a partitioning of the network into a small number of relatively large modules,average coherence dips below0,indicating that the modules have more external than internal connectivity.This is because a clustering algorithm is part of the module identification algorithm.Any clustering algorithm will identify“clusters”in whatever data it is given;whether or not they reflect real modules in the biological net-work.Only when the measured modular coherence is positive can confidence be placed in the biological re-ality of the modules.The spurious nature of the results at high threshold levels is also indicated by the sudden increase in standard deviation at the point where mod-ular coherence drops below zero.These modules are illusory.The yeast network has approximately equal coherence over most thresholds,indicating that the network has a strongly hierarchically modular organi-zation.In contrast to the yeast network,the random network (Fig.7)shows negative coherence over most of its range.Although the clustering algorithm is still identify-ing modules,as expected,they have no coherence, have more external than internal edges,and do notfit the definition of a module given earlier.At the higher threshold levels,coherence rises slightly above zero. Inspection of the clustered network reveals that the modules detected at the extreme of the graph are an artefact of the unequal lengths of the branches of theJ.Hallinan /BioSystems 74(2004)51–6259-0.8-0.6-0.4-0.200.20.40.6ThresholdFig.6.Coherence profile for the core S.cerevisiae protein–protein interaction network.The data is the mean of 100runs of the algorithm on the same network.Dashed lines represent ±1standarddeviation.-0.4-0.3-0.2-0.100.10.20.30.4ThresholdFig.7.Coherence profile for random networks with approximately the same number of nodes and edges as the yeast network.The data is the mean of 100runs of the algorithm for each of five randomly generated networks.Dashed lines represent ±1standard deviation.cluster tree.At the extremes of the tree,there tends to be one large cluster and a number of very small (one or two node)clusters.Within the large cluster most edges will lie between nodes in the same cluster.The number of tiny clusters,most of whose edges connect to nodes external to the cluster,is too small to drive the mean coherence below zero.The random network,therefore,shows no evidence of hierarchical modular-ity,or,indeed,of significant modularity at any level.The fact that the yeast network displays hierarchical modularity,while the random network does not,pro-vides a benchmark against which to assess the biolog-ical plausibility of the network evolution algorithms discussed in the introduction.The coherence profiles for the preferential attach-ment and gene duplication networks are shown in Figs.8and 9,respectively.The preferential attachment algorithm has been shown to produce scale-free networks (Albert and Barabasi,2000).Fig.8shows,however,that these networks exhibit no sign of modularity at any level of the hierarchy.In contrast,the net-works generated by the gene duplication algo-rithm have a coherence profile very similar to that of the yeast protein–protein interaction net-work,with significant modular coherence present at almost every level of the hierarchy.It appears that gene duplication is more likely to produce a。