Rigid Graph Alignment

TL;DR

This paper introduces a novel rigid graph alignment method that considers both network topology and edge lengths, significantly improving alignment quality in applications like brain networks and biomolecules.

Contribution

It formulates the rigid graph alignment problem and proposes a meta-algorithm based on expectation-maximization to enhance alignment accuracy by incorporating structural information.

Findings

Over two-fold improvement in brain network alignment quality

Effective on synthetic graphs with detailed performance analysis

Applicable to various metric space embedded networks

Abstract

Graph databases have been the subject of significant research and development. Problems such as modularity, centrality, alignment, and clustering have been formalized and solved in various application contexts. In this paper, we focus on databases for applications in which graphs have a spatial basis, which we refer to as rigid graphs. Nodes in such graphs have preferred positions relative to their graph neighbors. Examples of such graphs include abstractions of large biomolecules, functional connectomes of the human brain, and mobile device/ sensor communication logs. When analyzing such networks it is important to consider edge lengths; e.g., when identifying conserved patterns through graph alignment, it is important for conserved edges to have correlated lengths, in addition to topological similarity. In contrast to a large body of work on topological graph alignment, rigid graph alignment simultaneously aligns the network, as well as the underlying structure as characterized by edge lengths. We formulate the problem and present a meta-algorithm based on expectation-maximization that alternately aligns the network and the structure. We demonstrate that our meta-algorithm significantly improves the quality of alignments in target applications, compared to topological or structural aligners alone. We apply rigid graph alignment to functional brain networks derived from 20 subjects drawn from the Human Connectome Project (HCP) database, and show over a two-fold increase in quality of alignment over state of the art topological aligners. We evaluate the impact of various parameters associated with input datasets through a study on synthetic graphs, where we fully characterize the performance of our method. Our results are broadly applicable to other applications and abstracted networks that can be embedded in metric spaces -- e.g., through spectral embeddings.

Equations

$${\mathbf{X}}_{ij}=\{\begin{array}{ll}1& \text{if i∈VA is matched with j∈VB under M}\\ 0& \text{otherwise}\end{array}$$

$$\forall i\in [|V_A|],{\sum }_j{\mathbf{X}}_{ij}\le 1,\forall j\in [|V_B|],{\sum }_i{\mathbf{X}}_{ij}\le 1,$$

$$A\bullet \mathbf{X}B{\mathbf{X}}^T={\sum }_{ij}{A}_{ij}(\mathbf{X}B{\mathbf{X}}^T{)}_{ij}=\text{overlap of matching M}$$

$$C_A=C_B\hat{\text{R}}+\text{mathbbm}1\hat{t}$$

$$\mathbf{\Omega }=\left[\begin{array}{cc}\hat{\text{R}}& 0\\ \hat{t}& 1\end{array}\right]$$

$${\mu }_A=\frac{1}{n}\sum _{i=1}^n{C}_{A_i}{\overline{C}}_A=C_A-\text{mathbbm}1{\mu }_A$$

$${\mu }_B=\frac{1}{n}\sum _{i=1}^n{C}_{B_i}{\overline{C}}_B=C_B-\text{mathbbm}1{\mu }_B$$

$$\hat{\text{R}}=V{U}^T,$$

$$\hat{\text{t}}={\mu }_A-\hat{\text{R}}{\mu }_B$$

$${\Vert C_A-\mathbf{X}{C}_B\mathbf{\Omega }\Vert }_F^2$$

$$=\text{tr}[(C_A-\mathbf{X}{C}_B\mathbf{\Omega }{)}^T(C_A-\mathbf{X}{C}_B\mathbf{\Omega })]$$

$$=C_A\bullet C_A+\mathbf{X}{C}_B\mathbf{\Omega }\bullet \mathbf{X}{C}_B\mathbf{\Omega }-2C_A\bullet \mathbf{X}{C}_B\mathbf{\Omega }$$

$$=C_A\bullet C_A+C_B\mathbf{\Omega }\bullet C_B\mathbf{\Omega }-2C_A{\mathbf{\Omega }}^T{C}_B^T\bullet \mathbf{X}$$

$${\text{L}}_{i,j}=\{\begin{array}{ll}1& {\Vert C_{A_i}-C_{B_j}\Vert }_2^2\le ϵ\\ 0& \text{otherwise}\end{array}$$

$${\text{L}}_{i,j}=\{\begin{array}{ll}\text{exp}(-{\Vert C_{A_i}-C_{B_j}\Vert }_2^2)& {\Vert C_{A_i}-C_{B_j}\Vert }_2^2\le ϵ\\ 0& \text{otherwise}\end{array}$$

$${\text{L}}_{i,j}=\{\begin{array}{ll}\text{exp}(-{\Vert C_{A_i}-C_{B_j}\Vert }_2^2)& {\Vert C_{A_i}-C_{B_j}\Vert }_2^2\le d_i^k\\ 0& \text{otherwise}\end{array}$$

$${\mathbb{E}}_o$$

$$=\sum _{i,j,i^{\prime },j^{\prime }}\text{exp}(-{\Vert C_{A_i}-C_{B_{i^{\prime }}}\Vert }_2^2)\text{exp}(-{\Vert C_{A_j}-C_{B_{j^{\prime }}}\Vert }_2^2)A_{ij}{B}_{i^{\prime }{j}^{\prime }}$$

$${\mathbf{X}}^{(t)}=\underset{\mathbf{\Omega }}{\mathrm{argmax}}\mathbf{F}(\mathbf{X},{\mathbf{\Omega }}^{(t)})$$

$${\mathbf{\Omega }}^{(t+1)}=\underset{\mathbf{X}}{\mathrm{argmax}}\mathbf{F}({\mathbf{X}}^{(t)},\mathbf{\Omega })$$

Full text

In this paper, we focus on the problem of rigid graph alignment. In contrast to a large body of work on topological graph alignment, rigid graph alignment simultaneously aligns the network, as well as the underlying structure as characterized by edge lengths. We formulate the problem and present a meta-algorithm based on expectation maximization that alternately aligns the network (using any one of the large class of existing topological aligners) and the structure (using any one of the large class of existing rigid body transformations). We demonstrate that our meta-algorithm significantly improves the quality of alignments in target applications, compared to topological or structural aligners alone. We apply rigid graph alignment to functional brain networks derived from 20 subjects drawn from the Human Connectome Project (HCP) database, and show over a two-fold increase in quality of alignment over state of the art topological aligners. We evaluate impact of various parameters associated with input datasets through a study on synthetic graphs, where we fully characterize the performance of our method. Our results are broadly applicable to other applications and abstracted networks that can be embedded in metric spaces – e.g., through spectral embeddings.

1 Introduction and Motivation

Graph databases are commonly used to store data from a variety of complex systems such as financial transactions, social and communication networks, chemical reaction pathways, and biomolecular interactions. In an important subclass of such systems, the relative positions of nodes are fixed. While the network, or parts thereof, can be transformed from their preferred (relaxed) positions, their “potential energy” increases in the process, motivating the nodes to return to their relaxed state, i.e., the layout with the minimum energy. A natural way to think about such graphs is to visualize the edges as springs of different lengths. We refer to such graphs as rigid graphs. An intuitive example of such a graph is an abstraction of a molecule, with nodes corresponding to atoms, and edges to bonds. Bonds between specific atom pairs have well-defined lengths (e.g., a Hydrogen-Hydrogen bond has a preferred length of 74pm, a Carbon-Carbon bond, 154 pm, etc.). Bonds may be stretched or compressed, however, their potential energy increases (as determined by the Lennard Jones potential) in the process, and the molecule releases this energy in kinetic form to return to its native structure. Note that the precise interaction models (interaction potentials) are considerably more complex – however, this high-level description is illustrative of our concept.

To further motivate this problem and highlight its broad modeling scope, we consider common models for the human brain connectome. Brain activity can be modeled using graphs in which nodes represent regions of the brain that have specialized behavioural and/or structural properties. The coherent firing of pairs of regions is modeled by edges. In this application, the presence or absence of individual edges could vary across a population, whereas the physical positions of nodes are largely conserved in a reasonably constrained database, such as the Healthy Adult Dataset of the Human Connectome Project.

In this paper, we focus on the problem of aligning a pair of given rigid graphs. Informally stated, the goal is to find node/ edge correspondences across a pair of rigid graphs; e.g., across two given biomolecules or two brain connectomes. Traditional approaches to related problems use either structural embeddings that consider edge lengths (i.e., identifying structural transformations – translation, rotation, dilation, that maximize overlap of graphs) or connectivity (i.e., identifying node label permutations that maximize overlap in adjacency matrices) to find correspondences between the two graphs, but not both. However, in the case of rigid graphs, both of these sources of information are rich and complementary. We demonstrate that integrating these two sources into a single framework leads to markedly better correspondences. In this paper, we formulate the problem of rigid graph alignment, and propose a solution technique for this problem in which the rigidity of nodes as well as their connectivity are used to infer significantly better alignments that either structural or topological alignment alone.

We motivate rigid graph alignment in two distinct application contexts. The first one uses a growing and important database of 3D brain images from the Human Connectome Project. One of the major goals of functional brain studies is to understand spontaneous firings of neurons in absence of any stimulus – the “resting-state” activity. Resting state Magnetic Resonance Image (MRI) is useful in creation of a rough baseline network. A more interesting reason to study these images is that they record cognitive processes of individuals in absence of any activity. It is hypothesized that some of these cognitive processes are unique to individuals [16]. Thus, in principle, it must be possible to uniquely identify an individual from a population, using just their resting state functional brain image. This can be done by aligning the images structurally (through registration), topologically (through network alignment – Figure 1), or, as we demonstrate in this paper, using both structural and topological information, yielding superior results.

Our rigid graph alignment model incorporates both anatomical (structural) and functional (connectomic) information by iteratively improving on the quality of network match using a (topological) network aligner, and quality of structural match using rigid body transformations. This ensures that we: (i) do not suppress signals unique to an individual; and (ii) obtain anatomically consistent alignments. Owing to structural and functional uniqueness of individuals, we hypothesize and validate that the quality of alignment between two networks of the same individual (taken across two distinct imaging sessions) should be higher than that between two different individuals, and that “rigid graph alignment” is a critical methodological basis for such alignments.

We further motivate rigid graph alignment in the context of aligning large biomolecules drawn from existing databases of known drug molecules [4]. This problem is of significance, since it underlies biomolecular interactions, identification of drug targets, and design of drug molecules. Specifically, binding pockets in molecules correspond to local atomic structures (active sites) – two molecules that share active sites are likely to bind with the same ligand. This process of matching active sites can be formulated as a rigid graph alignment problem. This is an active area of research (please see Section 4), with commonly used methods relying on traditional graph or structural alignment, bootstrapped using techniques such as geometric hashing. These methods do not utilize constraints imposed by the bond lengths and bond angles into the graph alignment process, or the global structure of the molecules. Our novel formulation of rigid graph alignment factors the bond structure through rigid graph alignment, while simultaneously accounting for bond characteristics, using rigid body transformations.

Summary of Contributions The paper makes the following contributions: (i) it mathematically formulates the rigid graph alignment problem, and puts the problem in the context of existing graph alignment and structural alignment problem formulations; (ii) it proposes a solution to the rigid graph alignment problem based on expectation maximization; (iii) it comprehensively evaluates the performance of the proposed solution on real and synthetic graphs and demonstrates significant performance improvement over state-of-the-art aligners; and (iv) it characterizes the robustness of the solution and the impact of various problem parameters through a comprehensive study on synthetic graphs.

We note that our method does not compete with existing graph alignment and/or structural alignment methods. Rather, it leverages the distinct insights provided by both classes of methods to iteratively improve the quality of alignment. Our meta-algorithm can be realized with any pair of network and structural aligners to improve on the performance of either class of methods. The focus of this paper is on suitable formulation of the problem, derivation of the meta-algorithm, and demonstration of its superiority through detailed experiments for a specific choice of structural and network aligner.

2 Problem Formulation

2.1 Problem definition

We define the rigid graph alignment problem by first reviewing existing graph and structure alignment formulations, and use these to motivate our new problem.

Network alignment review. The literature on network alignment is vast – precluding a comprehensive review. As a representative of a broad class of techniques, we discuss the quadratic programming formulation of a network alignment objective [2]. Let $A=G(V_{A},E_{A})$ and $B=G(V_{B},E_{B})$ be two undirected graphs that we aim to align. For simplicity, we set $|V_{A}|=|V_{B}|=n$. Let $\mathbf{L}\in\mathbb{R}^{n\times n}$ be a matrix that denotes prior knowledge on likelihood of alignment (or known similarity) between vertices $V_{A}$ and $V_{B}$. If $\mathbf{L}_{ij}>0$, prior knowledge implies that node $i$ in $V_{A}$ potentially aligns with node $j$ in $V_{B}$, otherwise, node $i$ in $V_{A}$ cannot align with node $j$ in $V_{B}$. The aim of the alignment problem is to find a matching $M$ from $V_{A}$ to $V_{B}$ using only weights of $\mathbf{L}$ and adjacency information. We define $\mathbf{X}$ to be a matrix with the same dimensions as $L$ to encode matching $\mathcal{M}$:

[TABLE]

Under such a mapping, we say that an edge $(u,v)\in E_{A}$ overlaps with an edge $(u^{\prime},v^{\prime})\in E_{B}$ if $(u,u^{\prime})$ and $(v,v^{\prime})$ are in $\mathcal{M}$.

The goal in this graph alignment formulation is to find a matching $\mathcal{M}$ that maximizes a linear combination of overlap and matching weight. Note that, if $\mathbf{X}$ is constrained by a set of matching constraints,

[TABLE]

then

[TABLE]

where we use the adjacency matrices $A$ and $B$ for the graphs and $\bullet$ to represent a matrix inner product. The corresponding quadratic program is given by: {maxi}—s— XαL ∙X + βA ∙X B X^T \addConstraint∑i Xij ≤1 ∀j=1 …—V_B— \addConstraint∑j Xij ≤1 ∀i=1 …—V_A—,  X_ij ∈{ 0, 1}

Here, $\alpha,\beta$ are non-negative constants that allow for tradeoff between matching weights from the prior and the number of overlapping edges.

Review of structural alignment. We now summarize the Orthogonal Procrustes Problem, which aims to compute structural alignments. The problem, originally solved by Schönemann [35] aims to find an orthogonal transformation that reduces the distance between two matrices in the Frobenius norm. It is often used to find rigid body transformations that describe relationships between two objects. Formally, for any two matrices $Y$ and $Z$, the problem minimizes: {mini}—s— Ω——Y - Z Ω ——_F^2 \addConstraint Ω^T Ω = I

In our setting, the matrices $Y$ and $Z$ correspond to node coordinates $C_{A},C_{B}\in\mathbb{R}^{n\times d}$. Our objective is to measure the orthogonal transformation between the sets of corresponding nodes of the two graphs. In other words, we aim to find rotation $\hat{\textbf{R}}\in\mathbb{R}^{d\times d}$ and translation $\hat{t}\in\mathbb{R}^{1\times d}$ such that

[TABLE]

where, $\hat{\textbf{R}}^{T}\hat{\textbf{R}}=I$ and $\mathbbm{1}$ is the ones vector of length $n$. For convenience, we pad $C_{A}$ and $C_{B}$ by a vector of ones (for convenience in notation, we don’t explicitly show this detail), and combine the translation vector and the rotation matrix into $\bf{\Omega}$ as

[TABLE]

We rephrase the problem as: {mini}—s— Ω——C_A - C_B Ω ——_F^2

The algorithm due to Kabsch [21] was among the first solutions aimed at solving the problem for 2D and 3D coordinates. We use a related SVD-based method, due to Sabata et al. [34], since it is shown to stable by Eggert et al. [11].

[TABLE]

[TABLE]

Here, $C_{A_{i}}$ and $C_{B_{i}}$ refer to the coordinates of the $i$-th nodes of graphs $A$ and $B$. Hence, $\mu_{A},\mu_{B}\in\mathbb{R}^{n\times 1}$. Define $\textbf{H}=\overline{C}_{A}^{T}\overline{C}_{B}$, then the estimated rotation matrix is $\hat{\bf{R}}$ is given by:

[TABLE]

where $V$ and $U$ are orthonormal matrices that are obtained from SVD of H.

The optimal translation $\hat{\textbf{t}}$ is defined by:

[TABLE]

We note that we can estimate scaling between the two graphs using more generalized formulations of the orthogonal Procrustes problems. Instead, we assume that the coordinates are drawn from the same units (say meters). This assumption also circumvents the problem of having potentially different coordinate systems.

Rigid Graph Alignment. We now define our problem of rigid graph alignment. Let $C_{A}\in\mathbb{R}^{n\times d}$ and $C_{B}\in\mathbb{R}^{n\times d}$ represent coordinates of vertices from $V_{A}$ and $V_{B}$, respectively. Combining the network alignment objective function and rigidity metric, as described in Equations 2 and 2.1, our objective function for rigid graph alignment can be written as: {maxi}—s— X, ΩαL ∙X + βA ∙X B X^T - γ——C_A - X C_B Ω ——_F^2F = \addConstraint∑i Xij ≤1 ∀j=1 …—V_B— \addConstraint∑j Xij ≤1 ∀i=1 …—V_A—,  X_ij ∈{ 0, 1}

Here, we write $\mathbf{X}C_{B}$ to denote the conformally permuted set of coordinates of graph $B$ after alignment. The first term in this equation, $\mathbf{L}\bullet\mathbf{X}$, corresponds to the consistency between the prior $\mathbf{L}$ and the mapping of vertices across the two graphs $\mathbf{X}$. The second term, $A\bullet\mathbf{X}B\mathbf{X}^{T}$, corresponds to the alignment of the two networks, and the third term, $||C_{A}-\mathbf{X}C_{B}\mathbf{\Omega}||_{F}^{2}$, to the structural (mis)alignment. In rigid graph alignment, the prior $L$ is driven by spatial constraints. The weights $\alpha$, $\beta$ and $\gamma$ are parameters for the user to adjust relative importance of the prior, graph matching, and structural alignment. It can be seen that when $\alpha=1,\beta=0,\gamma=0$, the problem reduces to the maximum matching problem; when $\alpha=0,\beta=1,\gamma=0$, it is the solution to the problem of maximizing overlap; and when $\alpha=0,\beta=0,\gamma=1$, the problem reduces to rigid body transformation.

The structural alignment error can be minimized by solving for $\mathbf{\Omega}$ using generalized Procrustes method on the ordered set of vertices. To find a good network alignment solution, we need to find the optimal $\mathbf{X}$. For ease of analysis, we view $\mathbf{X}$ as a permutation matrix. Then, the structural error term can be expressed as

[TABLE]

For a given transformation matrix $\mathbf{\Omega}$, using Equation 10, we can rewrite the objective function as

{maxi}

—s— Xα⏟L ∙X_Update Prior + β⏞A ∙X B X^T^Network Alignment + γ⏟C_A Ω^T C_B^T ∙X_Structural AlignmentF = \addConstraint∑i Xij = 1 ∀j=1 …—V_B—\addConstraint∑j Xij = 1 ∀i=1 …—V_A—,  X_ij ∈{ 0, 1}

This form of the objective function suggests that the optimal $\mathbf{X}$ is one that maximizes the network alignment and structural alignment, while respecting the prior. Furthermore, it suggests that the prior $\mathbf{L}$ should be proportional to the similarity between the coordinates, i.e., $C_{A}\mathbf{\Omega}^{T}C_{B}^{T}$. However, this definition poses the potential problem of a dense prior, which can significantly increase the runtime of network alignment. To circumvent this issue, we constrain the number of non-zeros in the prior using a distance measure or number of neighbours. One such simple constraint is given by:

[TABLE]

where $\epsilon$ denotes a small neighborhood. Alternatively, we can set the $ij$-th entry of $\mathbf{L}$ to be inversely proportional to the distance between the $i$-th node of $A$ and the $j$-th node of $B$.

[TABLE]

Yet another option is to match node $i$ of graph $A$ with one of its $k$-nearest neighbours. The distance, $d^{k}_{i}$, of $i$ and the $k$-th closest neighbor is given by:

[TABLE]

The intuition for an rigid graph alignment technique is as follows: let $P_{ii^{\prime}}$ denote the likelihood of matching node $i$ of graph $A$ with $i^{\prime}$ of graph $B$. Since we use structural conformity as prior information to guide the network alignment, we say that $P_{ii^{\prime}}\propto\exp(-\left\lVert C_{A_{i}}-C_{B_{i^{\prime}}}\right\rVert_{2}^{2})$. We can similarly define $P_{jj^{\prime}}$. Assuming that the probabilities are independent, the expected edge overlap can be written as the product:

[TABLE]

We can see that the expected edge overlap increases as matched nodes move closer to each other. Hence, increased structural alignment results in increased edge overlap. In an analogous manner, accurate matches lead to more accurate transformation matrices, which minimizes the residual error after structural alignment. Premised on the observation that the two objectives depend on, and reinforce each other, we propose an algorithm that optimizes the two terms, alternately as follows:

[TABLE]

[TABLE]

This formulation directly motivates our rigid graph matching algorithm in the next section.

2.2 Rigid Graph Matching Algorithm

Our approach to rigid graph alignment splits the problem into two tasks and iterates to convergence: (i) align the graphs restricting the prior to pairs of nodes between graphs that are in spatial proximity (using definitions for L mentioned earlier) – i.e., maximizing the second term of the objective function, and (ii) align the coordinates using the current estimate of the alignment; i.e., maximizing the third term of the objective function. This requires an initial alignment to begin, which we discuss at the end of the section.

The vast majority of network alignment methods take as input $\alpha,\beta,L,A,B$ (or equivalent inputs). Hence, our goal in step (i) is to leverage the coordinates to estimate a matrix $L$ that constrains the set of alignments considered by the network alignment component. We devise a routine, get_prior that creates the matrix L based on the distance between vertices (using the current transformation $\Omega$). Our implementation of this routine assigns every node $i$ in $A$ with $k$ of its nearest neighbours as possible matches. Each candidate node is assigned a weight inversely proportional to the distance between that node and node $i$ (Equation 13). In terms of the objective function, note that this choice of $L$ ensures that the term $||C_{A}-\mathbf{X}C_{B}\mathbf{\Omega}||_{F}^{2}$ stays approximately the same after we use a network alignment method to optimize the permutation.

Step (ii) involves the use of a procedure for structural alignment. We used an SVD-based method to compute orthogonal transformation for structural alignment, as described in Section 2.1.

To recap, in step (i), we maximize the term $\alpha\mathbf{L}\bullet\mathbf{X}+\beta A\bullet\mathbf{X}B\mathbf{X}^{T}$ with the term $\gamma||C_{A}-\mathbf{X}C_{B}\mathbf{\Omega}||_{F}^{2}$ approximately fixed, whereas in step (ii), we maximize $-\gamma||C_{A}-\mathbf{X}C_{B}\mathbf{\Omega}||_{F}^{2}$ with the graph alignment term fixed.

The resulting rigid graph alignment procedure is given in Algorithm 1 111Code will be made available with final version of the paper.

The process is continued until convergence is achieved, in terms of preferred metrics such as edge/node overlap or in the solution to the Orthogonal Procrustes problem. The algorithm proposed (really a meta-algorithm) works with any pair of aligners (graph and structural) that use a prior. In our experiments, we present results where we used netalignmbp [2] as the graph aligner and the SVD-based structural aligner.

In the first iteration, when we have no prior knowledge of any correspondence between nodes, we use the following approach to bootstrap algorithm 1: we compute average pairwise distance for all nodes and populate the prior matrix with nodes that have similar distance profiles. We draw histograms for the distances of each node in both graphs. Then, we use correlation to measure the similarity of distance profiles of pairs of nodes across the two graphs. The prior matrix $\mathbf{L}$ is populated on the basis of this measure. We note that distance profiles are invariant to rotations and translations. Furthermore, correlation measures such as Pearson Correlation are invariant to scaling. In practice, we find that this heuristic works well both in the synthetic datasets, as well as the HCP dataset.

3 Experimental Evaluation

We present results from two sets of experiments – the first set of experiments is performed on a real dataset from the Human Connectome Project. These experiments demonstrate the superiority of our method on this important dataset, while validating a number of application hypotheses. The second set of experiments is performed on synthetic datasets. These experiments are used to assess the stability of our proposed method, comparison to state of the art techniques, and demonstrating runtime characteristics.

3.1 Analyzing the Human Functional Connectome

We present experimental results in the context of an important application in the analyses of brain connectomes. We demonstrate that our methods are capable of significant improvement in alignment quality over state of the art graph alignment techniques that do not consider structural rigidity of the graph.

3.1.1 Dataset

We construct human brain connectomes from the Young Adult Database. The Young Adult Database was created as part of the Washington University-University of Minnesota project of the Human Connectome Project (HCP) consortium, described in vanEssen et al. [13]. It contains both structural and behavioural data collected from over 1100 subjects, both male and female, between the ages of 21 and 35. The images include 3T and 7T structural Magnetic Resonance Imaging (MRI), diffusion MRI, Magnetoencephalography (MEG), Electroencephalography (EEG), functional MRI – at resting state and during tasks. A more detailed description of the dataset and acquisition protocol may be found in vanEssen et al. [14] The resting state functional MRI data of each subject is collected in two sessions, which are separated by a few days. Each session lasted 30 minutes: 15 minutes of left-to-right (LR) encoding, followed by 15 minutes of right-to-left (RL) encoding. Each voxel is isotropic, with dimensions $2mm\times 2mm\times 2mm$, and images were acquired once every 720 ms, as explained in Smith et al. [39]. In this paper, we use only LR scans of resting state functional MRI of 20 subjects.

3.1.2 Preprocessing Steps

Functional MRI data has four dimensions – three spatial dimensions, and one temporal dimension. In essence, it is a collection of time series data – one per voxel of the brain. The preprocessing steps remove spatial and temporal artifacts from the images. Head motion of subjects during acquisition is inevitable, and a potential source of errors. We use the the motion correction tool from FMRI Standard Library (FSL), called Motion Correction FSL’s Linear Image Registration Tool (MCFLIRT) [20]), and register each (volumetric) image in a session to the first time slice. Motion correction ensures consistent labeling of voxels between images in a session. This is followed by skull stripping, using FSL’s Brain Extraction Tool (BET) [38]. The images are then resampled to voxels of size $4mm$. This resampling is done to create networks whose dimensions are accessible to network alignment using state of the art methods. Non-brain voxels are masked out, along with voxels of low variance. The remaining, relevant voxels are vectorized, to create a $voxel\times time$ matrix. It has been shown that spontaneous firings observed in resting state functional MRI is best captured in low frequency fluctuations $<$ 0.1 Hz (Murphy et al. [29]), which is why we use a bandpass filter with limits of 0.001 Hz to 0.08 Hz.

A $voxel\times voxel$ similarity matrix is created by computing Pearson correlation between each pair of time series. We retain the top 5 percentile of correlation values and construct networks with voxels as nodes and high correlation valued edges between nodes. Prior to the creation of the similarity matrix, we regress out the global signal. Since we only use strongly positive values in subsequent steps we sidestep the issue of artificially induced negative correlations.

3.1.3 Rigid Graph Alignment Yields Higher Edge Overlap

In our first set of results, we show that Rigid Graph Alignment improves on a commonly used network alignment metric – edge overlap. For two adjacency matrices $A$ and $B$, edge overlap is defined as $A\bullet\mathbf{X}B\mathbf{X}^{T}$. We run Rigid Graph Alignment (Algorithm 1) for intra-subject (same subject across two sessions) and inter-subject (across subjects) analysis on 20 subjects, with equal weights given to the prior ($\alpha$), edge overlap ($\beta$), and structural alignment ($\gamma$), as discussed in Section 2.1 for a maximum of 20 iterations and a convergence threshold of $0.1\%$ in edge overlap. We find that the edge overlap at the end of first iteration is $20.18\pm 4.2\%$, whereas the edge overlap after rigid network alignment is $53.05\pm 12.5\%$. The improvement of edge overlap is due to increasingly accurate priors in the later iterations of the algorithm. The increase per iteration in typical intra-subject alignment is shown in Figure 3. We observe that the scores are largely stagnant in the first few iterations, which is attributed to the fact that the initial prior is not as informative as priors in subsequent iterations. As stated before, to bootstrap the algorithm, we use similarity of distance profiles to populate the initial prior matrix because the “drawings” of graphs could look very different. In the later iterations, successive transformations result in corresponding nodes of the two graphs being placed close together.

We also characterize the statistical significance of our performance improvement. To do this, we transform the coordinates by random transformation matrices (i.e., we randomly reorient the brains before alignment) in 100 trials. In each case, we found that the edge overlap by rigid graph alignment was higher than edge overlap by regular graph alignment. This implies that our performance improvements are significant and robust.

3.1.4 Residual Error in Structural Transformation as a Metric for Network Alignment

The problem of identifiability of the connectome, or the so-called brain fingerprint involves finding patterns in brain images (either structural or functional) that are unique to an individual. In principle, two networks belonging to the same subject should align “better” than two networks belonging to different subjects. In our context, we need to devise a scoring measure for alignments that can suitably distinguish intra-subject and inter-subject alignment.

Figure 4 shows the histogram of edge overlap at the end of the first iteration. The edge overlap for intra-subject networks is $20.73\pm 4.45\%$, whereas for inter-subject networks, it is $19.62\pm 4.28\%$. This suggests that edge overlap is not a strong discriminatory measure for the brain signature. In fact, the improvement in edge overlap as a consequence of our rigid graph alignment algorithm is observed both for inter-subject and intra-subject alignments, as can be seen in Figure 6. In such applications, where the position of the vertices are of importance, we show that the residual error in structural transformation, which we call “rigidity metric” (Equation 2.1) is a better indicator of quality of alignment.

Figure 5 shows the histogram of quality of inter-subject and intra-subject alignments. For intra-subject alignments, the rigidity metric, normalized by the number of vertices was found to be $1.52\pm 0.01\%$, whereas for inter-subject alignments, it was $1.99\pm 0.025\%$. This demonstrates significantly higher distinguishability from the rigidity metric.

3.2 Experiments on Synthetic Datasets

We test the performance of our formulation and method on synthetically generated preferential attachment and Erdos-Renyi ($G(n,p)$) graphs to assess the impact of various graph parameters.

Generation of Synthetic Graphs. Generation of synthetic rigid graphs involves two steps – spatial location of nodes and assignment of edges between them. To locate nodes, we first create a grid that is evenly spaced in all dimensions. Each grid point has a $d$-dimensional coordinate associated with it. A grid point is assigned a node on the basis of the outcome of an independent, (biased) coin toss.

To assign edges, we consider two strategies: Preferential Attachment [1] and Erdos Renyi ($G(n,p)$) models [17]. In Preferential Attachment or the Barabási-Albert model, we start with an initial graph of $n_{0}$. Then, an incoming node attaches itself to a subset of existing nodes with probability proportional to their degrees, i.e., $p_{i}=k_{i}/\sum_{j}k_{j}$, where $p_{i}$ is the probability of attaching to node $i$, which has degree $k_{i}$. In the $G(n,p)$ model, an edge exists between any pair of nodes $i$ and $j$ with a probability $p$, independent of other edges.

Perturbation Schemes For network alignment, we generate pairs of networks. The first graph is generated in accordance with the aforementioned procedure. The second network is created by adding noise both with respect to the spatial location of the nodes, as well connectivity profiles of the graphs. To add spatial noise, the coordinates are rotated along each axis by $\theta^{\circ}$. Then, each coordinate is translated by a random vector $t\in\mathbb{R}^{d}$. Following this, the coordinates of each node are independently perturbed by $C_{n}\in\mathbb{R}^{n\times d}$. Finally, edges are deleted independently with a probability $p_{d}$ and added with a probability $p_{a}$. This process yields two networks, with parametrizable differences, thereby allowing us to characterize the behaviour of both network alignment and structural alignment.

We quantify performance in terms of correctly identified edges, i.e., edge overlap (equation 2) and node overlap (or node correctness) – the fraction of correctly paired nodes, which corresponds to the precision for our problem. Note that this measure is applicable primarily to synthetic experiments, where the correct solution is known.

3.2.1 Impact of perturbation to the adjacency matrix

In this experiment, we add and delete edges, while holding the relative position of nodes. We allow the entire graph to be rotated along all axes. We set $p_{d}=p_{a}$ for these experiments. We repeat the experiment for $n=\{500,1000\}$, $d=3$, $\theta=\{1^{\circ},5^{\circ},30^{\circ},60^{\circ},90^{\circ},180^{\circ}\}$ for both preferential attachment and $G(n,p)$ graphs.

Figures 7 and 8 show the results of the experiment for $n=1000$, $d=3$, $\theta=60^{\circ}$ for preferential attachment graphs. The performance of both rigid graph alignment and a regular network alignment with respect to edge overlaps is comparable for low noise; however, rigid graph alignment is robust to higher levels of noise (Figure 7). In fact, the edge overlap in rigid graph alignment closely follows the pattern of true edge alignment. The drop in edge overlap for regular graph alignment corresponds to a drop in node alignment, as seen in Figure 8. The node overlap of rigid graph alignment is, on average, $97.67\pm 1.21\%$, whereas the same for regular network alignment is $65.85\pm 3.90\%$, across different rotations. In the absence of node noise, the spatial relationships between nodes are maintained, thereby leaving edge lengths of overlapped edges unchanged. This means that a correct rigid body transformation leaves the two graphs structurally identical, which explains the high values of node overlap when rigidity of nodes is accounted for.

The results of the same experiment on $G(n,p)$ graphs are presented in Figures 9 and 10. Similar patterns of behavior are observed in the case of rigid graph alignment. However, the absence of patterns in edges in $G(n,p)$, coupled with the fact that regular graph aligners do not leverage structural data, causes the edge and node overlap to be substantially lower.

3.2.2 Impact of perturbation by node movement

In this experiment, we perturb the physical position of nodes, while fixing the edge noise to zero, i.e., $p_{a}=p_{d}=0$. The nodes are perturbed along each coordinate by independently sampled, scaled random matrix $C_{n}$. The scaling is done so as to represent the perturbation as a proportion of the physical size of the network. The experiment was done for $n=\{500,1000\}$, $d=3$, $\theta=\{1^{\circ},5^{\circ},30^{\circ},60^{\circ},90^{\circ},180^{\circ}\}$ for both preferential attachment and $G(n,p)$ graphs.

The results for $n=1000$, $d=3$, $\theta=60^{\circ}$ in preferential attachment graphs is shown in Figures 11 and 12. In rigid graph alignment, the node overlap is $80.35\pm 1.21\%$, whereas the same for regular network alignment is $35.85\pm 3.90\%$. The decrease in node overlap for high node noise is due to the fact that large perturbations to individual nodes decrease the efficacy of the correction made by rigid body transformations. However, the transformations still make the graphs structurally similar, which explains the significant difference in node overlap between the two methods.

Figures 13 and 14 show results for the same experiment on $G(n,p)$ graphs. These figures show that an absence of patterns in edges of a $G(n,p)$ graph, coupled with poor priors due to perturbations in node positions make both algorithms more sensitive to noise.

3.2.3 Impact of perturbation of both nodes and edges

The effect on edge overlaps in preferential attachment networks of perturbing both nodes and edges is shown in Figures 15(b) and 15(a). In rigid graph alignment (Figure 15(b)), edge overlaps are more stable with moderate amounts of both edge and node noise. Edge overlap values track true edge overlap. Edge overlap in the case of regular network alignment, however shows a sharp decline for small amounts of noise. Note that the scales on the two figures are made identical to facilitate comparison.

The effect on node overlaps in the same setting is shown in Figures 16(b) and 16(a). Node alignment in rigid graph alignment is determined primarily by the amount of noise in the location of the nodes. However, regular alignment is far more sensitive to either noise.

3.2.4 Choice of graph aligners

In this section, we compare the performance of the previously used belief propagation algorithm (Figures 15(b) and 16(b)) with the method due to Klau et al. [23] in Figure 17 and Isorank [36] in Figure 18 on preferential attachment networks. As before, we vary edge and node noise and record the edge and node overlap in each case.

We see that the algorithm due to Klau et al. [23] is comparable with belief propagation in terms of edge overlap, however it does slightly worse in terms of node overlap when the edge noise is high ($>20\%$). Similarly, Isorank performs quite well in terms of edge overlap, however the node overlap degrades sharply with edge noise. Based on these observations, we conclude that our meta-algorithm can indeed be used with a wide variety readily available network aligners.

3.3 Runtime Considerations

Our rigid graph alignment method is an iterative process, where a network aligner is called in every iteration. As such, it can be used in conjunction with any aligner that accepts a prior. Hence, the runtime depends primarily on the base alignment technique, and is linear in the number of EM steps. Our choice of the network alignment algorithm – netalignmbp [2] takes $\mathcal{O}(nnz(A\otimes B)+|E_{L}|+matching)$, where $E_{L}$ is the number of edges in the prior $L$. Here, $\mathcal{O}(matching)$ is the complexity of bipartite matching (Khan et al.[22]), which depends on the matching algorithm used. In typical experiments, as we have shown, a small number of EM steps yield significant improvement in solution quality.

4 Related Literature

Network alignment has been an active area of research over several decades [12]. Early formulations of network alignment, in the form of exact graph matching (subgraph isomorphisms) in small chemical networks were analyzed by Sussenguth [40]. The notion of distance between graphs (or their similarity) was first quantified by Zelinka [43, 44] on the basis of isomorphisms, assuming that the graphs have the same number of nodes. An analytical solution to approximate weighted graph matching problem was explored by Umeyama [41]. This method used eigen-decompositions of weighted adjacency matrices to find qnear optimum matches using the Hungarian Method.

For larger graphs, inexact techniques are used for matching graphs, due to the high computational cost of isomorphism-based methods. These formulations largely fall into two classes – local aligners and global aligners. Local aligners (e.g., AlignNemo [9]) aim to find localized regions in graphs that align well. A localized network alignment approach was used by Lässig et al. [3] to align gene regulation networks in E. coli. The optimization criteria in this class of methods rewards local matches and does not penalize mismatches over entire graphs. Global aligners, on the other hand, use optimization functions that reward matches and penalize mismatches over the entire graph. In [37, 36], Singh et al. describe IsoRank, a pairwise global network alignment method for protein interaction networks. IsoRank is based on the intuition that a vertex in the first network should be matched with a vertex in the second network if and only if the neighbors of the two vertices are also well matched. A number of subsequent results have extended IsoRank, including approximations by Kollias et al. [24] and multiple graph alignment in IsoRankN by Liao et al. [26]. While IsoRank is formulated as a Quadratic Program (QP), there are linear relaxations used in network alignment by Klau et al. [23], which has quadratic memory requirements. In addition to the aforementioned IsoRankN, other examples for multiple network alignment include [18, 42, 30].

Many memory-efficient heuristic methods such as GRAAL [25] and GHOST [31] have been proposed, which require cubic runtime. A principled approach, based on low rank approximation of the similarity matrix used in Isorank was described by Mohammadi et al. [28], which does not require quadratic memory. Another method that relies on eigenvectors, called EigenAlign was introduced by Feizi et al. [15]. EigenAlign estimates the similarity matrix using the dominant eigenvector of a matrix related to the product-graph between the two networks.

Many alignment techniques incorporate prior knowledge of potential matches to restrict the search space of potential matches to a smaller subset of nodes, and to guide the search process. Examples of such algorithms include the previously described Isorank [37], where sequence similarity of proteins is often used as a prior for matching networks of interacting proteins. A message passing algorithm netalignmbp, proposed by Bayati et al. [2] is based on making greedy decisions constrained by a given prior, assuming that there are no cycles in the graph.

In this paper, we focus on global graph alignment techniques that incorporate prior knowledge to yield a one-to-one mapping of nodes. We use netalignmbp with an informative prior as our network alignment substrate. However, our algorithm can be used in conjunction with any network alignment technique that allows the use of a prior.

Network alignment in Graph databases The problem of graph alignment is also commonly applied to databases. Melnik et al. [27] introduce a meta algorithm to match different data structures, or models. These models are modeled as graphs, and network alignment algorithms are used to find matching nodes. In this application, network alignment is used to manipulate and maintain schemas and match results. However, human intervention is required to validate the matches. In [7], the problem of alignment between versions of evolving RDF databases was explored. In particular, the paper proposes a method to align nodes of two consecutive versions of an RDF database using a bisimulation-based approach. The approach is flexible to accommodate nodes without labels (called “blank” nodes). In another network alignment algorithm for databases, Zhang et al. [45] define constraints for consistency and align nodes accordingly.

In [19], the authors introduce a greedy approach to identify node matchings by aligning their latent feature representations. In the first step, graph attributes such as degree and neighbors are used to establish the identity of nodes. Following this, a similarity matrix is computed for pairs of nodes across the two graphs on the basis of the previously established attributes. Then, nodes are matched in a greedy fashion. The problem of alignment of social networks with location data was explored by Riederer et al. [32]. While these methods use node and/or edge attributes to construct similarity measures between nodes, they are not well suited for the applications that motivate our approach. We explicitly leverage the fact that positions of nodes are rich sources of information, and that their geometric transformations improve the quality of alignment.

Structural Alignment One of the commonly used approaches to structural alignment uses the iterative closest point (ICP) method of Besl et al. [5], and its variants [6, 8]. The generic steps of ICP described in [33] for matching two structures involve selection of source points from both bodies. ICP then matches the sets of points and computes weights corresponding to the matches appropriately. Points that do not match are rejected and the transformation that minimizes error is computed. This process is repeated to convergence. Our algorithm can be thought of as loosely following this framework, except that we use network alignment as the process to find suitable matching. The transformation itself, as stated in Equation 2.1, is called the Orthogonal Procrustes Problem [35]. Kabsch’s algorithm [21] provides an efficient solution to this problem in three dimensions. In fact, the graph matching problem as stated by [41] can be thought of as a two-sided Procrustes problem (Conroy et al. [10]), who approximate graph matching when vertices of the two graphs can be parcellated into groups.

5 Concluding Remarks

This paper formulates a novel and important problem of aligning a class of graphs called rigid graphs. The problem integrates topological and structural alignments into a single framework, and presents a meta-algorithm for computing rigid graph alignments. The method is versatile, in that it admits a range of topological and structural alignment techniques into an expectation maximization framework. The paper presents a detailed experimental study demonstrating significant performance improvements from rigid graph alignments over prior methods. Finally, in the context of human brain connectomes, it demonstrates significant improvements in distinguishability of connectomes and identification of individual signatures in brains. Our results open significant new avenues of research in the area of rigid graph aligners, akin to topological aligners two decades back.