Identification of co-evolving temporal networks.

BACKGROUND: Biological networks describes the mechanisms which govern cellular functions. Temporal networks show how these networks evolve over time. Studying the temporal progression of network topologies is of utmost importance since it uncovers how a network evolves and how it resists to external stimuli and internal variations. Two temporal networks have co-evolving subnetworks if the evolving topologies of these subnetworks remain similar to each other as the network topology evolves over a period of time. In this paper, we consider the problem of identifying co-evolving subnetworks given a pair of temporal networks, which aim to capture the evolution of molecules and their interactions over time. Although this problem shares some characteristics of the well-known network alignment problems, it differs from existing network alignment formulations as it seeks a mapping of the two network topologies that is invariant to temporal evolution of the given networks. This is a computationally challenging problem as it requires capturing not only similar topologies between two networks but also their similar evolution patterns.
RESULTS: We present an efficient algorithm, Tempo, for solving identifying co-evolving subnetworks with two given temporal networks. We formally prove the correctness of our method. We experimentally demonstrate that Tempo scales efficiently with the size of network as well as the number of time points, and generates statistically significant alignments-even when evolution rates of given networks are high. Our results on a human aging dataset demonstrate that Tempo identifies novel genes contributing to the progression of Alzheimer's, Huntington's and Type II diabetes, while existing methods fail to do so.
CONCLUSIONS: Studying temporal networks in general and human aging specifically using Tempo enables us to identify age related genes from non age related genes successfully. More importantly, Tempo takes the network alignment problem one huge step forward by moving beyond the classical static network models.

Background

Biological networks describe the interaction between molecules. They are frequently represented as graphs, where the nodes correspond to the molecules (e.g., proteins or genes) and the edges correspond to their interactions [1]. Formally, we denote a biolon class="Disease">gical network as G=(V,E) where V and E represent the set of nodes and the set of edges, respectively. Analysis of these networks enable the elucidation of cellular functions [2], the identification of variations in cancer networks [3], and the characterization of variations in drug resistance [4]. Studying biological networks led to numerous computational challenges as well as methods which address these challenges. Network alignment is one of the most important of these challenges [5] as it has a profound set of applications ranging from the detection of conserved motifs to the prediction of protein functions [6]. This problem aims to find a mapping of the nodes of two given networks in which nodes that are similar in terms of content (i.e. homology) and interaction structure (i.e. topology) are mapped to each other. Hence, we represent the alignment between two given networks G1=(V1,E1) and G2=(V2,E2) as a bijection function ψ:V1→V2, and the score resulting from alignment ψ as score(G1,G2|ψ). The network alignment problem seeks the function ψ that maximizes this score. We note that there are various ways to calculate the scoring function.

There are two categories of network alignment problem: local and global alignment. The former problem aims to find pairs of highly-conserved sub-networks in two given networks in which a sub-network of the query network is mapped to multiple sub-networks in the target network. Global network alignment aims to maximize the n class="Chemical">similarity in the networks in which all nodes in the query network are mapped to a set of nodes in the target network. Network alignment is a challenging task as the graph and subgraph isomorphism problems which are known to be GI and NP-hard [7], reduce to them. In “Related Work” section, we give a brief review of the methods addressing the global network alignment problem as the problem we consider in this paper is associated with that problem.

Biological networks have dynamic topolon class="Disease">gies [8]. There are various reasons behind this dynamic behavior. For example, genetic and epigenetic mutations can alter molecular interactions [9], and variation in gene copy number can affect the existence of interactions [10]. Due to this dynamic behavior, the topology of the network that models the molecular interaction evolve over time [11]. Majority of the previous work on alignment of biological networks assume the network topology is static [12]—an assumption that ignores the history of network evolution, and may lead to biased or incorrect analysis. For example, identifying causes and consequences of the influence of external stimuli is impossible when analyzing static topologies. To address this oversight, we define a biological network using a model that accounts for the evolution of the underlying network at consecutive time points. We refer to this model as a temporal network [13]. We view this model as containing a single snapshot of the network at each time point in a sequence of time points and thus, as a time series network. More formally, we denote a temporal network with t consecutive time points as , where G=(V,E) represents the topology of the network at the ith time point.

In this paper, we consider the problem of identifying co-evolving subnetworks in a given pair of temporal networks. We say that two subnetworks are co-evolving if their topologies remain similar to each other even though their topologies evolve (i.e. change) over time. We define this more formally as follows. We consider two input temporal networks and , where ∀i∈{1,2,…,t}, and represent and respectively at the time point i. Without losing generality, let be the query (smaller) network and be the target network, i.e., |V1|≤|V2|. An alignment of and maps to across all time points i. Thus, we represent the alignment of the two temporal networks and as a bijection of their nodes and denote it as a function ψ:V1→V2. We compute the score of the alignment ψ of and , denoted with , as the sum of the scores of the alignment at all time points. Hence, = . We assume is connected at all time points, but it maybe impossible to find an alignment that is connected in the target network at all time points.

It is worth emphasizing that the temporal network alignment problem described above is dramatically different than existing network alignment problems, which can be categorized as follows: (i) pairwise alignment, (ii) multiple network alignment, and (iii) dynamic network alignment. We illustrate these problems as well as the temporal one in Fig. 1. The pairwise network alignment problem (Fig. 1a) ignores that the network topology evolves. Although the multiple alignment problem (Fig. 1b) can conn class="Chemical">sider more than two networks at once, it lacks the ability to capture the temporal changes since it treats all networks as having static topologies. The dynamic network alignment problem (Fig. 1c) considers topological changes over time. It however, it seeks a different solution to the alignment problem at each time point. Thus, it can not identify co-evolving subnetwork. A new algorithm is needed to capture such evolving characteristics. Unlike these alignment problems, temporal network alignment (Fig. 1d) captures that network topologies co-evolve over time.

Fig. 1

This figure represents different network alignment problems in different types of biological networks. a This represents the alignment between two input static networks. b This represents the alignment between multiple time points where each network represent a different organism. c This represents the alignment between two input networks where one of them is static and one of them is dynamic. Here, there exist different alignment between the static network and each version of the dynamic network. d This represents the alignment between two input temporal networks where each have time specific snapshots that was taken at three specific time points. Here, the alignment is persist across all time points

Contributions in this paper. We develop an efficient algorithm, Tempo, to identify co-evolving subnetworks in a given pair of the temporal networks. More specifically, we aim to find subnetworks of given networks which have similar evolving topologies over time. Briefly, our algorithm first finds an initial alignment between the input networks and using the similarity score between pairs of aligned nodes across all time points. It then performs a dynamic programming strategy that maximizes the alignment quality (i.e. score) by repeatedly altering the aligned nodes in the target network. We demonstrate the efficiency and accuracy of Tempo using both real and synthetic data. We compare the running time and the quality of the alignments found by Tempo against those of three existing alignment algorithms, IsoRank [12], MAGNA++ [14] and GHOST [15]. Note that all these networks are tailored towards optimizing alignment at a single time point. To have a fair comparison, we allow each of these methods to consider each time point independently then apply the resulting alignments to all other time points and took the average. We show Tempo has competitive running time and generates significantly better alignments. We use a human brain aging [16] dataset, and integrate this dataset to analyze three phenotypes—two age related diseases (Alzheimer’s and Huntington’s) and one disease that is less prone to aging (Type II diabetes). We perform gene ontology analysis on the aligned genes reported by our algorithm and compared algorithms. Our algorithm could successfully align genes of the phenotype query (i.e. the underlying disease) to strongly related genes in the target network despite their evolving topologies unlike other algorithms. Consequently, we could predict disease-related genes based on the generated alignment using tempo which suggests that Tempo generates alignments that reflect the evolution of nodes topologies through time as well as their homological similarities while other methods only focuses on static and independent topologies. Lastly, we observe that alignments of age related phenotype is significantly higher than alignment of non age phenotype which reflects their high evolution rates and shows that Tempo could identify between different queries.

Related Work

One of the key studies on pairwise global network alignment is IsoRank [12], which is based on the conjecture that two nodes should be matched if their respective neighbors can also be matched. It formulates the alignment as an eigenvalue problem and computes the similarity between pairs of nodes from two given networks as a combination of their homological and topological similarities. It obtains the global alignment of the two given networks using their maximum weight bipartite match with the scores as the weights. The GRAAL (GRAph ALigner) family [17] of global network alignment methods use the graphlet degree similarity to align two networks. Briefly, the graphlet-degree of a node counts the number of graphlets (i.e. induced subgraph) that this node touches, for all graphlets on 2 to 5 nodes. GRAAL [18] first selects a pair of nodes (one from each of the two given networks) with high graphlet degree signature similarity as the seed of the alignment, and greedily expands the alignment by iteratively including new pairs of similar nodes. The SPINAL algorithm [19] iteratively grows the alignment based on apriori computed node similarity score. MAGNA [20] optimizes the edge conservation between two networks using a genetic algorithm. There are several other methods for pairwise network alignment [15, 21–25]. Although the underlying algorithms of these methods vary, the end goal is similar to those discussed above.

Several algorithms address the multiple network alignment [26-28]. IsoRankN [29] extends IsoRank. It adopts spectral clustering on the induced graph of pairwise alignment scores. The algorithm developed by Shih et al. [30] is a seed-expansion heuristic that first selects a set of node pairs with high n class="Chemical">similarity scores using a clustering algorithm, and then expands these pairs by aligning nodes that maximizes the number of the total conserved edges of aligned nodes.

INQ [31] aligns a dynamically evolving query network with one static target network. It uses ColT [32] to find an initial alignment of the initial query, then it observes the differences between the topologies of the already aligned query network and the new query network, and finally, uses these differences to refine the alignment found for the previous query and generate alignment of the current query network. n class="Chemical">DynaMAGNA++ [33] aligns two dynamic networks. It assigns a value to each node based on how the incident edges and graphlets change through dynamic events. It assigns each node a value based on dynamic graphlet degree vector (DGDV) of graphlets up to size four. It considers a pair of nodes from two networks similar if their DGDVs are similar.

Problem Formulation

In this section, we develop a new scoring function, , that integrates the similarities of the aligned nodes and their evolving topolon class="Disease">gies, and includes a penalty for each disconnected component in the aligned subnetworks of the target network at each time point. Next, we introduce the terminology and discuss how we drive our scoring function.

Given a network G=(V,E) and a subset of nodes , we define the induced subnetwork of in G as the nodes in and all incident edges (i.e., ). We denote this induced network as . We say two nodes u and v in G are connected if there exists a path between u and v in G. We say a subset of nodes in G form a connected component if all pairs of nodes in that subset are connected in G. We define a subset of nodes in G as a maximum connected component if the following conditions hold: (i) is a connected component in G, and (ii) there is no node in which is connected to a node in . In the rest of the paper, we use the term “connected component” instead of “maximum connected component”. We denote the number of connected components of a given network G with NCC(G).

Given two temporal networks with t time points, and , we denote the n class="Chemical">similarity between a pair of nodes u∈V1 and v∈V2 at time point i (1≤i≤t) with S(u,v). We use an existing pairwise alignment method to calculate S(u,v). The alignment function ψ maps all nodes in V1 to a subset of the nodes in V2. We denote this subset with Ψ(V1) (i.e. Ψ(V1) = {ψ(u)|∀u∈V1}). We note that ψ yields an induced subnetwork of for each time point i, and each induced subnetwork forms one or more connected components. See Figure 2a for an illustration of this latter point. We denote the number of connected components of the induced subnetwork at time point i as . If the number of connected components at time point i is greater than one then the corresponding induced subnetwork is disconnected. We incur a penalty to account for the missing edges which would connect the disconnected components, and apply this penalty to each disconnected component.

Fig. 2

This figure represents an alignment between two networks G1 and G2. Each node in the query network G1 has a one-to-one mapping with a node in the network G2. The dashed line between two nodes emphasizes that they are mapped to each other. a This represents a hypothetical alignment where a is aligned with b for all 1≤i≤11. The induced subnetwork of the aligned nodes in G2 forms three connected components; C1={b1,b2,b3,b4},C2={b5,b6,b7}, and C3={b8,b9,b10,b11}. Gap nodes are {b12,b13,b14}. b After swapping b11 with b14. This swapping results in two connected components in G2. c After swapping b8 with b14. The aligned nodes in G2 form four connected components

This figure represents an alignment between two networks G1 and G2. Each node in the query network G1 has a one-to-one mapping with a node in the network G2. The dashed line between two nodes emphasizes that they are mapped to each other. a This represents a hypothetical alignment where a is aligned with b for all 1≤i≤11. The induced subnetwork of the aligned nodes in G2 forms three connected components; C1={n class="Gene">b1,b2,b3,b4},C2={b5,b6,b7}, and C3={b8,b9,b10,b11}. Gap nodes are {b12,b13,b14}. b After swapping b11 with b14. This swapping results in two connected components in G2. c After swapping b8 with b14. The aligned nodes in G2 form four connected components

The minimum number of edges needed to join connected components is . We penalize each edge insertion with a constant value denoted with δ, where δ≥S(u, v), ∀u∈V1, v∈V2 and i∈{1,2,…,t}. We define the score of the alignment ψ() at time point i as: We define the temporal network alignment as

Methods

Overview. Our algorithm for solving the temporal network alignment problem has two phases. The first phase finds an initial alignment between the input networks and using the similarity score between pairs of aligned nodes across all time points. We discuss how we calculate the similarity score later in this section. The induced subnetwork of obtained by this alignment may be disconnected since this phase ignores the penalty incurred by edge insertions. The second phase reduces the number of connected components, improving the alignment score. In the second phase, we improve the alignment between the input networks by swapping a subset of the nodes in that are aligned with nodes in with other nodes in . In order to swap a node v∈Ψ(V1) with v∈V2−Ψ(V1), we update the alignment function ψ() to ψ′() such that ∀u one of the two conditions is satisfied: (i) ψ′(u)=v if ψ(u)=v; and (ii) ψ′(u)=ψ(u) if ψ(u)≠v. Figure 2 illustrates this. Here, initially b11 is aligned to a11 (Figure 2a). Swapping b11 with b14 updates the alignment function so that b14 is aligned to a11 (Figure 2b). We observe that this swapping reduces the number of connected components in the induced subnetwork of G2 by one. Notice that if we swap b8 with b14 (instead of b11 with b14) then the number of connected components increases (Figure 2c).

We note that the number of connected components may simultaneously decrease at one time point and increase at other time points when we swap two nodes. We prove that the problem of finding the subset of node swaps that minimizes the number of connected components across all time points is NP-hard. We give a reduction from the Maximum Coverage problem [34] to this problem later in this section.

Algorithm details. Tempo takes two networks ( and ) and the maximum number of allowed swaps (denoted as k) as input. In the following, we explain the two phases of our method in detail.

PHASE I (INITIALIZATION). Here, we construct an initial alignment of and . There exists several algorithms to perform pairwise alignment of two static networks at a single time point. Each of these methods asn class="Chemical">sign similarity scores to all node pairs (one from the first network and one from the second) and then choose the alignment that maximizes the total score of all aligned node pairs. We adopt one of these methods to obtain the similarity scores of each network pairs at each time point i, and use the outputted scores to calculate an initial alignment. We denote the similarity of the node pair (u,v), u∈V1 and v∈V2 generated by such method at the ith time point with S(u,v).

Following, we describe how we adopt pairwise alignment methods to generate similarity scores in temporal networks that are needed to calculate an initial alignment. For that purpose, we conn class="Chemical">sider adopting IsoRank. We note that our choice of such method has no impact on our method. Recall that IsoRank perform pairwise network alignment. Thus, our modifications of IsoRank are meant to adopt it to temporal networks. First, we calculate the homology score between all pairs of nodes (u,v) where u∈V1 and v∈V2 as the similarity score of their sequences using BLAST [35]. We denote the homology score between u and v as H[u,v]. Next, we calculate the topological similarity matrix at the ith time point, denoted as A, as follows. First, we initialize A to be the zero matrix. Next, for u,w∈V1 and v,z∈V2 we let if , where N(v|G) denotes the neighbors of v in network G. Conceptually, A[(u,v),(w,z)] models the topological support that the node pair (u,v) gives to the alignment of their neighboring pair (w,z) at the ith time point. We integrate the homology and the topology scores for and at the ith time point iteratively using a mixing parameter α. We initialize . We then update the similarity between node pairs at iteration r as . We stop this iterative process when = .

We note that in subsequent iterations of the above formulation, the homological similarity of each node pair (w,z) propagates their neighboring pairs (u,v) by a function governed by the topology matrix and the mixing parameter α. We explain three issues arising from these iterations. First, as the number of neighbors of w and z increases, the similarity propagating to each neighbor pair decreases because the number of ways to align nodes w and z without altering the topological similarity grows with increasing number of their neighbors. Secondly, as the value of α decreases, the contribution of the homological similarity to the final similarity value between each node pair grows and the contribution of the topological similarity decreases. In the extreme case when α = 0, the topological similarity has no contribution. Lastly, the iterations above are guaranteed to converge since A is a column stochastic matrix (i.e., the values at each column add up to one). We denote the converged vector at the ith time point with S and call it a score vector. Each entry S[u,v] in this vector shows the similarity (homology and topology combined) between nodes u and v.

We generate an initial alignment ψ0 as follows. We first construct a weighted bipartite network G = as follows: we insert an edge in G between each pair of nodes (u,v) such that u ∈V1 and v ∈V2. We set the weight of the edge (u,v) as the similarity between nodes u and v aggregated over all time points. We denote the similarity as S(u,v) = . The maximum-weight bipartite matching algorithm maps each node in V1 to a node in V2 [36]. This mapping represents the initial alignment, ψ0. We call the nodes in V2 that are not mapped to any node in V1 as gap nodes and denote with F=V2−Ψ(V1).

PHASE II (SELECT k SWAPPING PAIRS). Here, we describe our dynamic programming algorithm that selects a set of k swaps which maximize the alignment score by reducing the number of connected components in the induced alignment across all time points of (see Eq. 1).

We denote a set of r swaps with Δ={(u1,v1),(u2,v2),…,(u,v)} with ∀i≠j, u≠u and v≠v. We denote the alignment after applying the swaps in a given set Δ as ψ. Let us denote the optimal set of r swaps for the alignment ψ with . Also, for a given u∈Ψ(V1), we denote the optimal set of r swaps for the alignment ψ which contains the swap pair (u, v), ∃v∈F, with .

Our algorithm works iteratively. In the first iteration, our algorithm selects one swapping pair for each aligned node u∈Ψ(V1) as

At each subsequent iteration r where 2 ≤r ≤k, for each aligned node u∈Ψ(V1), our algorithm selects a set of r swapping pairs denoted with by adding one swapping pair (u, v), ∃v∈F, to the previously selected r−1 pairs as follows.

Here Θ represents the necessary conditions to include the (u,v) swap pair with a set of r−1 swap pairs as

The first condition above ensures that node u is swapped with a gap node and the second ensures the dynamic programming iterates over all size r−1 swap sets for all aligned nodes of . The third condition ensures that the aligned node u has not already been swapped in the r−1 sized swap set. The final condition is the dual of the previous one, as it ensures that the gap node v has not already been swapped in the r−1 sized swap set. When these conditions hold, the two nodes u and v can be swapped and included into the existing set of r−1 swaps without conflicting with any of the existing swaps.

We report the output of the algorithm at end of the kth iteration as set of k swaps with the highest alignment score using equation

Complexity Analysis. We represent the set cardinalities |V1|, |V2|, and |F| with m, n, l, respectively. The complexity of our algorithm is . We note that . This value is either given as input or we set it to . Next, we provide the derivation of this complexity.

Here we analyze the complexity of our method. Recall that we represent |V1|, |V2|, and |F| with m, n, l respectively. We refer to “Related Work” section as we discuss the phases of our method. For each phase, we explain its complexity. We then summarize the complexity of all phases to denote the overall complexity of our method. These phases are;

(1) PHASE I (CONSTRUCT INITIAL ALIGNMENT). In this phase, we calculate the similarity score between node pairs of the input two networks based on their homology and their topology. First to calculate the topology vector A, we need to trace neighbors of all node pairs which is performed in . Thus, the complexity of calculating the topology score for all time points is . We then integrate the homology and topology score by multiplying the topology and the homology vectors in . The algorithm repeat the previous step, let us say for c times to converge . We select the initilat alignment un class="Chemical">sing the weighted-bipartite matching algorithm in . Thus, in this scenario, the complexity of this phase becomes .

(2) PHASE II (SELECT k SWAPPING PAIRS). This phase is performed in two steps. The first step performs the initialization process of the dynamic programming algorithm, in which we calculate the profit of swapping a gap node f with an aligned node v. In order to to this, we calculate the number of components that f can connect if swapped with v using depth first search through all time points in . The second step performs the iterative process of selecting k swapping pairs where the maximum number of iterations is (k−1). The process combines a gap node f (i.e. 1≤l≤|F|) with a set from swapping pairs from the previous iteration where the maximum number of sets is l. Due to resolving the conflict nodes issue, each combination may trace all profits of all n class="Gene">gab nodes in the current combination. This process is performed in . Thus, the complexity of the second step of phase II is = . Hence, the complexity of phase II is .

In summary, the complexity of our method considering all the previous phases is .

Proof of correctness. Here, we formally proof the correctness of our algorithm. We say that swapping the pair of nodes (u, v) is proper if that the swapping does not increase the number of connected components of the aligned nodes. We first prove that our algorithm will always find a proper swapping node u from the set of aligned node for each gap node v. We first present a lemma which is necessary for the proof of our first theorem. Let us denote the degree of a node v (i.e. number of edges connected to this node) within a component C=(V,E) of the induced subnetwork at time point i with the function deg(v|C).

Lemma 1

Given an undirected subnetwork of , where |V| = z and is acyclic network (has no cycle) within its topology, then .

Proof

Since C is a connected subnetwork with no cycles, the number of edges in C equals z−1 edges. Each edge belongs to an undirected network increases the sum of the network nodes degrees by two. Thus, . _

Lemma 2

Given a gap node v that connects at least two connected components, there exist at least one aligned node u which we can swap with v without increan class="Chemical">sing the number of connected component.

We formally prove this by induction on the size of connected components that u belongs to.

BASE CASE. We consider a component C=(V,E) where |V|=2 and v is connected to C through u, and assume u belongs C. If we swap v with u, then C will contain u and v which corresponds to one component. Thus, the number of connected components of C is still one after swapping.

INDUCTION HYPOTHESIS. We assume there exists a node u for all components of size q nodes that can be swapped without disconnecting its component. We consider two cases of one component C where v is connected to through u. The first case is when C contains at least one cycle with the set of nodes, V={v1,v2,…,v}. It follows that for each node u∈V and u≠u, u can be swapped with v without disconnecting C. In the second case, C represents acyclic network with no cycles. Next, we prove our theorem in this case by contradiction. First, we assume that the number of nodes in C with degree equal to 1 is less than 2. Consequently, , which contradicts Lemma 1. Thus, the number of nodes in C with degree equal to 1 is at least 2 nodes and thus, ∃v,w∈C st. deg(v|C) = 1 and deg(w|C) = 1 and v≠w. Therefore, we can swap v with either v or w. __

Next, we prove that swapping a gap node v with an aligned node u at each iteration will increase the alignment score , showing that the alignment score will always improve by our dynamic programming algorithm.

Theorem 1

Given a value of δ where δ is greater than or equal to S(ψ(u), u) for all u∈V2. At each iteration of our algorithm, monotonically increases.

We assume that our algorithm chooses one pair of nodes to swap; a gap node v and aligned node u which will connect x number of components. We note that the condition x ≥2 must be satisfied for v to be considered for swapping. Also, it follows from Lemma 2 that if we swap v and u then the number of connected components will not increase. Thus, the difference in the score equals D=δ(x−1)−p where p is the difference in pairwise score from swapping (i.e. p = S(ψ(u), u) - S(ψ(u), v)). n class="Chemical">Since δ is greater than or equal to S(u,v)∀u∈V1 and ∈V2, then δ(x−1)≥p. Consequently, D ≥ 0 and will not decrease. __

Proof of NP-hardness. Here, we prove that our problem is NP-hard. To do that, we reduce the Maximum Covern class="Gene">age Problem (MCP), which is known to be NP-hard [37], to our problem. Given a positive integer and a collection of sets, S={S1,S2,…,S}, MCP seeks the subset S such that and the number of covered elements is maximized.

We reduce MCP to an instance of our problem. Let U={x1,x2,…,x} be the union of elements in S (i.e. ). We construct a target temporal network with one time point G2=(V2,E2) as follows. We initialize G2 as V2=∅ and E2=∅. Next, we add a node a in G2 for each element x∈U. Also, for each set S∈S, we add two nodes f and b in V2. Formally, V2={a1,a2,…,a}∪{b1,b2,…,b}∪{f1,f2,…,f}. Next, we populate the set of edges E2. To do that, for all S∈S and x∈S, we insert the edge (f, a) in E2. In addition, for all pair of sets S,S∈S, where in E2. Finally, for a given query network G1=(V1,E1), we construct the set of nodes in G2 aligned to those in G1 as Ψ(V1)={a1,a2,…,a}∪{b1,b2,…,b}. Thus, the set of gap nodes is {f1,f2,…,f}. Notice that, the subnetwork of G2 induced by Ψ(V1) has m+n nodes but it contains no edges as all the edges in G2 are connected to a gap node by our construction. Thus, the alignment yields n+m connected components as each node in Ψ(V1) represents a component.

Recall that each swapping operation swaps an aligned node with a gap node. Also, n class="Disease">recall that the optimization problem we solve for aligning temporal networks aims to find at most k swaps, such that after applying those swaps, the number of connected components NCC(Ψ(V1)∣G2) is minimized (see “Problem Formulation” section). We call this optimization problem minimum Connected Component Problem (mCCP) in the rest of this proof. Next, we prove that MCP is maximized if and only if mCCP is minimized.

First, we prove that if there exists a solution to mCCP, then there exists a solution to MCP. In other words, we prove that minimizing mCCP maximizes MCP. Let us denote the nodes corresponding to the elements in a set S with . In our problem instance, a swap operation swaps f with a node in the set V2−A−{f}. This is because all nodes in A are connected to f, and thus swapping f with a node not in A ensures that all nodes in S∪{f} form one connected component. Therefore, to minimize the number of connected components, we swap f with one of the nodes which is not a part of this connected component. To ensure that, we swap f with a node in the set {b1,b2,…,b}. n class="Chemical">Since all nodes in this set are disconnected, swapping f with any node in this set will yield the same number of connected components. Let us assume that the solution to mCCP performs k swaps. Following from the discussion above, without losing generality, we assume that these swaps are {(f1,b1),(f2,b2),…,(f,b)}. Notice that after these swaps, the nodes in forms one connected component, and all remaining nodes are isolated. Let us denote the number of connected components after these swaps with β. Let us denote the number of nodes in with τ. Notice that τ also reflects the number of elements covered in . We have β=(m−k)+(n−τ)+1.

In the formulation above, the first term (m−k) is the number of nodes b which are not swapped with a gap node. Since all those nodes are isolated, each one forms a connected component by itself. The second term (n−τ) is the number of nodes a which are not included in the set . These nodes remain isolated even after swapping of nodes. The last term (i.e., 1) is the connected component containing the nodes in . After minor algebraic manipulation, we rewrite the equation above as β=(m+n−k+1)−τ. In this equation, the parameters m, n, and k are input to the given mCCP problem, and thus we denote the first term above with the constant c=m+n−k+1. Therefore, we have β=c−τ. In this equality the smaller the value of β is, the larger τ gets. Thus, minimizing the number of connected components β in mCCP maximizes the nuumber of elements covered in MCP.

Second, we prove that if there exists a solution to MCP, then there exists a solution to mCCP. In other words, we prove that maximizing MCP minimizes mCCP. Let us assume that the solution to MCP is . The number of elements covered by this solution is = . By constructing an instance of mCCP as described above, we have swaps denoted with the set . Consequently, after performing these swaps, the nodes in forms one connected component, and all the remaning nodes are isolated. Let us denote the number of connected components with . We have

After minor algebraic manipulation, we rewrite the equation above as Since m, n, and are input parameters, we have , where c is a constant . In this equality, the larger the value of is, the smaller gets. Thus, maximizing in MCP results in maximizing in mCCP.

Lastly, the proof we describe above reduces an instance of MCP to an instance of mCCP in polynomial time and space as it requires only building a network with nodes and edges. Thus, we conclude that the mCCP problem is NP-hard.

Results and Discussion

We evaluate the performance of our algorithm on synthetic and real data. Next, we describe both datasets in detail.

Real Dataset. We obtain our real dataset from two sources. The first one is the human brain an class="Disease">ging dataset [16]. This dataset contains microarray human brain gene expressions profiles obtained from 55 individuals spanning 37 ages from 20 to 99 years. Data from each individual is collected in at least two of the four brain regions namely, the hippocampus, entorhinal cortex, superior-frontal gyrus, and postcentral gyrus. These samples were preferentially selected where tissue was available, thus the number of tissues vary across different individuals. In total, transcription values for 173 samples are collected. Overall, the dataset contains 9426 genes with different expression across ages. The ages in this dataset are not uniformly spaced. In order to bring consecutive time gaps to a more uniform values, we remove two data points which have an age gap of more than 5 years from their successive age values, leading to 35 ages. The samples from each age group were found to be correlated [16]. Thus, to construct different correlated temporal networks from these dataset, we form temporal networks that each has interleaved age groups. We select five temporal networks each having seven time points. Next, we explain how we do that for the first temporal network. We start with the first (i.e., youngest) time point in the aging data. We then skip the next four time points and take the sixth time point in aging data iteratively until we have seven time points. Similarly, for 1sion data to form age-specific PPI networks, we set a cut-off on the gene-expression value. All the interactions that have a lower transcription value for either or both the proteins are removed from the corresponding age-specific network. We use the protein-protein interaction (PPI) network data from BioGRID [38]. For the second source, we select phenotype specific query temporal networks from this dataset. We use two neurodegenrative disorders which are conjectured to be age-related (Alzheimer’s and Huntington’s) and a third one which we expect to be less prone to aging (Type II diabetes). We retrieve the gene sets specific to these three diseases from KEGG database [39]. We form three query PPI temporal networks by keeping only the interactions where both the interactors are from each of the three phenotype-specific (Alzheimer’s, Huntington’s or Type II Diabetes) gene set. We form temporal networks of phenotype disease by taking the intersection of phenotype genes and temporal networks of aging dataset.

Synthetic dataset. We generate synthetic networks to observe the performance of our method under a wide spectrum of parameters classified under two categories; (i) network size and (ii) temporal model parameters, namely number of time points, temporal rate, and cold rate. We vary the target network size to take values from {100, 250, 500, 750, 1000}. We fix the network density to two edges per node on the average (i.e., mean node degree is set to four). We randomly select as a connected subnetwork of . We set the size of the query network to 50 nodes. We generate target network using Barabási-Albert (BA) [40] model as this model produces scale-free networks. In order to explain the parameters in the second category, we describe how we generate the query and target networks and at the first time point. We then explain how we use the parameters in this category to build the query and target networks at the remaining time points.

We generate the subsequent networks for the remaining time points using the three parameters in the second category above as follows. The first parameter is the number of time points t in and . We use 5, 10, 15, and 20 time points in our experiments. n class="Disease">Recall that we select a subnetwork of the target network as the first query network . We mark all nodes and edges in within this subnetwork as cold nodes and edges respectively. We mark all other nodes and edges in as hot. Next, we iteratively generate the networks and at the ith time point (i>1) from and respectively as follows. Let us denote temporal and cold rates (two real numbers) with ε and ε respectively such that 0≤ε≤ε≤1. Let us denote the ratio of cold edges to the total number of edges in the target network with γ. We calculate the hot rate, denoted with ε, from temporal rate and cold rate as ε = (ε−εγ)/(1−γ). Conceptually, hot and cold rates model the rate of evolution of hot and cold edges between two consecutive time points respectively. More specifically, for each subsequent time point i, we generate by randomizing as follows. We iterate over all edges in . For each edge e, if it is a cold edge we remove it with probability ε and insert a new edge between two randomly chosen cold nodes. If e is a hot edge, we remove it with probability ε and insert a new hot edge between two random nodes (with at least one being a hot node). We generate query networks at subsequent time points using almost the same procedure with the only difference being that all edges are cold. We generate datasets by varying ε and ε to take the values {0.05, 0.1, 0.2, 0.4, 0.8} and {0.05, 0.1, 0.2} respectively. For each parameter setting we generate 10 target and query temporal networks.

Recall that, we generate the scoring matrix based on both homology and topology n class="Chemical">similarities. We generate the homology score between two pair of nodes u∈V1 and v∈V2 as follows. If v was originally selected as cold node and u is the same as v, then we generate a homology score between u and v from log-normal distribution [41] with mean 2μ and standard deviation σ. Otherwise, we randomly generate the homology score between u and v from log-normal distribution with mean μ and standard deviation σ. In this way, we allow nodes in query network to be likely to align to nodes in the target network that were originally extracted from. In this paper, we set μ and σ to be 2 and 0.25 respectively. Notice that the homology scores do not change through time points, although topology scores do. Thus, evolution through time points of query and target networks may affect how the query is aligned to the cold region in the target network. We set the edge insertion penalty δ to be .

We compare the accuracy and running time of our algorithm against IsoRank, MAGNA++ and GHOST. IsoRank, n class="Chemical">MAGNA++ and GHOST are designed to align two networks at a single time point. We therefore find the alignment using each of these methods at each time point, impose the alignment to all the time points and report the average. We analyze the biological significance of our results on real data by performing gene ontology analysis and exploring publication evidence. We implemented Tempo in C++, performed all experiments on a computer equipped with AMD FX(tm)-8320 Eight-core Processor 1.4 GHz CPU, 32 GB of RAM running Linux operating system, and used α=0.7 unless otherwise stated.

Evaluation on real data

In this section, we evaluate Tempo on the real data. We first evaluate the significance of alignment score un class="Chemical">sing Tempo. We calculate the z-score by comparing the score of aligned nodes to the score of 1000 randomly selected alignments of the same number of nodes. We compare our results to those of IsoRank. We repeat this experiment for three different disease network queries: Alzheimer’s, Huntigton’s and Type-II diabetes. Figure 3 shows the results. Our results demonstrate that Tempo yields highly significant alignments, and outperforms IsoRank in terms of z-score. We also observe that z-scores of non-age related disease (diabetes) is lower than those of age-related diseases (i.e. Alzheimer and Huntington’s). Although there are some fluctuations in the z-score with growing time gap between query and target networks, we observe that the z-score tends to increase for Alzheimer’s and Huntington’s disease unlike the Type-II diabetes. This suggests that age-related pathways have higher evolution rate than other pathways. Thus, we conjecture that Tempo, which takes all time points into consideration, is suitable for capturing evolving topologies.

Fig. 3

The average z-score of our method using real data of three different diseases; Alzheimer’s, Huntington’s and Type-II diabetes. The x-axis shows which time points was selected to represent the target network. The y-axis shows the z-score of IsoRank (white bars) against our method (black bars)

The average z-score of our method un class="Chemical">sing real data of three different diseases; Alzheimer’s, Huntington’s and Type-II diabetes. The x-axis shows which time points was selected to represent the target network. The y-axis shows the z-score of IsoRank (white bars) against our method (black bars)

Next, we consider the biolon class="Disease">gical significance of our results by identifying aligned gene pairs in which the aligned genes are different, and determining prior evidence that these gene pairs are biologically relevant. We use Tempo to identify 4, 4 and 6 such pairs for Alzheimer’s, Huntington’s and Type-II diabetes, respectively. We note that Alzheimer’s, Huntington’s and Type-II diabetes query sizes are 39, 36, and 23. Thus, the percentages of the different genes found to all the genes in the alignment are 10% to 26%. IsoRank only mapped genes to themselves, suggesting that IsoRank only considers static topologies while our algorithm could map genes based on homological similarities as well as evolving topologies. MAGNA++ and GHOST could only map few genes to themselves while other mapped genes were poorly related.

For each combination of disease and differently mapped gene pairs identified by Tempo, we first search PubMed for publication evidence specific to that disease. For instance, in case of Alzheimer’s disease, the gene n class="Gene">DAB1 that was selected by Tempo and was identified as a potential gene that encode proteins related to functions in biological pathways relevant to the disease [42]. Genes found by Tempo for type II diabetes, for example gene ACTA1, has remarkable change in gene expression value that was observed for the in diabetic samples compared to non-diabetic samples [43]. Moreover, significant up-regulation of GRB2 is observed in transgenic samples compared to controls [44].

Evaluating signaling pathways. In order to determine the biolon class="Disease">gical processes of the aligned genes found by Tempo in gene aging dataset, we perform the gene ontology analysis of the aligned genes in target network using Gene Ontology Consortium [45]. We identify the biological processes or signaling pathways that play significant roles in the disorder. Notice that, while the aligned genes represent many pathways, we only focus on pathways that are related to the underlying query phenotype. We compare all significant pathways later in this section. We calculate how many related pathways found by our method (Tempo) against MAGNA and GHOST and their significance. We also counted the frequency of those pathways when used different range of time points. Table 1 present the results. We find references of certain pathways that are related to specific neurodegenerative disorders (Alzheimer’s and Huntigton’s diseases). For genes we identify when we use Alzheimer’s disease as a query network, we find two pathways, namely Alzheimer disease-amyloid secretase and Alzheimer disease-presenilin are related to Alzheimer’s disease [46]. Various growth factors alter the brain development process at younger age, that manifest as a variety of risk factors at an older age and eventually results in aging-related diseases such as Alzheimer’s and Huntigton’s diseases [47]. For the genes we identify when we use type II diabetes phenotype as a query, we find two pathways that they are commonly associated with type II diabetes [48] namely Insulin/IGF pathway-protein kinase B signaling cascade and Insulin/IGF pathway-mitogen activated protein kinase kinase/MAP kinase cascade. On the other hand, MAGNA or GHOST found at most one pathway with very low significance and did not appear through all tested target networks (see Table 1). In conclusion, studying temporal networks in general and human aging specifically using Tempo enables us to identify age related genes from non age related genes successfully. More importantly, Tempo takes the network alignment problem one huge step forward by moving beyond the classical static network models.

Table 1

Number and significance of functional pathways associated with the underlying disease observed among the aligned genes of target network

Disease	Tempo	MAGNA++	GHOST
Alzheimer	2 / 4 / 2.29E-14	1 / 2 / 2.14E-03	1 / 2 / 3.32E-04
Huntigton’s	1 / 4 / 1.15E-22	0	0
Diabetes	2 / 4 / 2.29E-09	1 / 1 / 2.2E-01	0

Each cell lists the results in the form x/y/z. Here, x represents number of pathways identified, y denotes the number of time points at which these pathways are observed, and z is the statistical significance (p-value) of the least significant of these pathways. The cell with the value 0 implies that no pathways were found

Each cell lists the results in the form x/y/z. Here, x represents number of pathways identified, y denotes the number of time points at which these pathways are observed, and z is the statistical significance (p-value) of the least n class="Chemical">significant of these pathways. The cell with the value 0 implies that no pathways were found

Next, we compare significant pathways which are related to query phenotype to the rest of the pathways of aligned genes un class="Chemical">sing Tempo as well as MAGNA and GHOST. In order to perform this comparison, we present the percentage of genes that contributes to the significant pathways which are related to the query disease. We show the results for Alzheimer disease. Results are similar for the other two queries. Recall that using our algorithm we could find two pathways that are related to Alzheimer disease while MAGNA and GHOST find only one (see Table 1). Figure 4 presents the results. The results demonstrate that the aligned genes result from our method have two pathways that are associated with Alzheimer while MAGNA and GHOST results in only one. In addition our method finds alignments in target network with substantial fraction of genes that contributes to the pathways which are associated with the query disease (15.6% and 18.7% of the genes in the two pathways). On the other hand, resulting alignments of MAGNA and GHOST contributes with a very small fraction to pathways associated with Alzheimer; more precisely 2.2% and 1.64% of the genes for MAGNA and GHOST respectively.

Fig. 4

This figure represents the percentage of genes that contributes to each pathway of the resulting aligned genes in the target network. We point to the significant related pathways of the query disease (Alzheimer). a Tempo b MAGNA c GHOST

This figure represents the percentage of genes that contributes to each pathway of the resulting aligned genes in the target network. We point to the n class="Chemical">significant related pathways of the query disease (Alzheimer). a Tempo b MAGNA c GHOST

Evaluation of recovered query. In this experiment, we evaluate the recovered query region from gene an class="Disease">ging dataset by our algorithm, Tempo, against MAGNA++ and GHOST on real dataset. The recovered region computes the percentage of genes in the query network that were mapped to themselves in the target network despite their evolving topologies. Tables 2, 3, and 4 present the results for Alzheimer’s, Huntington’s, and Type II diabetes respectively. The results show that our algorithm significantly outperform both MAGNA++ and GHOST by aligning similar genes despite their evolving topologies. On the other hand, MAGNA++ and GHOST could poorly align small portion of the query genes to themselves. This suggests that our algorithm could successfully capture the evolving topologies of the genes through time points while other algorithms fail to do so since they consider aligning each time point independently.

Table 2

Percentage of recovered query genes from gene aging dataset when using Alzheimer’s phenotype as query

Target time points	Tempo	MAGNA++	GHOST
First 7	94.87	2.56	0
Second 7	97.43	5.13	0.36
Third 7	97.43	2.56	0
Forth 7	97.43	2.56	0

Table 3

Percentage of recovered query genes from gene aging dataset when using Huntington’s phenotype as query

Target time points	Tempo	MAGNA++	GHOST
First 7	90.9	0.36	0
Second 7	86.36	0	0
Third 7	95.45	0.73	0
Forth 7	95.45	0.73	0

Table 4

Percentage of recovered query genes from gene aging dataset when using Type II diabetes phenotype as query

Target time points	Tempo	MAGNA++	GHOST
First 7	97.22	2.56	0
Second 7	97.22	2.56	0
Third 7	97.22	5.12	0
Forth 7	97.22	2.56	0

Percentage of recovered query genes from gene an class="Disease">ging dataset when using Alzheimer’s phenotype as query

Percentage of recovered query genes from gene an class="Disease">ging dataset when using Huntington’s phenotype as query

Percentage of recovered query genes from gene an class="Disease">ging dataset when using Type II diabetes phenotype as query

Evaluation on synthetic dataset

Evaluation of recovered region. In this experiment, we compare the accuracy of the alignment generated by Tempo against that of IsoRank, n class="Chemical">MAGNA++, and GHOST. We recall that we select the original query network from a subset of nodes and their edges from the target network, and then evolve the query through time points. Here, we evaluate the accuracy by calculating the percentage of the aligned nodes from query network that are paired with the same nodes of the target network that they were originally selected from. We refer to this percentage as recovered region. We illustrate the results in Fig. 5, which demonstrate that Tempo recovers high percentage of the query networks compared to other methods. As the temporal rate increases, the accuracy of Tempo improves dramatically while that of IsoRank remains nearly stagnant and while MAGNA++ and GHOST continue to generate alignments with low recovery rates. Growing the temporal rate while keeping the cold rate unchanged means that the topology of the query network (i.e., cold edges) is evolving slower than the rest of the temporal network (i.e., hot edges). This implies that Tempo can capture the variation in such evolutionary rate while competing alignment strategies which fail to do so.

Fig. 5

The percentage of recovered query in the resulting alignment varying ε and ε to take the values {0.05, 0.1, 0.2, 0.4, 0.8} and {0.05, 0.1, 0.2} respectively. The x-axis shows temporal rate, ε and cold rate, ε (these are the parameters used for constructing synthetic temporal network, with varying evolution rates. The y-axis shows the percentage of recovered query of IsoRank, MAGNA++, and GHOST against Tempo. The error bars show the 80-percentile of the recovered query based on the 10 repetitions of each parameters setting

The percentage of recovered query in the resulting alignment varying ε and ε to take the values {0.05, 0.1, 0.2, 0.4, 0.8} and {0.05, 0.1, 0.2} respectively. The x-axis shows temporal rate, ε and cold rate, ε (these are the parameters used for constructing synthetic temporal network, with varying evolution rates. The y-axis shows the percentn class="Gene">age of recovered query of IsoRank, MAGNA++, and GHOST against Tempo. The error bars show the 80-percentile of the recovered query based on the 10 repetitions of each parameters setting

Evaluation of induced conserved structure. Next, we evaluate the topological quality of the alignment generated by Tempo through comparison with IsoRank, n class="Chemical">MAGNA++, and GHOST. For this purpose, we measure the shared topological structure between and which is preserved under the alignment function ψ through all time points i. Induced conserved structure (ICS) measures the percentage of edges from that are aligned to edges in to the total edges of the induced subnetwork , and is one of the most common measures of topological quality [14]. Formally, = . Figure 6 presents the results, which demonstrate that Tempo generates alignments with high quality based on ICS compared to other algorithms. We note that GHOST was created to optimize ICS, however, Tempo outperforms GHOST on this measure—especially when the temporal rate is high since the performance of GHOST degrades.

Fig. 6

The induced conserved structure (ICS) score of the resulting alignment varying ε and ε to take the values {0.05, 0.1, 0.2, 0.4, 0.8} and {0.05, 0.1, 0.2} respectively. The x-axis shows temporal rate, ε and cold rate, ε. The y-axis shows the ICS score of GHOST, MAGNA++, and IsoRank against our method (Tempo)

Evaluation of edge correctness. In this experiment, we evaluate the topological quality of the alignment generated by our method against IsoRank, n class="Chemical">MAGNA++, and GHOST. For this purpose, we measure the shared topological structure between and which is preserved under the alignment function, ψ through all time points i. Edge correctness (EC) is one of the most common measures of topological quality [14, 15]. It has a similar computations to ICS. Basically, it measures the percentage of edges from that are aligned to edges in to the total edges of smaller network. More specifically, . Figure 7 presents the results. The results demonstrate that our algorithm generates alignments with high quality based on EC compared to other algorithms.

Fig. 7

The Edge correctness (EC) score of the resulting alignment varying ε and ε to take the values {0.05, 0.1, 0.2, 0.4, 0.8} and {0.05, 0.1, 0.2} respectively. The x-axis shows temporal rate, ε and cold rate, ε. The y-axis shows the EC score of GHOST, MAGNA++, and IsoRank against our method (Tempo)

Evaluation of statistical significance of the alignment. We compare the statistical n class="Chemical">significance of the alignments generated by Tempo against that of existing methods. In order to ensure that our experiments do not give any advantage to our algorithm, we use IsoRank to generate initial alignments for Tempo and thus, compare the statistical significance against IsoRank only.

(I) Varying evolution rate. In this experiment, we evaluate the effect of varying the temporal rate (ε) and cold rate (ε) on the significance of the score of the alignments produced by Tempo and that of IsoRank. We generate synthetic networks of sizes {100, 250, 500, 750, 1000} and 20 time points. We fix the network density to two edges per node on the average, and vary ε and ε (ε≤ε) to take the values {0.05, 0.1, 0.2, 0.4, 0.8} and {0.05, 0.1, 0.2}, respectively. Next, we randomly selected 50 nodes from target network 1000 times, and calculate the alignment score of each, i.e., each random selection corresponds to an alignment. We calculate the mean and standard deviation of these 1000 scores and generate the z-score of the alignment generated by Tempo using this mean and standard deviation. Hence, we denote the score generated from our method by S∗, and denote the mean and standard deviation of 1000 scores generated from the random selections with S and σ, respectively. We calculate the z-score of our method as (S∗−S)/σ. We calculate the z-score of the IsoRank method in a similar manner. Figure 8 presents the average z-score values across all target network sizes. The results show that as we increase the temporal rate, the z-score of Tempo significantly increases while the z-score of IsoRank increases by small amount. As the evolution rate increases, the topology of the alignment found by Tempo differs significantly from the topology of rest of the network, and thus, it becomes more challenging to find the correct alignment. However, Tempo continues to generate accurate and significant results especially for large evolution rates unlike IsoRank which considers each single time point independently. We observe the same pattern as we increase cold rate.

Fig. 8

The average z-score of Tempo across network sizes {100, 250, 500, 750, 1000} varying ε and ε to take the values {0.05, 0.1, 0.2, 0.4, 0.8} and {0.05, 0.1, 0.2} respectively. The x-axis shows temporal rate, ε and cold rate, ε. The y-axis shows the z-score of IsoRank (white) against Tempo (black)

The average z-score of Tempo across network n class="Chemical">sizes {100, 250, 500, 750, 1000} varying ε and ε to take the values {0.05, 0.1, 0.2, 0.4, 0.8} and {0.05, 0.1, 0.2} respectively. The x-axis shows temporal rate, ε and cold rate, ε. The y-axis shows the z-score of IsoRank (white) against Tempo (black)

(II) Varying time points. In this experiment, we evaluate how the z-scores of Tempo and IsoRank differ as the input networks evolve and deviate from each other. More specifically, we consider aligning the query network with each of the four target sets we have which have evolving time points (i.e. older n class="Gene">ages) as we move to later target sets. First, we measure the z-score of aligning the query to the first target set (i.e., containing time points 2, 7, 12, …) then we measure the z-score of aligning the query to the second target set (i.e., containing time points 3, 8, 13, …) and so on. We present the average z-score across all temporal and cold rates. Figure 9a presents the results. The results show that Tempo continues to generate alignment with high score significance as we evolve the network. We observe the same pattern for IsoRank, however, Tempo outperforms IsoRank—especially when the time points are distant. This confirms the fact that as the target and query networks evolve and deviate from each other, Tempo is able to take into account the evolution through consecutive time points and generate accurate alignments that persist.

Fig. 9

The average z-score of Tempo (black) against IsoRank (white) a varying target time points, the x-axis shows time point selected, and b varying network size, the x-axis shows network sizes in terms of number of nodes

(III) Varying network size. In this experiment, we compare the n class="Chemical">significance of the alignment generated by Tempo against IsoRank as the target network size increases and the query becomes small with respect to the target. We average the z-score across all evolution rates and vary target network size to take values {100, 250, 500, 750, 1000}. Figure 9b presents the results, which show that the significance of the alignment (best alignment) increases as we increase the size of the underlying target network. We expect this behavior since we compare the aligned nodes (50 nodes) to a random selection of 50 nodes from the underlying target network. Thus, the chance of selecting the best alignment decreases. That said, Tempo was able to identify the accurate alignment which results in high significant values.

Evaluation of running time. In this experiment, we evaluate the running time of our algorithm using synthetic dataset for network sizes as well as number of time points (t). We report the average running time over all values of ε and ε with each parameter combination tested 10 times. We also report the running time for IsoRank, MAGNA++, and GHOST for aligning two networks at a single time point. Figure 10 presents the results. The results demonstrate that Tempo successfully scales to large target networks. The running times of both Tempo and IsoRank grow linearly with increasing target network size and the number of time points (t). We notice that MAGNA++ has similar behavior than IsoRank, while GHOST has an exponential running time. The running time of Tempo is more than that of IsoRank, which is unsurprising since Tempo computes alignment across multiple time points. That said, Tempo has practical running time even for large networks with many time points. More importantly, unlike IsoRank, Tempo considers the network topology at all time points while aligning networks. As we present later in this section, as a natural consequence of the extra effort our method puts to consider all time points, the alignment it finds is significantly more accurate than that of IsoRank which considers only one time point at a time.

Fig. 10

The total running time of IsoRank and Tempo for synthetic networks varying target network size from {100, 250, 500, 750, 1000}, and varying t from 5 to 20. The x-axis shows the input network sizes. The y-axis shows the total running time in seconds

Conclusion

In this paper, we modeled the problem of network alignment between two given temporal networks and proposed a new alignment score function. We developed a novel method to solve this problem by optimizing the alignment score and generating a persist alignment through all time points. Our algorithm incorporates a dynamic programming approach which iteratively refines the alignment to monotonically increase the alignment score. We adapted IsoRank, MAGNA++, and GHOST which are used for pairwise static network alignment, to align two temporal networks by aligning snapshots at each time point independently. We compare the quality and significance of the resulting alignment of both our method and other methods as well as their running time. We observed that the running time of our algorithm is reasonable, compared to other methods, with growing the size of the target network and number of time points, t. In addition, the results showed that our method generates significantly more accurate alignments than that of IsoRank, MAGNA++, and GHOST especially for large evolution rates where finding the correct alignment becomes hard which indicated that our algorithm could capture temporal evolution of the two input networks unlike existing methods. Our experimental results on human aging dataset suggests that age-related pathways (i.e. Alzheimer and Huntington’s) have higher evolution rate than other pathways (i.e. diabetes) and thus, our method could capture such evolving topologies. Furthermore, we performed gene ontology analysis on aligned gene pairs and found that our method could successfully align genes from target network that are similar to genes of the query or significantly related to the underlying query phenotype unlike existing methods which failed to do so.