Link Prediction with Continuous-Time Classical and Quantum Walks

Protein–protein interaction (PPI) networks consist of the physical and/or functional interactions between the proteins of an organism, and they form the basis for the field of network medicine. Since the biophysical and high-throughput methods used to form PPI networks are expensive, time-consuming, and often contain inaccuracies, the resulting networks are usually incomplete. In order to infer missing interactions in these networks, we propose a novel class of link prediction methods based on continuous-time classical and quantum walks. In the case of quantum walks, we examine the usage of both the network adjacency and Laplacian matrices for specifying the walk dynamics. We define a score function based on the corresponding transition probabilities and perform tests on six real-world PPI datasets. Our results show that continuous-time classical random walks and quantum walks using the network adjacency matrix can successfully predict missing protein–protein interactions, with performance rivalling the state-of-the-art.

partners for proteins without known links, self-interacting proteins, or links between proteins that have longer paths between them.Given the low coverage of the current PPI databases, this can be a significant drawback.It is therefore highly desirable to complement the existing frameworks with methods relying on the exploration of the whole network, and consequently be able to incorporate information provided by longer paths.This is one of the the main goals of our paper as we propose two novel random-walk based link prediction methods.

II. THE METHOD
Consider a network modelled by an undirected graph G = (V, E), where V is the set of nodes of size n and E is the set of edges.We allow for the existence of self-edges, so that for any node i, the edge (i, i) may or may not be present in E. The adjacency matrix of G is the n × n matrix defined by The graph Laplacian is defined as L = D−A, where D is the degree matrix defined by D = diag j A 1j , . . ., j A nj .The precise details of the random walks we employ will be described in the next subsections.For now, it suffices to consider the notion of a probability transition matrix that evolves over time, denoted by P (t); for a graph G, the probability of the random walker being at node v at time t, given that it began at node u, is thus P uv (t).For a fixed time t, we define the score S(i, j; t) between two non-adjacent nodes i and j at time t to be where N (v) denotes the neighbourhood of node v, and k v = j A vj is the degree of node v. Equations ( 1) and ( 2) handle the cases of distinct nodes and self-edges, respectively.While the score in Equation ( 1) is superficially similar to the one proposed in [20], the fact that we are using continuous-time random walks leads to several key differences: the continuous-time nature of our method allows for a wider range of time parameters t to use; in the continuous-time setting there is symmetry in the transition probabilities, i.e.P ij (t) = P ji (t) for all nodes i, j; and finally, there is a close relationship in the implementation of classical and quantum walks in the continuous-time setting.Regardless of which type of random walk is used, we must choose a value t, representing the time-duration of the walk.We start the walk at time t 0 = 0, and let it run for a time t, at which point we extract the scores for the target edges from the probability distributions.In the case of a continuous-time classical random walk, the expected time it takes for a random walker to leave a node i is 1/k i .This motivates the idea that the amount of time we let the random walk run should be related to the degree distribution of the network.Our experiments show that setting t = 1/ k , where k is the average node degree in the graph, yields good results for both classical and quantum walks.
Next, we describe the different type of walks we use in more detail.

A. Continuous-time random walks
A continuous-time (classical) random walk (CRW) is a Markov process with state space V characterized by a rate matrix Q and initial distribution p(0) over the set of nodes.Here, we consider edge-based random walks [23] (as opposed to node-based), which are characterized by setting Q = −L, where L is the Laplacian of the underlying graph.In this case, the dynamics are governed by the equation Intuitively, the random walker operates as follows.Every edge of the graph is associated with an independent Poisson process with unit intensity.When the walker is at some node, it will remain there until one of the Poisson processes at an incident edge jumps, at which point the walker follows that edge to the corresponding neighbour, and the process repeats.Note that this implies that, on average, a random walker will spend less time waiting at a higher degree node than at a lower degree node.

B. Continuous-time quantum walks
In contrast to a classical random walk, a quantum walk on a network evolves according to the laws of quantum physics.A major implication of this is that the trajectories of the walker across the network can interfere constructively or destructively.This interference causes the evolution of the quantum walker to sometimes be significantly different from the classical one [2,8].
A continuous-time quantum walk (QRW) [12] on a graph G is defined by considering the Hilbert space H spanned by the orthonormal vectors {|i } n i=1 , corresponding to the n nodes of the graph and the unitary transformation e −itH .This transformation implies that the state vector in H at a time t after starting from initial time t 0 = 0 is given by the evolution where U (t) = e −itH is the unitary evolution operator, and H is the Hamiltonian.In general, the Hamiltonian H can be almost any Hermitian matrix related to G as long as it describes the structure of the network [36], but the most common choices are the Laplacian L or the graph adjacency matrix A [34].In this paper we will use both the Laplacian and the adjacency matrix as the Hamiltonian separately, and therefore can compare different realizations of quantum walks for the link prediction task.
In order to obtain a probability transition matrix analogous to the one in Eq. ( 3), we must take the square of the modulus of the entries of U (t).The entries of the probability transition matrix are given by ( These transition probabilities can then be used to compute scores for missing edges as described in Equations ( 1) and ( 2) above.
We emphasize that the usage of continuous-time quantum walks for link prediction is a new direction of research, with very few studies conducted so far.The method proposed in [25], in particular, appears to be competitive with some state-of-the-art link prediction methods in certain real networks.While some aspects of their algorithm are similar to the quantum version of our random walk algorithm, the implementation details and calculation of the link prediction scores are very different.Moreover, their algorithm requires entanglement with an additional ancilla.While this would be feasible in a hypothetical implementation on a quantum computer, the typical sizes of relevant real networks are way beyond the capabilities of current and near-term quantum hardware.Simulations on classical computers are required, but the presence of the extra ancilla increases the complexity of simulations.

III. RESULTS
In this section we first describe the datasets and networks that were used for evaluation, then we compare our methods to some well-known and state-of-the-art link prediction algorithms.

A. Datasets and evaluation
We tested our link prediction methods on four different PPI networks.Three of them are high-quality PPI networks from the HINT database [10]: M. musculus and S. cerevisiae are networks resulting from combining data from high-throughput screenings with literature-curated small-scale experiments, and H. sapiens consists of only the latter, representing a subset of the human interactome.The fourth network we tested on, HI-AP-MS, is another large-scale map of the human interactome, but in this case was produced by affinity purification followed by mass spectrometry [15].Some statistics of these networks are listed below in Table I, and their degree distributions are shown in Figure 1.H. sapiens -lit.
FIG. 1. Complementary cumulative degree distributions.For each degree value k (x-axis), the proportion of nodes with degree greater than or equal to k (y-axis) is shown, each on a logarithmic scale.
For each dataset, we randomly removed P % of the edges in the original network, for P ∈ {10, 20, 30, 40, 50} and reserved these edges as positive test cases.All of the nonexistent edges were used as negative testing data.These positive and negative edges were used for evaluation, and the remaining (100 − P )% existing edges were used for running the models in question.This process was repeated 20 times for each P , and the results were averaged (see results in the next subsection).

B. Comparison to other methods
In order to test our methods, we selected 5 other popular link prediction methods to compare against: L3 relies on a weighted counting of paths of length three, and was designed specifically to predict links in PPI networks [17]; preferential attachment (PA) defines a score between two disconnected nodes by multiplying their degrees [4,19]; common neighbours (CN) is a straightforward heuristic that assigns a score to edge (u, v) defined by the number of neighbours that u and v have in common; Adamic-Adar (AA) is an adaptation of the common neighbours idea, but adds more weight to less connected neighbours [1]; SPM is based on perturbations of the adjacency matrix of the graph [22].
The following tables show average precision [32] and area under the receiver operating characteristic (ROC) curve [14] values for the 4 different networks described in Section III A. Each value is averaged over 20 runs (20 randomly selected edges removals).We compare three variations of our proposed methods, labelled as 'QRW-A', 'QRW-L', and 'CRW', referring to quantum random walks using the network adjacency matrix as Hamiltonian, quantum random walks using the network Laplacian matrix as Hamiltonian, and classical random walks, respectively.For completeness, in Figures 2-5 we also include plots showing the relationship of average precision and area under the ROC curve as a function of edge removal fraction.

IV. DISCUSSION
The experimental results in the previous section show that our methods perform well on a variety of PPI networks.In particular, we see that our classical random walk method yields the best performance of all algorithms tested with respect to the area under the receiver operating characteristic curve when 50% of the edges are removed (Table III).The quantum walks using the adjacency Hamiltonian yield the best average precision on one of the human proteome datasets when 10% of edges are removed (Table IV).Furthermore, the adjacency Hamiltonian almost always beats out the Laplacian as the better choice when comparing the results of quantum walks.However, as is usually the case with link prediction methods, there is no single method that wins in all cases.One must pay careful attention to the goals of the application in question and decide on which metric to rely on accordingly.
Finally, we mention a few points about the computational complexity of our algorithm and its implementation.The bottleneck of our algorithm, in either the classical or quantum case, is the computation of the matrix exponential  appearing in Equation (3) and Equation ( 5), which is a very well-studied problem with a long history [24].Our experiments were done using the 'matrix exp' function in PyTorch [29], which is an implementation of the Taylor polynomial approximation algorithm described in [3].The problem is thus reduced to a constant number of matrix multiplications, another well-studied problem that can be solved more quickly than the naive O(n 3 ) method; for example, in O(n 2.376 ) time using the Coppersmith-Winograd algorithm [9].It is also worth noting that in this implementation of matrix exponentiation, and many others, the norm the matrix being exponentiated has an impact on running time, so that using a small t, as tends to be the case in our algorithm, may help in this regard.

V. CONCLUSIONS
Although experimental methods have greatly improved in the past ten years, most interactomes remain far from being complete.It is therefore important to discover new computational methods for inferring interactions from incomplete datasets.We have described a class of algorithms based on continuous-time random walks that rank among the best link prediction methods tested on PPI networks.
Furthermore, the quantum continuous-time random walks described here are among the first successful quantuminspired link prediction methods.Although we have found that using the reciprocal of the average degree provides a good time length for which to run the random walks, many further options can still be explored: using cross-validation to choose a more optimal value, or using times that depend on the walker's location are immediate candidates.Another open direction of research involves the choice of Hamiltonian.Our experimental results demonstrate a strong sensitivity on the Hamiltonian used for controlling the quantum walks.While the adjacency matrix yields better results than the Laplacian on most of the networks we tested, it would be beneficial to understand why this is the case.This also indicates the potential for improvement if better Hamiltonians can be found for the purpose of link prediction.Further investigations in this direction may yield better methods and insights into both the networks being studied, and the quantum walks being employed.
Competing interests: The authors declare no competing interests.Authors contributions: MG, HS conceived the algorithm.GGP, MACR, SM designed and directed the research.MG, JM implemented the algorithms and ran simulations.MG, JM, HS wrote the first version of the manuscript.All authors contributed to scientific discussions and to the writing of the manuscript.

FIG. 2 .FractionFIG. 3 .FIG. 4 .
FIG.2.Average precisions (left) and area under the ROC curve (right) as a function of the fraction of true links that have been removed from the literature-curated Homo sapiens network found in[10].Plotted values are the averages over 20 trials.

FIG. 5 .
FIG.5.Average precisions (left) and area under the ROC curve (right) as a function of the fraction of true links that have been removed from the Saccharomyces cerevisiae network found in[10].Plotted values are the averages over 20 trials.

TABLE I .
Some properties of the networks that were tested.|V | : number of nodes, |E| : number of edges, k : average degree, ρ : network density, C : average clustering, A : assortativity, SIPs: number of self-interacting proteins (self-edges).

TABLE III .
Area under the ROC curve for 50% edge removals, averaged over 20 trials.