Bayesian Inference for Duplication–Mutation with Complementarity Network Models

Abstract We observe an undirected graph G without multiple edges and self-loops, which is to represent a protein–protein interaction (PPI) network. We assume that G evolved under the duplication–mutation with complementarity (DMC) model from a seed graph, G0, and we also observe the binary forest Γ that represents the duplication history of G. A posterior density for the DMC model parameters is established, and we outline a sampling strategy by which one can perform Bayesian inference; that sampling strategy employs a particle marginal Metropolis–Hastings (PMMH) algorithm. We test our methodology on numerical examples to demonstrate a high accuracy and precision in the inference of the DMC model's mutation and homodimerization parameters.


Introduction
As a result of breakthroughs in biotechnology and high-throughput experiments thousands of regulatory and protein-protein interactions have been revealed and genome-wide proteinprotein interaction (PPI) data are now available.Protein-protein interactions are one of the most important components of biological networks as they are fundamental to the functioning of cells.To gain a better understanding of why these interactions take place, it is necessary to view them from an evolutionary perspective.The evolutionary history of PPI networks can help answer many questions about how present-day networks have evolved and 1 arXiv:1504.01794v1[stat.CO] 8 Apr 2015 provide valuable insight into molecular mechanisms of network growth [12,9].However, inferring network evolution history is a statistical and computational challenging problem as PPI networks of extant organisms provide only snapshots in time of the network evolution.
There has been recent work on reconstructing ancestral interactions (e.g.[5,7,11]).The main growth mechanism of PPI network is gene duplication and divergence (mutations) [15]: all proteins in a family evolve from a common ancestor through gene duplications and mutations and the protein network reflects the entire history of the genome evolution [14].
In this paper we follow [10] and we develop computational methods to infer the growth history and the parameters under the given model incorporating not only the topology of observed networks, but also the duplication history of the proteins contained in the networks.In their article [10] propose a maximum likelihood approach.The authors establish a neat representation of the likelihood function and it is this representation which is used in this article.The duplication history of the proteins can be inferred independently by phylogenetic analysis [11,13].
The approach we adopt here is first to obtain a numerically stable estimate of the likelihood function, under fixed parameters; this achieved via the sequential Monte Carlo (SMC) method (see [4] and [8]).This approach can then be used to infer the parameters of the model, from a Bayesian perspective, as well as the growth history, via a Markov chain Monte Carlo (MCMC) method.To the best of our knowledge, this has not been considered in the literature, although related ideas have appeared for simpler models in [16].Our computational strategy not only improves on likelihood estimation in comparison to [10], but also provides a natural setup to perform posterior inference on the parameters of interest.This article is structured as follows.In Section 2 we detail the model and associated computational method for statistical inference.In Section 3 our numerical results are presented.
In Section 4 the article is concluded with some discussion of future work.

Model and Methods
We follow similar notation and exposition as in [10] to introduce the protein-protein interaction network, its duplication history, and the duplication-mutation with complementarity (DMC) model [14].In particular, the notions of adjacency and duplication are made concrete there.We also introduce the associated Bayesian inference problem with which this work is primarily concerned (i.e., that of inferring the parameters of the DMC model).We then describe a particle marginal Metropolis-Hastings (PMMH) algorithm [1] that can be used to perform such inference.

PPI network and DMC model
Consider an undirected graph G without multiple edges and self-loops, where the nodes represent proteins and the edges represent interactions between those proteins.Such a graph is called a PPI network, and as in [10], we denote the vertex set by V (G), the edge set by E(G), and the number of nodes in G by |V (G)|.All nodes which are adjacent to a node v (not including v itself) comprise the neighbourhood of v, and that neighbourhood is We assume that G evolved from a seed graph G 0 via a series of duplication, mutation, and homodimerization steps under a DMC model.Under the DMC model, at each time step t, the graph G t evolves from G t−1 by the following processes in order: 1.The anchor node u t is chosen uniformly at random from V (G t−1 ), and a duplicate node v t is added to G t−1 and connected to every member of N Gt−1 (u t ).This is the duplication step, and it yields an intermediary graph denoted G * t−1 .
2. For each w ∈ N G * t−1 (u t ), we uniformly choose one of the two edges in {(u t , w), (v t , w)} ⊆ E(G * t−1 ) at random and delete it with probability (1 − p).This is the mutation step, and the parameter p is henceforth referred to as the mutation parameter.
3. The anchor node u t and the duplicate node v t are connected with probability p c to finally obtain G t .This is the homodimerization step, and the parameter p c is henceforth known as the homodimerization parameter.
The DMC model is Markovian, and we denote the transition density at time t (which encompasses the three afore-mentioned steps) by p M (G t | G t−1 ), where M := (p, p c ).If we assign to a seed graph some prior density p M (G 0 ), then the density of the observed graph where

Bayesian inference
In practice, one will not have access to the parameters (p, p c ) and they must be inferred given G. Thus, in the Bayesian setting, our objective is to consider the posterior density where p(M) is some proper prior for (p, p c ) that we assume can easily be computed (at least pointwise up to a normalising constant).
The total number of growth histories grows exponentially with n [10], and so any computations involving (1), and thus ( 2), (e.g. the evaluation of p M (G)) could potentially become very expensive.In the following sections, we reformulate the inference problem in the same manner as in [10] to alleviate this issue.

Duplication history
As in [10], let Γ be a binary forest, i.e. a collection of rooted binary trees.The authors of [10] describe a scheme that encodes the duplication history of a growth history H within a series of duplication forests, (Γ 0 , Γ 1 , . . ., Γ n ), where each forest Γ t corresponds to a graph G t .We describe that scheme here.
Consider a trivial forest Γ 0 , whose only two isolated trees each consist of a single node.
Each of those isolated nodes will correspond to a node within the seed graph G 0 .To build Γ 1 from Γ 0 , one replaces an anchor node u 1 from Γ 0 with a subtree, {u 1 , v 1 }, consisting of two leaves (v 1 is the duplicate node included G 1 but not G 0 ).This process continues until one builds the series of forests (Γ 0 , Γ 1 , . . ., Γ n = Γ) to correspond to H.
As highlighted in [10], the duplication forest Γ (corresponding to G) is uniquely determined by H and a list of anchor nodes, π = (u 1 , . . ., u n ).The important thing to emphasize here, is that given the duplication forest Γ and G, one now only needs to infer the duplication nodes sequence θ(H) = (v 1 , v 2 , . . ., v n ) to reconstruct the complete growth history H.
For instance, at the first step backwards, knowledge of Γ n = Γ together with v n uniquely identifies the anchor node u n , thus one can reconstruct G n−1 and Γ n−1 ; this is then repeated for the remaining backward steps.Thus, given θ(H), one can construct the growth history H backward-in-time using the backward operators defined in [10, Section 2.4], which

Bayesian inference given the duplication history
Now suppose that in addition to G, a practitioner is given Γ corresponding to G. Our new objective -and the primary inference problem with which this work is concerned -is to consider the posterior density π(M|G, Γ).Notice that we have the joint distribution: where θ = (v 1 , . . ., v n ) is a sequence of duplication nodes compatible with the observed G, Γ, and G θ 0 , Γ θ 0 , . . ., G θ n , Γ θ n the corresponding reconstructed history.We are thus interested in the parameter posterior: The density p M (G 0 , Γ 0 ) is typically a trivial term which could be ignored in practice.As the duplication forest Γ limits the number of allowable anchor-and-duplicate node pairs, one can see that the number of possible growth histories is reduced.

Methods
We will now present an SMC algorithm that can sample the latent growth histories from the DMC model given the fixed parameters M := (p, p c ).We then show that this algorithm can be employed within a PMMH algorithm, as in [1], to sample from the posterior (3) and infer M (and even θ|G, Γ).
An SMC algorithm simulates a collection of N samples (or, particles) sequentially along the index t via importance sampling and resampling techniques to approximate a sequence of probability distributions of increasing state-space, which are known pointwise up to their normalising constants.In this work, we use the SMC methodology to sample from the posterior distribution of the latent duplication history: backwards along the index t via Algorithm 1 in the Appendix.The technique provides an unbiased estimate of the normalising constant [3, Theorem 7.4.2],p M (G, Γ): where each W i t is an un-normalised importance weight computed in Algorithm 1.Note that under assumptions on the model, if N > cn for some c < ∞, then the relative variance of the estimate is O(n/N ); see [2].It is remarked that, as in [16], one could also use the discrete particle filter [6], with a possible improvement over the SMC method detailed in Algorithm 1; see [16] for some details.This SMC can be employed within a PMMH algorithm to target the posterior of M in (3).One can think of the deduced method as an MCMC algorithm running on the marginal M-space, but with the SMC unbiased estimate pM (G, Γ) replacing the unknown likelihood pM (G, Γ).More analytically, we can consider all random variables involved in the method and write down the equilibrium distribution in the enlarged state space, with M-marginal the target posterior p M (G, Γ).Following [1] and letting φ i t denote a sample (G t , Γ t ) at time t, the extended equilibrium distribution is written as: where Ψ M is the probability of all the variables associated to Algorithm 1, with a j k , l ∈ {1, . . ., N } and the φ's being the simulated variables at each step of Algorithm 1.
A PMMH algorithm (see Algorithm 2) samples from (5), and one can remove the auxiliary variables from the samples to obtain draws for the parameters from (3).Furthermore, one could even save the sampled growth histories with particle index l to obtain draws from the joint posterior π(M, θ | G, Γ).However, in this work, we are primarily interested in the inference of M.

Results
The variance of the estimate (4) plays a crucial role in the performance of Algorithm 2, as ( 4) is used to compute the acceptance probability within the PMMH algorithm.Thus, we first tested the variability of (4) as computed by Algorithm 1 to understand how the variance changes with |V (G)|.We then ran Algorithm 2 to sample from the posterior (3) and infer M for a given pair of observations (G, Γ).We present the details of those experiments below.

Variance of pM (G, Γ)
We simulated a graph G and a forest Γ from the DMC model [14]  As remarked above, if N > cn for some c < ∞, then the relative variance of pM (G t , Γ t ) is O(n/N ). Figure 2 confirms that the variance increases linearly and that increasing the value of N with |V (G)| (at least linearly) will help to control the variance.However, the plots show that the relative variance is still high, which means that N will have to be large to ensure satisfactory performance of the PMMH in practice.

Parameter inference
We separately simulated a graph G and a forest Γ from the DMC model with the parameters again set as (p = 0.7, p c = 0.7), and we set |V (G)| = 15.Given only (G, Γ), we inferred (p, p c ) with each parameter having a uniform prior on the interval [0.1, 0.9].We set the number of particles within Algorithm 1 to be N = 2000, and we ran the PMMH algorithm to obtain 10,000 samples from the extended target.Figure 3 in Section A illustrates good mixing of the PMMH algorithm and accurate inference of the parameters (p, p c ).The trace plots show that the algorithm is not sticky, and the autocorrelation functions give evidence to an approximate independence between samples.The posterior densities are also interesting, in that they are clearly different from the uniform priors and they show that the PMMH algorithm spends a majority of the computational time sampling the true parameter values.

Discussion
We have introduced a Bayesian inferential framework for the DMC model [14], where, as in [10], one assumes the pair (G, Γ) is observed and the parameters (p, p c ) are unknown.
We then described how an SMC algorithm can be used to simulate growth histories and ultimately be employed within PMMH to target the posterior distribution of the parameters (3), thereby opening up the possibility of performing Bayesian inference on the DMC model.

Numerical tests demonstrated that Algorithm 1 can have a high variability when |V (G)|
is large and N is not sufficiently high, and this limits the scope of the inference problem which can be tackled using the complete Algorithm 2. However, the proposals used in the example experiments within Algorithms 1 and 2 are naive, as the method simply chooses a candidate duplicate node v i t at random from all permitted nodes given the current G i t , Γ i t .It is reasonable to assume that more sophisticated proposal densities could reduce the variance of the SMC and/or improve the mixing of the PMMH, thereby allowing one to perform inference when |V (G t )| is large and N is smaller.This could be explored in a future work.

A Figures
Variability of pM (G, Γ) Recall that the seed graph, G 0 , has two nodes, and note that we did not compute pM (G 0 , Γ 0 ) because it is trivial.
denotes the collection of growth histories.In this work, a seed graph will always be the graph consisting of two connected nodes; thus, |V (G 0 )| = 2 and p M (G 0 ) = 1.Note that we are summing over all possible growth histories by which G can evolve from a seed graph.Also note that a growth history H induces a unique sequence of duplicate nodes, θ(H) = (v 1 , . . ., v n )[10].
with the parameters set as (p = 0.7, p c = 0.7), where |V (G)| = 40.We saved each pair (G t , Γ t ) for 1 ≤ t ≤ 40, and we ran Algorithm 1 fifty times per pair (with N = |V (G t )| * 5) to compute fifty unbiased estimates of p M (G t , Γ t ) for 1 ≤ t ≤ 40.In the top of Figure 2 in Section A, we plot the relative variance of the estimate (or, the variance divided by the square of the expected value) per each value of |V (G t )|.We repeated the experiment two more times, with N = |V (G t )| * 10 and N = |V (G t )| * 20, and the associated output is also presented in Figure 2.

Figure 2 :
Figure 2: All plots illustrate the relative variance of pM (G t , Γ t ), per |V (G t )| on the horizontal axis; the relative variance is the variance divided by the square of the expected value.In the top plot, the number of SMC particles used to compute each pM (G t , Γ t ) is |V (G t )| * 5.In the middle and bottom plots, that number is |V (G t )| * 10 and |V (G t )| * 20, respectively.

Figure 3 :
Figure 3: Plots associated with p (p c ) are at left (right).The top figures are trace plots, with PMMH iteration running along the horizontal axes and parameter value running along the verticals.The middle figures are plots of the autocorrelation functions (with lag running along the horizontal axes), and at the bottom we present the parameter posterior densities.