Mean Field Analysis of Algorithms for Scale-free Networks in Molecular Biology

TL;DR

This paper provides a mean field analysis of algorithms for generating scale-free networks in molecular biology, deriving expressions for the degree distribution exponent and validating them through simulations.

Contribution

It extends mean field analysis to various algorithms, including biological-inspired ones, and assesses their ability to produce scale-free networks.

Findings

Derived mean field expressions for the power-law exponent γ.

Validated theoretical predictions with numerical simulations.

Identified certain algorithms that do not generate scale-free networks.

Abstract

The sampling of scale-free networks in Molecular Biology is usually achieved by growing networks from a seed using recursive algorithms with elementary moves which include the addition and deletion of nodes and bonds. These algorithms include the Barabasi-Albert algorithm. Later algorithms, such as the Duplication-Divergence algorithm, the Sol\'e algorithm and the iSite algorithm, were inspired by biological processes underlying the evolution of protein networks, and the networks they produce differ essentially from networks grown by the Barabasi-Albert algorithm. In this paper the mean field analysis of these algorithms is reconsidered, and extended to variant and modified implementations of the algorithms. The degree sequences of scale-free networks decay according to a powerlaw distribution, namely $P(k) \sim k^{-\gamma}$, where $\gamma$ is a scaling exponent. We derive mean field expressions for $\gamma$, and test these by numerical simulations. Generally, good agreement is obtained. We also found that some algorithms do not produce scale-free networks (for example some variant Barabasi-Albert and Sol\'e networks).

Tables

Table 1. Table 1: Connectivity data for Duplication-Divergence Networks.

$n$	$q=0.4$	$q=0.5$	$q=0.6$
$3125$	$25.9$	$11.4$	$5.93$
$6250$	$31.3$	$12.6$	$6.14$
$12500$	$37.3$	$13.6$	$6.33$
$25000$	$44.8$	$14.4$	$6.55$
$50000$	$51.9$	$15.5$	$6.64$
$100000$	$60.1$	$16.8$	$6.75$
$200000$	$70.3$	$17.7$	$6.88$

Table 2. Table 2: Connectivity data for Solé Networks.

$n$	$\delta =0.25$	$\delta =0.75$
$100$	$2.95$	$1.50$
$200$	$4.46$	$1.94$
$400$	$7.59$	$2.94$
$800$	$14.75$	$5.36$
$1600$	$30.46$	$10.64$
$3200$	$59.94$	$21.26$
$6400$	$122.78$	$45.57$
$12800$	$245.35$	$85.18$
$25600$	$496.87$	$170.35$
$51200$	$994.54$	$340.76$

Table 3. Table 3: Connectivity data for iSite Networks.

$n$	Column 2	Column 3	Column 4	Column 5
$3125$	$22.385$	$20.701$	$4.756$	$6.648$
$6250$	$26.524$	$25.752$	$4.770$	$6.556$
$12500$	$31.395$	$29.137$	$4.677$	$6.579$
$25000$	$37.808$	$35.308$	$4.733$	$6.358$
$50000$	$45.931$	$42.244$	$4.579$	$6.299$
$100000$	$54.830$	$50.035$	$4.584$	$6.204$
$200000$	$64.668$	$59.284$	$4.649$	$6.071$
Column 2:	$I=3$, $p=0.5$	$q=0.4$, $r=0.3$
Column 3:	$I=5$, $p=0.5$	$q=0.4$, $r=0.3$
Column 4:	$I=3$, $p=0.5$	$q=0.05$, $r=0.8$
Column 5:	$I=5$, $p=0.5$	$q=0.05$, $r=0.8$

Table 4. Table 4: Computational Time Complexity of Implemented Algorithms.

Algorithm	$n=6250$	$n=12500$	$n=25000$	$n=50000$	$\tau$
Bar-Alb ($p=0$)	$0.602$	$2.51$	$9.03$	$38.0$	$1.97$
Mod Bar-Alb ($\lambda =2$, $p=A=0$)	$0.618$	$2.55$	$10.1$	$36.3$	$1.96$
Var Bar-Alb ($\alpha =2$, $a=0$)	$1.35$	$4.46$	$16.4$	$--$	$--$
Dupl-Div ($p=1$, $q=0.4$)	$0.349$	$0.862$	$2.04$	$5.01$	$1.28$
Dupl-Div ($p=1$, $q=0.6$)	$0.155$	$0.319$	$0.635$	$1.31$	$1.02$
Mod Dupl-Div ($p=1$, $q=0.4$)	$0.340$	$0.891$	$2.45$	$7.09$	$1.46$
Mod Dupl-Div ($p=1$, $q=0.6$)	$0.165$	$0.338$	$0.699$	$1.44$	$1.04$
Solé ($\delta =0.25$, $\alpha =0.005$)	$4.84$	$20.5$	$91.0$	$436.0$	$2.16$
Solé ($\delta =0.75$, $\alpha =0.005$)	$6.10$	$20.0$	$79.5$	$323.2$	$1.92$
iSite ($p=0.5$, $q=0.01$, $r=0.8$, $I=1$)	$0.114$	$0.234$	$0.454$	$0.925$	$1.00$
iSite ($p=0.5$, $q=0.01$, $r=0.8$, $I=2$)	$0.110$	$0.216$	$0.458$	$0.878$	$1.01$
iSite ($p=0.5$, $q=0.01$, $r=0.8$, $I=3$)	$0.106$	$0.217$	$0.432$	$0.857$	$1.00$
iSite ($p=0.5$, $q=0.01$, $r=0.8$, $I=4$)	$0.107$	$0.231$	$0.422$	$0.848$	$0.98$
iSite ($p=0.25$, $q=0.01$, $r=0.8$, $I=4$)	$0.104$	$0.249$	$0.415$	$0.844$	$0.98$
iSite ($p=0.75$, $q=0.01$, $r=0.8$, $I=4$)	$0.108$	$0.216$	$0.437$	$0.867$	$1.00$
Mod iSite ($p=0.5$, $q=0.1$, $r=0.8$, $s=0.1$, $I=4$)	$0.288$	$0.560$	$1.102$	$2.53$	$1.04$

Equations

$$P(k)\simeq C_o^{-1}{k}^{-\gamma }.$$

$$C_o=\sum _{k=1}^{n}P(k).$$

$$P(k)\sim \frac{1}{\log n}{k}^{-1}$$

$$\zeta (k)=\log P(\alpha k)-\log P(k)=-\gamma \log \alpha +o(1).$$

$$⟨k{⟩}_n$$

$$\langle k\rangle_{n}\simeq\cases{\hbox{\large$\frac{n}{\log n}$},&\hbox{if $\gamma=1$};\\ \log n,&\hbox{if $\gamma=2$}.}$$

$$\displaystyle E_{n}=\cases{\hbox{\large$\frac{n^{2}}{2\log n}$},&\hbox{if $\gamma=1$};\\ \left(\hbox{\large$\frac{\gamma-1}{2(\gamma-2)}$}\right)n^{3-\gamma},&\hbox{if $1<\gamma<2$};\\ \frac{1}{2}n\log n,&\hbox{if $\gamma=2$};\\ \left(\hbox{\large$\frac{\gamma-1}{2(\gamma-2)}$}\right)n,&\hbox{if $\gamma>2$} .}$$

$$2E_n=\sum _j{k}_j(n)$$

$$\Delta E_n=p+(1-p)\lambda +(1-p)\frac{nA}{2E_n}.$$

$$2E_n\frac{d}{dn}{E}_n=2(p+(1-p)\lambda )E_n+(1-p)nA.$$

$$E_n=\frac{n}{2}((p+(1-p)\lambda )+\sqrt{(p+(1-p)\lambda {)}^2+2(1-p)A})=Cn,$$

$$k_j(n+1)=k_j(n)+\frac{p}{n}+\frac{(1-p)(\lambda k_j(n)+A)}{2E_n}$$

$$\frac{d}{dt}{k}_j(t)=\frac{p}{t}+\frac{(1-p)(\lambda k_j(t)+A)}{2Ct}.$$

$$\frac{d}{dt}{k}_j(t)=\frac{p}{t}+\frac{(1-p)k_j(t)}{2t}$$

$$\frac{d}{dt}{k}_j(t)=\frac{Q}{t}+\frac{P}{t}{k}_j(t)$$

$$k_j(t)=\left(1+\frac{Q}{P}\right)(t/t_j{)}^P-\frac{Q}{P}$$

$$\left(1+\frac{Q}{P}\right)(t/t_j{)}^P<\kappa or,\; equivalently,t_j>t{\left(\frac{Q/P+\kappa }{1+Q/P}\right)}^{-1/P},$$

$$P[k_j(t)<\kappa ]=P\left[(t_j/t)>{\left(\frac{Q/P+\kappa }{1+Q/P}\right)}^{-1/P}\right]=1-{\left(\frac{Q/P+\kappa }{1+Q/P}\right)}^{-1/P}.$$

$$P(\kappa )=P[k_j(t)=\kappa ]=\frac{\partial }{\partial \kappa }P[k_n(t)<\kappa ]=\frac{(P+Q{)}^{1/P}}{(P\kappa +Q{)}^{1+1/P}}.$$

$$\gamma$$

$$\gamma =3+\frac{2p}{(1-p)\lambda }.$$

$$\gamma =1+\frac{1}{1-p}+\frac{\sqrt{1+2(1-p)A}}{1-p}.$$

$$\gamma =\frac{3-p}{1-p}.$$

$$⟨k{⟩}_n\simeq \frac{{\int }_1^{n}kP(k)dk}{{\int }_1^{n}P(k)dk}\simeq \frac{2C}{2C-(1-p)\lambda },$$

$$\log P(k)-\log P(\frac{1}{2}k)=-\gamma \log 2+o(1)$$

$$\Delta E_n=p+\frac{(1-p)({\sum }_j(k_j(n){)}^{\alpha }+A)}{2E_n}.$$

$$2E_n\frac{d}{dn}{E}_n=2p{E}_n+(1-p)nA+(1-p)\sum _j(k_j(n){)}^{\alpha }.$$

$$2E_n\frac{d}{dn}{E}_n\approx 2p{E}_n+(1-p)nA+(1-p)n^{1-\alpha }(2E_n{)}^{\alpha }.$$

$$k_j(t+1)=k_j(t)+\frac{p}{t}+\frac{(1-p)((k_j(t){)}^{\alpha }+A)}{2D{t}^{\beta }}.$$

$$\frac{d}{dt}{k}_j(t)=\frac{p}{t}+\frac{(1-p)((k_j(t){)}^{\alpha }+A)}{2D{t}^{\beta }}.$$

Full text

Keywords: Biological networks, Network algorithms, Scale-free networks

1 Introduction

Many systems in nature and society are described by means of complex networks [9]. Some of these systems include the cell [18], chemical reactions [16], the world wide web [6], social interactions [7], etc. It is generally found that many system, though different in nature, produce networks which are scale-free and exhibit similar properties [2, 3].

The main property of scale-free networks is that their degree distribution decays as a power law [2, 4] – this shows that there is no characteristic scale for the degrees, which is why the networks are called scale-free. The average degree of a scale-free network offers little insight into the real topology of the network [3] since most nodes have degrees which are far away from the average degree of the network. Nodes of high degree are called hubs and though small in number for realistic networks, they are over-represented compared to the number of hubs in random networks. These hubs play an important role in dynamical processes which occur in scale-free networks.

Scale-free networks also exhibit an unexpected degree of robustness – this is the property that such networks maintain their dynamic properties even when many nodes and bonds fail to transmit signals (suffer high failure rates) [9]. However, these networks remain vulnerable to failure of hub nodes, since these nodes play a significant role in maintaining the network’s connectivity.

In this paper the mean field approach to the analysis of algorithms for sampling scale-free networks inspired by processes in molecular biology is presented. In addition, numerical testing and, in some cases, verification, of the mean field approach will be examined. The focus will be on four widely used and discussed algorithms in the literature, nameley, the Barabasi-Albert algorithm [4, 5], the Duplication-Divergence algorithm [26, 27], the Solé algorithm [24] and the iSite algorithm [14, 15].

The Duplication-Divergence, Solé and iSite algorithms are inspired by modelling networks in biological models of protein-protein interaction evolution, and all these algorithms are based in one way or another on two ideas: growth by preferential attachment [11], and growth and changes (mutations) in networks induced by the duplication, deletion or replacement of nodes or bonds (these are elementary moves which mutate the network by adding, deleting or moving some of its bonds or nodes).

Growth by preferential attachment is implemented by adding bonds preferentially to nodes of high degree. This increases the probability that a node will grow to be a hub in the network, and the resulting network has an increased probability that it will contain hubs [4]. The Barabasi-Albert algorithm uses preferential attachment to grow scale-free networks by attaching bonds to nodes with a probability which is proportional to the degrees of nodes [2]. A mean field analysis of the Barabasi-Albert algorithm was done in reference [5].

The Duplication-Divergence algorithm [27, 26] generates scale-free networks by implementing elementary moves which mutate and grow the network. These are duplication (the duplication of existing nodes and bonds) and divergence (local changes made to existing bonds and nodes) elementary moves. These moves model processes which are thought to underlie the evolutionary mechanisms by which protein interaction networks evolve [27, 23, 26]: The duplication of genes is a mechanism which generates genes coding for new proteins during evolution and the divergence step is a model for the mutation of duplicated genes. After a duplication of a gene, two genes (one the progenitor gene, the other the progeny gene) coding for the same protein are obtained, and these mutate over time to drift away from one another in gene space, giving rise to modified proteins when translated by cellular machinery [23]. Biologically, the duplication step may result in a new protein interaction between two mutating copies of the same gene (this is called heteromerization), and the divergence step is a model of subfunctionalization (a process where interactions between proteins are lost).

Closely related to the Duplication-Divergence algorithm is the Solé algorithm [24, 23]. This algorithm grows networks by duplication of nodes, and mutates the network by rewiring it (this algorithm does not implement the heteromerization of the duplicated genes) [6]. It then implements a process of deleting some bonds on the duplicated nodes (modelling evolutionary changes due to subfunctionalization).

The iSite algorithm [14, 15] is a refinement of the Duplication-Divergence and Solé algorithms. This algorithm introduces more complex nodes which each contains interaction sites as models of protein and protein complexes with localized interaction sites where the interactions with other proteins take place. These localized interaction sites are iSites. Such iSites may be involved in many interactions, but each interaction is related to only two iSites, one on each of the proteins involved. That is, iSites are models of the concept of domains on protein surfaces where the actual interactions take place between two proteins. The implementation of the algorithm on nodes containing iSites proceeds by duplication of nodes, and the mutation of iSites through subfunctionalization and heteromerization (namely, the subfunctionalization of iSites leading to loss of protein interactions, and heteromerization where new interactions are introduced between existing iSites). In this model the subfunctionalization is of iSites, leading to the loss of all bonds incident with the iSite (contrary to the situation in other algorithms, for example the Duplication-Divergence algorithm, where subfunctionalization leads to the loss of bonds, rather than nodes).

This paper is organised as follows. We first consider the general properties of scale-free networks, including their scaling and connectivity properties. These ideas are then applied to the analysis of particular algorithms. The Barabasi-Albert model is considered first together with a modified version of the algorithm, and a variant of the algorithm. Mean field theory for the modified and variant algorithms is developed, giving mean field values for the scaling exponent $\gamma$. These results are compared to numerical results obtained by generating networks using implementations of the algorithms.

The Duplication-Divergence algorithm and networks generated by it are considered next. The algorithm is also modified, and mean field theory is developed to find mean field values for the scaling exponent. The mean field predictions are then compared to numerical results generated by implementing the algorithm and sampling networks.

A similar approach is followed for the Solé algorithm. However, in this model the degree distribution is not integrable, and our results indicate that the networks generated by this algorithm are not scale-free. Instead, the degree distribution must be modified. This gives a testable scaling hypothesis for Solé networks, which is tested numerically by generating networks and examining their scaling, as well as by computing the connectivity of Solé networks and comparing it to the mean field predictions. This shows that the size of Solé networks of order $n$ is $O(n^{2})$, while the connectivity is $O(n)$ – this implies that Solé networks are rich in bonds.

Finally, the iSite algorithm is presented and examined developing a mean field approach to determine its scaling properties. The algorithm is also modified, and the resulting mean field results are tested numerically.

The paper is completed in the conclusion section, where our main results are briefly considered and reviewed.

2 Scale-free networks

Scale-free networks of order $n$ are characterised by degree sequences $\{d_{k}\}$ which follow a power law distribution (where $d_{k}$ is the number of nodes of degree $k$ and \hbox{\normalsize\frac{1}{n}}d_{k} is the fraction of nodes of degree $k$).

If $\langle d_{k}\rangle$ is the average degree distribution, then \hbox{\normalsize\frac{1}{n}}\langle d_{k}\rangle is proportional to the probability $P(k)$ that a node has degree $k$. In scale-free networks, the probability $P(k)$ decays like a powerlaw with exponent $\gamma$:

[TABLE]

Here, $\gamma$ is the scale-free network exponent. The constant $C_{o}$ is a normalisation constant given by

[TABLE]

As $n\to\infty$, it is necessary that $\gamma>1$ for $P(k)$ to be summable (and $C_{o}<\infty$). In this case $C_{o}$ converges to a constant as $n\to\infty$. Thus, if $\gamma>1$ then the network is said to be integrable with scaling exponent $\gamma$ (in this event equation (1) is the scaling of the limiting degree distribution with $C_{o}>0$ finite and $P(k)\to 0$ as $k\to\infty$).

The case that $\gamma=1$ gives rise to a logarithmic correction. Since $\sum_{k=1}^{n}k^{-1}\sim\log n$, this gives the distribution

[TABLE]

for networks of (large) order $n$. This network is said to be not integrable, but for asymptotic values and fixed values of $n$ the decay of $P(k)$ will appear to be proportional to $k^{-1}$.

Since $P(k)$ is the probability that a node in a network has degree $k$, the average degree sequence $\{\langle d_{k}\rangle_{n}\}$ over randomly generated networks of order $n$ is given approximately by $\langle d_{k}\rangle\sim n\>\!P(k)$, for $n$ large. It is not known that the degree sequence is self-averaging (that is, that the degree sequence $\{d_{k}\}$ has distribution $d_{k}\sim n\>\!P(k)$ as $n\to\infty$ for a single randomly generated scale-free network).

This powerlaw decay of degree sequences shows that nodes of large degree (that is, for large $k$) are more common in scale-free networks (compared to randomly generated networks, where they are exponentially rare). These nodes of large degree are called hubs. A precise definition of a hub in a network is somewhat arbitrary, but for the purpose of this paper, a “hub” in a network of order $n$ is defined as a node of degree exceeding $\lfloor\sqrt{n}\rfloor$.

The exponent $\gamma$ can be estimated from numerical data by computing the average degree sequence $\{\langle d_{k}\rangle\}$ and then plotting $\log P(k)/\log k$ against $1/\log k$ (for networks of order $n\gg k$). Extrapolating the data to $k=\infty$ using a linear or a quadratic regression gives the value of $\gamma$ as the $y$-intercept of the graph. This method works well if $P(k)$ scales with $k$ as in equation (1). However, strong corrections to the powerlaw behaviour may make the extrapolation difficult or inaccurate.

A second method to estimate $\gamma$ is to note that if $\gamma>1$, then for a fixed value of $\alpha>0$,

[TABLE]

Experimentation with numerical data shows that by plotting $\zeta(k)$ against \hbox{\normalsize\frac{1}{k}}\log k good results are obtained, and linear or quadratic regressions of $\zeta(k)$ against \hbox{\normalsize\frac{1}{k}}\log k can be used to estimate $\gamma$.

If it is assumed that $P(k)$ is well approximated by equation (1) for all $k\geq 1$, then the average connectivity of a network of order $n$ with average degree distribution proportional to $P(k)=C_{o}n^{-\gamma}$ is given by

[TABLE]

Observe that the asymptotic estimate is very poor if $\gamma\approx 2$, and if $n$ is small.

The cases $\gamma=1$ and $\gamma=2$ can also be determined; this gives

[TABLE]

The coefficient $\frac{\gamma-1}{\gamma-2}$ may be modifed if $P(k)$ is not well approximated by the powerlaw decay for smaller values of $k$ in equation (1). These results, however, do show that the connectivity is a constant independent of $n$ (for large $n$) if $\gamma>2$.

The expected number of bonds in the network is given by E_{n}=\hbox{\normalsize\frac{1}{2}}n\>\!\langle k\rangle_{n}. Assuming the powerlaw relation in equation (1), it follows that

[TABLE]

Of course, if $\gamma<1$, then $E_{n}=\Theta(n^{2})$ and since a complete graph has \hbox{\normalsize\frac{1}{2}}n(n-1) bonds, this implies that these graphs are dense in the sense that \liminf_{n\to\infty}\hbox{\normalsize\frac{1}{n^{2}}}E_{n}>0. For all values of $\gamma\geq 1$ the above shows that \limsup_{n\to\infty}\hbox{\normalsize\frac{1}{n^{2}}}E_{n}=0, and the graphs are sparse.

These results are useful in examining numerical data for scale-free networks. For example, $\gamma$ can be estimated by examining degree sequences averaged over randomly sampled networks (from equation (1)), or alternatively by using equation (4). The connectivity $\langle k\rangle_{n}$ approaches a constant if $\gamma>2$ (as in equation (5)) or grows as a powerlaw with $n$ if $\gamma<2$, and with logarithmic corrections if $\gamma=1$ or $\gamma=2$ (as in equation (6)). Alternatively, the average size $E_{n}$ (the number of bonds in a network of order $n$) can be considered, using the results in equation (7).

3 Barabasi-Albert networks and the Barabasi-Albert algorithm

The Barabasi-Albert algorithm is a recursive algorithm which grows networks (or clusters of nodes and bonds) from a seed node. This algorithm was introduced in reference [4] and reviewed in 2002 in a seminal paper [2], and its elementary move was inspired by processes underlying the (presumed) evolution of scale-free networks seen in the physical world. The elementary move is a preferential attachment of new nodes (and bonds) to hubs (nodes of high degree) in the network. The algorithm is initiated by a single node, and then new nodes and bonds are recursively attached, with new bonds preferentially attached to existing nodes of large degree.

A Barabasi-Albert network of order $N$ nodes is grown as follows:

Barabasi-Albert algorithm:

1. 1.

Initiate the network with one node $x_{0}$; 2. 2.

Suppose that the network consists of nodes $\{x_{0},x_{1},\ldots,x_{n-1}\}$ of degrees $\{k_{0},k_{1},\ldots,k_{n-1}\}$; 3. 3.

Append a new node $x_{n}$ by executing step (a) or step (b):

1. (a)

With probability $p$: Select $x_{j}$ uniformly and attach $x_{n}$ to it by inserting the bond ${\langle x_{j}{\sim}x_{n}\rangle}$; 2. (b)

With default probability $1-p$: Attach $x_{n}$ by adding bonds ${\langle x_{j}{\sim}x_{n}\rangle}$ independently with probability $\frac{k_{j}}{\sum_{j}k_{j}}$; 4. 4.

Repeat step 3 until a network of order $N$ is grown. ∎

Step 3(a) is a random attachment of a node and bond, and step 3(b) attaches a node with bonds preferentially to existing nodes of high degree. The algorithm has a single parameter $p$. If $p=1$ then the algorithm grows acyclic (and connected) networks of order $N$ (these are random trees).

On the other hand, if $p=0$, then step 3(b) is executed on each iteration. New bonds are created with probabilities q_{j}=\hbox{\normalsize\frac{k_{j}}{\sum_{j}k_{j}}} for $j=0,1,\ldots,n-1$ when the $n$-th node is added. This shows that the expected number of bonds added in this step is on average $\sum_{j}q_{j}=1$. That is, on average $1$ bond is added in each iteration, and the average sum of degrees $\sum_{j}k_{j}$ should be equal to $2n$ by handshaking after $n$ iterations. This suggests that the algorithm grows a sparse graph with increasing $n$. However, since bonds are appended preferentially on growing hubs, the largest clusters in the network should be dominated by growing hubs.

For values of $p\in(0,1)$ the algorithm adds either (wih probability $p$) a single bond randomly, or it adds a collection of bonds (on average one bond) preferentially. This grows simple networks of order $N$ and size $N-1$, typically not connected unless acyclic.

In figure 1 an example of a Barbasi-Albert network of order $122$ with $p=0$ is shown (left) and the right is a network of size $380$. The appearance of hubs in these networks is clearly seen: In the network on the left there are $5$ nodes of degrees exceeding $\sqrt{122}$, the largest of degree $31$, and in the network on the right there are $3$ hubs of degrees exceeding $\sqrt{380}$, the largest of degree $63$.

3.1 Modified Barbasi-Albert networks

Barabasi-Albert networks are relatively sparse networks. A modification of the algorithm can be introduced to grow denser networks. For example, one may replace step 3(b) by

1. 1. 3(b).

With default probability $1-p$: Attach $x_{n}$ by adding bonds ${\langle x_{j}{\sim}x_{n}\rangle}$ with probability q_{j}=\min\{\hbox{\normalsize\frac{\lambda,k_{j}+A}{\sum_{j}k_{j}}},1\} (where $\lambda$ and $A$ are non-negative parameters of the algorithm);

Since $k_{j}\ll\sum_{j}k_{j}$ in Barabasi-Albert networks, one may assume that $\lambda k_{j}+A\leq\sum_{j}k_{j}$ for values of $\lambda$ and $A$ which are not too large (and so $q_{j}\leq 1$).

In figure 2 two examples of Modified Barabasi-Albert networks are shown, one a sparse network with $\lambda=0.5$, $A=0$ and $p=0$, and the second a denser network with $\lambda=2.0$, $A=0$ and $p=0$. In both cases the algorithm was iterated $200$ times; the sparse network has order $203$ and two hubs of degrees $\{15,17\}$, and the dense network has order $172$ with seven hubs of degrees $\{15,15,16,17,19,27,33\}$.

3.2 Variant Barbasi-Albert networks

A variant Barbasi-Albert algorithm can be introduced by changing step 3(b) in the Barbasi-Albert algorithm to

1. 1. 3(b).

With default probability $1-p$: Attach $x_{n}$ by adding bonds ${\langle x_{j}{\sim}x_{n}\rangle}$ with probability q_{j}=\min\{\hbox{\normalsize\frac{k_{j}^{\alpha}+A}{\sum_{j}k_{j}}},1\}, (where $\alpha$ and $A$ are non-negative parameters of the algorithm);

The effect of the parameter $\alpha$ is to increase the probability of adding bonds to the hubs of the network if $\alpha>1$, and to decrease this probability if $\alpha<1$. In the case that $\alpha>1$ networks dominated by a single very large hub are obtained (see figure 3 (right network)), while networks with $\alpha<1$ are more sparse and not dominated by a few hubs (see figure 3 (left network)). The left network in figure 3 was grown by putting $\alpha=0.15$ and $A=0$ and has order $327$. None of the nodes in this network has degree which exceeds $\sqrt{327}$, and so none qualify as hubs. A denser network is obtained if $\alpha=1.15$ and $A=0$, as shown in figure 3 on the right. This network is dominated by hubs of degrees $\{22,24,26,42,43,116\}$ and has order $351$.

3.3 Mean field theory for Modified Barabasi-Albert networks

Let $k_{j}(n)$ be the degree of node $j$ after $n$ iterations. A mean field calculation of $k_{j}(n)$ is done by assuming that $k_{j}(n)$ is equal to its expected value for each $n$; that is, $k_{j}(n)=\langle k_{j}(n)\rangle$ for each $j$ and $n$.

The modified Barabasi-Albert algorithm appends bonds to a network of order $n$ as follows: Step 3(a) is executed with probability $p$, and a bond (and the $(n+1)$-th node) is appended with uniform probability on one of the $n$ existing nodes. The probability that node $j$ gets a bond in this way is $\frac{p}{n}$ and on average one bond is attached with probability $p$.

If step 3(b) is done instead, then the expected number of bonds added in the mean field is approximately \sum_{j}\hbox{\normalsize\frac{\lambda,k_{j}(n)+A}{\sum_{j}k_{j}(n)}}=\lambda+\hbox{\normalsize\frac{nA}{\sum_{j}k_{j}(n)}}. The total number of bonds in the network is

[TABLE]

by handshaking. Thus, the increment in the number of bonds when the next node is appended is

[TABLE]

Approximate this by a differential equation

[TABLE]

This can be solved to obtain

[TABLE]

where $C$ is a function of $(p,\lambda,A)$ defined by this expression. Notice that $E_{n}$ grows linearly in $n$, so that Barabasi-Albert graphs will be necessarily sparse as $n\to\infty$ (and by equation (7) the scaling exponent is $\gamma>2$).

With each iteration the mean field value of $k_{j}(n)$ (the degree of the $j$-th node after $n$ iterations) increments by

[TABLE]

since $2E_{n}=\sum_{j}k_{j}(n)=2C\>\!n$, and since the probabilty of adding a bond to node $j$ is $\frac{\lambda k_{j}(n)+A}{\sum_{j}k_{j}(n)}$. This can again be approximated by a differential equation: Take $n\to t$, a continuous time variable, and let $k_{j}(n)\to k_{j}(t)$, the continuous mean field degree of node $j$. Then

[TABLE]

The initial condition is to assume that node $j$ is added at time $t_{j}$. Putting $A=0$ and $\lambda=1$ gives $C=1$ and the equation

[TABLE]

which was also derived in reference [5]. In this event the solution is k_{j}(t)=\hbox{\normalsize\frac{1+p}{1-p}}(t/t_{j})^{(1-p)/2}-\hbox{\normalsize\frac{2p}{1-p}} (assuming the initial condition $k_{j}(t_{j})=1$).

More generally, equation (13) can be cast in the general form

[TABLE]

where Q=p+\hbox{\normalsize\frac{(1-p)A}{2C}} and P=\hbox{\normalsize\frac{(1-p)\lambda}{2C}}, with solution

[TABLE]

using again the intial condition $k_{j}(t_{j})=1$.

The mean field degree distribution can be determined from this solution. The probability that node $j$ has degree $k_{j}(t)$ smaller than $\kappa$ at time $t$ is denoted by $P[k_{j}(t)<\kappa]$. Since $k_{j}(t)<\kappa$ if

[TABLE]

this is also the probability P\left[(t_{j}/t)>\left(\hbox{\normalsize\frac{Q/P+\kappa}{1+Q/P}}\right)^{-1/P}\right]. If the node $t_{j}$ is chosen uniformly from the $n$ available, then

[TABLE]

The mean field degree distribution is the derivative of this to $\kappa$:

[TABLE]

For large $\kappa$ this shows that the modified Barbasi-Albert network is scale-free with exponent

[TABLE]

Putting $A=0$ gives the exponent

[TABLE]

This is the mean field exponent of a modified Barabasi-Albert network. For small $\lambda<1$ the exponent is large, indicating a network with few nodes (if any) of high degree. For large $\lambda>1$, $\gamma\searrow 3^{+}$. This is a lower bound on $\gamma$ for modified Barabasi-Albert networks.

If $\lambda=1$, then the exponent $\gamma$ is given by

[TABLE]

In this model one similarly finds that $\gamma\geq 3$, and in fact, if $p=0$, then $\gamma=2+\sqrt{1+2A}$. The parameter $A$ may be used to tune the exponent $\gamma$ for any given $p$.

If both $\lambda=1$ and $A=0$, then the known expression for $\gamma$ for Barabasi-Albert networks is recovered, namely

[TABLE]

Notice that $\gamma\geq 3$ and that $\gamma=3$ if $p=0$ [5].

The connectivity of modified Barabasi-Albert networks is given by

[TABLE]

where $2C=((p+(1-p)\lambda)+\sqrt{(p+(1-p)\lambda)^{2}+2(1-p)A})$. Since 2-\gamma=1-\hbox{\normalsize\frac{1}{P}}, equation (5) gives \langle k\rangle_{n}\simeq\hbox{\normalsize\frac{1}{1-P}}. Inserting the value of $P$ gives the result above as well.

In the figure 4 the probability $P(k)$, that the degree of a Barabasi-Albert network is equal to $k$, is examined by plotting $\log P(k)/\log(k+1)$ against $1/\log(k+1)$ where $P(k)$ was estimated for values $n\in\{6250,12500,25000,50000,100000,200000\}$ and for $p=0$. The curves should intersect the vertical axis at $-\gamma$. Least squares fit of the data to quadratic curves gives $6$ estimates for $\gamma$, which average to $\gamma=3.026\pm 0.076$, very close to the theoretical value $\gamma=3$ from equation (20) (for $p=0$ and $\lambda=1$).

Data collected for the same values of $n$ and for $p=0.5$ cannot be successfully analysed by regressions with quadratic curves, but cubic curves give the average value $\gamma=5.161\pm 0.068$, which are not equal to but still fairly well approximated by $\gamma=5$ prediced by equation (20) for $p=0.5$ and $\gamma=1$.

When $p=0.8$ the plots are strongly curved and extrapolation to estimate $\gamma$ is more difficult. In this case a different approach is needed. Putting \alpha=\hbox{\normalsize\frac{1}{2}} in equation (4) gives

[TABLE]

so that a plot of \zeta(k)=(\log P(k)-\log P(\hbox{\normalsize\frac{1}{2}}k))/\log 2\to-\gamma as $k\to\infty$. That is, plotting $\zeta(k)$ against $\frac{1}{k}$ gives a curve with $y$-intercept equal to $-\gamma$. Better results are obtained when plotting against \hbox{\normalsize\frac{1}{k}}\log k. In this case a linear extrapolation gives $\gamma=11.67\pm 0.41$ and a quadratic extrapolation gives $\gamma=11.6\pm 2.6$. These results are close to the mean field prediction $\gamma=11$ for $p=0.8$. Incidently, if $p=0.5$ then this kind of analysis show that $\gamma=5.47\pm 0.14$ (linear extrapolation) or $\gamma=4.4\pm 1.0$ (quadratic extrapolation), and if $p=0$, then the results are $\gamma=3.088\pm 0.022$ (linear extrapolation) and $\gamma=2.86\pm 0.18$ (quadratic extrapolation).

If $\lambda=2$ and $p=A=0$ then the algorithm grows modified Barabasi-Albert networks with $\gamma=3$ (the mean field estimate given by equation (19)). Estimating $\gamma$ by plotting $\zeta(k)$ against \hbox{\normalsize\frac{1}{k}}\log k gives the estimate $\gamma=3.019\pm 0.098$ (linear extrapolation) and $\gamma=2.62\pm 0.33$ (quadratic extrapolation).

The connectivity of Modified Barabasi-Albert networks should converge quickly to a constant with increasing $n$ (by equation (5)) since $\gamma>2$. Computing it for Barabasi-Albert networks (with $\lambda=1$ and $A=0$) gives $\langle k\rangle_{n}\approx 3.16$ for $p=0$, $\langle k\rangle_{n}\approx 2.28$ for $p=0.5$ and $\langle k\rangle_{n}\approx 2.08$ for $p=0.8$, and for $n=12500$. Increasing $n$ does not change these results.

3.4 Mean field theory for Variant Barabasi-Albert networks

In this model the increment in the number of bonds when the $(n+1)$-th node is appended is given by

[TABLE]

Approximating this with a differential equation gives

[TABLE]

The right hand side can be approximated as follows: For $\alpha>1$ the algorithm should grow dense networks with nodes of high degree. Assuming that $k_{j}(n)\approx k_{\ell}(n)$ for all $\ell$ shows that \sum_{j}(k_{j}(n))^{\alpha}\approx n(k_{j}(n))^{\alpha}\approx n\left(\hbox{\normalsize\frac{1}{n}}\sum_{j}k_{j}(n)\right)^{\alpha}=n^{1-\alpha}(2E_{n})^{\alpha}. Using this approximation gives

[TABLE]

If $A=p=0$, then the differential equation can be solved directly to obtain $E_{n}\simeq 2^{(\alpha-1)/(2-\alpha)}n$, provided that $\alpha>1$. This shows that $E_{n}$ is linear in $n$, which may be expected if $\alpha$ is not too much larger than $1$.

Numerical experimentation shows that $E_{n}$ grows linearly in $n$ for values of $\alpha$ not too much larger than $1$. For example, if $p=0.5$, $A=1$ and $\alpha=1$ then \hbox{\normalsize\frac{1}{n}}E_{n}\to 1.207\ldots, if $\alpha=1.5$ then \hbox{\normalsize\frac{1}{n}}E_{n}\to 1.539\ldots, but if $\alpha=2$ then \hbox{\normalsize\frac{1}{n}}E_{n} increases slowly with $n$. Similarly, if $p=0$, and $A=1$, then, if $\alpha=1$, \hbox{\normalsize\frac{1}{n}}E_{n}\to 1.366\ldots, and if $\alpha=1.5$, \hbox{\normalsize\frac{1}{n}}E_{n}\to 2.399\ldots, but if $\alpha=2$ then \hbox{\normalsize\frac{1}{n}}E_{n} increases slowly with $n$ and for even larger values of $n$ this growth accelerates.

The recurrence for the degree of the $j$-th node may be approximated by a differential equation similar to equation (13): Assuming that $E_{n}=Dn^{\beta}$, replacing $n\to t$ (a continuous time variable), gives the recurrence

[TABLE]

This can be approximated by the differential equation

[TABLE]

If $\alpha=1$ and $\beta=1$ then the solution of this equation gives the Barabasi-Albert case with $\gamma=3$. Proceed by considering the case $A=p=0$ and the initial condition $k_{j}(t_{j})=1$. Assume that $\alpha=1+\epsilon$. Then the equation becomes

[TABLE]

A perturbative approach for small $\epsilon$ can be done by expanding (k_{j}(t))^{\epsilon}=\hbox{exp}(\epsilon\log k_{j}(t))=1+\epsilon\log k_{j}(t)+\hbox{\normalsize\frac{1}{2}}\epsilon^{2}\log^{2}k_{j}(t)+\cdots. Truncating this at $O(\epsilon^{2})$ and putting $g(t)=\log k_{j}(t)$ gives the differential equation

[TABLE]

Using the initial condition $g(t_{j})=\log k_{j}(t_{j})=0$ the solution of this equation is

[TABLE]

In the case $\beta>1$ suppose that $\delta=\beta-1$ and that $\delta$ is small. Then approximate

[TABLE]

With this approximation the solution for $g(t)$ above can be expanded in $\epsilon$ and $\delta$ to give the first order approximations

[TABLE]

Proceed by solving the above quadratics for \log(\hbox{\normalsize\frac{t}{t_{j}}}) in terms of $g(t)$. Expand the solution in $\epsilon$ and $\delta$ and keep only the first few terms. In the case that $\beta=1$ this gives

[TABLE]

Since $g(t)=\log k_{j}(t)$, the probability that $k_{j}(t)<\kappa$ is given by

[TABLE]

Taking the derivative to $\kappa$ gives the distribution function in the case that $\beta=1$:

[TABLE]

These networks are thus not scale-free. For small values of $k$ the $\log k$ terms are slowly varying, and the networks will appear to be scale-free with $\gamma=1+2D$. However, with increasing $k$ the exponent reduces in value and the connectivity of the network will become dependent on $k$ in the way seen in equation (5) for small values of $\gamma$.

Notice that if $D=1$ and $\epsilon=0$ (or $\alpha=1$), then the above reduces to $P(k)\sim k^{-3}$, as expected for Barabasi-Albert networks.

If $\beta>1$, then a similar approach to the above may be considered. Solving the expression for $g(t)$ above for \log(\hbox{\normalsize\frac{t}{t_{j}}}) and keeping only terms to $O(\epsilon)$ and $O(\delta)$ gives

[TABLE]

This shows that

[TABLE]

This shows that

[TABLE]

This gives an effective exponent $\gamma_{k}=1+2D+\delta\log t_{j}+D(2D\delta-\epsilon)\log k$ which decreases in size if $2D\delta-\epsilon<0$ and increases in size if $2D\delta-\epsilon>0$. Since $\delta=\beta-1$ and $\epsilon=\alpha-1$, and for small $\alpha$ numerical simulations show that $\beta\approx 1$, it is normally the case that $2D\delta-\epsilon<0$. This means that the networks will first appear scale-free with constant connectivity until $k$ becomes large enough in which case the connectivity will increase with $k$, as seen above.

3.5 Numerical results on Variant Barabasi-Albert networks

In figure 5 data for networks with $p=0$ and $\alpha=1.1$ and $\alpha=0.5$ is shown. Since $\alpha=1.1$ is still very close to $1$, the results above show that these networks should still appear scale-free, and with connectivity a constant. This is indeed the case. For $n=6250$ the data gives $\langle k\rangle_{n}=3.149$, and increasing $n$ to $n=200000$ gives $\langle k\rangle_{n}=3.176$. That is, the connectivity of the networks are insensitive to $n$ over this range. Least squares fits to the curves with quadratic polynomials in order to determinate the value of $\gamma$ give the average $\gamma=2.857\pm 0.068$. This result is consistent with a constant value of the connectivity of networks of these size ranges. With increasing $n$, it is expected that $\gamma$ will decrease in value (that is, the value given here is an effective value), and eventually, the connectivity will start to increase.

Networks generated with $p=0$ and $\alpha=0.5$ turned out to be sparse with low connectivity. For example, for $n=100000$, the connectivity is $\langle k\rangle_{n}=1.036$ and this decreases even further for $n=200000$, where $\langle k\rangle_{n}=1.020$. Attemps to extract an exponent $\gamma$ from the data for these networks were not succesful, the regressions did not settle on a value, but are strongly dependent on $n$. Notice that the mean field analysis above does not apply to networks with $\alpha<1$.

Putting $\alpha=2$ gives networks with average connectivity which increases with $n$. For example, if $n=100$, then $\langle k\rangle_{n}=43$, for $n=500$, $\langle k\rangle_{n}=260$ and for $n=1000$, $\langle k\rangle_{n}=527$. On the other hand, if \alpha=\hbox{\normalsize\frac{3}{2}}, then $\langle k\rangle_{n}=3.08$ if $n=100$, $\langle k\rangle_{n}=3.27$ if $n=500$, and $\langle k\rangle_{n}=3.31$ if $n=1000$, and it appears that for small values of $n$ the connectivity does not change quickly with increasing $n$.

4 Duplication-Divergence networks

Biological models of protein evolution are usually presented in terms of two processes, namely (1) a duplication event involving a gene sequence in DNA, and (2) a (random) mutation of duplicated genes which then drift from one another in genetic space [28, 8, 17]. The mutations of duplicated and mutated genes change the proteome and the network of protein interactions: If the protein is self-interacting, then the duplicated proteins interact, and the mutated genes code for proteins with altered interactions (some gained, others weakened or lost) with other proteins.

The Duplication-Divergence algorithm models these processes in order to grow a network, and was used in order to estimate the rates of duplication and mutation in the protein interaction networks [27]. There is a rich and large literature reporting on modeling protein interaction networks using models which include processes of duplication and divergence [13, 19, 25, 22].

Since proteomic networks appear to be scale-free [12, 20], it seems likely that duplication and divergence processes should grow scale-free networks and that this should also be seen in computer algorithms which grow networks using duplication and divergence elementary moves. Duplication can be implemented by selecting nodes and duplicating them, and their incident bonds, in a network. Divergence is implemented by altering the bonds incident on particular nodes, namely either by deleting, adding or moving bonds. In the Duplication-Divergence algorithm these moves are implemented by selecting nodes uniformly for duplication to progenitor-progeny pairs, and by deleting bonds incident to either the progenitor node or its progeny. Notice that since nodes of high degree have a larger probability of being adjacent to a node selected for duplication, these nodes have a larger probability of receiving new bonds in the duplication process – in this way there are events of preferential attachment in this algorithm [11, 23].

The basic elementary move of the Duplication-Divergence algorithm is illustrated in figure 6. The algorithm is implemented as follows.

Duplication-Divergence algorithm:

1. 1.

Initiate the network with one node $x_{0}$ and apply the following steps iteratively; 2. 2.

Duplication: Choose a node $\upsilon$ uniformly and duplicate by creating node $\upsilon^{\prime}$; 3. 3.

For all bonds ${\langle w{\sim}\upsilon\rangle}$ incident with $\upsilon$, add the bonds ${\langle w{\sim}\upsilon^{\prime}\rangle}$; 4. 4.

With probability $p$ add the bond ${\langle\upsilon{\sim}\upsilon^{\prime}\rangle}$; 5. 5.

Divergence: delete one bond of the pair $\{{\langle w{\sim}\upsilon\rangle},{\langle w{\sim}\upsilon^{\prime}\rangle}\}$ incident with $\upsilon$ or with its duplicated node $\upsilon^{\prime}$ with probability $q$ (for each $w$ adjacent to both $\upsilon$ and $\upsilon^{\prime}$ independently); 6. 6.

Stop the algorithm when a network of order $N$ is grown. ∎

The algorithm has two parameters $(p,q)$.

The parameter $p$ is the probability that the protein corresponding to the progenitor node $\upsilon$ is self-interacting. If it is (with probability $p$) then the bond ${\langle\upsilon{\sim}\upsilon^{\prime}\rangle}$ is added to the network and it represents the interaction between $\upsilon$ and $\upsilon^{\prime}$.

The parameter $q$ controls the model of divergence in this algorithm. As $\upsilon$ and $\upsilon^{\prime}$ diverge from one another, one bond in each pair of bonds incident with $\upsilon$ and $\upsilon^{\prime}$ is lost independently, with probability $q$. The result is that the network mutates as bonds (interactions) are lost (while they are created by the duplication process).

A slightly modified algorithm is found by changing step 5 in the algorithm to find a modification of the Duplication-Divergence algorithm which assumes that one of the duplicated pair mutates, while the other remains stable.

1. 1. 5.

Divergence: Consider all bonds ${\langle w{\sim}\upsilon^{\prime}\rangle}$ incident with the duplicated node $\upsilon^{\prime}$ and delete these independently with probability $q$.

The Duplication-Divergence algorithm tends to grow disconnected networks, while the Modifed Duplication-Divergence algorithm is more likely to grow networks with a single component (that it, connected networks).

4.1 Mean field theory for Duplication-Divergence networks

Let $k_{j}(n)$ be the degree of node $j$ after $n$ iterations. The algorithm appends nodes by duplicating them (the probability that a node $\upsilon$ is duplicated in a network of order $n$ is $\frac{1}{n}$), adds bonds by inserting a bond between a node and its duplicate with probability $p$, and remove bonds by selecting one bond between node-duplicate pairs and other nodes independently and deleting it with probability $q$. Let $2\>\!E_{n}=\sum_{j}k(n)$ be twice the total number of bonds after $n$ iterations. Then, if $k_{j}(n)$ is the degree of node $j$ at time $n$, and node $j$ is duplicated, the number of bonds in the network $E_{n}$ increases in the mean field by

[TABLE]

This follows since $k_{j}(n)$ bonds are created in the duplication move in the mean field, and another bond is created between the $j$-th node and its duplicate with probability $p$. The number of deleted bonds in the mean field is $q\,k_{j}(n)$.

Notice that $2\>\!E_{n}=\sum_{j}k_{j}(n)=n\>\!a_{n}$ where $a_{n}=\langle k_{j}(n)\rangle$ is the average degree. In the mean field approximation one substitutes $k_{j}(n)$ in the recurrence (38) by its network average $a_{n}$. Then equation (38) can be casted as a recurrance for $a_{n}$:

[TABLE]

Let $n\to t$, where $t$ is a continuous time variable, and approximate this recurrence by the differential equation

[TABLE]

The initial condition is $a_{1}=1$, and this has solution

[TABLE]

Since E_{n}\simeq\hbox{\normalsize\frac{1}{2}}n\,a_{n}, it follows that

[TABLE]

Comparison to equation (7) shows that, if q<\hbox{\normalsize\frac{1}{2}},

[TABLE]

In this case $E_{n}=O(n^{2(1-q)})+O(n)$ and that while $2(1-q)>1$, the term $O(n)$ is a strong correction to the growth in $E_{n}$ for even large values of $n$. In other words, the degree distribution $P(k)$ of the network will be strongly corrected from the powerlaw distribution in equation (1).

If q=\hbox{\normalsize\frac{1}{2}}, then by solving equation (40), $a_{t}=1+2p\log t$ (so that $a_{1}=1$). Since E_{n}=\hbox{\normalsize\frac{1}{2}}n\,a_{n}, this shows that

[TABLE]

In this case $\gamma=2$ by equation (7), but notice the subtle domination of the $n\log n$ term. In numerical work this will be very hard to see.

The case q>\hbox{\normalsize\frac{1}{2}} is considered by noting that a_{t}\simeq\hbox{\normalsize\frac{2p}{2q-1}} as $t\to\infty$. This shows that

[TABLE]

This shows that $\gamma\geq 2$ by equation (7).

Putting the above together gives

[TABLE]

with a logarithmic correction if q=\hbox{\normalsize\frac{1}{2}}.

Comparing the coefficient in equation (7) with equation (45) gives a refined estimate \gamma=1+\hbox{\normalsize\frac{2p}{1+2p-2q}}\geq 2, provided that $2q<1+2p$. For example, if $q=0.75$ then $p>0.25$. However, numerical work shows this estimate to be too small, and estimating $\gamma$ in this regime for this model remains an open question.

The power law decrease in $P(k)$ in equation (1) is only asymptotic for this algorithm; and there should be corrections in particular for q<\hbox{\normalsize\frac{1}{2}}. From the results above the average connectivity can be computed: Since E_{n}=\hbox{\normalsize\frac{1}{2}}n\,\langle k_{j}(n)\rangle,

[TABLE]

From these results $P(k)$ can be calculated. Since $\langle k\rangle_{n}\simeq\int_{1}^{n}k\,P(k)\,dk$, it follows that \hbox{\normalsize\frac{d}{dn}}\langle k\rangle_{n}=n\,P(n). Thus, using this approach gives

[TABLE]

where the case q>\hbox{\normalsize\frac{1}{2}} is unknown since the dependence of the exponent $\gamma$ on the parameters $(p,q)$ is not known. Notice the change in behaviour at the critical value q=\hbox{\normalsize\frac{1}{2}}; this was already observed numerically in reference [27].

The modified Duplication-Divergence algorithm has the same recurrence (41), and so the values for $\gamma$ and relations for $\langle k\rangle_{n}$ and $P(k)$ remain unchanged for this algorithm. Notice that this implementation preserves the degree of the selected node, and tends to give a duplicated node with lower degree (while the (unmodified) implementation tends to lower the degrees of both the selected and duplicated nodes). As a result, networks generated with the modified algorithm have, on average, more nodes of degree equal to one (and so appear more tree-like).

4.2 Numerical results on Duplication-Divergence networks

In figure 7 two networks grown with the Duplication-Divergence algorithm are shown. Both networks were grown with $p=1$ and have order $300$. The network on the left was grown with divergence parameter $q=0.4$, and that on the right, with the higher mutation rate $q=0.6$.

In figure 8 data for networks grown with $p=0.75$ and $q=0.4$ are shown. The curves on the right were obtained by plotting $(\log P(k))/\log(k+1)$ averaged over $100$ networks of sizes $\{3125,6250,12500,25000,50000,100000,200000\}$ against $1/\log(k+1)$. The mean field value of $\gamma$ is denoted by the bullet on the left-hand axis. These data show that convergence to this value is very slow – this indicates strong corrections to scaling arising in equation (42).

An alternative approach is to estimate $\gamma$ by plotting $\zeta(k)=(\log P(2k)-\log P(k))/\log 2$ as a function of $\log(k+1)/k$ (see equation (4) with $\alpha=2$). The results are also strongly curved data (left in figure 8), and while the results are not inconsistent with the mean field value $\gamma\approx 1.9$ in this model, however, it seems difficult to extrapolate these curves to a limiting value of $\gamma$.

If q=0.60>\hbox{\normalsize\frac{1}{2}} then the results in figure 9 are seen. The curves of $\zeta(k)=(\log P(2k)-\log P(k))/\log 2$ as a function of $\log(k+1)/k$ have straightened considerably, and each can be extrapolated by a quadratic least squares to obtain an estimate $\gamma_{n}$ for each value of $n=3125\times 2^{\ell}$ (for $\ell=0,1,2,\ldots,6$). This gives estimates $\{9.68,8.52,7.99,7.95,7.82,7.58,7.05\}$ which can be extrapolated by a least squares fit of $\gamma_{n}=\gamma+A/\log n$, giving the estimate $\gamma\approx 2.87$, which is slightly larger than the value predicted by the mean field formula \gamma=1+\hbox{\normalsize\frac{2p}{1+2p-2q}} (see the paragraph following equation (46)). This suggests that the approach to limiting behaviour in this model is quite slow, consistent with the remarks after equation (46) in the previous section.

The average connectivity $\langle k\rangle_{n}$ is expected to behave according to equation (47). In table 1 $\langle k\rangle_{n}$ is listed for $p=0.75$ and $q=0.40$, $q=0.50$ and $q=0.60$. If $q=0.4$, then equation (47) suggests that $\langle k\rangle_{n}\simeq 8.5\,n^{0.2}$. Computing $\langle k\rangle n\times n^{-0.2}$ from the data in table 1 gives $\{5.18,5.45,5.65,5.91,5.96,6.01,6.12\}$. Plotting these results against $1/\log n$ and then linearly exptrapolating as $n\to\infty$ gives $7.98$, close to the value of $8.5$ predicted in equation (47).

If $q=0.5$, then equation (47) suggests that $\langle k\rangle_{n}\simeq 1.5\log n$ since $p=0.75$. Dividing the results in table 1 by $\log n$ for each value of $n$ gives the results $\{1.42,1.44,1.44,1.42,1.43,1.46,1,45\}$. The average of this is close to the predicted value of $1.5$.

Finally, if $q=0.6$ then the data appear to approach a constant. Extrapolating these results using the model $A+B/\log(n)$ gives the estimated limiting value $8.72$. By equation (5) this indicates that $\gamma=2.13$, a value which is quite close to $2.15$, the value predicted by the formula \gamma=1+\hbox{\normalsize\frac{2p}{1+2p-2q}} in the paragraph following equation (46).

5 Solé evolutionary networks

The Solé model [24, 23] modifies Duplication-Divergence model by using duplication and network rewiring as the basic elementary moves. As before, the duplication of nodes is an implementation of gene duplication, and the network rewiring is based on the loss and gain of protein interactions in the bulk of the network [6]. Thus, the algorithm grows networks based on a model of gene duplication and the rewiring of protein interactions; both these processes drive the evolution of the interactome.

The elementary move of the algorithm is as follows: A node in the network is chosen uniformly and randomly, and duplicated to form a progenitor-progeny pair. The progeny will have the same interactions as the progenitor. This network is updated in the rewiring step which has two parts: Bonds incident with the progeny protein are deleted with probability $\delta$, and new bonds are added in the network between nodes (excluding the progenitor protein) are created with probability $\alpha$. This implementation differs in two ways from the Duplication-Divergence algorithm. In the Solé model there are no self-interacting nodes, and the formation of new bonds in the rewiring steps only occurs in the Solé model.

The basic iterative step of the Solé algorithm is shown in figure 10 and a Solé evolutionary network of order $N$ nodes is grown as follows:

Solé evolutionary algorithm:

1. 1.

Initiate the network with one node $x_{0}$ and apply the following steps iteratively; 2. 2.

Choose a node $\upsilon$ uniformly and duplicate it to a new node $\upsilon^{\prime}$; 3. 3.

For each bond ${\langle w{\sim}\upsilon\rangle}$ incident with the chosen node $\upsilon$, add the bond ${\langle w{\sim}\upsilon^{\prime}\rangle}$ incident with the duplicated node $\upsilon^{\prime}$; 4. 4.

Delete each bond ${\langle w{\sim}\upsilon^{\prime}\rangle}$ added in step 3 with probability $\delta$ independently; 5. 5.

For all nodes $u$ not adjacent to the chosen node $\upsilon$, create the bond ${\langle u{\sim}\upsilon^{\prime}\rangle}$ with probability $\alpha$; 6. 6.

Stop the algorithm when a network of order $N$ is grown. ∎

The algorithm has two parameters $(\delta,\alpha)$. If $\delta=0$ and $\alpha=1$ then the algorithm grows complete simple networks. More generally, if $\alpha>0$ then on average roughly $\alpha N$ bonds are added to a network of order $N$. This shows that the algorithm grows networks of size $O(N^{2})$ – that is, Solé networks are rich in bonds.

5.1 Mean field theory for Solé networks

Let $E_{n}$ be the total number of bonds in a Solé network after $n$ iterations of the algorithm, and let $\langle k\rangle_{n}$ be the connectivity of the network (that is, the average degree of nodes) after $n$ iterations (so that $2\,E_{n}=n\langle k\rangle_{n}$). In the mean field approximation the node in step 2 of the algorithm has degree $\langle k\rangle_{n}$ and this number of bonds is added in step 3, while, in a similar way, $\delta\langle k\rangle_{n}$ bonds are removed in step 4. In step 5 there are $n-\langle k\rangle_{n}$ choices in the mean field for the node $u$ not adjacent to $\upsilon^{\prime}$ and each bond ${\langle u{\sim}\upsilon^{\prime}\rangle}$ is added with probability $\alpha$. This shows that the number of bonds after $n+1$ iterations is given by the recurrance relation

[TABLE]

Since $2\,E_{n}=n\langle k\rangle_{n}$ this becomes

[TABLE]

which is a mean field recurrance relation for $E_{n}$.

Taking $n\to t$, a continuous time variable, and approximating $E_{n}$ by $E_{t}$, and approximating the finite difference as a derivative, gives the following differential equation for $E_{n}$:

[TABLE]

Solving this equation and letting $t\to n$ again gives the approximate mean field solution for $E_{n}$:

[TABLE]

Equation (52) shows that the number of bonds is proportional to $n^{2}$, so that networks created by this algorithm are dense, except when $\alpha=0$. Comparison to equation (7) suggests that $\gamma\leq 1$ in this model. Notice that there is no logarithmic factor in the denominator, and that $E_{n}=\Theta(n^{2})$. This is consistent with a mean field value $\gamma<1$ (and this requires that $P_{n}(k)$ be modified so that it is a normalisable probability distribution). With these results, it is reasonable to expect that, in the mean field,

[TABLE]

If $\alpha=0$ then equation (52) gives $E_{n}\sim n^{2-2\delta}$ and comparison to equation (7) gives

[TABLE]

5.2 Numerical results for Solé networks

Similar to Barabasi-Albert and Duplication-Divergence networks, Solé networks can be grown numerically by implementing the algorithm as given above, using sparse matrix routines to efficiently store the adjacency matrix of the network. The larger size of networks makes these more difficult to grow, and our algorithms sampled efficiently to networks of size $51200$ bonds.

Solé networks are rich in bonds. This is seen, for example, in equation (52), which shows that $E_{n}\propto n^{2}$ if $\alpha>0$. In figure 11 two examples of networks generated by the Solé algorithm are shown. If $\delta<0.5$, then the networks have a dense appearance dominated by a few hubs. If $\delta>0.5$, then the networks appear more extended, often with no nodes qualifying as hubs under the definition that the degree of a hub in a network of order $n$ is at least $\lfloor\sqrt{n}\rfloor$. The networks in figure 11 were generated with $\alpha=0.005$, and increasing the value of $\alpha$ quickly increases the number of bonds.

The mean field result that $\gamma\leq 1$ has implications for the scaling of Solé networks. In particular, $P_{N}(k)$ in equation (1) is not normalisable for infinite networks if $\gamma\leq 1$ and so is not a valid candidate degree distribution in this model. The degree distribution can be modified to

[TABLE]

where $D(x)$ is a function of the combined (or scaled) variable $x=n^{-\phi}k$. That is, as $n\to\infty$, $k$ is rescaled by $n^{-\phi}$ and $k^{\gamma}P(k)$ approaches a limiting distribution proportional to $D(x)$.

This can be tested numerically by plotting $n^{\gamma}P(k)$ as a function of $x=n^{-\phi}k$. For the proper choices of $\gamma$ and $\phi$ it is expected that $n^{\gamma}P(k)\simeq C_{o}D(x)$ for a wide range of values of $n$ (that is, the data should approach a limiting curve as $n\to\infty$). The result is shown in figure 12 for $(\delta=0.25,\alpha=0.005)$ and $(\delta=0.75,\alpha=0.005)$. These are plots on the same graph for $n=100\times 2^{n}$ for $n\in\{6,7,8,9\}$ (other curves at smaller values of $N$ are left away to give a clearer picture).

The data for $\delta=0.75$ are the cluster of peaks to the left, rescaled by choosing $\phi=1$ and \gamma=\hbox{\normalsize\frac{1}{2}}, while the cluster of peaks to the right is for $\delta=0.25$ with $\phi=1$ and \gamma=\hbox{\normalsize\frac{2}{3}}. With increasing $n$ the data appear to approach a single underlying curve if \gamma=\hbox{\normalsize\frac{1}{2}} in the one instance, and \gamma=\hbox{\normalsize\frac{2}{3}} in the other instance. Both these values are consistent with the mean field expectation that $\gamma\leq 1$ in this model. Further refinements in this scaling assumption may be necessary, since the curves are still becoming narrower with increasing $n$. It is not clear that these approach a limiting curve as $n\to\infty$, although the data for $\delta=0.75$ suggest this to be the case. In these cases the curves are sharply peaked with a mean of about $0.02$ if $\delta=0.25$ and about $0.007$ if $\delta=0.75$.

Since the curve $D(x)$ is sharply peaked at a constant value $c_{o}$ of the rescaled variable $x$, the connectivity of Solé networks is estimated by treating $D(x)$ as concentrated at $c_{o}$ and then (assuming that $\phi=1$ and approximating the connectivity)

[TABLE]

In other words, the connectivity of Solé networks should increase linearly with $n^{\phi}$ (and since $\phi=1$, linearly with $n$). In table 2 the connectivities of Solé networks for $\delta=0.25$ and $\delta=0.75$ (with $\alpha=0.005$) are listed. Non-linear least squares fits to the data show that $\phi=1.01$ when $\delta=0.25$ and $\phi=0.99$ when $\delta=0.75$. That is, these results are consistent with the value $\phi=1$ seen above.

6 The iSite model of network evolution

Protein interaction networks evolve by mutations in proteins which change the interactions of the proteins in the network. In the Duplication-Divergence algorithm, a mutated protein loses its interactions randomly. This random deletion of interactions is a good first order approximation to the evolution of networks. The iSite model refines this by giving structure to nodes in the network by introducing iSites on nodes as localities of the interaction sites on a protein [14, 15]. Subfunctionalization of interaction sites in the iSite model is implemented by silencing iSites, and adding interactions with reduced probability if the iSite is not silenced.

The implementation of the iSite algorithm relies in the first place on duplication of nodes, and then subfunctionalization of iSites on the nodes. The subfunctionalization of iSites is implemented by randomly deleting of bonds incident to duplicated iSites, and by the silencing of iSites by turning them off. These processes are models of random mutations which cause the loss of information in the genome (and leave behind non-coding remnants of genes). A process of spontaneously creating new iSites is not in the iSites algorithm, although this is a possible refinement which may be introduced in the algorithm.

The elementary move of the iSite algorithm is illustrated schematically in figure 13. A uniformly chosen node is duplicated into a progenitor-progeny pair (and so also duplicating the iSites of the progenitor onto the progeny). If the duplicated iSite is self-interacting, then bonds are added between the iSite on the progenitor and the duplicated iSite on the progeny with probability $p$ – this allows for subfunctionalization of the duplicate iSites. Bonds incident with the iSites on the progenitor are duplicated with reduced probability $r$, and iSites on the progenitor or progeny nodes are silenced with probability $q$. If an iSite is silenced, then all bonds incident with it are deleted. Notice that subfunctionalization enters in several ways, both in the duplication of self-interacting iSites, in the duplication of bonds, and in the silencing of iSites.

iSite evolutionary algorithm:

1. 1.

Initiate the network with one node $x_{0}$ with $I$ active iSites (each of which is self-interacting with probability $p$) and iterate the following steps; 2. 2.

Choose a progenitor protein $\upsilon$ uniformly in the network and duplicate it, and its associated iSites $A$, to a successor protein $\upsilon^{\prime}$ with duplicated iSites $A^{\prime}$;

• (a)

A duplicated iSite $A^{\prime}\in\upsilon^{\prime}$ is active with probability $q$ if it is duplicated from an active iSite on $A\in\upsilon$, and silenced otherwise;

• (b)

An active duplicated iSite $A^{\prime}\in\upsilon^{\prime}$ is self-interacting with probability $p$ if it is duplicated from a self-interacting iSite on $A\in\upsilon$, and not self-interacting otherwise;

• (c)

If a silenced iSite $A$ is duplicated to iSite $A^{\prime}$, then $A^{\prime}$ is also silenced; 3. 3.

Add bonds as follows:

• (a)

If iSite $A\in\upsilon$ is self-interacting and $A$ is duplicated to iSite $A^{\prime}\in\upsilon^{\prime}$, then add the bond ${\langle A{\sim}A^{\prime}\rangle}$ if $A^{\prime}$ is not silenced;

• (b)

If ${\langle A{\sim}B\rangle}$ is a bond incident with iSite $A$ on the progenitor $\upsilon$, and $A$ is duplicated to iSite $A^{\prime}$ on the duplicate $\upsilon^{\prime}$, then ${\langle A{\sim}B\rangle}$ is duplicated to ${\langle A^{\prime}{\sim}B\rangle}$ with probability $1-r$ provided that $A^{\prime}$ is not silenced; 4. 4.

Iterate the algorithm from step (2) and stop the iterations when a network of order $N$ is grown. ∎

6.1 Mean field theory for the iSite model

Let nodes in the network correspond to proteins, and let $i_{j}(n)$ be the number of active iSites on node $j$ after $n$ iterations of the algorithm. Denote the degree of node $j$ by $k_{j}(n)$ (that is the total number of bonds with one end-point in node $j$), and let $E_{n}$ be the number of bonds of the network (this is the size of the network). Then $2\>\!E_{n}=\sum_{j}k_{j}(n)$.

The average number of active iSites per node is i(n)=\hbox{\normalsize\frac{1}{n}}\sum_{j}i_{j}(n). With each iteration $i(n)$ iSites are created, of which $q\>\!i(n)$ are silenced, in the mean field. This gives the following recurrance relation for $i(n)$:

[TABLE]

The exact solution of this recurrance is

[TABLE]

where $\Gamma$ is the gamma function with the property that $\Gamma(x+1)=x\,\Gamma(x)$ and $\Gamma(1)=1$. Notice that $i(0)=I$, where $I$ is the number of iSites on the source node $x_{0}$.

For large $n$ the $\Gamma$-function and the factorial have well known asymptotics (namely the Stirling approximation [29]), so that

[TABLE]

This shows that with increasing $n$ the total number of iSites grows proportionally to $n^{1-q}$. If $q=0$, then this is linear in $n$ since no iSites become silenced, and if $q=1$, then the number approaches a constant.

The total number of bonds in the network increases after $n$ iterations by the recurrance

[TABLE]

since there are on average \hbox{\normalsize\frac{2}{n}}E_{n} bonds incident to each node, and the probability that each one of them is duplicated is $1-r$, and there are on average $i(n)$ iSites per node, and the probability that each of these is self-interacting is $p$.

Using the asymptotic solution for $i(n)$ and approximating this recurrence by a differential equation gives

[TABLE]

This equation can be solved, and using the initial condition $E_{1}=0$, the result is

[TABLE]

Thus, the average degree of a node is equal to \hbox{\normalsize\frac{2}{n}}E_{n}, so that the connectivity of iSite evolutionary networks is given by

[TABLE]

in the mean field. This shows that the large $n$ value of $\langle k\rangle_{n}$ is dominated by the larger of $-q$ and $1-2r$. In particular, if r<\hbox{\normalsize\frac{1}{2}}(1+q), then $\langle k\rangle_{n}\sim n^{1-2r}$.

By equation (7) one may determine $\gamma$ for this model:

[TABLE]

If $2r=1+q$, then a different solution is obtained, namely

[TABLE]

This shows that $\gamma=2+q$ in this case as well, but there is also a logarithmic correction to the growth of $E(t)$, and so there is a logarithmic factor in the expression for $\langle k\rangle_{n}$.

6.2 Modified iSite evolutionary algorithm

The subfunctionalization of proteins can be refined by introducing in the iSite algorithm the probability of creating new iSites on the progeny node with a probability $s$. This changes the algorithm as follows.

Modified iSite evolutionary algorithm: Implement the algorithm as above but introduce the parameter $s$ and create new active iSites by replacing step 2 in the iSite evolutionary algorithm by

• 2.

Choose a progenitor node $\upsilon$ uniformly in the network and duplicate it, and its associated iSites $A$, to a progeny node $\upsilon^{\prime}$ with duplicated iSites $A^{\prime}$;

• (a)

A duplicated iSite $A^{\prime}\in\upsilon^{\prime}$ is active with probability $q$ if it is duplicated from an active iSite on $A\in\upsilon$, and silenced otherwise;

• (b)

An active duplicated iSite $A^{\prime}\in\upsilon^{\prime}$ is self-interacting with probability $p$ if it is duplicated from a self-interacting iSite on $A\in\upsilon$, and not self-interacting otherwise;

• (c)

If a silenced iSite $A$ is duplicated to iSite $A^{\prime}$, then $A^{\prime}$ is also silenced;

• (d)

With probability $s$ create an active iSite $C$ on the progeny node $\upsilon^{\prime}$, where $C$ is self-interacting with probability $p$.

The recurrance for the average number of active iSites per node $i(n)$ (see equation (58)) is modified to

[TABLE]

in the Modified iSite evolutionary algorithm. The exact solution is obtained by replacing $q$ by $q-s$ in equation (58), and the asymptotic approximation of the solution is given by

[TABLE]

as seen in equation (59).

The total number of bonds in the network, $E_{n}$, still satisfies equation (60), and so it follows from equations (62), (63) and (64), that for the modified iSite evolutionary algorithm (notice the condition that $q<r+s$):

[TABLE]

This shows that the connectivity of Modified iSite networks is given by

[TABLE]

The value of the scaling exponent is seen from above to be given by

[TABLE]

with a correction factor in the expression for $\langle k\rangle_{n}$ if $2r=(1+q-s)$.

6.3 Numerical results for iSite networks

The iSite algorithm was coded and networks were grown to compute averaged statistics. Examples of iSite networks generated by the algorithm are shown in figure 14. The algorithm was then used to sample networks of size up to $200,000$.

The connectivity $\langle k\rangle_{n}$ of iSite networks for $I=3$ iSites per node, and with $p=0.5$, $q=0.4$ and $r=0.3$, is shown in table 3. By equation (5), \log\langle k\rangle\simeq\log\hbox{\normalsize\frac{\gamma-1}{2-\gamma}}+(2-\gamma)\log n. Least squares fit to the data in Column 2 gives \log\hbox{\normalsize\frac{\gamma-1}{2-\gamma}}\approx 1.0211, and $(2-\gamma)=0.258$. Solving for $\gamma$ gives in the first instance $\gamma=1.735$ and in the second $\gamma=1.742$. Since $2r<1+q$ in this case, the mean field value of $\gamma$ is $\gamma=1+2r=1.6$, close to these estimated values.

Data for $I=5$ and with the same values of $(p,q,r)=(0.5,0.4,0.3)$ are shown in table 3 as well. Changing the value of $I$ (the number of iSites per node) should not change the value of $\gamma$, and this appears to be the case here. A least squares fit to the data in Column 3 and determining $\gamma$ as above gives $\gamma=1.737$ and $\gamma=0.7498$, very close to the values above.

If $p=0.5$, $q=0.05$ and $r=0.8$, then $2r>1+q$, and in this case $\gamma=2+q$. If the number of iSites per node is $I=3$, then the data in table 3 gives a constant value for $\langle k\rangle$, and for $I=5$ a slightly decreasing numerical estimate. The mean field value of $\gamma$ in these cases is $2.05$, and a least squares fit gives $\gamma\approx 2.009$ if $I=3$ and $\gamma\approx 2.022$ if $I=5$ (where the coefficient of $\log n$ in the least squares fit is $2-\gamma$). These results are consistent with the mean field results obtained above, since it shows that the value of $\gamma$ is close to $2+q$.

7 Conclusions

In this paper a number of algorithms used for generating networks in molecular biology were examined. Mean field theory for the algorithms was in some cases reviewed, and in other cases newly presented, and also refined. The algorithms include the Barabasi-Albert [2], Duplication-Divergence [26], Solé [24] and iSite algorithms [14, 15], and these were in some cases modified by the introduction of more general elementary moves.

The efficient implementation of these algorithms was also examined, and sparse matrix routines (or, more general, hash-coding; see for example reference [21]) were used to optimize the implementation. This gives computer algorithms which can generate very large networks efficiently, and networks of order $200,000$ nodes were routinely sampled. We also explored even larger networks, up to order 3 million, but did not use those in our data analysis.

The adjacency matrix of a network of size $E$ bonds can be stored (using sparse matrix routines) in an array of size $O(E)$. This means that the implementation of these network growth algorithms has average case space complexity $O(E)$.

Hash coding allows for the efficient implementation of routines which search, insert or delete entries in arrays storing the networks. These routines have average time complexity $O(1)$ [10], (and worst case time complexity $O(E)$ for searches, inserting and deleting bonds, due to collisions if a hash table is densely populated).

Generally, the time complexity of algorithms should grow as $O(E^{\tau})$ if networks of size $E$ are grown (where $\tau$ is an exponent dependent on the particular algorithm). For example, networks of size $E$ bonds can be generated using $O(E)$ computer memory, and the Duplication-Divergence and iSite algorithms can be implemented with $O(n^{\tau})$ time complexity to grow networks of order $n$ nodes (and where $n\leq E$). An examination of these algorithms (the Duplication-Divergence and iSite algorithms) suggests that an optimal implementation will have $\tau\approx 1$ (if the size of the hash tables is much larger than $n$).

The Barabasi-Albert and Solé algorithms (with their modified and variant implementations) should have average time complexity of $O(n^{2})$ for growing networks of order $n$ nodes. This follows because each iteration of the algorithms has to explore all nodes in the current network for the possible insertion of new bonds.

Data on the time complexity of the algorithms are shown in table 4. The data displayed are the average time $T$ to grow one network of order $n$. Assuming that $T=C_{0}n^{\tau}$ and fitting $\log T$ to $\log n$, least squares estimates of $\tau$ can be obtained. For example, it is expected that $\tau=2$ for the Barabasi-Albert algorithm, while the estimate obtained in the table is $\tau\approx 1.97$. This is consistent with the expectation that the time complexity of the algorithm is $O(n^{2})$ in an optimal implementation. This is similarly seen for the modified and variant implementation of the Barabasi-Albert algorithm, and for the Solé algorithm.

The time complexity of the remaining algorithms is $O(n)$, and this is found consistently, except for the Duplication-Divergence algorithm for $q=1$ and $q=0.4$ (and also for the modified implementation of this algorithm). In these cases the algorithm samples denser networks (see figure 7) which takes up larger amounts of memory, making the implementation less efficient.

The results in this paper raise some questions about the sampling of scale-free networks by random iterative growth algorithms:

• •

In some cases, see for example reference [27], the parameters of the algorithms were set to grow networks with properties similar to that of real protein interaction networks. The values of the parameters are then used to estimate the rate of subfunctionalization (or mutation) in the genome. The results are dependent on the algorithm, and so further refinement of algorithms may be needed before useful estimates can be made.

• •

The mean field approaches are useful in some models (for example the Barabasi-Albert algorithm, and the iSite algorithm), but are poorer approximations in other models (the variant Barabasi-Albert algorithm, the Duplication-Divergence algorithm and its modification, and the Solé algorithm). Can the mean field approach be improved to give a better approximation to these algorithms?

• •

Investigation of some numerical properties of the networks (for example the connectivity) suggests that the algorithms may be self-averaging. That is, networks are generated with properties which converge to the statistical averages of these properties over a sample of networks generated by the algorithm. This is, for example, illustrated in figure 16 for the connectivity of Barabasi-Albert networks. As the network is grown, its connectivity appears to approach the average connectivity over a large sample of networks.

• •

In this paper some algorithms were modified in ways not done before in the literature (this includes the modified Barabasi-Albert, the Duplication-Divergence, the Solé and iSite models). Exploring the properties of these modified algorithms, including their usefulness as models of networks in molecular biology, will be the subject of future investigation.

Lastly, these algorithms grow networks using a probabilistic set of rules to implement an elementary move. Each realised network $N_{n}$ of order $n$ is obtained with some probability $p(N_{n})$, so that the function $p(N_{n})$ is a probability distribution over networks of order $n$. Determining $p(N_{n})$ for any of the algorithms presented here seems difficult, and general properties of $p(N_{n})$ remain unknown (other than averages of network properties over $p(N_{n})$ are scale-free if the algorithm grows scale-free networks).

Acknowledgements

EJJvR is grateful to NSERC (Canada) for support in the form of a Discovery Grant