Bayesian statistical modelling of human protein interaction network incorporating protein disorder information

The advent of high-throughput technologies allows the large-scale identification of cellular components and their interactions. This wealth of experimental data is assembled in biological networks (transcriptional regulatory, protein-protein interaction, metabolic, signaling, phosphorylation networks etc.). The question of systems biology araises: do the structural, architectural principles of biological networks can reveal functional properties of cellular systems. Analysis of biological networks has become one of the emergent topics. It was suggested that biological networks are scale-free with the majority of nodes having a small number of connections but with relatively fewer hubs possessing a high degree of connectivity [1, 2]. Starting from the work of Jeong et al [3], which stated that highly connected proteins are more likely to be lethal when knocked out, many works showed that hubs represent vulnerable points in the cell. Previous approaches, however, used different, rather subjective definitions of hubs. Some of them defined hubs as the top 20% of high-degree nodes [4], while other defined them as those interacting with ≥ 10 partners [5].

Beyond a purely local notion of centrality such as the node degree, more global centrality measures incorporating network-wide properties have attracted much interest. E.g. betweenness centrality which measures the total number of shortest pathes going through a node [6, 7]. Further centrality measures like eccentricity, centroid etc. are described in [8, 9]. Recently, the pairwise diconnectivity index of a network's element was proposed, defined as the fraction of initially connected pairs of nodes which become disconnected if the element is removed from the network [10]. The new measure quantifies the importance of an individual element in sustaining communication in the network. Hence, different measures result in different scoring of nodes in a network and lead to different hypotheses about the importance of biological entities.

Previous approaches for the analysis of biological networks were descriptive, graph-based ones. In this paper, a statistical approach is employed, originating from social network analysis field [11]. It is important to recognize that biological processes underlying the network data as well as experimental measurements are of stochastic nature: some of the interactions might be unreliable and some of the real relations might be missed or not measured. This concern leads to the employment of the statistical approach based on probability distributions. A probabilistic model is defined which tries to predict edge probabilities and captures network properties summarizing them in form of statistical parameters. By learning the model from available network data, estimates of the parameters are obtained. I.e. the empirical network is treated as an outcome of a statistical process that generated the network data. Thus, generic questions can be addressed: How the biological network emerged? Are there local biological processes that generated the network? In the statistical approach, the global network structure emerges as an agglomeration of local structural configurations, which rely on local interdependencies of edges. In fact, it is the dependency that may explain the deviation of real biological networks from random ones.

At the heart of the statistical modelling of networks is a general family of graph probability distributions called exponential random graph models (ERGMs) or the p*-model [12]. Each possible edge in a network is represented by a probabilistic variable and the network model expresses the value of each variable as a stochastic function of network structural properties. Various dependencies among the network variables are considered reflecting the interactive nature of processes from which a given network emerges. The dependencies give rise to different local structural configurations and correspond to different classes of models. The simplest dependency structure assumes no dependencies among the network variables (all edges are independent). Allowing the probabilities of edges to be equil implies the class of Bernoulli graphs or Erdős-Rényi graphs (p0-model). The structural configuration considered in this model is just a single edge. Further on, the p1-model (dyadic independence model) assumes the dependency between reciprocated edges connecting two nodes in a directed graph [13]. The Markov random graph model considers Markov dependencies, whereby two possible edges are assumed to be dependent if they share a node. (Such type of dependency resembles a Markov spatial process). The assumption leads to the consideration in the network model of three-nodal network structures like stars (2-path, 2-star, 3-star etc.) and triangles. The Hammersley-Clifford theorem [14] yields the expression of the graph probability distribution in terms of network statistics, being counts of local network configurations. In the formulation of the model, each statistic comes with its own parameter reflecting the importance of the corresponding configuration in the network under study. They are a kind of 'regression coefficients' which, when positive, indicate that a particular network configuration is relevant for the emergence of edges in the network. The choice of statistics present in the model embodies different assumptions about relevant local processes that might generate the network. Recently, new statistics were proposed including higher-order triangulations, alternating k-stars, geometrically weighted degree (GWD) [15, 16].

Exponential random graph models for biological networks were first explored in [17]. There, several networks (RegulonDB E.coli, ChIP-chip S.cerevisiae etc.) were fitted using 2-5-degree statistics, 2-star, edges and GWD, while iteratively increasing the complexity of the models. The results show that only basic statistics, namely edges and 2-star, have a positive effect, determining the network.

In this paper, we leave apart complex structural configurations and, instead, turn to a node-oriented model, namely to p2-model, which is a generalization of the p1-model [18]. The p1-model considers dyads (a pair of directed edges between two nodes) and represents the probability of any dyad as the function of global features of the graph (density and tendency towards reciprocity) and of individual nodal features (tendencies to send edges and to receive them). The model specification thus includes two global parameters - density (also called edges) and reciprocity, as well as local, node-specific parameters: expansiveness and attractiveness. The density determines the probability of an edge between a pair of nodes irrespective of individual characteristics of the nodes and captures the overall formation of edges in the network. The attractiveness and expansiveness describe the abilities of a node to attract and to produce edges beyond the overall density, hence individual nodal contributions to the network emergence. Estimates of the expansiveness and attractiveness parameters can be used to rank the nodes and to reveal the most important nodes in the network from the perspective of the statistical modelling.

Given the node-specific parameters, the p1-model implies the conditional independence between dyads. The p2-model is an extension of the p1-model in that the nodal parameters are modelled as random effects, a formulation that makes it possible to include into the model node- and dyad-specific covariates, as will be shown below.

Although the p1-model was defined for directed graphs, it can be modified to be applicable for undirected graphs. Then, the local, node-specific parameters are called sociality and reflect the propensity of an individual node to be connected to other nodes in the network. We term social nodes those with positive sociality parameters i.e. positively influencing the formation of edges in the network.

In this paper, we apply the undirected graph p2-model to a human protein interaction network. We assess the node-specific parameters of the model and use them to infer social nodes, which are important for the emergence of the network. Thus, our view of essential nodes in biological networks is based not on the descriptive measures, but, alternatively, on the parameters of the statistical model learned from network data.

An apparent advantage of the present approach is the possibility to introduce exogeneous biological knowledge into the network model and to study it in the light of structural properties of the network. This permits to investigate important biological hypotheses and to examine, what biological properties have effect on the connectivity of proteins in the network. In this paper, we introduce into the protein interaction network model the information on protein disorder and study, if this property influences the protein interactivity. Further benefit of the node-specific sociality parameters is that they provide a basis for partitioning the nodes into groups of structural similarity. One important question is: How many such structural groups are contained in the network? The connectivity pattern between the groups displays organizational properties of the network, in abstraction from the level of individual proteins. We analyze the protein interaction network to reveal its organizational principles.

It was well established in biology that proteins fold to their unique native conformations as determined by their amino acid sequences, and that each protein's function originates from a specific three-dimentional (3D) strucuture. In the last decade, however, numerous biologicaly active proteins were found that fail to maintain stable ordered 3D-structures under physiological conditions [19–22]. These proteins are called natively unfolded or intrinsically disordered. Bioinformatics analysis has revealed that about 25-30% of eukaryotic proteins are mostly disordered, and that more than half of eukaryotic proteins have long regions of disorder [23, 24]. It was predicted that 70% of signaling proteins and vast majority of cancer- and cardiovascular disease-associated proteins have long disordered regions [25, 26]. It was shown that disordered proteins play a number of crucial roles in regulation, signaling and control processes, many post-translational modifications (ubiquitination, methylation, phosphorylation etc.) occur within the regions of intrinsic disorder [27]. In [28–30], the authors revealed 238 Swiss-Prot functional keywords positively correlated with long disordered regions in proteins. The common view now is that disordered regions carry out a large variety of functions [31, 32].

It is postulated that the intrinsic disorder plays an important role in the interactivity of proteins [33]. The intrinsically disordered proteins can bind to large numbers of diverse targets or can facilitate other proteins binding to many targets. Haynes et al demonstrate that intrinsic disorder is a distinctive and common characteristic of eukaryotic hub proteins [5]. Dosztanyi et al [34] and Patil and Nakamura [35] investigate the structural properties of hubs that enable them to interact with multiple proteins and conclude that global flexibility and extended interaction surface provided by disordered regions play a significant role in the binding ability of hubs. In this paper, we study the effect of disorder on the connectivity of proteins.

The Bayesian inference of the p2-model for an undirected network was applied to the human protein interaction network (see Data). The code for the software OpenBUGS is available on [51]. We used 50000 MCMC iterations for the 'burn-in' and 100000 iterations for the estimation of parameters. ('Burn-in' period of the MCMC are the first iterations which are used for the chain to get to convergence sufficiently away from its initial values, but are discarded in the estimations of the posterior values). The classification of proteins into ordered and disordered based on their sequences was performed with the classifier POODLE-W (see Methods). We obtained 130 proteins (out of total 401) predicted to be disordered (see Additional Files 2 and 3).

The mean (standard deviation) of the density parameter θ, estimated from the Markov chain, is -5.68 (0.13). This means that, by absence of the sociality parameters α, the probability p of each edge in the network is 0.0034 (calculated by log( ) = -5.68). This probability can produce only 287 edges out of = 80200 possible. The probability of edges further increases by positive effect of the nodes with high sociality. There were 177 nodes (from total 401) found with positive estimated means of sociality parameter α, whereas 69 of them had 2.5% quantile greater than 0. We propose to consider these 69 nodes as having the true positive influence on the formation of edges in the protein interaction network and call the respective proteins social. The identifiers of the proteins with the estimated means of their sociality parameters sorted in the decreasing order, as well as the estimated standard deviations, 2.5% and 97.5% quantiles, the proteins' degrees of connectivity in the network and the disorder class labels are demonstrated in Table 1.

Protein	Symbol	Degree	alpha mean	sd	2.5%quantile	97.5%quantile	Disorder
NP_473376.1	UNC119	56	3.48	0.13	3.14	3.70	1
CAD39125.1	RIF1	53	3.41	0.13	3.06	3.63	1
AAH12509.1	EEF1A1	49	3.30	0.14	3.03	3.57	0
AAH13918.1	EEF1G	41	3.06	0.14	2.77	3.34	0
NP_000537.2	TP53	36	2.89	0.15	2.51	3.13	1
AAK55500.1	CRMP1	36	2.89	0.15	2.59	3.18	0
BAA92615.1	KIAA1377	35	2.85	0.15	2.47	3.10	1
NP_002816.1	PTN	34	2.81	0.15	2.42	3.06	1
NP_036564.1	SETDB1	32	2.73	0.16	2.34	2.99	1
XP_351098.1	CHD3	30	2.65	0.16	2.25	2.91	1
NP_005068.2	TLE1	29	2.61	0.16	2.21	2.87	1
NP_057144.1	CGI125	28	2.56	0.17	2.23	2.88	0
NP_009107.1	C14orf1	26	2.47	0.17	2.13	2.80	0
NP_874368.1	HTATIP	25	2.42	0.17	2.08	2.74	0
NP_002961.1	SAT	23	2.31	0.18	1.96	2.65	0
NP_003371.1	VIM	22	2.26	0.18	1.82	2.56	1
NP_004146.1	SERPINB9	20	2.14	0.19	1.76	2.49	0
AAH33561.1	DKFZP564O0523	19	2.08	0.19	1.69	2.44	0
NP_006824.2	COPS6	18	2.01	0.20	1.62	2.39	0
O00231	PSMD11	18	2.01	0.19	1.62	2.38	0
NP_009153.2	ZHX1	17	1.95	0.20	1.48	2.28	1
NP_002564.1	PAFAH1B3	17	1.94	0.20	1.55	2.32	0
NP_008896.1	ZNF24	16	1.87	0.20	1.39	2.22	1
NP_001060.1	TUBB	15	1.79	0.21	1.37	2.19	0
NP_000994.1	RPLP1	14	1.72	0.21	1.21	2.08	1
NP_006640.2	SDCCAG16	13	1.63	0.22	1.11	2.00	1
NP_777576.1	UBR1	13	1.63	0.22	1.18	2.04	0
P04183	TK1	13	1.63	0.22	1.17	2.05	0
NP_002937.1	RPA2	13	1.62	0.22	1.18	2.05	0
NP_006312.1	ARIH2	12	1.53	0.23	1.07	1.97	0
NP_002281.1	LAMA4	12	1.53	0.23	1.06	1.96	0
NP_060267.2	BTBD2	11	1.44	0.24	0.96	1.88	0
NP_003940.2	HAP1	11	1.43	0.24	0.88	1.83	1
Q9Y2X7	GIT1	11	1.43	0.23	0.96	1.88	0
NP_002037.2	GAPD	11	1.43	0.24	0.95	1.88	0
NP_004630.2	BAT3	10	1.33	0.25	0.75	1.74	1
Q13332	PTPRS	10	1.32	0.24	0.83	1.78	0
NP_005878.1	DLEU1	10	1.32	0.25	0.81	1.79	0
BAB14293.1	ASC1p100	10	1.32	0.24	0.82	1.79	0
NP_001316.1	CSTF2	9	1.21	0.26	0.62	1.65	1
NP_002148.1	HSPE1	9	1.21	0.25	0.68	1.69	0
NP_061730.1	PCDHA4	9	1.20	0.26	0.68	1.69	0
CAD97612.1	IMMT	8	1.08	0.27	0.45	1.54	1
NP_057103.1	LUC7L2	8	1.08	0.27	0.45	1.53	1
CAB72445.1	BRD7	8	1.08	0.27	0.45	1.53	1
NP_061960.1	ARFRP2	8	1.08	0.27	0.53	1.58	0
NP_002203.1	ITGB4BP	8	1.07	0.27	0.52	1.57	0
NP_005997.2	ZNF145	8	1.07	0.27	0.52	1.58	0
NP_008998.1	MYST2	7	0.93	0.29	0.28	1.42	1
NP_060719.3	CDK5RAP2	7	0.93	0.29	0.28	1.41	1
AAB96331.1	APLP1	7	0.93	0.29	0.27	1.41	1
NP_001253.1	CDKN2C	7	0.93	0.29	0.34	1.46	0
NP_036569.1	SNAPAP	7	0.93	0.28	0.34	1.46	0
NP_000362.1	TTR	7	0.93	0.29	0.34	1.46	0
Q96RU7	C20orf97	7	0.93	0.28	0.34	1.46	0
AAH33094.1	IKBKAP	7	0.93	0.29	0.34	1.46	0
NP_004441.1	ERH	7	0.93	0.29	0.33	1.46	0
NP_005774.2	APACD	7	0.93	0.28	0.35	1.45	0
NP_001782.1	CDC42	7	0.93	0.29	0.33	1.46	0
NP_002613.2	PFDN1	6	0.76	0.31	0.07	1.28	1
NP_004302.1	ARL3	6	0.76	0.30	0.14	1.33	0
NP_002680.2	POLA2	6	0.76	0.30	0.13	1.32	0
NP_001460.1	G22P1	6	0.76	0.30	0.13	1.33	0
NP_444252.1	PFN2	6	0.76	0.30	0.13	1.33	0
NP_006077.1	TUBB4	6	0.76	0.31	0.13	1.33	0
NP_005251.1	GDF9	6	0.76	0.30	0.13	1.32	0
NP_002936.1	RPA1	6	0.76	0.30	0.13	1.32	0
NP_071921.1	FTS	6	0.76	0.30	0.13	1.32	0
AAH08720.1	CRELD1	6	0.76	0.30	0.13	1.31	0

Table 1 contains the list of 69 social proteins revealed by the statistical analysis. NCBI identifiers (RefSeqIds) and symbols (as used by [36]) of the proteins are displayed. For each protein, the degree of the respective node in the network is showed. Columns 4-7 present the estimated mean, standard deviation, 2.5% and 97.5% quantiles of the sociality parameter α. The last column contains the disorder group of each protein. (The proteins are sorted in the decreasing order of the mean of the sociality parameter.)

The degrees of connectivity of the social proteins are greater or equal 6. The definition of hubs as the top 20% of high-degree proteins would suggest 80 proteins, hence also with degree 5.

From the 69 social proteins, 23 are disordered. The corresponding functional keywords obtained from Swiss-Prot [52] are presented in the Additional File 4. It can be seen, that these proteins are involved in biological processes like transcription, transcription regulation, cell cycle, mRNA processing, ribosome biogenesis, Wnt signaling. The proteins possess coiled coil, Zinc-finger, repeat domains, as well as homeobox, EGF-like and bromodomain. The majority of proteins are associated with alternative splicing and polymorphism i.e. coding sequence diversity that potentially relies on protein intrinsic disorder. The fact that alternatively splicing events map to regions of disorder much more often than to regions of structure was emphasized in [53]. The combination of disorder and alternative splicing was proposed to provide a mechanism for generating signaling diversity and enabling the evolution of cell differentiation and multicellularity. Alternative splicing and polymorphism are proposed to increase proteome size and regulatory and signaling network complexity, thus generating organism complexity.

We used the p2-model in order to incorporate in it exogeneous nodal attributes, namely the disorder groups of the respective proteins. The parameter γ reflecting the effect of the disorder on the sociality of the nodes in the network was estimated 0.06 (0.095), thus revealing that protein disorder has the positive effect on the nodal sociality, and hence on the probability of edges in the network. However, the 2.5% quantile of γ was estimated as -0.105 and the 97.5% quantile was 0.275. This indicates that the effect of the disorder is still small and cannot fully explain the sociality of nodes in the network. The parameters a, a kind of 'residuals' after considering the effect of disorder, appeared to have high estimated values. The present finding supports the hypotheses that proteins folding property disorder plays a role in the protein network achitecture, however, other reasons enabling proteins to interact with large numbers of proteins should also be considered.

The goodness-of-fit plot for the degree statistics, which compares the degree distributions of the observed and simulated networks, is presented in the Additional File 5. There were 100 networks generated from the model. It can be seen, that the model slightly under-predicts the degrees of the nodes, and hence, the number of edges in the network.

Further on, we calculated the Euclidian distance matrix between the sociality parameters of all the nodes in the network, and used it to cluster the nodes by hierarchical complete-linkage clustering. The heatmap with the cluster diagram on top is presented in Figure 1. The clustering suggests five groups of nodes, whereby the nodes in each group possess similar patterns of connetivity. Table 2 lists the clusters, their sizes and the range of degrees of the nodes contained in each cluster, respectively. Moreover, to investigate the common functional properties of proteins contained in each cluster, we performed the gene ontology analysis of the respective genes. Namely, we obtained the GO terms significantly enriched in each cluster as compared to other proteins in the set (we employed the R package GOSim [54]). The taxonomy molecular function ("MF") was used. The enriched GO terms for each cluster are displayed in the last column of Table 2. The inspection of the GO terms suggests that "very high"-connected proteins (Cluster 1) are involved in remodeling of chromatin structure, modifications of histones, translation elongation, DNA methylation, DNA replication, alternative splicing, thus related to transcriptional machinery i.e. to fundamental genetic processes. The proteins are associated with development, differentiation, response to diverse cellular stresses etc. Cluster 2 ("high"-degree proteins) probably contains proteins with a broad diversity of functions - only one enriched GO term "zinc ion binding" was found. Cluster 3 with "moderate"-degree proteins is associated with signal transduction and, via hormone binding, with metabolic processes. Clusters 4 and 5 ("low"- and "very low"-connected proteins) are involved in the specific regulation of transcription. The proteins exhibit transcription factor activity. One additional function is the structural molecule activity i.e. contribution to the structural integrity of a complex or assembly.

Cluster	Size	Degrees	GO Terms enriched (Molecular Function, "MF")
Cluster 1	16	22-56	GO:0000739 DNA strand annealing activity
			GO:0003682 chromatin binding
			GO:0003746 translation elongation factor activity
			GO:0004145 diamine N-acetyltransferase activity
			GO:0004157 dihydropyrimidinase
			GO:0004864 protein phosphatase inhibitor activity
			GO:0016455 RNA polymerase II transcription mediator activity
			GO:0018024 histone-lysine N-methyltransferase activity
			GO:0050681 androgen receptor binding
Cluster 2	32	8-20	GO:0008270 zinc ion binding
Cluster 3	86	4-7	GO:0003899 DNA-directed RNA polymerase activity
			GO:0005125 cytokine activity
			GO:0019900 kinase binding
			GO:0042562 hormone binding
Cluster 4	131	2-3	GO:0003777 microtubule motor activity
			GO:0005088 Ras guanyl-nucleotide exchange factor activity
Cluster 5	136	1	GO:0003700 transcription factor activity
			GO:0003723 RNA binding
			GO:0003779 actin binding
			GO:0005198 structural molecule activity

Table 2 presents the groups of proteins obtained by clustering the respective sociality parameters. The number of proteins in each cluster, their degrees and molecular functions found by GO-significance analysis are displayed.

We summarized the whole protein network as a diagram representing clusters (i.e. groups of structurally similar proteins) and relations between them (Figure 2). The size and grey color intensity of the nodes in the diagram reflect the degrees of the network nodes contained in each cluster. The width of the ties in the diagram reflects the number of network interactions between the nodes belonging to two clusters (presented as percentage of the total 911 interactions). Interactions between nodes of the same cluster are displayed with the arrows. All protein interactions are now distributed as flows in the diagram. The diagram reveals the global organization of the network, abstracted from the level of individual nodes. The flows indicate that Cluster 1 is tightly connected to all other clusters, thus playing a role important to all proteins in the network. Indeed, it is associated with the basal, generally important process - transcriptional control. Interestingly, Cluster 2 is much higher connected to itself than other clusters, depicting the presence of a tight middle-layer of proteins in the network. Longer path "moderate"-"very high"-"low"-degree proteins might indicate a certain sequence of processes in the proteins interactivity. It can be postulated that the connectivity of proteins is related to the generality (or specificity) of their functions.

Although we used a subpart of all available protein interactions, the most reliable one, we believe that the tendencies revealed here will remain valid, as the biological knowledge will accumulate. Our method can be readily applied to larger networks.

We proposed a method for the analysis of biological networks based on the statistical modelling with the exponential random graph model, namely p2-model. We applied the method to a human protein interaction network. Bayesian approach enabled flexible model definitions and facilitated the model inference.

The statistical modelling handles the observed network as a probabilistic outcome of an underlying process of network formation. It allows to make explicit assumptions on the model that generated the network data and precise statements on the validity of hypotheses about the network.

In contrast to the descriptive methods of the analysis of biological networks, the ERGM approach attempts to address the origin of global network architecture and to explain the emergence of the network from local interaction patterns. It is more biologically plausible, that biological phenomena are reflected in the local network structural configurations, than that some proteins are responsible for the organization of the global network structure. It is a good challenge to reveal, what biological phenomena or aspects of cellular life are essential from the local structural point of view.

Previously, the problem of network formation was regarded by the evolutionary approaches. The duplication-divergence model considered gene duplications and gain and loss of protein interactions through mutations as the driving force shaping the structure of the network ([55, 56]). In [57], the domain-based evolutionary model was introduced. The goal was mainly to explain (or reproduce) the scale-free property of the protein network.

The exponential random graph model tries to describe the structural properties of the network with a limited number of parameters. Here, we used one global parameter representing the density of interactions in a network, and local, node-specific parameters reflecting the propensity of each individual node to connect to other nodes. i.e. the model is node-oriented. When positive, the node-specific parameters indicate that the respective nodes contribute to the formation of edges in the network. This feature suggests a novel definition of essential nodes in the network, namely social, as revealed based on the parameters of the fitted model. Previously, hubs, as highly interactive nodes occupying central positions in the network, were considered to be essential as they are important for the maintenance of the global network architecture - the scale-free architecture. Removing hubs increases the proportion of unreachable nodes and the shortest path lengthes in the network more than removing non-hubs [58]. The previous definitions of essential nodes in biological networks are based on descriptive measures, in contrast to the generic approach proposed here.

The main advantage of our approach is that it provides a possibility for biological hypothesis testing. Attributes of nodes (as well as of edges) can be inserted into the model to explore the role of these properties in the emergence of the network. Here, we examined the influence of protein disorder on the nodal connectivity and revealed the positive effect of the disorder. We confirmed the recent idea that protein disorder provides a mechanism facilitating proteins binding diversity [33]. However, protein disorder exhibits an effect, but do not fully explains the proteins connectivity, henceforth other information on proteins should be examined in the context of the protein interaction network.

Notably, previous approaches were based on the paradigm "from the structural importance of nodes to the functional". Our work challenges an opposed way of thinking, and enables to deduce, which properties of biological entities are important for the structure of the biological network.

In [59], the role of intrinsic disorder was studied with respect to the "date" and "party" classifications of proteins. Date proteins are those building transient interactions with their partners, while party proteins interact with their partners simultaneously (this classification of proteins was done based on the correlation of mRNA expression). It was suggested that intrinsic disorder in date proteins may facilitate transient interactions. The date and party classification of human proteins can be incorporated in our framework together with the disorder information, to explore the role of both properties in the organization of the protein interaction network.

The results of [60] indicate that multi-domain proteins are overrepresented among hubs in the interaction network of yeast. Also, repeated structural motifs are enriched in hub proteins. E.g. the WD40-repeat protein acts in the formation of a multi-protein corepressor complex, where the repeating units serve as a rigid scaffold for protein interactions. The information about domain content and repeats would be valuable to consider in our modelling.

Whereas previous works firstly determined hubs by setting an arbitrary cutoff and then related particular biological properties to them (e.g. "date" and "party" hubs, singlish- and multi-interface hubs, see [59, 61, 62]), our approach allows each protein in the network to exhibit an impact according to its biological feature.

Further examples of biological properties, which could be incorporated in our model, are genes biological essentiality, as well as proteins functional groups. The relation between two proteins is actually explained by a biological function they are involved in (regulatory, signalling, metabolic etc.). Inference of this function remains a great challenge. The proposed statistical modelling approach allows to check, what qualities of proteins stimulate the interactions between them, what functional characteristics determine the emergence and the structure of the protein network.

We provided a novel way of revealing the organizational principles of a biological network as being constructed from a few building blocks comprizing individuals with equivalent connectivity.

The ERGM approach provides a link between the network motifs research emergent in the systems biology and the generic perspective. Various types of network motifs significantly enriched as compared to randomized networks were revealed by the study of transcription regulatory network of E.coli: feed-forward loops (FFL), single-input modules (SIM), dense overlapping regulons (DOR) [63]. The ERGM modelling of networks offers a natural way of assessing importance of the network motifs. Moreover, it should explain the biological reasons for the formation of motifs. Using the terms of the statistical network analysis field, the feed-forward loop is the triangle, the SIM is the out-star. The work of [17] was a first attempt to learn the higher-order exponential random graph models from biological networks and to obtain the importance of these structural configurations. However, only basic configurations (edges and 2-pathes) showed a positive propensity to be present in the networks. This suggests that further configurations and their respective statistics, characteristic to the biological networks, must be developed (e.g. longer paths). Novel, four-element motifs Bi-fan, Bi-FFL and Bi-parallel were found recently in the combined transcription factor-binding and phosphorylation network of yeast [64].

To summarize, the statistical modelling of networks presented here provides a framework combining three research directions: 1) the generic analysis of networks, 2) the structural analysis of networks and 3) the biological analysis of individual entities, all aiming at understanding the functional principles of complex biological systems.

The approach described here is very promising for the comparative analysis of networks involving the same nodes, since it provides the parameters on which the networks can be compared. An important application will be the cross-species comparative analysis of biological networks (proteins interaction, metabolic etc.) involving conserved nodes. Also, the approach provides new perspectives on the dynamic network analysis that investigates how network structure changes over time.