A Markov Chain Approach to Randomly Grown Graphs

Abstract

A Markov chain approach to the study of randomly grown graphs is proposed and applied to some popular models that have found use in biology and elsewhere. For most randomly grown graphs used in biology, it is not known whether the graph or properties of the graph converge (in some sense) as the number of vertices becomes large. Particularly, we study the behaviour of the degree sequence, that is, the number of vertices with degree 0,1,…, in large graphs, and apply our results to the partial duplication model. We further illustrate the results by application to real data.

1. Introduction

Over the past decade, networks have played a prominent role in many different disciplines including theoretical physics, technology, biology, and sociology [1–9]. Particularly in biology, networks have become fundamental for the description of complex data structures. The appeal of networks may, at least partly, be due to the fact that in addition to being based on a rigorous mathematical base [10–14], they also provide a convenient graphical representation of the data which allows for visual interpretation. Examples of complex data structures that can be described by networks include food webs in ecology, sexual partner networks in sociology, and protein interaction networks in biology.

The canonical model in random graph theory has been Erdös-Renyi random graphs, where each of a fixed number of vertices has a Poisson-distributed number of links to other vertices. A Poisson number of links have turned out barely to be realistic for many empirically observed networks, and other models have been suggested to accomodate the discrepancies between theory and observation. Barabási and Albert [2] proposed a simple stochastic model, the preferential attachment (PA) model, whereby the network gradually is built up by adding one vertex at a time until the network reaches the desired size. This model is able to account for the scale-free degree distribution that is observed in some empirical networks, but not many of the other features and motifs that are found in real networks (e.g., [15–18]). Therefore, for mathematical and statistical analysis of network data, many other stochastic models have been proposed, in particular models that fall in the class of randomly grown graphs (RGGs; see next section for a definition) which share the property of the PA model of gradual growth. Overviews of different models and their properties can be found in [13, 16, 19, 20].

While the PA model has been under close mathematical scrutiny (e.g., [20]), other RGGs have been treated less extensively (e.g., [19, 21]) and mostly in the context of considering a continuous time approximation to the original discrete time Markov process (e.g., [13, 22, 23]). In this paper, we specifically address the question of the behavior of the vertex degrees as the number of vertices grows large. For a class of RGGs (including the PA model), the existence of a limiting degree distribution has been proven and its analytical form has been derived [21]. However, for most RGGs applied in biology, it is not known whether a limiting distribution exists, letting alone its form.

Biologically, it is of great interest to know whether the network stabilizes as it grows, or whether the degree distribution is a function of the size of the network, even for large network sizes. It relates to the question whether the network in an evolutionary perspective reaches an equilibrium, such that adding new vertices does not change the overall connectivity of the network. For example, in relation to protein interaction networks where vertices represent proteins and edges represent physical interactions between proteins, both scenarios seem a priori possible. Proteins may be able to engage in an unlimited number of interactions, or the number of interactions may be restricted by a number of factors such as space, time, and protein production rates. With the increasing statistical interest in analyzing complex biological networks with respect to evolutionary and functional properties [1, 5, 9, 13, 14, 24], it is also becoming of interest to understand the mathematical properties of the models.

We study a large class of RGGs that allows the construction of a simple, but time-inhomogeneous, Markov chain. For a given RGG, the corresponding Markov chain can be used to address questions about the RGG, for example, questions about the degree distribution. In particular, we focus on a special RGG, the partial duplication model, which has recently been used in the study of biological protein interaction networks [16, 18, 25, 26] and has formed the basis for new and more biologically realistic models (e.g., [16, 27]). The partial duplication model has two parameters (𝑝 and 𝑞) and we give conditions under which the chain is ergodic or transient. Further, based on the time-inhomogeneous Markov chain, we define a time-homogeneous Markov chain and a continuous time, time-homogeneous Markov process, and demonstrate that these, in general, are easier to study and apply than the original chain. Proofs rely on general theory of discrete Markov processes, which makes it easy to prove similar results for other RGGs.

Finally, we apply our results to a collection of real protein interaction data.

2. RGGs

An RGG is a Markov chain {𝐺_𝑡}_𝑡≥𝑠 on undirected graphs such that 𝐺_𝑡 has exactly 𝑡 vertices, and the set of vertices of 𝐺_𝑡 is a subset of the set of vertices of 𝐺_𝑡+1 for all 𝑡≥𝑠. Note that we do not require the set of edges of 𝐺_𝑡 to be a subset of the set of edges of 𝐺_𝑡+1.

Denote by 𝑛_𝑡(𝑘) the expected number of vertices of degree 𝑘 at time 𝑡. Since, by assumption, the graph 𝐺_𝑡 has exactly 𝑡 vertices, the expected relative frequency of vertices of degree 𝑘 at time 𝑡 is 𝑓_𝑡(𝑘)=𝑛_𝑡(𝑘)/𝑡. For many RGGs, one can derive a recursive formula for 𝑛_𝑡(𝑘), often referred to as the master equation [13]. Here, we consider a very general class of master equations given by𝑛_𝑡+1(𝑘)=_𝑗≥0𝐴(𝑡)_𝑘,𝑗𝑛_𝑡(𝑗), (2.1)where 𝐴(𝑡) for all 𝑡≥𝑠 is an infinite real matrix with 𝐴(𝑡)_𝑘,𝑗=0 for 𝑘>𝑗+1, and such that all columns sum to the same number 𝑎(𝑡). Furthermore, assume for suitable real numbers 𝑏_𝑘,𝑗 that𝐴(𝑡)_𝑘,𝑗=𝑏1−^𝑘,𝑘𝑡𝑏for𝑘=𝑗,_𝑘,𝑗𝑡for𝑘≠𝑗 (2.2)with 𝑏_𝑘,𝑗=0 for 𝑘>𝑗+1. The latter condition guarantees that the vertex degree can increase by at most one. By construction, 𝑛_𝑡(𝑘)=0 for 𝑘+1>𝑡, and hence 𝐴(𝑡) is effectively a 𝑡×(𝑡−1) matrix. We assume that the entries (2.2) in this submatrix are positive.

One particular example of a model fulfilling the conditions above is the partial duplication model (details are found in Section 4). The master equation is given by𝑛_𝑡+1(𝑘)=(1−𝑞+𝑘𝑝𝑡)𝑛_𝑡(𝑘)+(1−𝑞)_𝑗≥𝑘(𝑗𝑘)𝑝^𝑘(1−𝑝)^𝑗−𝑘𝑛_𝑡(𝑗)𝑡+𝑞+(𝑘−1)𝑝𝑡𝑛_𝑡(𝑘−1)+𝑞_{𝑗≥𝑘−1}(𝑗𝑘−1)𝑝^𝑘−1(1−𝑝)^{𝑗−𝑘+1}𝑛_𝑡(𝑗)𝑡. (2.3)For several other models, the master equation takes a similar form. Among these models are the duplication-divergence model [16], an approximation to the duplication-mutation model [22, 23], and the models discussed in [21] after a suitable modification (see Section 5.2). Generally, (2.1) is fulfilled whenever the expected degree change in a vertex depends on the degree only, and not on the degrees of the other vertices.

It follows immediately from (2.1) that𝑓_𝑡+1(𝑘)=_𝑗≥0𝐵(𝑡)_𝑗,𝑘𝑓_𝑡(𝑗), (2.4)where 𝐵(𝑡) is the transpose of (𝑡/(𝑡+1))𝐴(𝑡), and by assumption all rows of 𝐵(𝑡) sum to 𝑏(𝑡)=(𝑡/(𝑡+1))𝑎(𝑡). It follows that1=_𝑘≥0𝑓_𝑡+1(𝑘)=_{𝑘≥0𝑗≥0}𝐵(𝑡)_𝑗,𝑘𝑓_𝑡(𝑗)=𝑏(𝑡)_𝑗≥0𝑓_𝑡(𝑗)=𝑏(𝑡), (2.5)that is, 𝑏(𝑡)=1 and the matrices {𝐵(𝑡)}_𝑡≥𝑠 describe a Markov chain with time-dependent transition probabilities.

Proposition 2.1. Assume that ∑_𝑘≥0𝑏_𝑗,𝑘<∞ for all 𝑗≥0. If 𝑓_𝑡(𝑗)→𝑓(𝑗) pointwise for all 𝑗≥0, then {𝑓(𝑗)}_𝑗≥0 satisfies0=_𝑘≠𝑗𝑏_𝑗,𝑘𝑓(𝑘)−(1+𝑏_𝑗,𝑗)𝑓(𝑗),_𝑗≥0𝑓(𝑗)≤1. (2.6)

Proof. The second part of the proposition is a simple application of Fatou's lemma. By using (2.4), the definition of 𝐵(𝑡), and ∑_𝑘≥0𝑏_𝑗,𝑘<∞, it follows that(𝑡+1)(𝑓_𝑡+1(𝑗)−𝑓_𝑡(𝑗))=_𝑘≠𝑗𝑏_𝑗,𝑘𝑓_𝑡(𝑘)−(1+𝑏_𝑗,𝑗)𝑓_𝑡(𝑗)⟶𝑑_𝑗for𝑡⟶∞ (2.7)for some real number 𝑑_𝑗, and it remains to prove that 𝑑_𝑗=0. Note that^𝑡_𝑛=𝑠𝑑_𝑗(𝑛)=(𝑡+1)𝑓_𝑡+1(𝑗)−𝑠𝑓_𝑠(𝑗)−^𝑡_𝑛=𝑠𝑓_𝑛(𝑗), (2.8)and by using Cesaro's lemma, we get𝑑_𝑗=lim_𝑡→∞1𝑡^𝑡_𝑛=𝑠𝑑_𝑗(𝑛)=𝑓(𝑗)−𝑓(𝑗)=0. (2.9)

Consider the jump chain corresponding to the Markov chain {𝐵(𝑡)}_𝑡≥𝑠, that is, the Markov chain with transition probabilities 𝐵(𝑡)_𝑗,𝑘/(1−𝐵(𝑡)_𝑗,𝑗) for 𝑗≠𝑘, unless 𝐵(𝑡)_𝑗,𝑗=1 in which case the probability is put to 0. The jump chain has time-independent transition probabilities given by𝑝_𝑗,𝑘=𝑏_𝑘,𝑗1+𝑏_𝑗,𝑗for𝑘≠𝑗, (2.10)and 𝑝_𝑗,𝑗=0 for all 𝑗≥0. If 1+𝑏_𝑗,𝑗=0, then 𝑝_𝑗,𝑘=0. Occasionally, we consider a slightly modified jump chain (still with time-independent transition probabilities) which is allowed to stay in the same state with positive probability.

If a stationary distribution {𝜋_𝑗}_𝑗≥0 for the jump chain exists, it fulfills𝜋_𝑗=_𝑘≠𝑗𝜋_𝑘𝑏_𝑗,𝑘1+𝑏_𝑘,𝑘∀𝑗≥0. (2.11)Assume that inf_𝑗≥0(1+𝑏_𝑗,𝑗)>0 and put 𝜋^′_𝑗=𝜋_𝑗/(1+𝑏_𝑗,𝑗). Then we obtain that0=_𝑘≠𝑗𝑏_𝑗,𝑘𝜋^′_𝑘−(1+𝑏_𝑗,𝑗)𝜋^′_𝑗∀𝑗≥0, (2.12)and hence {𝜋^′_𝑗}_𝑗≥0 is a solution to the equation in Proposition 2.1. Furthermore, we may normalize {𝜋^′_𝑗}_𝑗≥0 to get a distribution, and hence (2.11) and (2.12) may be used to transfer a stationary distribution for the jump chain to the limit of the time-inhomogeneous Markov chain and vice versa.

In our main example, the partial duplication model (see Section 4 for details), we have 𝑏_0,0=2𝑞−1 and𝑏_𝑗,𝑗=𝑞+𝑗𝑝−(1−𝑞)𝑝^𝑗−𝑞𝑗𝑝^𝑗−1(1−𝑝)for𝑗≥1, (2.13)and hence the assumption inf_𝑗≥0(1+𝑏_𝑗,𝑗)>0 is fulfilled if 𝑞>0.

3. A Continuous Time Approximation

In this section, we show that the time-inhomogeneous Markov chain converges to a continuous time, time-homogeneous Markov process after a suitable time transformation.

Denote by 𝑇_𝑖 the time of the 𝑖th jump in the time-inhomogeneous chain after a given time 𝑡₀, and let 𝐽_𝑖 be the state to which it jumps. Set 𝑇₀=𝑡₀ and 𝐽₀=𝑗₀, the state of the chain at time 𝑡₀. To simplify notation further, introduce 𝑠_𝑖=(𝑡_𝑖,𝑗_𝑖) and 𝑆_𝑖=(𝑇_𝑖,𝐽_𝑖).

Note that at time 𝑡, the probability of staying in state 𝑗 is 𝐵(𝑡)_𝑗,𝑗=1−(𝑏_𝑗,𝑗+1)/(𝑡+1). In particular, if we let 𝛼_𝑖=𝑏_{𝑗𝑖−1,𝑗𝑖−1}+1, then𝑃(𝑇_𝑖>𝑡_𝑖∣𝑆_𝑖−1=𝑠_𝑖−1)=^𝑡𝑖_{𝑢=𝑡𝑖−1+1}𝛼(1−^𝑖𝑡𝑢+1)≈(^𝑖𝑡_𝑖−1)^−𝛼𝑖 (3.1)for large 𝑡_𝑖−1 and 𝑡_𝑖. Now consider the transformation 𝑍_𝑖=log𝑇_𝑖−log𝑇_𝑖−1=log𝑇_𝑖/𝑇_𝑖−1. It follows that𝑃(𝑍_𝑖>𝑧∣𝑆_𝑖−1=𝑠_𝑖−1)=𝑃(𝑇_𝑖>𝑡_𝑖−1𝑒^𝑧∣𝑆_𝑖−1=𝑠_𝑖−1)⟶𝑒^{−𝛼𝑖𝑧} (3.2)as 𝑡_𝑖−1→∞. That is, in the limit, the transformed waiting time is exponentially distributed with parameter 𝛼_𝑖.

Proposition 3.1. Let 𝑋_𝑠0(𝑧), 𝑧≥0, take the value of the time-inhomogeneous Markov chain at time 𝑡, where 𝑡=⌊𝑡₀𝑒^𝑧⌋ and ⌊𝑥⌋ denotes the integer part of 𝑥. At time 0, 𝑋_𝑠0(0)=𝑗₀. For fixed 𝑗₀, the process 𝑋_𝑠0(𝑧) converges to a continuous time, time-homogeneous Markov process as 𝑡₀→∞.

Proof. Clearly, the process 𝑋_𝑠0(𝑧), 𝑧≥0, is Markovian by definition. Let 𝑍^⋆_𝑖 be the time of the 𝑖th jump, that is, 𝑇_𝑖=⌊𝑡₀𝑒^𝑍⋆𝑖⌋=𝑡₀𝑒^𝑍⋆𝑖 and 𝑍_𝑖=𝑍^⋆_𝑖−𝑍^⋆_𝑖−1 in the notation above. It follows from (3.2) that𝑃_𝑠0(𝑍^⋆_𝑖>𝑧∣𝑍^⋆_𝑖−1=𝑧^⋆_𝑖−1,𝐽_𝑖−1=𝑗_𝑖−1)⟶𝑒^{−𝛼𝑖(𝑧−𝑧⋆𝑖−1)} (3.3)for 𝑡₀→∞. (Subscript 𝑠₀ in 𝑃_𝑠0 is used to underline the implicit dependency of 𝑠₀=(𝑡₀,𝑗₀).) Recall the transition probabilities (2.10) in the original jump chain. It follows immediately that𝑃_𝑠0(𝐽_𝑖=𝑗_𝑖∣𝑍^⋆_𝑖−1=𝑧^⋆_𝑖−1,𝐽_𝑖−1=𝑗_𝑖−1)=𝑃_𝑠0(𝐽_𝑖=𝑗_𝑖∣𝐽_𝑖−1=𝑗_𝑖−1𝛽)=^𝑖𝛼_𝑖, (3.4)where 𝛽_𝑖=𝑏_{𝑗𝑖−1,𝑗𝑖}. Combined with (3.3) this shows that, in the limit as 𝑡₀→∞, the rate of jumping to 𝑗_𝑖 from 𝑗_𝑖−1 is 𝛽_𝑖. More precisely, it demonstrates that 𝑋_𝑠0(𝑧), 𝑧≥0, converges to a continuous time, time-homogeneous Markov process with transition rate matrix 𝑄={𝑞_𝑗,𝑘}_{𝑗,𝑘≥0} given by 𝑞_𝑗,𝑘=𝑏_𝑘,𝑗 for 𝑗≠𝑘, and 𝑞_𝑗,𝑗=−𝑞_𝑗=∑_𝑘≠𝑗𝑞_𝑗,𝑘. This sum is indeed finite because by assumption 𝑏_𝑘,𝑗=0 for 𝑘>𝑗+1 (see Section 2).

Note that a stationary equation {𝜋_𝑗}_𝑗≥0 for the continuous-time Markov chain fulfills the equation in Proposition 2.1 with 𝑓(𝑗) replaced by 𝜋_𝑗.

4. The Partial Duplication Model

Consider the model {𝐺_𝑡}_𝑡≥𝑠, where 𝐺_𝑠 is a simple graph with 𝑠 vertices, and where 𝐺_𝑡+1 is obtained from 𝐺_𝑡 in the following way: introduce a new vertex 𝑣 and choose 𝑢∈𝐺_𝑡 uniformly. With probability 𝑞, connect 𝑣 and 𝑢. Independently of each other, connect each neighbor of 𝑢 to 𝑣 with probability 𝑝.

In this section, we follow the path outlined in the previous section. That is, we first find the jump chain corresponding to the partial duplication model. As already stated in Section 1, the master equation is given by𝑛_𝑡+1(𝑘)=(1−𝑞+𝑘𝑝𝑡)𝑛_𝑡(𝑘)+(1−𝑞)_𝑗≥𝑘(𝑗𝑘)𝑝^𝑘(1−𝑝)^𝑗−𝑘𝑛_𝑡(𝑗)𝑡+𝑞+(𝑘−1)𝑝𝑡𝑛_𝑡(𝑘−1)+𝑞_{𝑗≥𝑘−1}(𝑗𝑘−1)𝑝^𝑘−1(1−𝑝)^{𝑗−𝑘+1}𝑛_𝑡(𝑗)𝑡. (4.1)

It can be seen in the following way: the first term corresponds to the case where a vertex of degree 𝑘 keeps its degree, and this is the case unless one of two things happens: (i) the vertex is copied and receives a link to the new vertex, or (ii) it receives a link because one of its 𝑘 neighbors is copied. The probabilities of these two events are 𝑞/𝑡 and 𝑘𝑝/𝑡, respectively. Similarly, the third term corresponds to the case where a vertex of degree 𝑘−1 gets a new link in one of the above-mentioned ways. The two remaining terms correspond to the cases where the new vertex has degree 𝑘. The new vertex has degree 𝑘 when a vertex of degree ≥𝑘 is copied and receives exactly 𝑘 links to the neighbors of the copied vertex and no link to the copied vertex, or if a vertex of degree ≥𝑘−1 is copied and receives a link to the copied vertex and exactly 𝑘−1 links to the neighbors of the copied vertex.

The cases 𝑞=0 and 𝑞=1 have been studied in [19, 26], respectively. Note, however, that the master equation given in [26] is incorrect. For general 𝑞, the model has been discussed in [18]. It follows immediately that𝑓_𝑡+1(𝑘)=(1−𝑞+𝑘𝑝𝑡)𝑡𝑓𝑡+1^𝑡(𝑘)+1−𝑞𝑡+1_𝑗≥𝑘𝑏_𝑗(𝑘)𝑓_𝑡+(𝑗)𝑞+(𝑘−1)𝑝𝑓𝑡+1^𝑡𝑞(𝑘−1)+𝑡+1_{𝑗≥𝑘−1}𝑏_𝑗(𝑘−1)𝑓_𝑡(𝑗), (4.2) where we, in order to simplify notation, define𝑏_𝑗𝑗𝑘(𝑘)=()𝑝^𝑘(1−𝑝)^𝑗−𝑘. (4.3)From (4.2), we may read off the description of the matrix 𝐵(𝑡). Its entries satisfy that(𝑡+1)𝐵(𝑡)_𝑗,𝑘=(1−𝑞)𝑏_𝑗(𝑘)+𝑞𝑏_𝑗(𝑘−1)for𝑘<𝑗,𝑡−(𝑞+𝑗𝑝)+(1−𝑞)𝑏_𝑗(𝑗)+𝑞𝑏_𝑗(𝑗−1)for𝑘=𝑗,𝑞+𝑗𝑝+𝑞𝑏_𝑗(𝑗)for𝑘=𝑗+1, (4.4)and 𝐵(𝑡)_𝑗,𝑘=0 otherwise. An easy calculation shows that(𝑡+1)_𝑘≠𝑗𝐵(𝑡)_𝑗,𝑘=1+𝑞+𝑗𝑝−[(1−𝑞)𝑏_𝑗(𝑗)+𝑞𝑏_𝑗(𝑗−1)] (4.5)from which it follows that the probability of jumping from state 𝑗 is1−𝐵(𝑡)_𝑗,𝑗=1+𝑞+𝑗𝑝+𝑡+1(1−𝑞)𝑏_𝑗(𝑗)+𝑞𝑏_𝑗(𝑗−1)𝑡+1. (4.6)Motivated by this formula, we allow the jump chain to stay in state 𝑗 with probability (1−𝑞)𝑏_𝑗(𝑗)+𝑞𝑏_𝑗(𝑗−1), and it follows that the transition probabilities 𝑝_𝑗,𝑘 in the modified jump chain satisfy that(1+𝑞+𝑗𝑝)𝑝_𝑗,𝑘=(1−𝑞)𝑏_𝑗(𝑘)+𝑞𝑏_𝑗(𝑘−1)for𝑘≤𝑗,𝑞+𝑗𝑝+𝑞𝑏_𝑗(𝑗)for𝑘=𝑗+1, (4.7)and 𝑝_𝑗,𝑘=0 otherwise.

In particular, the chain is irreducible if and only if 0<𝑞<1. If 𝑞=0, the state 0 is absorbing, and if 𝑞=1, the state 0 is not reachable from any other state. If state 0 is ignored, the resulting chain is irreducible for 𝑞=1.

4.1. Classification of States

We first recall a theorem from [28]. The theorem is reformulated in [29], and we will use that formulation. If 𝑞=1, then we ignore the state 0, and since in this case all 𝑝_𝑗,0 are zero, the conditions stated in theorems below stay the same.

Theorem 4.1. Let {𝑝_𝑗,𝑘}_{𝑗,𝑘≥0} be a Markov chain. If there exist a sequence of non-negative real numbers {𝑥_𝑗}_𝑗≥0 and an integer 𝑁≥1 such that^∞_𝑘=0𝑝_𝑗,𝑘𝑥_𝑘≤𝑥_𝑗∀𝑗≥𝑁,𝑥_𝑗⟶∞for𝑗⟶∞, (4.8)then the chain is ultimately recurrent.

Applied to the partial duplication model the theorem states that if there is a sequence {𝑥_𝑗}_𝑗≥0 of nonnegative real numbers with 𝑥_𝑗→∞ such that^𝑗+1_𝑘=0𝑥_𝑘𝑝_𝑗,𝑘≤𝑥_𝑗∀𝑗≫0, (4.9)then, if 𝑞=0, the probability of ultimate absorption in 0 is 1. If 𝑞≠0, the conclusion of the theorem is that all states are persistent.

The solution 𝑝 of log(𝑝)+𝑝=0, where log denotes the natural logarithm, is known as the omega constant, and we denote it by Ω. We have Ω≈0.5671.

Proposition 4.2. Let 𝑝<Ω in the partial duplication model. If 𝑞=0, the probability of ultimate absorption in 0 is 1, and if 𝑞>0, the Markov chain is persistent.

In [26] it is claimed that for 𝑞=0 there exists a limiting distribution different from the one we find.

Proof. Let 𝑥_𝑗=log(𝑗+1). Then {𝑥_𝑗}_𝑗≥0 is a nonnegative sequence of real numbers with 𝑥_𝑗→∞, and hence it suffices to show that, for the choices of 𝑝 and 𝑞 in the proposition, the sequence satisfies (4.9). Since log is a concave function, Jensen's inequality implies that 𝐸(log(𝑋+1))≤log(𝐸(𝑋)+1) for a positive random variable 𝑋. In particular, using this for binomially distributed random variables, we get^𝑗+1_𝑘=0𝑥_𝑘𝑝_𝑗,𝑘≤(1−𝑞)log(𝑗𝑝+1)+𝑞log(𝑗𝑝+2)+(𝑞+𝑗𝑝)log(𝑗+2)1+𝑞+𝑗𝑝, (4.10)and hence we need only prove that the right-hand side of this inequality is less than or equal to log(𝑗+1) for 𝑗≫0. This may, for 𝑗≫0, be rewritten as𝑞log(𝑗𝑝+2𝑗𝑝+1)+𝑞log(𝑗+2𝑗+1)+log(𝑗𝑝+1𝑗+1)+𝑗𝑝log(𝑗+2𝑗+1)≤0, (4.11)and here the first two log-terms converge to 0, while the two remaining terms converge to log(𝑝) and 𝑝, respectively. Here we have used that𝑗𝑝log(𝑗+2𝑗𝑗+1)=𝑝𝑗+1log(1+1/(𝑗+1))1/(𝑗+1)⟶𝑝for𝑗⟶∞. (4.12)Note that since 𝑝<Ω by assumption, we have log(𝑝)+𝑝<0, and hence the inequality in (4.11) holds for all 𝑗≫0.

Since zero is the only absorbing state, it follows that for 𝑝≥Ω, a limiting distribution takes the form (𝑎₀,0,0,…), with 𝑎₀≤1. To infer the behaviour of the Markov chain for other values of 𝑞, we first recall a result proved in [30].

Theorem 4.3. Let {𝑝_𝑗,𝑘}_{𝑗,𝑘≥0} be an irreducible, aperiodic Markov chain. If there exist a sequence of positive real numbers {𝑥_𝑗}_𝑗≥0 and an integer 𝑁≥1 with^∞_𝑘=0𝑝_𝑗,𝑘𝑥_𝑘≤𝑥_𝑗∀𝑗≥𝑁,𝑥_𝑗⟶0for𝑗⟶∞, (4.13)then the chain is transient.

Let Φ denote the golden ratio conjugate. That is, Φ is the unique positive real number 𝑝 satisfying that 1/𝑝=𝑝+1. We have Φ≈0.6180.

Proposition 4.4. Let 𝑞>0 in the partial duplication model. Then the Markov chain is transient for all 𝑝>Φ.

Proof. Put 𝑥_𝑗=1/(𝑗+1) for all 𝑗≥0. Then 𝑥_𝑗>0 for all 𝑗≥0, and 𝑥_𝑗→0. Thus, in order to apply Theorem 4.3, we only need to verify that {𝑥_𝑗}_𝑗≥0 is a solution to the inequalities in (4.9).
It follows from a straightforward calculation that(1+𝑞+𝑗𝑝)^𝑗+1_𝑘=0𝑝_𝑗,𝑘𝑥_𝑘≤1+(𝑗+1)𝑝𝑞+𝑗𝑝𝑗+2 (4.14)such that {𝑥_𝑗}_𝑗≥0 is a solution if the right-hand side of this inequality is less than or equal to (1+𝑞+𝑗𝑝)/(𝑗+1) for 𝑗≫0. This is equivalent to1𝑝−𝑞+𝑗𝑝𝑗+2≤1for𝑗≫0, (4.15)and the left-hand side converges to 1/𝑝−𝑝 as 𝑗→∞. Since 𝑝>Φ, it follows that 1/𝑝−𝑝<1, and hence the inequality in (4.15) holds for all 𝑗≫0.

Let 𝑞>0 such that the chain is irreducible. One may ask for which 𝑝 the chain is ergodic. By Proposition 4.4, a necessary condition is 𝑝<Φ. However, as we will see, this may not be sufficient. To see this, we first recall another theorem from [28].

Theorem 4.5. Let {𝑝_𝑗,𝑘}_{𝑗,𝑘≥0} be an irreducible, aperiodic Markov chain. If there exist an 𝑁≥1 and a nonnegative sequence {𝑥_𝑗}_𝑗≥0 of real numbers such that^∞_𝑘=0𝑝_𝑗,𝑘𝑥_𝑘≤𝑥_𝑗−1for𝑗≥𝑁,^∞_𝑘=0𝑝_𝑗,𝑘𝑥_𝑘<∞for𝑗<𝑁, (4.16)then the chain is ergodic.

In the partial duplication model, the second condition in the theorem is always fulfilled since 𝑝_𝑗,𝑘=0 for 𝑘>𝑗+1. Let 𝑋_𝑡 denote the state of the chain at time 𝑡. If there exists 𝑁≥0 and 𝜀>0 such that𝐸[𝑋_𝑡∣𝑋_𝑡−1=𝑗]≤𝑗−𝜀∀𝑗≥𝑁, (4.17)then this 𝑁, together with 𝑥_𝑗=𝑗/𝜀, will work in Theorem 4.5. This is pointed out in [28].

Proposition 4.6. Let 𝑞>0. Then the Markov chain is ergodic for all 𝑝<1/2.

Proof. We find𝐸[𝑋_𝑡∣𝑋_𝑡−1=𝑗]−𝑗=𝑗(2𝑝−1)+2𝑞11+𝑞+𝑗𝑝⟶2−𝑝for𝑗⟶∞. (4.18)Note that 𝑝<1/2 implies 2−1/𝑝<0, and hence 2−1/𝑝≤−𝜀 for all sufficiently small 𝜀>0. That is, for a large 𝑁, (4.17) is fulfilled.

In general, it is not an easy task to actually find the stationary distribution of the jump chain or the time-inhomogeneous Markov chain. For 𝑞=1, an attempt to solve (2.4) has been made in [19]. They assume that {𝑓_𝑡(𝑗)}_𝑗≥0 converges and show that, under this assumption, the limit (for 𝑝>0) has a power-law tail. However, this does not establish the existence of a stationary distribution. Further, the power-law they provide for 𝑝>Ω is in fact not a distribution. In the special case 𝑝=0, the stationary distribution is 𝜋_𝑗=(1/2)^𝑗 for 𝑗≥1.

It is natural to ask what happens for the values of 𝑝 not covered in the propositions above. In general, this is difficult. However, if Ω is not the maximal upper bound in Proposition 4.2, the culprit must be the particular choice of {𝑥_𝑗}_𝑗≥0. Indeed, the damage provided by the use of Jensen's inequality is not severe. This may be seen in the following way: denote by 𝜇_𝑘(𝑗) the 𝑘th central moment of a binomially distributed random variable 𝑋 with parameters 𝑗 and 𝑝. From [31], we get 𝜇_𝑘(𝑗)=𝑂(𝑗^−𝑘/2), and by expanding log(𝑋+1) as a Taylor series around 𝑗𝑝, it follows that 𝐸[log(𝑋+1)]=log(𝑗𝑝+1)+𝑂(𝑗⁻¹).

4.2. Application to Protein Interaction Networks

We used the computer program developed for [18] to estimate the parameters under the partial duplication model for different protein interaction networks. The Plasmodium falciparum (P. falciparum) dataset is obtained from [32], and the remaining datasets are downloaded from the Database of Interacting Proteins (http://dip.doe-mbi.ucla.edu). Curiously, we note that according to Proposition 4.6, all pairs of 𝑝 and 𝑞 correspond to ergodic Markov chains, indicating that the networks stabilize as the number of vertices becomes large.

For one of the networks, P. falciparum, we conducted some further experiments where 50 networks were simulated with the same number of vertices as in P. falciparum (1271) and the degree distribution was computed. All simulations were started from an initial network of two vertices connected by an edge. Furthermore, 1271 runs of the corresponding Markov chain were performed, and the degree distribution was calculated and compared to the degree distribution obtained from the simulated networks. Here, the initial state of the Markov chain is 1. The length of the runs was varied, as shown in Figure 1.

(a)

(b)

Figure 1 

Shown is the distribution of vertex degrees of 50 simulated networks (solid) and that of numbers simulated from the corresponding Markov chain (dashed), using parameters estimated for the P. falciparum dataset. In addition, the observed degree distribution for P. falciparum is shown (dot-dashed).

The simulations indicate that the Markov chain approach may be used to approximate the degree distribution. This is particularly useful for simulation of large networks in terms of memory usage; storing the connections between vertices requires memory capacity proportional to the number of vertices times the average number of connections. Simulation of the corresponding Markov chain requires memory capacity proportional to the current value of the chain.

The empirical degree distribution for P. falciparum shows that the partial duplication model does not provide a perfect fit. For example, no zero-degree vertices are included in the dataset (experimenter's choice), and this needs to be incorporated into the model.

5. Other Models

We have applied the Markov chain approach to other models, and in this section we briefly present some of the results.

5.1. The Duplication-Divergence Model

The duplication-divergence model is an extension of the partial duplication model, and it has been used for analysis of protein interaction networks as well [15, 16, 27, 33]. However, the model is slightly more complicated than the partial duplication model, and it has three parameters 𝑝, 𝑞, and 𝑟. A step in this model is as follows: pick a vertex 𝑢 in the graph uniformly, and add a new vertex 𝑣. Connect 𝑢 and 𝑣 with probability 𝑞, and create an edge 𝑒_𝑤 between 𝑣 and 𝑤 whenever there is an edge 𝑒^′_𝑤 between the vertices 𝑢 and 𝑤. Now modify the pairs (𝑒_𝑤,𝑒^′_𝑤) independently of each other in the following way: with probability 𝑝, keep both edges; otherwise, with probability 𝑟, keep 𝑒_𝑤 and delete 𝑒^′_𝑤, and with probability 1−𝑟, keep 𝑒^′_𝑤 and delete 𝑒_𝑤.

One can derive a master equation and go through the construction of the modified jump chain. In this case, the transition probabilities 𝑝_𝑗,𝑘 satisfy that(𝑗𝑝+2)𝑝_𝑗,𝑘=(1−𝑞)𝑏_𝑗(𝑘,1−𝜓)+𝑞𝑏_𝑗(𝑘−1,1−𝜓)+(1−𝑞)𝑏_𝑗(𝑘,𝑝+𝜓)+𝑞𝑏_𝑗(𝑘−1,𝑝+𝜓)for𝑘≤𝑗,𝑗𝑝+𝑞𝑏_𝑗(𝑗,1−𝜓)+𝑞𝑏_𝑗(𝑗,𝑝+𝜓)for𝑘=𝑗+1, (5.1)and 𝑝_𝑗,𝑘=0 otherwise. Here 𝜓=(1−𝑝)(1−𝑟), and 𝑏_𝑗(𝑘,𝑠) is the binomial probability from (4.3) with 𝑝 replaced by 𝑠.

In order to apply Theorem 4.1, we put 𝑥_𝑗=log(𝑗+1). It follows from simple calculations, again using Jensen's inequality, that {𝑥_𝑗}_𝑗≥0 is a solution to (4.9) if 𝑝 and 𝑟 satisfy thatlog(1−𝜓)+log(𝑝+𝜓)+𝑝<0. (5.2)Note that in the special cases 𝑟=0 and 𝑟=1, the left-hand side of the inequality reduces to log(𝑝)+𝑝, the same inequality as seen earlier. Actually, for 𝑟=0 the model is the partial duplication model. It follows that if 𝑟=0 or 𝑟=1, a solution 𝑝 of (5.2) must satisfty that 𝑝<Ω. For 0<𝑟<1 an exact upper bound on 𝑝 is harder to derive. For these values of 𝑟, the solution 𝑝 is less than Ω and attains a minimum 𝑝≈0.5235 for 𝑟=1/2.

5.2. Another Class of Models

We believe that the Markov chain approach presented in this paper may be used to infer the behavior of other classes of models. In [21], simple models with master equations on the form𝑛_𝑡+1𝑎(𝑘)=(1−^𝑘𝑡)𝑛_𝑡𝑎(𝑘)+^𝑘−1𝑡𝑛_𝑡(𝑘−1)+𝑐_𝑘, (5.3)where 𝑎_𝑘 and 𝑐_𝑘 are nonnegative numbers, are studied. The resulting master equation for the relative frequencies 𝑓_𝑡(𝑘) may be written in matrix form as1(𝑡𝑡+1𝟏𝐜1𝐴(𝑡)𝐟)(_𝑡1𝐟)=(_𝑡+1), (5.4)where 𝟏=(111…), and 𝐜 and 𝐟_𝑡 are the column vectors consisting of all the numbers 𝑐_𝑘 and 𝑓_𝑡(𝑘), respectively. The matrix 𝐴(𝑡) is given by𝐴(𝑡)=𝑡−𝑎₀𝑎000…₀𝑡−𝑎₁00…0𝑎₁𝑡−𝑎₂0…00𝑎₂𝑡−𝑎₃…⋮⋮⋮⋮⋱. (5.5)

Note that columns of the partitioned matrix in (5.4) sum to 𝑡+1. That is, when divided by 𝑡+1, the transpose of this matrix represents a Markov chain with time-dependent transition probabilities. We identify the countable set of states with 𝐍∪{−∞} where the artificial state −∞ accounts for the first row and the first column in the partitioned matrix.

We may compute the corresponding jump process, and again it turns out that its transition probabilities 𝑝_𝑗,𝑘 are time-independent. We may get rid of the state −∞ by simply forgetting the time we spend there. That is, for 𝑗,𝑘≥0, we replace 𝑝_𝑗,𝑘 by the sum 𝑝_𝑗,𝑘+𝑝_𝑗,−∞𝑝_−∞,𝑘, and this leads to a Markov chain with transition probabilities given by𝑝_𝑗,𝑘=𝑎_𝑗+𝑐_𝑗+11+𝑎_𝑗𝑐for𝑘=𝑗+1,_𝑘1+𝑎_𝑗otherwise. (5.6)

These jump chains are in fact all ergodic, and the stationary distribution of the time-inhomogeneous Markov chains has been derived in [21].

5.3. Other Extentions

Still other models do not fall under the conditions and assumptions introduced in this paper. For example, the master equation of the most general form of the duplication-mutation model [22, 23] depends on terms 𝑂(1/𝑡²), and the columns of 𝐴(𝑡) do not sum to the same number 𝑎(𝑡) because of 𝑂(1/𝑡²) terms, and because the requirement 𝐴(𝑡)_𝑘,𝑗=0 for 𝑘>𝑗+1 is not fulfilled.

Some of these problems may be circumvented at the cost of a more technical and elaborate exposition, but often the results need to be stated as limiting results. For example, if the columns of 𝐴(𝑡) do not sum to the same number, the jump chain in (2.10) should be considered as emerging in the limit 𝑡→∞.

Furthermore, one may choose to ignore terms of order 𝑂(1/𝑡²) in the master equation. As 𝑡→∞, the influence from higher-order terms often becomes insignificant, justifying such an approximation. This is, for example, the case for the duplication-mutation model.

Table 1 

Parameters estimated from protein interaction data.

Acknowledgments

M. Knudsen is supported by the Centre for Theory in the Natural Sciences, University of Aarhus. C. Wiuf is supported by the Danish Cancer Society and the Danish Research Councils. They would like to thank an anonymous reviewer for valuable suggestions that improved the clarity of the paper.