Structure- and Function-Aware Substitution Matrices
via Differentiable Graph Matching

Paolo Pellizzoni 1 Carlos Oliver 1 Karsten Borgwardt 1

Abstract

for alignments, such as Blast3D and Foldseek’s 3Di [8],[9]. These matrices focus on structural fragments and target
specific protein families [10], organisms [11], and phylogenetic knowledge [12]. Alignment-free methods predict
the function of biomolecules using neural networks, as seenin DeepFRI, RNAmigos, and ChemBERT [13]–[16], butthey lack interpretability compared to substitution matrices. Ideally, substitution matrices would use prior knowledge without needing pre-curated alignments. The DEDAL
model learns, in a data-driven fashion, substitution costs
using Pfam annotations, but relies on sequence data [17].

Substitution matrices, which are crafted to quantify the functional impact of substitutions or deletions in biomolecules, are central component of
remote homology detection, functional element
discovery, and structure prediction algorithms.
However, they are often limited to sequence data
and the conditioning on external priors can only
be given implicitly through the curation of theground-truth alignments they are crafted on. Here
we propose an algorithmic framework, basedon regularized optimal transport, for learning
graph-based substitution matrices from data, conditioned on any functional knowledge. In particular, our graph-neural-network-based model learns
to produce substitution matrices and graph matchings such that the resulting metric correlates withthe function at hand. Our method shows promising performance in functional similarity classification and shows potential for interpreting thefunctional importance of molecular substructures.

The graph edit distance [18], [19] models the similarity
between graphs based on an optimal matching between
their nodes, generalizing sequence alignment. The problem
is NP-hard, so several attempts have been made towards
developing heuristics, including machine-learning-based
approaches [20]–[25].
In a recent paper [26], the problem of automatically obtaining task-specific substitution matrices for biological structures has been tackled taking inspiration from the graph edit
distance literature. The substitution matrices are learnt in a
metric learning framework based on a functional prior, andthey provide an interpretable explanation of which biochemical substructure drive the function at hand.

1. Introduction
Alignment algorithms, a crucial tool in bioinformatics, use
substitution matrices to measure the impact of changes inthe sequence and structure of biomolecules, with alignment
quality depending on the choice of these matrices [1].

The optimal transport problem [27], [28], which we use
in this paper, is in general concerned with finding theminimum-cost transport plan between two distributions, andit has been used in tasks like domain adaptation [29], graphkernels [30], and generative adversarial networks [31].

Substitution matrices, like BLOSUM62 and PAM for aminoacids [2], [3], BLASTN for DNA [4], and specific matrices for RNA and chemicals [5], [6], are based on the frequency of modifications in homologous molecules [7]. Lowfrequency changes indicate functional significance, while
high-frequency changes suggest neutrality.

1.1. Contributions
This work extends our recent paper: Structure- andFunction-Aware Substitution Matrices via Learnable Graph
Matching. In: RECOMB 2024 [26]. In particular, we:

Conservation of structural data led to structural alphabets
1. generalize the GMSM (Graph Matching Substitution
Matrices) model to work with entropy-regularized optimal transport. This allows the model to have lower
time complexity and to be fully differentiable.

1

Department of Machine Learning and Systems Biology, Max
Planck Institute of Biochemistry, Germany.. Correspondence to:
Paolo Pellizzoni <pellizzoni@biochem.mpg.de>.
Published at the 2 nd Differentiable Almost Everything Workshop
at the 41 st International Conference on Machine Learning, Vienna, Austria. July 2024. Copyright 2024 by the author(s).

2. provide experimental evidence that the entropy regularization helps in finding better substitution matrices.
1

Structure- and Function-Aware Substitution Matrices via Differentiable Graph Matching

Figure 1. Architecture of GMSM. (a) Biochemical structures are transformed into graphs. (b) For each graph, its nodes are represented asa structure-aware embeddings using the same GNN. (c) The model computes the substitution matrix from node embeddings and obtains
the graph alignment with respect to the learned substitution matrix.

where Π denotes the set of bijections from V1+ to V2+ , andcv denotes the cost for node edits and ce for edge edits.

2. Preliminaries
In what follows, we denote as a graph a tuple G = (V, E),with V the set of nodes and E ⊆ V2 the set of undirectededges. Both nodes and edges can have discrete labels.

The graph edit distance, which is closely related to the
Fused Gromov-Wasserstein distance [34], is known to be
NP-hard [18]. The bipartite graph matching heuristic [32]to the graph edit distance is a heuristic technique to getapproximate GEDs in polynomial time. The main idea isto transform the problem into a linear assignment problem, which is well-known to be polynomial-time solvable (e.g. with the Hungarian algorithm), by disregarding edge edit operations. Let n = |V1 | + |V2 |. Let Π
n×nbe the set of matrices π ∈ R+
such that π1n = 1n ,T
π 1n = 1n . Then, the linear assignment problem to besolved is dC (G1 , G2 ) := minπ∈Π ⟨C, π⟩ with ⟨·, ·⟩ being
the Frobenius inner product, and with node-edit cost matrix
C. Many approximation algorithms to the GED modify
the node assignment costs in order to account for edge edit
operations [35]. The common technique is to represent
nodes as rooted substructures, such as neighborhoods [36]or subgraphs [37], and to compute the assignment costs accordingly. In fact, in this work we generalize this approachby adding an entropic regularization to the linear assignment
problem in order to make it differentiable.

2.1. Graph Neural Networks
Message passing graph neural networks (GNNs), given
a graph G, iteratively produce for each node v ∈ VG ,at each level k = 1, . . . , K, the embeddings hkv ∈ Rdk
as hkv = fupd hvk−1 , fagg {huk−1 : u ∈ N (v)} , where
fagg and fupd are the aggregate and the update operations,
respectively. The first layer of the GNN is fed with the initialnode embeddings h0v , e.g. one-hot encodings of the nodelabels. Finally, one can get a graph-level readout hG byaggregating the last-level embeddings via a function fout .
2.2. Graph Edit DistanceThe graph edit distance [18], [19] is a distance that assesses
the similarity between two graphs. In particular, it is computed as the minimum cumulative cost of the edit operations
required to transform one graph into another via node andedge insertion, deletion, and substitution, each with an associated cost. When the edit costs form a metric [32], thegraph edit distance can be equivalently defined as a graphmatching problem:

3. Methods
3.1. GMSM Architecture

Definition 2.1 (Graph Edit Distance [33]). Let G1 =(V1 , E1 ) and G2 = (V2 , E2 ) be the source and the target
graphs respectively. Let V1+ = V1 ∪ {ε1 , . . . , ε|V2 | } be thevertex set of G1 enriched with |V2 | dummy nodes ε to allow
for insertion and deletions. Let the same hold for V2+ , with
|V1 | dummy nodes. The graph edit distance (GED) between
G1 and G2 is defined byX
cv (vi , π(vi ))+GED(G1 , G2 ) = minπ∈Π

We describe the architecture of the method GMSM, which
is designed to learn expressive and interpretable substitution
matrices for biological structures and, at the same time,
when given two such graph-represented structures, to output
an interpretable alignment of the two graphs based on such
substitution matrices.
Taking inspiration from the bipartite graph matching heuristic for graph edit distance, we represent our graphs as bag oflearnable node features, and compute the edit distance between graphs by computing the optimal assignment between
such features. More formally, given a graph G we repre-vi ∈V1+Xce (vi , vj , π(vi ), π(vj )),vi ,vj ∈V12

Structure- and Function-Aware Substitution Matrices via Differentiable Graph Matching

Proposition 3.1. Let π ∗ = argminπ∈Π ⟨C, π⟩ + εΩ(π) andd(C) = ⟨C, π ∗ ⟩. Then d(C) is differentiable with gradientsent each node v ∈ V as a parametric function ψθ (v, G),
which is implemented by a graph neural network (GNN)parametrized by θ. For a pair of graphs G1 , G2 , the setof node features is computed using the same parametricfunction, a siamese network [38], as shown in Figure 1.

∇C d(C) = argminπ∈Π ⟨C, π⟩ + εΩ(π).

Then, the cost of substituting v with u is given by cu,v =∥ψθ (v, G1 ) − ψθ (u, G2 )∥2 . Node insertion and deletion
costs are obtained by computing the distance to a learnableembedding for a dummy isolated node ε. In particular, given
two graphs G1 , G2 , the model computes the corresponding
edit cost matrix CG1 ,G2 (ψθ ) = (cu,v )u∈V + ,v∈V + , which is1
2
a submatrix of a global matrix C(ψθ ).

∗When ε = 0, we define πA(C) the solution returned by theassignment algorithm A (e.g. Hungarian algorithm). We
˜ C d(C) = π ∗ (C) and apply gradient descent todefine ∇A
GMSM with this re-defined gradient. When ε > 0, we
instead use the true gradient.

4. Experimental Evaluation

Finally, the matching between the node features is computed by solving an optimal transport problem
P [27], [28],regularized with the entropic term Ω(π) = ij πij log πij :In this section, we provide experimental evidence that thesubstitution matrices learned by our method indeed distill
useful information about the conditioning priors on which it
was trained, and in particular that the entropic regularization
can help the learning task. The code is publicly available1 .

π ∗ = argminπ∈Π ⟨CG1 ,G2 (ψθ ), π⟩ + εΩ(π),
dψθ ,ε (G1 , G2 ) := ⟨CG1 ,G2 (ψθ ), π ∗ ⟩.

As shown in [26], the function dψθ ,0 (·, ·) is symmetric andsatisfies the triangle inequality.

4.1. Experimental setup
We tackle the similarity-based classification task. Given two
graphs, the task is to predict whether or not they belong to
the same class, solely as a function of their learned distance.
In particular, we evaluate the triplet accuracy and the pair
AUROC, smilarly to [23]. The goal of this task is to evaluate
whether the learned dissimilarity between graphs correlates,
at both the short range and long range scale, with the conditioning priors on which it was trained on. We focus here
on datasets of small molecules (see Section C.2), but themethod can be applied also to, e.g., proteins and RNA [26].

Unlike the original version of the model, which
 was limited to ε = 0 and had O max(n1 , n2 )3 time complexity, the regularizedversion of GMSM can run inO max(n1 , n2 )2 time [28], allowing for faster training.
Note that the function dψθ ,ε depends on the size of thegraphs that are being compared. To avoid biasing the model
into considering smaller graphs as more similar, we normalize the distance by the sum of the number of nodesof the two graphs [22] and use the graph dissimilarity
ˆ 1 , G2 ) = (|V1 | + |V2 |)−1 dψ ,ε (G1 , G2 ) as the outputd(Gθof the GMSM model, rather than the unnormalized distance.

We ask how the substitution matrices learned by GMSM
compare to alternative architectures, isolating the effects
of graph-matching, structure-awareness and prior information. To this end, we report the result of three baselinearchitectures: WL kernel, which is an approach that only
captures structural information, Siamese-GNN, which has
the same GNN architecture, but replaces the graph matching
step by computing a graph-level pooled embedding for eachgraph and by using the distance between such embeddings,
and GMSM with fixed uniform costs between all nodes.
Moreover, we report the results for GMSM both with andwithout entropic regularization during training. Since thegoal is to obtain interpretable cost matrices between simple
substructures, we fix the number of message passing layers
in GMSM to 2. For the sake of a fair comparison, we do thesame for the baselines. In order to allow the model to output
one-to-one graph matchings, we set ε = 0 at inference time.

3.2. Training GMSM
Informally, we want that the graph dissimilarity induced by
the learned substitution matrix C(ψθ ) correlates with thefunctional labels at hand. In particular, we would like forgraphs belonging to the same class to have low distance,
and for graphs belonging to different classes to have higher
distance, as common in the metric learning setting [39]. In
particular, we use the margin loss proposed in [40].
Crucially, since one wants the cost matrix C := CG1 ,G2 (ψθ )to be learnt based on the data at hand, the graph dissimilarity should have a well-defined gradient with respect to
C, to allow for the optimization of the loss function. The
original model GMSM, which sets ε = 0, relied on thefact that when the solution to the optimal transport problem
is unique, the dissimilarity is differentiable with respect to
C by Danskin’s theorem. In fact, by adding a non-zero
entropic regularization, the assignment problem becomes
convex and everywhere differentiable [28], [41].

Note that this experiment is meant to isolate the effect ofdifferent algorithmic components on the task rather than tomaximize molecule property prediction performance.
1

3

github.com/BorgwardtLab/GraphMatchingSubstitutionMatrices

Structure- and Function-Aware Substitution Matrices via Differentiable Graph Matching

Table 1. Similarity-based classification on small molecules

Mutagenicity

MethodNCI1

AIDSTrip. acc.

AUROC

Trip. acc.

AUROC

Trip. acc.

AUROC

WL kernel
Siamese-GNN
GMSM (uniform costs)

0.535 ±0.004
0.720 ±0.005
0.518 ±0.004

0.531 ±0.003
0.729 ±0.005
0.520 ±0.004

0.519 ±0.005
0.693 ±0.005
0.550 ±0.004

0.516 ±0.004
0.708 ±0.004
0.543 ±0.004

0.411 ±0.008
0.981 ±0.001
0.944 ±0.003

0.427 ±0.005
0.993 ±0.000
0.933 ±0.002

GMSM (ε = 0)
GMSM (ε ≥ 0)0.655 ±0.005
0.663 ±0.005

0.652 ±0.005
0.659 ±0.005

0.677 ±0.003
0.689 ±0.003

0.675 ±0.004
0.686 ±0.004

0.980 ±0.001
0.981 ±0.002

0.988 ±0.000
0.986 ±0.000

4.2. Similarity-based classification
We see in Table 1 that GMSM without entropic regularization (ε = 0.0) outperforms the simple WL kernel anduniform-cost baselines, showing that training on class labels
yields indeed task-specific edit costs, as expected. Moreover,
the Siamese-GNN baseline, which is not forced to assign
similarities based on a graph matching, is more expressive
and yields better results, at the expense of not providing theinterpretable edit costs, which are the focus of this work.
Allowing for entropic regularization during training provides a performance boost, suggesting that the smoothnessof the objective leads to better convergence.
Figure 2. Two sample chemicals from the Mutagenicity dataset,shown on the left, and the edit cost matrix calculated by GMSM.

4.3. Regularization at inference time
If one is willing to sacrifice getting one-to-one graph matchings in output, the model can be run with entropic regularization also at inference time. As shown in Section D.1, thiscan yield better classification performance.

this is due to the known mutagenic effect of halogens [42]and thus substituting halogens would not be likely to affectthe class label.

4.4. Retrieval

In general, these interpretable edit cost matrices offer insights into the substructures driving the function at hand,showcasing the potential of GMSM as tool for hypothesis
generation in biochemistry.

We also investigate, in Section D.2, the retrieval task, where
the model returns the graphs in a database that are mostsimilar to a query graph, testing the quality of the learned
graph dissimilarity at short scales. Results show that entropic regularization can yield better results, although the
improvement is not as marked as in the classification task.

At the same time, we notice that these matrices reflect only
the way the particular instance of the model classifies pairsof graphs, and is therefore subject to change dependingon the model architecture, the training data, and the initial
parameters of the model.

5. Analysis of the edit cost matrices
Figure 2 shows an example of two chemicals from the Mutagenicity dataset that have a very different graph topology,
and therefore would be classified as dissimilar by method
such as the WL kernel, but that have the same label. In thiscase, they are both known to be mutagenic. In particular,
GMSM (trained with ε = 1.0) learns to assign a low graphdissimilarity between the two chemicals.

6. Conclusion
In this paper, we have extended the GMSM model, whichallows to obtain task-specific substitution matrices for biochemical structures, to use entropy-regularized optimaltransport in its graph-matching module. We show that thismodification allows to lower the time complexity of themodel and at the same time to make its output differentiablewith respect to the model parameters. This in turn allows
the model to converge more easily to good solutions.

This example showcases the explainability of GMSM. Indeed, we can observe that the substructures containing halogens (Br and Cl) have very low edit costs. We hypothesize
4Structure- and Function-Aware Substitution Matrices via Differentiable Graph Matching

References
[1]

[2]

[3]

[4]

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

O. Gotoh, “Multiple sequence alignment: Algorithms andapplications,” Advances in biophysics, vol. 36, pp. 159–
206, 1999.
D. Song, J. Chen, G. Chen, et al., “Parameterized blosum
matrices for protein alignment,” IEEE/ACM Transactions
on Computational Biology and Bioinformatics, vol. 12,
no. 3, pp. 686–694, 2014.
W. J. Wilbur, “On the PAM matrix model of protein evolution.,” Molecular biology and evolution, vol. 2, no. 5,
pp. 434–447, 1985.
A. C. Eklund, P. Friis, R. Wernersson, and Z. Szallasi, “Optimization of the blastn substitution matrix for predictionof non-specific DNA microarray hybridization,” Nucleic
acids research, vol. 38, no. 4, e27–e27, 2010.
I. L. Hofacker, S. H. Bernhart, and P. F. Stadler, “Alignmentof RNA base pairing probability matrices,” Bioinformatics,
vol. 20, no. 14, pp. 2222–2227, 2004.
R. P. Sheridan, “The most common chemical replacements
in drug-like compounds,” Journal of chemical informationand computer sciences, vol. 42, no. 1, pp. 103–108, 2002.
S. Henikoff and J. G. Henikoff, “Amino acid substitution
matrices,” Advances in protein chemistry, vol. 54, pp. 73–
98, 2000.
C.-H. Tung, J.-W. Huang, and J.-M. Yang, “Kappa-alphaplot derived structural alphabet and blosum-like substitution matrix for rapid search of protein structure database,”
Genome biology, vol. 8, no. 3, pp. 1–16, 2007.
M. van Kempen, S. S. Kim, C. Tumescheit, et al., “Fast andaccurate protein structure search with Foldseek,” Nature
Biotechnology, pp. 1–4, 2023.
F. Keul, M. Hess, M. Goesele, and K. Hamacher, “Pfasum:
A substitution matrix from Pfam structural alignments,”
BMC bioinformatics, vol. 18, pp. 1–14, 2017.
R. A. Sutormin, A. B. Rakhmaninova, and M. S. Gelfand,“Batmas30: Amino acid substitution matrix for alignmentof bacterial transporters,” Proteins: Structure, Function,
and Bioinformatics, vol. 51, no. 1, pp. 85–95, 2003.
J. M. Koshi and R. A. Goldstein, “Context-dependent optimal substitution matrices,” Protein Engineering, Design
and Selection, vol. 8, no. 7, pp. 641–645, 1995.
V. Gligorijević, P. D. Renfrew, T. Kosciolek, et al.,“Structure-based protein function prediction using graphconvolutional networks,” Nature communications, vol. 12,
no. 1, p. 3168, 2021.
M. Kulmanov and R. Hoehndorf, “DeepGOPlus: Improved
protein function prediction from sequence,” Bioinformatics,
vol. 36, no. 2, pp. 422–429, 2020.
C. Oliver, V. Mallet, R. S. Gendron, et al., “Augmentedbase pairing networks encode RNA-small molecule binding preferences,” Nucleic acids research, vol. 48, no. 14,
pp. 7690–7699, 2020.
S. Chithrananda, G. Grand, and B. Ramsundar, “Chemberta: Large-scale self-supervised pretraining for molecular property prediction,” arXiv preprint arXiv:2010.09885,
2020.
F. Llinares-López, Q. Berthet, M. Blondel, O. Teboul, and
J.-P. Vert, “Deep embedding and alignment of protein sequences,” Nature Methods, vol. 20, no. 1, pp. 104–111,
2023.

[19]

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]
[29]

[30]

[31]

[32]

5

H. Bunke and K. Riesen, “Graph edit distance–optimaland suboptimal algorithms with applications,” Analysis ofComplex Networks: From Biology to Linguistics, pp. 113–
143, 2009.
M. Neuhaus and H. Bunke, “A probabilistic approach tolearning costs for graph edit distance,” in Proceedings ofthe 17th International Conference on Pattern Recognition,
2004. ICPR 2004., IEEE, vol. 3, 2004, pp. 389–393.
Y. Bai, H. Ding, S. Bian, T. Chen, Y. Sun, and W. Wang,“Simgnn: A neural network approach to fast graph similarity computation,” in Proceedings of the twelfth ACM
international conference on web search and data mining,
2019, pp. 384–392.
R. Ranjan, S. Grover, S. Medya, V. Chakaravarthy, Y. Sabharwal, and S. Ranu, “Greed: A neural framework for learning graph distance functions,” Advances in Neural Information Processing Systems, vol. 35, pp. 22 518–22 530,
2022.
K. D. Doan, S. Manchanda, S. Mahapatra, and C. K. Reddy,“Interpretable graph similarity computation via differentiable optimal alignment of node embeddings,” in Proceedings of the 44th International ACM SIGIR Conference on
Research and Development in Information Retrieval, 2021,
pp. 665–674.
Y. Li, C. Gu, T. Dullien, O. Vinyals, and P. Kohli, “Graph
matching networks for learning the similarity of graphstructured objects,” in International conference on machinelearning, PMLR, 2019, pp. 3835–3845.
P. Riba, A. Fischer, J. Lladós, and A. Fornés, “Learning
graph edit distance by graph neural networks,” Pattern
Recognition, vol. 120, p. 108 132, 2021.
J. Heo, S. Lee, S. Ahn, and D. Kim, “EPIC: Graph augmentation with edit path interpolation via learnable cost,”
arXiv preprint arXiv:2306.01310, 2023.
P. Pellizzoni, C. Oliver, and K. Borgwardt, “Structure- andfunction-aware substitution matrices via learnable graphmatching,” in Research in Computational Molecular Biology, 2024.
M. Cuturi, “Sinkhorn distances: Lightspeed computationof optimal transport,” in Advances in Neural Information
Processing Systems, vol. 26, 2013.
G. Peyré and M. Cuturi, Computational optimal transport,
2020. arXiv: 1803.00567 [stat.ML].
N. Courty, R. Flamary, D. Tuia, and A. Rakotomamonjy,“Optimal transport for domain adaptation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 39,
no. 9, pp. 1853–1865, 2017.
M. Togninalli, E. Ghisu, F. Llinares-López, B. Rieck, and K.
Borgwardt, “Wasserstein weisfeiler-lehman graph kernels,”
in Advances in Neural Information Processing Systems,
2019.
M. Arjovsky, S. Chintala, and L. Bottou, “Wasserstein
generative adversarial networks,” in Proceedings of the34th International Conference on Machine Learning, D.
Precup and Y. W. Teh, Eds., ser. Proceedings of Machine
Learning Research, vol. 70, PMLR, Jun. 2017, pp. 214–
223.
S. Bougleux, L. Brun, V. Carletti, P. Foggia, B. Gaüzère,and M. Vento, “Graph edit distance as a quadratic assignment problem,” Pattern Recognition Letters, vol. 87,
pp. 38–46, 2017, Advances in Graph-based Pattern Recognition, ISSN: 0167-8655.

Structure- and Function-Aware Substitution Matrices via Differentiable Graph Matching
[33]

[34]

[35]

[36]

[37]

[38]
[39]

[40]

[41]

[42]

[43]

[44]
[45]

[46]

[47]

[48]

H. Bunke and K. Riesen, “Graph classification based on dissimilarity space embedding,” in Structural, Syntactic, andStatistical Pattern Recognition: Joint IAPR International
Workshop, SSPR & SPR 2008, Orlando, USA, December
4-6, 2008. Proceedings, Springer, 2008, pp. 996–1007.
V. Titouan, N. Courty, R. Tavenard, C. Laetitia, and R.
Flamary, “Optimal transport for structured data with application on graphs,” in Proceedings of the 36th International
Conference on Machine Learning.
D. B. Blumenthal, N. Boria, J. Gamper, S. Bougleux, and
L. Brun, “Comparing heuristics for graph edit distancecomputation,” The VLDB journal, vol. 29, no. 1, pp. 419–
458, 2020.
D. B. Blumenthal and J. Gamper, “Improved lower bounds
for graph edit distance,” IEEE Transactions on Knowledge
and Data Engineering, vol. 30, no. 3, pp. 503–516, 2017.
V. Carletti, B. Gaüzere, L. Brun, and M. Vento, “Approximate graph edit distance computation combining bipartite
matching and exact neighborhood substructure distance,”in Graph-Based Representations in Pattern Recognition:
10th IAPR-TC-15 International Workshop, GbRPR 2015,
Beijing, China, May 13-15, 2015. Proceedings 10, Springer,
2015, pp. 188–197.
D. Chicco, “Siamese neural networks: An overview,” Artificial neural networks, pp. 73–94, 2021.
R. Hadsell, S. Chopra, and Y. LeCun, “Dimensionality reduction by learning an invariant mapping,” in 2006 IEEE
computer society conference on computer vision and pattern recognition (CVPR’06), IEEE, vol. 2, 2006, pp. 1735–
1742.
C.-Y. Wu, R. Manmatha, A. J. Smola, and P. Krahenbuhl,“Sampling matters in deep embedding learning,” in Proceedings of the IEEE international conference on computer
vision, 2017, pp. 2840–2848.
M. Cuturi, O. Teboul, and J.-P. Vert, “Differentiable ranking
and sorting using optimal transport,” in Advances in Neural
Information Processing Systems, vol. 32, 2019.
H. Brem, A. B. Stein, and H. S. Rosenkranz, “The mutagenicity and dna-modifying effect of haloalkanes,” Cancer
research, vol. 34, no. 10, pp. 2576–2579, 1974.
K. Xu, W. Hu, J. Leskovec, and S. Jegelka, “How powerful
are graph neural networks?” In International Conference
on Learning Representations, 2018.
C. Morris, M. Ritzert, M. Fey, et al., “Weisfeiler and leman
go neural: Higher-order graph neural networks,” 2019.
F. Schroff, D. Kalenichenko, and J. Philbin, “Facenet: A
unified embedding for face recognition and clustering,” inProceedings of the IEEE conference on computer visionand pattern recognition, 2015, pp. 815–823.
P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Lio,and Y. Bengio, “Graph attention networks,” arXiv preprint
arXiv:1710.10903, 2017.
C. Morris, N. M. Kriege, F. Bause, K. Kersting, P. Mutzel,
and M. Neumann, “Tudataset: A collection of benchmark datasets for learning with graphs,” arXiv preprint
arXiv:2007.08663, 2020.
J. Kazius, R. McGuire, and R. Bursi, “Derivation and validation of toxicophores for mutagenicity prediction,” Journal of medicinal chemistry, vol. 48 1, pp. 312–20, 2005.

[49]

[50]

6

K. Riesen and H. Bunke, “IAM graph database repository
for graph based pattern recognition and machine learning,”
in Structural, Syntactic, and Statistical Pattern Recognition,
Joint IAPR International Workshop, SSPR & SPR 2008,
Orlando, USA, December 4-6, 2008. Proceedings, N. da
Vitoria Lobo, T. Kasparis, F. Roli, et al., Eds., ser. Lecture
Notes in Computer Science, vol. 5342, Springer, 2008,
pp. 287–297.
N. Wale and G. Karypis, “Comparison of descriptor spaces
for chemical compound retrieval and classification,” in Proceedings of the 6th IEEE International Conference on Data
Mining (ICDM), IEEE Computer Society, 2006, pp. 678–
689.

Structure- and Function-Aware Substitution Matrices via Differentiable Graph Matching

A. Additional backgroundA.1. Optimal transportn
In the discrete version of optimal transport, we consider two finite sets of with associated weights. Let {ai }mi=1 and {bj }j=1
be two such sets with corresponding weight vectors, or probability distributions, u ∈ Rm and v ∈ Rn such that u⊤ 1n = 1
and v ⊤ 1m = 1. Here, 1n and 1m are vectors of ones of appropriate dimensions. The goal is to find a transportation planm×n
π ∈ R+
that transfers the mass from {ai } to {bj } while minimizing the total transportation cost. This can be formulated
as the linear program:dC (u, v) =min ⟨π, C⟩,

m×n
π∈R+

subject to the constraints π1n = u and π ⊤ 1m = v, where cij represents the cost of transporting mass from ai to bj andC = (cij )i=1,...,m;j=1,...,n . In [27] it was proposed an entropy-regularized version of the optimal transport problem, whichallows the resulting distance to be differentiable with respect to C and to solve the problem in quadratic time using the
Sinkhorn algorithm.
A closely related problem is the linear assignment problem, where the sets have the same size, i.e. n = m, and the weightn×n
vectors are defined as u = v = 1n . Then, the feasible assignments are the set Π of matrices π ∈ R+
such that π1n = 1n ,T
π 1n = 1n . The task is then to find minπ∈Π ⟨π, C⟩. The problem can be solved in cubic time, e.g., with the Hungarianalgorithm.
A.2. Graph neural networks
Message passing graph neural networks (GNNs), given a graph G, iteratively produce for each node v ∈ VG , at each levelk = 1, . . . , K, the embeddings hkv ∈ Rdk as hkv = fupd hvk−1 , fagg {huk−1 : u ∈ N (v)} , where fagg and fupd are theaggregate and the update operations, respectively. The first layer of the GNN is fed with the initial node embeddings h0v , e.g.
one-hot encodings of the node labels.
In [43] it was shown that there exist injective functions fagg , fupd yielding GNNs that are provably as expressive as color
refinement.
An example of such functions that leads to models that are provably as expressive as color refinement [44], denoting ∥ asconcatenation, isX
hkv = mlp hvk−1
huk−1 ∈ Rdk .
u∈N (v)

In fact, provided that different node labels are encoded to linearly independent h0v ’s, even the following simpler architecture,
denoting with σ a nonlinear function such as ReLU, is as expressive as color refinement [44]:Xhuk−1 ∈ Rdk .
hkv = σ W1k hvk−1 + W2k(1)u∈N (v)

We then let the final node embedding be a function of the per-layer node embeddings as ψ(v, G) = f ({hkv : k =0, . . . , K}). We denote (v, G1 ) =ψ (u, G2 ) if ψ(v, G1 ) = ψ(u, G2 ), and denote by Vψθ the set of equivalence classes
induced by =ψ . Moreover, we will denote, with abuse of notation, G1 =ψ G2 if {{ψ(v, G1 = (V1 , E1 )) : v ∈ V1 }} ={{ψ(v, G2 = (V2 , E2 )) : v ∈ V2 }}, and denote for brevity with Gψ the set G\ =ψ of equivalence classes induced by =ψ ongraphs.
A.3. Heuristics for the graph edit distanceThe bipartite graph matching heuristic [32] to the graph edit distance is an heuristic technique to get approximate GEDs inpolynomial time. The main idea is to transform the problem into a linear assignment problem, which is well-known to bepolynomial-time solvable (e.g. with the Hungarian algorithm), by disregarding edge edit operations and only considering
node edit ones.
7

Structure- and Function-Aware Substitution Matrices via Differentiable Graph Matching

Then, the linear assignment problem to be solved isX

dC (G1 , G2 ) := min⟨C, π⟩ = min
π∈Π

π∈Π

cvi ,π(vi ) ,vi ∈V1+

with ⟨·, ·⟩ being the Frobenius inner product, and with cost matrix C defined as
 c
c
···
∞ 
···
c
1,1

 ...

 cn1 ,1C=
 cε,1

 ..
.
∞

..

.
···
···
..
.
···

1,n2

1,ε

..
.

..
.
∞
0
..
.
0

cn1 ,n2
∞
..
.
cε,n2

..
.

..

.
···
···
..
.
···



cn1 ,ε 

0 

.. 
.
0

where ci,j denotes the cost of a node substitution vi → vj , ci,ε denotes the cost of a node deletion ui → ε, and cε,j denotes
the costs of a node insertion ε → vj with vi ∈ V1 and vj ∈ V2 ; and where Π is the set of permutation matrices of size
|V1 | + |V2 |, representing bijections from V1+ to V2+ .
If one chooses ci,j = cv (vi , vj ), i.e. setting the cost of assigning node vi to node vj as the node edit cost in the original
problem, then dC (G1 , G2 ) lower bounds GED(G1 , G2 ).
Many approximation algorithms to the GED modify the node assignment costs in order to account for edge edit operations[35]. The common technique is to represent the nodes as some local substructure around them, such as neighborhoods [36]or subgraphs [37], and to compute the assignment costs accordingly. In many practical applications, representing the graphas such a bag of local structures is enough to compute a good approximation to the true edit distance [35].
In fact, in this work we generalize this approach by representing nodes as their GNN embeddings and by adding an entropic
regularization to the linear assignment problem in order to make it differentiable.
A.4. Metric learning
Common losses for metric learning are the contrastive loss [39] and the triplet loss [45], which act on pairs and triplets ofgraphs, respectively. We use the margin loss proposed in [40], which reportedly yields better results than the contrastive one.
Moreover, it can work with randomly sampled pairs, while the triplet loss usually requires hard or semi-hard sample mining
to work properly [40].
In particular, for a pair of graphs G1 , G2 , the loss is defined asℓmargin (G1 , G2 ) = max 0, α + y(dˆC (G1 , G2 ) − β) ,with y = 1 for positive pairs, i.e. both graphs belonging to the same class, and y = −1 for negative pairs. This loss strives
to push the graph dissimilarities of graphs belonging to the same class to be less than β − α and the dissimilarities of graphs
belonging to different classes to be more than β + α. In our experiments, we set β = 0.5 and α = 0.1.

B. Architecture of GMSM
Our model, GMSM, is implemented as follows. The GNN message passing layers are realized by Equation 1, with thenonlinearity realized by a ReLU. We fix the number of layers to K = 2. The final node embeddings ψ(v, G) are obtainedK

by concatenating the embeddings for each layer and by normalizing them, i.e. hv =k=0

hkv and ψ(v, G) = hv /∥hv ∥2 .

Then, the cost of substituting the structure rooted at v with the one rooted at u is given by cu,v = ∥ψθ (v, G1 ) − ψθ (u, G2 )∥2 .
Node insertion and deletion costs are obtained by computing the distance to a learnable embedding for a dummy isolated
node ε, that is cv,ε = cε,v = ∥ψθ (v, G1 ) − ψθ (ε)∥2 . Note that the normalization of the embeddings makes it so that the edit
costs are bounded in [0, 2].
This yields a global learnable substitution matrix C(ψθ ) = (cu,v )u,v∈Vψθ of costs between elements of Vψθ ∪ {ε}, whichcorrespond to nodes and their rooted subgraph explored by the GNN. In particular, given two graphs G1 , G2 , the modelcomputes the corresponding edit cost matrix CG1 ,G2 (ψθ ) = (cu,v )u∈V + ,v∈V + , which is a submatrix of C(ψθ ).
1

82

Structure- and Function-Aware Substitution Matrices via Differentiable Graph Matching

Finally, the matching between the
Pnode features is computed by solving an optimal transport problem [27], [28], regularizedwith the entropic term Ω(π) = ij πij log πij :π ∗ = argminπ∈Π ⟨CG1 ,G2 (ψθ ), π⟩ + εΩ(π),

dψθ ,ε (G1 , G2 ) := ⟨CG1 ,G2 (ψθ ), π ∗ ⟩.
ˆ 1 , G2 ) = (|V1 | + |V2 |)−1 dψ ,ε (G1 , G2 ), as well as the edit cost matrix
Finally, the model outputs the graph dissimilarity d(Gθand the transport plan.
With respect ot the original work [26], we substituted the GNN layers from GAT layers [46] to layers given by Equation 1[44] and we concatenate the embedding of all the layers. This yields a provably more expressive model [44].

C. Experimental setup
The model is trained taking the same number of random positive and negative pairs from the training set. In particular,
for each anchor graph we sample positive pairs by selecting 10 random graphs with the same label, and for each each ofsuch positive pairs we sample negative pairs by selecting 10 random graphs with different label. The model parameters are
optimized with respect to the margin loss, using β = 0.5 and α = 0.1, with the Adam optimizer. We set the learning rate to
0.001. The best model is selected using the validation set pair AUROC metric on the classification task, using early stopping.
For the GMSM with regularization, we choose ε ∈ {0.1, 1.0} based on the validation set pair AUROC metric. In particular,
for the Mutagenicity and NCI1 datsets, the highest validation AUROC is obtained for ε = 1.0, while for the AIDS
datset it is obtained for ε = 0.1. These models are then used at inference time on both the similarity-based classification
task and the retrieval task. The models are tested using the same sampling strategy for random positive and negative pairs
from the validation and test set.
C.1. Tasks
The first task we tackle is the similarity-based classification task. Given two graphs, the task is to predict whether or not they
belong to the same class, solely as a function of their learned distance. In particular, we evaluate two metrics. The firstone is the triplet accuracy. Namely, given an anchor graph, one graph from the same class of the anchor (positive pair) and
one graph from another class (negative pair), the triplet of graphs is considered a successful prediction if the distance ofthe positive pair is lower than the one of the negative pair. The second is pair AUROC, the area under the ROC curve forclassifying pairs of graphs as similar or not based on a distance threshold. The goal of this task is to evaluate whether thelearned dissimilarity between graphs correlates, at both the short range and long range scale, with the conditioning priors on
which it was trained on. We report the means and standard deviations over 5 different random samples from the training set.
The second task we evaluate is the retrieval task. In particular, given a query graph, the task is to return a set of graphs thatare the most similar to the query one. The goal of this task is to evaluate the quality of the learned graph similarity at very
short scales. Indeed, for the retrieval task it does not matter if some positive pairs are at a high distance as long as there areenough positive pairs at a very short distance, which will be returned as hits by the retrieval procedure. In particular, we takeas queries the graphs of the test set, and search for the hits in the training set, which serves then as the searchable database.
We report the precision@k (P@k), for k ∈ {10, 50}. In particular, we report the mean and the standard deviation over all
queries.
C.2. Datasets
In our experimental evaluation we three small molecule datasets. All datasets are split into training, validation and test sets
at random and with ratios {0.8, 0.1, 0.1}.
The datasets are obtained from the TUDataset [47] and contain small molecules annotated with a label on mutagen activity(Mutagenicity [48], [49]), on activity against non-small cell lung cancer (NCI1 [50]) and on evidence of anti-HIV
activity (AIDS [49]).
9Structure- and Function-Aware Substitution Matrices via Differentiable Graph Matching

D. Additional experimental results
D.1. Sensitivity to entropic regularization
Table 2 report the test set metrics for GMSM models that are trained and run at inference time with different levels ofentropic regularization.
As shown by the results, adding an entropic regularization at inference time can boost the classification performance. In
particular, for models trained with high regularization (ε = 1.0), using regularization also at inference time provides asignificant increase in the metrics. This could be due to the fact that when one removes the regularization at inference time,
there is shift in the distributions of the predicted graph dissimilarities. Maintaining the regularization, this shift is removed.
For models trained with smaller or no regularization, this clear behavior does not appear. Adding regularization at inferencetime can nonetheless prove to be an effective hyperparameter to select with a validation set.
Table 2. Effect of entropic regularization at inference time

Mutagenicity

NCI1

AIDS

Training

Inference

Trip. acc.

AUROC

Trip. acc.

AUROC

Trip. acc.

AUROC

GMSM (ε = 0.0)

ε = 0.0
ε = 0.1
ε = 1.00.655 ±0.005
0.663 ±0.003
0.651 ±0.004

0.652 ±0.005
0.661 ±0.004
0.648 ±0.005

0.677 ±0.003
0.664 ±0.004
0.660 ±0.004

0.675 ±0.004
0.661 ±0.004
0.658 ±0.005

0.980 ±0.001
0.977 ±0.001
0.984 ±0.001

0.988 ±0.000
0.985 ±0.000
0.988 ±0.000

GMSM (ε = 0.1)

ε = 0.0
ε = 0.1
ε = 1.0

0.651 ±0.005
0.660 ±0.003
0.646 ±0.005

0.648 ±0.005
0.659 ±0.004
0.643 ±0.005

0.679 ±0.003
0.676 ±0.002
0.662 ±0.003

0.676 ±0.004
0.674 ±0.004
0.659 ±0.005

0.981 ±0.002
0.979 ±0.001
0.984 ±0.002

0.986 ±0.000
0.985 ±0.000
0.987 ±0.000

GMSM (ε = 1.0)

ε = 0.0
ε = 0.1
ε = 1.0

0.663 ±0.005
0.677 ±0.004
0.680 ±0.004

0.659 ±0.005
0.673 ±0.004
0.682 ±0.005

0.689 ±0.003
0.652 ±0.003
0.690 ±0.003

0.686 ±0.004
0.643 ±0.004
0.692 ±0.004

0.973 ±0.001
0.972 ±0.001
0.980 ±0.001

0.980 ±0.000
0.980 ±0.000
0.984 ±0.000

D.2. Retrieval
In the retrieval task, we observe as in the classification task, that GMSM models are in general more expressive that the
WL kernel and the baseline with uniform costs, although the latter performs surprisingly well on AIDS. Moreover, theSiamese-GNN model in general outperforms GMSM. Interestingly, the GMSM model trained with entropic regularizationand selected based on the validation performance on the classification task, as described in Section C, preform very similarly
to GMSM trained with no entropic regularization.
Table 3. Retrieval performance on small molecules

MethodMutagenicity

NCI1

AIDS

P@10

P@50

P@10

P@50

P@10

P@50

WL kernel
Siamese-GNN
GMSM (uniform costs)

0.700 ±0.276
0.785 ±0.297
0.646 ±0.262

0.643 ±0.248
0.779 ±0.293
0.613 ±0.215

0.701 ±0.256
0.777 ±0.338
0.597 ±0.229

0.627 ±0.203
0.771 ±0.332
0.575 ±0.179

0.892 ±0.180
0.988 ±0.105
0.998 ±0.021

0.802 ±0.193
0.985 ±0.116
0.991 ±0.073

GMSM (ε = 0)
GMSM (ε ≥ 0)0.773 ±0.294
0.774 ±0.284

0.751 ±0.271
0.757 ±0.268

0.773 ±0.292
0.772 ±0.289

0.744 ±0.278
0.747 ±0.268

0.991 ±0.076
0.988 ±0.094

0.988 ±0.090
0.985 ±0.111

D.3. Stability of the learnt substitution matrices
The learnt substitution matrices depend on the model parameters that are learnt from the training data at hand, and therefore
correspond in local minima in the loss function landscape. Because of this, they can depend on several factors, including thespecific model architecture, the training data and the initialization of the parameters. For example, we observe that switchingfrom the GAT-based architecture in the original GMSM paper [26] to the more expressive architecture we use here leads todifferent edit costs, e.g., for the graphs of Figure 2.

10