Expressive Graph Informer Networks
Jaak Simm Adam Arany Edward De Brouwer Yves Moreau
stat.ML cs.LG




























































  
  
  Applying machine learning to molecules is challenging because of their natural representation as graphs rather than vectors. Several architectures have been recently proposed for deep learning from molecular graphs, but they suffer from information bottlenecks because they only pass information from a graph node to its direct neighbors.
    
  
  
  
  Here, we introduce a more expressive route-based multi-attention mechanism that incorporates features from routes between node pairs.
  We call the resulting method .
  A single network layer can therefore attend to nodes several steps away.
  
  
  We show empirically that the proposed method compares favorably against existing approaches in two prediction tasks: (1) 13C Nuclear Magnetic Resonance (NMR) spectra, improving the state-of-the-art with an MAE of 1.35 ppm and (2) predicting drug bioactivity and toxicity.
  Additionally, we develop a variant called injective  that is provably as powerful as the Weisfeiler-Lehman test for graph isomorphism.
  Furthermore, we demonstrate that the route information allows the method to be informed about the nonlocal topology of the graph and, thus, even go beyond the capabilities of the Weisfeiler-Lehman test.
  




Graph Neural Networks, Graph Convolution, NMR, bioactivity prediction, toxicity prediction, molecules



Expressive Graph Informer Networks
    Jaak Simm12,
Adam Arany12,
Edward De Brouwer1, and
Yves Moreau1
1ESAT-STADIUS,
KU Leuven, 
Leuven, 3001, Belgium
2Equal contribution as first authors.
Manuscript received x; revised y. 
Corresponding author: J. Simm (email: jaak.simm@esat.kuleuven.be).

2020 IEEE. Personal use of this material is permitted.
Permission from IEEE must be obtained for all other uses,
in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.

    
================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================
































§ INTRODUCTION



Graphs are used as a natural representation for objects in many domains, such as compounds in computational chemistry or protein–protein interaction networks in systems biology.
Machine learning approaches for graphs fall into two main lines of research: the spectral and the spatial approach.
The spectral approach relies on the eigenvalue decomposition of the Laplacian of the graph and is well suited for problems involving a single fixed graph structure.


By contrast, spatial graph methods work directly on the nodes and edges of the graph.

Convolutional neural networks (CNN) have inspired many spatial methods.
Similarly to a local convolution filter running on the 2D grid of an image, spatial approaches update the hidden vectors of a graph node by aggregating the hidden vectors of its neighbors.
Several spatial graph methods have been proposed, which vary in how to carry out this update.
The most straightforward approach is to sum (or average) the hidden vectors of neighbors, and then transform the results with a linear layer and a nonlinearity (e.g., as proposed in Neural Fingerprints <cit.> and Graph Convolutional Networks <cit.>).
Instead of a simple dense feedforward layer, some methods use GRU-based gating <cit.>, edge hidden vectors <cit.>, and attention to the neighbors <cit.>.
One of the advantages of spatial methods over spectral methods is that they can be straightforwardly applied to problems involving multiple graphs.

However, the update step of a node in spatial methods only has access to its own neighborhood, which limits the information flow throughout the whole graph. Increasing the accessible neighborhood of each node requires the stacking of several layers.
The dilated filters in CNNs <cit.> are an example of solving this neighborhood limitation in images, by providing an efficient way to gather information over longer distances (see Figure <ref>, (A) and (B)).

In this work, we propose a method that flexibly aggregates information over longer graph distances in one step, analogously to dilated filters.
Our  approach is inspired by the sequence-to-sequence transformer  model <cit.>.
The core contribution of our work is the introduction of route-based multi-head self-attention (RouteMHSA), which allows the attention mechanism to also access route information and, thus, base its attention scores both on the features of the nodes and the route between them.
This enables  to gather information from nodes that are not just direct neighbors.
Furthermore, the proposed approach can straightforwardly use edge features, as they can be incorporated into the route features. In the case of chemical compounds, this allows using the annotation of the bonds between the atoms (single, double, aromatic, etc.).
The central idea is illustrated in Figure <ref>.
Additionally, we investigate the expressiveness of , by describing a variant that is provably as powerful as the Weisfeiler-Lehman (WL) test and showing empirically that it can even go beyond the WL test. 



Our  compares favorably with state-of-the-art graph-based neural networks on 13C NMR spectrum and drug bioactivity prediction tasks, which makes it suitable for both node-level and graph-level prediction tasks.




§ PROPOSED APPROACH



 §.§ Setup

In this work, we consider learning problems on graphs where both nodes and edges are labeled.
For a graph G_i with N_i nodes,  we denote its node features by a matrix X_i ∈ℝ^N_i × F_nodes and the route information by tensor P_i ∈ℝ^N_i × N_i × F_route, F_nodes and F_route are the dimensions of the nodes and route features, respectively and the route information X_i is computed using the adjacency matrix of the graph and the edge labels.
For example, for chemical compound graphs, where nodes are atoms, the vector between a pair of atoms in the route information tensor (e.g., vector P_i[k,l] for atoms k and l) can contain information on the type of route that connects the two atoms (rigid, flexible, aromatic) and how far they are from each other (shortest path).

We propose a method that works for both (1) graph-level tasks, such as graph classification where the whole graph G_i is assigned a label, and (2) node-level tasks where one should predict a value or a label for nodes (or subset of nodes) of the graph.



 §.§ Dot-product self-attention

Our approach is inspired by the encoder used in the transformer network <cit.>, which proposed dot-product attention for sequence models formulated as

    Attn(Q, K, V) = σ(
            1/√(d_k)Q K^⊤) V,

where Q, K ∈ℝ^N_i × d_k and V ∈ℝ^N_i × d_v are the query, key, and value matrices and σ is the Softmax function, which makes sure that the attention probabilities for each node sum to one.
This attention can be used as self-attention by projecting the input d-dimensional hidden vectors H ∈ℝ^N_i × d into queries, keys and values. Q, K and V are computed from the input hidden vectors with corresponding weight matrices W_Q ∈ℝ^d_k × d, W_K ∈ℝ^d_k × d and W_V ∈ℝ^d_v × d (i.e., Q = W_Q H^T, K = W_KH^T and V = W_V H^T).
Intuitively, the k-th row of Q is the query vector of the k-th node (attending node), the l-th row of K is the key vector of node l (attended node) and the l-th row of V is the value vector for node l (value of the attended node).



To encode sequences with dot-product attention, the original transformer network <cit.> adds position encodings to the hidden vectors H.
In particular, sine waves with different frequencies and phases were used in the original paper to encode the position of the sequence elements.
This allows the network to change attention probabilities based on the positions.



 §.§ Route-based dot-product self-attention

One aspect to note is that in the case of graphs, there is no analogue for the global position of nodes.

Therefore, we propose a novel stable relative addressing mechanism for graphs.
The proposed mechanism adds a new component to the dot-attention that depends the features of the route between the two nodes in the graph.
This route component is made up of a query and key part.
For the route query Q_R ∈ℝ^N_i × d_r we map the hidden vectors of the nodes (H) by linear transformation

    Q_R = W_Q^route H,

where d_r is the dimension of the route query and W_Q is the query weight matrix.
Similarly, for the route key K_R ∈ℝ^N_i × N_i × d_r we map the route features P (for readability we drop the sample index i) by linear transformation

    K_R = W_K^route P.

The route information allows the attention mechanism to access the topology of the graph.
Note that while Q_R is a matrix of route queries, one for each node of the graph, K_R is a tensor of rank 3, with each slice K_R[k,l] corresponding to the route from node k to node l.

The route query and key can be now contracted into a N_i × N_i matrix of (logit) attention scores, which can be added to the original node-based scores QK^⊤.
Using the einsum convention[We follow the convention introduced by Numpy and used also in Pytorch.], let

    Q_R ⊗ K_R = einsum_kf,klf→kl(Q_R, K_R)

or equivalently

    (Q_R ⊗ K_R)[k,l] = Q_R[k]^⊤ K_R[k,l].

The intuition behind this tensor product is straight-forward: we compute the scalar product between the query vector Q_R[k] of the node k and the key vector K_R[k,l] of the route from k to l.


Now the route-based attention probabilities A ∈ [0,1]^N_i × N_i are given by

    A = σ(
            1/√(d_k + d_r)(Q K^⊤ + Q_R ⊗ K_R)
        ),

where we also added the size of the route key d_r to the normalization for keeping the values within trainable (non-plateau) ranges of the softmax function.
This enables the network to use both the information on the routes (Q_R ⊗ K_R), as well as the node hidden vectors (Q K^⊤), to decide the attention probabilities.

Analogously, we project the route information into the value component V of the dot-product attention.
Firstly, we map the route data P into

    V_R = W_V^route P,

where V_R ∈ℝ^N_i × N_i × d_v is the tensor of rank 3 with slice V_R[k,l] corresponding to the value vector for the route from node k to node l.
The route value tensor V_R can be then weighted according to the attention probabilities A from Eq. <ref> giving as the final route-based graph attention:

    RouteAttn(Q,K,V,Q_R,K_R,V_R) 
    
      = AV + einsum_kl,klv→kv(A, V_R).





 §.§ Locality-constrained attention

The route-based attention defined in Eq. <ref> is general and allows unconstrained access throughout the graph.
In some applications, localized access is preferred and this can be easily forced by adding a masking matrix to the attention scores:

    A_LC = σ(
            1/√(d_k + d_r)(Q K^⊤ + Q_R ⊗ K_R) + M_route),

where M_route∈{-∞,0}^N_i × N_i is 0 for unmasked routes and minus infinity for masked connections (in practice, we use a large negative value).
This allows us, for example, to mask all nodes whose shortest distance from the current node is larger than 3 or some other number, thereby creating an attention ball around the current node.
Similarly, it is possible to use the mask M_route to create an attention shell (i.e., attending nodes within a given range of shortest distances).

In our experiments, we only used the attention ball with different values for the radius.



 §.§ Multi-head and mini-batching

Similarly to the transformer architecture, we use multiple heads.
This means several route-based self-attentions are executed in parallel and their results are concatenated.
This allows the network to attend to several neighbors at the same times.



  
Pool node Additionally, we introduce a pool node that has a different embedding vector than the graph nodes.
The pool node has no edges to the graph nodes, but is unmasked by M_route, so that the attention mechanism can always attend and read, if required.
This idea has two motivations.
First, it allows the information to be easily shared across the graph.
Second, the pool node can be used as a “not found” answer when the attention mechanism does not find the queried node within the graph nodes.



  
Mini-batching The proposed algorithm supports mini-batching via a few simple modifications.
To introduce the batch dimension, the input tensors H and P must first be padded with zeros to have the same size, after which they are stacked.
The linear projections Q,K,V,Q_R,K_R,V_R are now indexed by both batch (sample) and head.
We replace all matrix multiplications by batched matrix products and in einsum we introduce leading batch and head dimensions.
For the attention mechanism, we add an additional mask over the nodes, M_node, with large negative values for non-existing nodes (padded) for each sample:

    A_LCB = σ(
            1/√(d_k + d_r)(Q K^⊤ + Q_R ⊗ K_R) + M_LCB),

where M_LCB = M_route + M_node.



  
Learnable module Finally, we use the mini-batched attention probabilities A_LCB for RouteAttn from Eq. <ref>, which takes in (1) hidden vectors of the nodes H, with size (B, N_nodes, N_hidden), (2) route information P, with size (B, N_nodes, N_nodes, F_route), and (3) masks M_node and M_route.
We use multiple heads and concatenate their results, giving us a route-based multi-head self-attention module (RouteMHSA):

    RouteMHSA
           = Concat(X_1, …, X_N_head)
        
    with 
        X_i    = RouteAttn(
            W_i,Q H,
            W_i,K H,
            W_i,V H,
           
            W_i,Q^route H,
            W_i,K^route P,
            W_i,V^route P),

where W are the learnable parameters of the module.

The values for hidden size N_hidden, key size d_k, value size d_v, route key size d_r, and the number of heads N_heads can be chosen independently of each other.
However, in our experiments we set N_hidden equal to d_k N_heads, with N_heads∈{6, 8}.
We also use d_v = d_k = d_r.















 §.§ Computational complexity

The computational complexity of the RouteMHSA is quadratic with respect to the number of nodes in the graph, because the attention mechanism computes all pairwise attention probabilities.
In practice, this means that the method can scale to graphs of a few hundred, maximally 1,000 nodes.
This is more than enough for handling drug-like compounds, which typically have less than 100 (non-hydrogen) atoms.










§ ARCHITECTURE OF THE NETWORK











The architecture of Graph Informer is inspired by transformer networks <cit.>, where the output of the multi-head attention is fed through a linear layer and added to its input, then a subsequent layer normalization <cit.> is applied. Additionally, a feedforward network (FFN) was applied to each hidden vector separately.
However, in our graph experiments, for both node-level and graph-level tasks, we found this architecture to be quite difficult to train (for more details, see Appendix <ref>).

Instead, we found that creating a residual-style network <cit.> with LayerNorm solved the gradient flow issue and was easy to train.
The layer for this architecture can be expressed as

    T     = H + LayerNorm(Linear(RouteMHSA(H)))
        
    
        H'    = T + LayerNorm(FFN(T)),

with H' being the output of the block (i.e., updated hidden vectors). The architecture is depicted in Figure <ref>.
 uses both neuron-level and channel-level dropout, with dropout rate 0.1.









§ EXPRESSIVENESS OF 

Next, we investigate the general ability of  to distinguish between distinct but similar graphs, also known as the graph isomorphism problem.



 §.§ The Weisfeiler-Lehman test


We propose a variant of our approach that is provably as powerful as the Weisfeiler-Lehman test of isomorphism of dimension 1 (WL dim=1). Generalizing the notation of <cit.> to include route information, P, we write the update function of our model as

    h_v^k = ϕ( h_v^k-1,
        {(h_u^k-1, P_v, u) : u ∈𝒩(v) }),

where  { (h_u^k-1, P_v, u) : u ∈𝒩(v) } is the multiset of node vectors, h_u^k-1, and their route information vectors P_v,u, which can be attended by node v.

Such a network is provably as powerful as the Weisfeiler-Lehman test of isomorphism if the function ϕ(.) is injective[We assume here that the output function is also injective.] with respect to the hidden vectors h_v^k-1 and h_u^k-1 arguments, as follows straightforwardly from the theorem by <cit.>.

RouteMHSA satisfies this condition with respect to the first argument (h_v^k-1) because one of the heads can attend to itself. However, to make it provably injective with respect to the second argument, we can replace the softmax mapping in Eq. <ref> by an elementwise sigmoid. We call this version the Injective .



 §.§ Beyond the Weisfeiler-Lehman test

Several graph neural networks, including GIN <cit.> and NeuralFP <cit.>, are limited in their expressiveness to the Weisfeiler-Lehman test of dim=1.
This means that if the WL test cannot distinguish two graphs, then also these methods will not be able to separate them.
By contrast, the  can aggregate route information and, therefore, has the capability to go beyond WL dim=1 as the routes between the nodes reveal nonlocal topology.
There are several types of graphs known to be indistinguishable by the WL dim=1 test. For illustrative purposes, we examine the family of regular graphs.
For an example of two non-isomorphic 6-node regular graphs, see Figure <ref>.






  §.§.§ Example: regular graphs

The WL dim=1 test fails to distinguish any regular graphs of the same size and degree because the nodes of the regular graph have all the same degree.

In the case of node-level prediction tasks, this also implies that all nodes of a regular graph get the same prediction, while their neighboring topology can be radically different.

We illustrate the capability of the  to distinguish between all regular 8-node 3-degree graphs (RegN8D3).
As a route feature, we use a histogram of the lengths of all routes between two nodes up to length k.
In other words, for a given node pair (a,b) the route feature is a vector of size k, containing the numbers of (possibly cyclic) routes connecting a and b for lengths 1, … k:

    P = Stack(A, A^2, …, A^k),

where A is the adjacency matrix and Stack creates a tensor of size N × N × k.
Computing the histogram takes k multiplications of the adjacency matrix, which gives complexity O(N^2 k).

There is a total of five RegN8D3 graphs, depicted in Appendix <ref>.
It should be noted that one of the RegN8D3 graphs corresponds to a chemically well-known compound called cubane (G4 in the appendix). While none of them are separable by the WL test or neural networks such as GIN <cit.>, GAT <cit.>, and NeuralFP <cit.>, a randomly initialized and untrained 1-layer Injective  with k=4 can separate them all, for details see the Appendix <ref>.

Additionally, we ran the same network on all regular graphs up to size 10[Graphs were taken from <http://www.mathe2.uni-bayreuth.de/markus/reggraphs.html>].
As Table <ref> shows, Graph Informer (both injective and ordinary) can distinguish them all.
This contrasts with the WL test and standard aggregating graph neural networks, which are provably unable to separate any of them (for details, see <cit.>).

We can show that with these features  has the following capability:

  Injective  using route features from Eq. <ref> with k=N can distinguish any pair of graphs with different spectra.

The proof of this lemma follows from the fact that injective  is injective with respect to route features P and that graphs are cospectral (in terms of their adjacency matrix) if and only if Tr(A^r) are equal for both matrices for all values of r <cit.>.
Notably, our method offers a more detailed view because it can also distinguish cospectral graphs (e.g., the Tesseract (Q4) and Hoffman graphs, which are both 16-node and 4-regular).




We hypothesize that in practice this capability of  gives strong additional expressiveness and makes it difficult to construct counterexamples because the route information allows the method to capture the nonlocal topology of the graph in detail.



§ RELATED RESEARCH

There are several neural network approaches that can directly work on the graph inputs.
The key component in many of the works is to propagate hidden vectors of the nodes H using the adjacency matrix[The adjacency matrix here is assumed to include the self-connection (i.e., A_ii = 1).] A, either by summing the neighbor vectors, A H, or by averaging them as Ã H, where Ã = A (D)^-1 is the normalized adjacency matrix with D being the diagonal matrix of degrees.
Neural Fingerprint <cit.> proposed to use the propagation AH, combined with a node-wise linear layer and a nonlinearity to compute fingerprints for chemical compounds.
Based on the degree of the node, Neural Fingerprint uses a different set of linear weights.
Similarly, Graph Convolutional Networks (GCN, <cit.>) propose to combine AH with a linear layer and a nonlinearity for semi-supervised learning, but by contrast with Neural Fingerprint, they use the same linear layer for all nodes, irrespective of their degree.
Kearnes et al. <cit.> proposed the Weave architecture, where they use embedding vectors also for the edges of the graph and the propagation step updates nodes and edges in a coupled way.

In a separate line of research, <cit.> proposed Graph Neural Networks (GNN), which also use AH to propagate the hidden vectors between the nodes, but instead of passing the output to the next layer, the same layer is executed until convergence, which is known as the Almeida-Pineda algorithm <cit.>.
A recent GNN-based approach is called Gated Graph Neural Networks (GGNN) <cit.>, where at each iteration the neighborhood information AH is fed into a GRU <cit.>.
An iterative graph network approach, motivated by the Message Passing algorithm in graphical model inference <cit.> was proposed by <cit.>.
Gilmer et al. <cit.> generalized the existing approaches to the Neural Message Passing (NeuralMP) algorithm and showed the effectiveness of the approach for predicting simulated quantum properties based on 3D structures of compounds.

Hamilton et al. <cit.> propose a representation learning method, GraphSAGE, that does not require supervised signals and inductively embeds nodes in large graphs.

The method closest to ours is Graph Attention Networks (GAT) <cit.>.
Similarly to us, they proposed to use multi-head attention to transform the hidden vectors of the graph nodes. However, GAT cannot integrate route features and its attention is only limited to the neighbors of every node.
For graph isomorphism testing, this means it fails to separate any regular graphs, similar to WL dim=1.
















In contrast to the spatial methods reviewed above, spectral methods use the graph Fourier basis.
One of the first spectral approaches to introduce this concept  in the context of neural networks was SCNN <cit.>.
To avoid the costly eigenvalue decomposition,  Chebyshev polynomials were proposed to define localized filters <cit.>.
Based on this work,  proposed a neural network that combines     spectral methods with metric learning for graph classification was proposed <cit.>.

For sake of completeness, we note that another loosely related line of research has been working on the problem of transforming an input graph to a target graph and are known as graph transformers <cit.>. Those works usually do not use the attention mechanism used in the transformer network <cit.>. Yet, in a recent communication, also called graph transformers, <cit.> proposed a node-level attention mechanism for the same graph-to-graph task. In contrast, our model tackles node- or graph-level classification and regression tasks, and uses route-level attention mechanism.

In graph isomorphism research, as already mentioned, <cit.> showed that injective update rule and readout are crucial to guarantee expressiveness at the level of WL dim=1.
In a recent publication, <cit.> proposed a graph neural network as powerful as WL dim=3.

In one of our main application areas, Nuclear Magnetic Resonance (NMR) chemical shift prediction (covered in Section <ref>), is an important element in structure elucidation.
Recently, <cit.> proposed an imitation learning approach to solve this NMR inverse problem.



§ EVALUATION

We compare our method against several baselines on two tasks. We first consider a node-level task where, given the graph of the compound, the goal is to predict 13C NMR peaks for all of the carbon atoms, see Section <ref>.
For the graph-level tasks, we consider the drug–protein activity data set ChEMBL <cit.> and compound toxicity prediction from MoleculeNet <cit.>, see Section <ref>.
In all tasks and for all methods, we used the same feature set (see Appendix <ref>). For methods that could not use route features but supported edge features, we used the bond type (single, double, triple, aromatic).


Each method was executed using one NVIDIA P100 with 16GB RAM with the data fed using 9 Intel(R) Xeon(R) 6140 Gold cores.
For the implementation of our algorithm, we used Pytorch.



 §.§ Baselines

Next, we overview the baselines we used for the experiments.
The code for these methods were available online, for the exact URLs, see Appendix <ref>

GGNN <cit.>:  We used the Pytorch implementation of Gated Graph Neural Networks, which was designed for node-level prediction on graphs. We also modified it to work on graph-level tasks by adding the same graph pooling layer used in our network (global average pooling).
For both tasks, we used 5 propagation steps.

NeuralFP <cit.>:  We implemented the Neural Fingerprint in Pytorch. The implementation follows the original method using softmax graph pooling.
We also extended the method to work on node-level tasks by removing the pooling and feeding the last layer hidden vectors to the regression output layer (see Section <ref>).

GCN <cit.>:  For the Graph Convolutional Networks, we adapted the Pytorch implementation to support mini-batch training.
We also adapted the code to support graph classification by incorporating the same graph pooling layer used in our network (global average pooling).

Weave <cit.>:  For the Weave network, we used the Deepchem implementation.
The implementation only supports graph-level predictions.

GAT <cit.>:  We used the Tensorflow implementation of Graph Attention Networks.
The implementation only supports node-level prediction. 

NeuralMP <cit.>:   We used the Tensorflow implementation of the neural message passing model where we used the edge network message function (best performing architecture) and added a node-level readout function for the NMR spectrum prediction (same as in all other compared methods).

GIN <cit.>   We used the Graph Isomorphism Network with the epsilon parameter set to 0, as suggested for the best performance of this architecture in their paper.



 §.§ Model selection

All methods used the validation set to choose the optimal hidden size from  {96, 144, 192, 384} and dropout from {0.0, 0.1} <cit.>. For our method, attention radius from {1,2,3,4} and head count from {6,8} were also selected.
All methods were trained with the Adam optimizer <cit.>. 
The validation set was used for early stopping (i.e., we executed all methods for 100 epochs and picked the best model).
At epochs 40 and 70, the learning rate was decreased by a factor of 0.3.
Finally, the best early stopped model was used to measure the performance on the test set.



 §.§ Node-level task: NMR 13C spectrum prediction



13C NMR is a cost-efficient spectroscopic technique in organic chemistry for molecular structure identification and structure elucidation.
If 13C NMR measurements can be combined with a highly accurate prediction, based on the graph (structure) of the molecule, it is possible to validate the structure of the sample.
It is crucial to have highly accurate predictions, because in many cases several atoms have peaks that are very close to each other.
If such predictions are available, the chemist can avoid running additional expensive and time-consuming experiments for structure validation.







To create the data set, we used 13C NMR spectra from NMRShiftDB2[Available here: <https://nmrshiftdb.nmr.uni-koeln.de/>]. After basic quality filtering (see Appendix <ref>) 25,645 molecules remained with peak positions raging from -10 to 250 ppm (parts per million) with mean and standard deviation of 95 and 52 ppm.



In accordance to the most common Lorentzian peak shape assumption in NMR spectroscopy, we minimize the Mean Absolute Error (MAE) <cit.>.
For evaluation, we split the data into train-validation-test sets containing 20,000, 2,500, and 3,145 molecules, respectively, and report results over 3 repeats with different splits.
All methods use the same node-level output head, consisting of a tanh layer followed by a linear output layer.
The MAE loss is minimized only for the atoms that have measurements for the 13C peaks (i.e., the carbon atoms).

For the results, see Table <ref>.

Our proposed  reaches 1.35 MAE, which is also better than the state-of-the-art results from computational chemistry literature, which have reported 1.43 MAE for 13C for the NMRShiftDB2 <cit.>.
At these levels of accuracy, many molecules with densely packed 13C NMR spectra can be resolved.
This class of molecules with densely packed spectra are common for aromatic ring systems, which are prevalent in drug-like compounds.
As an example, see Figure <ref> depicting 2-chloro-4-fluorophenol, which has three carbons, labelled 1, 3, 4 with peaks within the 115.2 to 116.6 ppm range.
Figure <ref> displays the predictions from the molecule when it was in the test set. The predictions match closely the true spectrum and, critically, predict the true ordering, which allows us to perfectly match peaks.
For more details and the table with exact prediction values, see Appendix <ref>.

For a  with 3 layers and 6 heads all having radius 2, we depict the attention probabilities in Figure <ref> for the bottom carbon (number 1 in Figure <ref>) of same 2-chloro-4-fluorophenol compound.
We can see that all heads in Layer 1 strongly use the pool node as a reference.
Also, we can see that in Layer 1 the first head focuses on the atom itself, the second head attends to all neighbors within radius 2, and the third head focuses on the halogen atoms.







 §.§ Results for graph-level tasks


We benchmark the methods on two drug bioactivity data sets: ChEMBL v23 by <cit.> (94 tasks, 87,415 molecules, and 349,192 activity points) and Tox21 from MoleculeNet by <cit.> (12 tasks, 7,831 molecules, and 77,946 activity points).
For ChEMBL, we used a repeated holdout-validation, and for Tox21 we used the single train-validation-test split given by MoleculeNet (see Appendix <ref> for the details).

In the case of ChEMBL, because of well-known series effect in drug design, allocating molecules randomly to the training and test sets results in overly optimistic and inconsistent performance estimates <cit.> because highly similar molecules end up in both training and test sets. Therefore, we apply clustered hold-out validation (See Appendix <ref>).  
We allocate the clusters (rather than individual molecules) to the training (60%), validation (20%), and test sets (20%). We repeat this split 5 times.



Our and all adapted methods (GCN, GIN and GGNN) use a graph pooling layer consisting of a linear layer, ReLU, mean pooling, and a linear layer.
We minimize the cross-entropy loss and Table <ref> reports the AUC-ROC averages over tasks for each data set.
 outperforms all baseline methods in both ChEMBL and Tox21 data sets.



§ CONCLUSION

In this work, we proposed a route-based multi-head attention mechanism that allows the attention to use relative position and the type of connection between the pair of nodes in the graph.
Because of its route features, our approach can incorporate information from nondirect neighbors efficiently—similarly to the use of dilated filters in CNNs.
In our theoretical analysis, we showed that a variant of  is provably as powerful as the WL dim=1 test.
Our empirical evaluation demonstrated that the proposed method is suitable both for node-level and graph-level prediction tasks, and delivers significant improvements over existing approaches in 13C NMR spectrum and drug bioactivity prediction.





























































































































§ NODE AND ROUTE FEATURES

In addition to the following features, we fed the one-hot encoding of the atom types by concatenating it to the node feature vector, see Table <ref>.
The atom types that occurred in the two data sets are {C, N, O, Cl, F, S, I, Br, P, B, Zn, Si, Li, Na, Mg, K}.
As the route features (Table <ref>) contain all edge labels (single, double, triple and aromatic bond type), all information about the graph topology is retained.







§ BASELINES USED IN EXPERIMENTS







    

    
    
    Method     Paper     Source code  

    GGNN     <cit.>     <https://github.com/JamesChuanggg/ggnn.pytorch/> 

    NeuralFP     <cit.>     Our implementation based on the paper 

    GCN     <cit.>     <https://github.com/tkipf/pygcn.git>  

    Weave     <cit.>     <https://github.com/deepchem/deepchem/> 

    GAT     <cit.>     <https://github.com/PetarV-/GAT> 

    NeuralMP     <cit.>     <https://github.com/brain-research/mpnn> 

    GIN     <cit.>     Adapted GCN code to use sum instead of 

             average and added epsilon 

    
    
    





§ PREPROCESSING OF THE CHEMBL DATA SET





 §.§ Extracting human drug–protein interaction data

We extracted IC50 values for Homo sapiens protein targets from ChEMBL version 23 SQLite release <cit.>  using the following SQL query, and computed the negative log10 of IC50 in nM, which we refer to as pIC50. Extreme values that cannot correspond to meaningful real measurement results are treated as database errors and filtered out (IC50 < 10^-5 or 10^9 ≤ IC50). The ChEMBL SQL query is shown in Table <ref>.

 











 §.§ Filtering

First, we sanitized the molecules by removing inactive ions (Cl-, Na+, K+, Br-, I-) from the structure.
Then, we removed all non-drug-like molecules containing heavy metals or more than 80 non-hydrogen atoms. We thresholded the activity for each protein with four pIC50 thresholds (5.5, 6.5, 7.5, and 8.5) taking into account the relation column.
We then kept protein–threshold pairs that had measurements at least on 2,500 compounds containing at least 200 positive and 200 negative samples.
This filtering process resulted in 94 protein–threshold pairs as binary classification tasks on 87,415 molecules. The matrix contains 349,192 activity points.



 §.§ Clustered hold-out validation


We clustered the molecules with similarity threshold 0.6 Tanimoto on ECFP6 features computed using RDKit <cit.>. For clustering, we used the leader-follower method <cit.>. We then assigned each cluster to the training, test, or validation sets randomly. We adapted the clustered validation pipeline from <cit.>.



§ TOX21 EXPERIMENTAL SETUP


We used the Tox21 split by MoleculeNet, with a total of 7,831 molecules divided into 6,264 train, 783 validation, and 784 test compounds.
Each method and each hyperparameter was trained 3 times and the best performing configuration was chosen using its validation AUC-ROC. All runs also used early stopping with validation AUC-ROC.



§ NMR DATA SET


To create the data set, we used 13C NMR spectra from NMRShiftDB2, available from <https://nmrshiftdb.nmr.uni-koeln.de/>.
We kept only molecules where all carbon atoms where assigned a 13C peak position and removed any compounds that failed RDKit <cit.> sanitization or Gasteiger partial charge calculation.
This resulted in 25,645 molecules.
The peak positions are measured in parts per million (ppm) and ranges from approximately -10 to 250 ppm. Note that the ppm value measures the shift of the peak relative to a specific reference carbon atom.
In our data set, the mean and standard deviation are 95 and 52 ppm.



§ DESCRIPTION OF THE NMR APPLICATION



In molecules, the chemical environment of the atomic nuclei has an influence on their measured resonance frequency, the NMR peak position.
This influence is called the chemical shift, and measured in parts per million (ppm).
The nominal value of this shift is independent from the measurement equipment and can be used to compare chemical environments between molecules measured in different laboratories.

We are specifically interested in the following scenario. An organic chemist intends to synthesize a compound. After the synthesis and purification steps, the identity of the product needs to be validated.
The chemist can measure the 13C NMR spectrum of the product, and compare it to the predicted spectrum from the structure of this compound.
To compare the two spectra, the correct peaks corresponding to the same atom need to be aligned.
There is no way to tell from a 13C NMR spectrum alone which peak corresponds to which atom. Therefore, wrongly ordered peaks in the prediction would result in an erroneous assignment.

As in the example in the main text, we use the 13C NMR spectrum of the compound 2-chloro-4-fluorophenol to illustrate the assignment problem.
The measured and predicted chemical shifts are reported in Table <ref>.
As can be seen in Figure <ref>, the three unsubstituted aromatic carbons (1,3,4) have very close peaks in the 115.2 to 116.6 ppm region. 
For the sake of completeness note that the intensity of the peaks (value of the y-axis) in the plot are directly related to the number of direct hydrogen neighbors of the carbon, and therefore does not need to be predicted <cit.>. 

 




§ IMPROVING THE CONVERGENCE OF 


Precisely following the strategy used for transformer networks <cit.>, we get the following architecture for graphs: 

    T     = LayerNorm(H + Linear(RouteMHSA(H)))
        
    
        H'    = LayerNorm(T + FFN(T)),
        
    FFN(x)    = W_2 ReLU(W_1 x + b_1) + b_2,

and H' is the output of the block (i.e., updated hidden vectors).

In our experiments, we observed that the output of RouteMHSA(H) dominates the input H in Eq. <ref>.
This causes convergence issues as normalizing this sum by the LayerNorm causes the gradient through the input H branch to vanish.
This was the main motivation for introducing the pure-residual style network described in Section <ref>.



§ ISOMORPHISM TEST


Randomly initialized untrained 1-layer injective  was used to separate regular graphs when using histogram of route counts as route features.
For the histogram, we used routes up to length (k=4) (i.e., the feature vector for each pair of nodes contains the counts of routes of length 1, 2, 3, and 4).
After the embeddings for all nodes have been computed the embedding for the graph is their sum. 
Finally the two graphs were considered different if their embeddings differed more than 10^-4.
We used the following setup:







 §.§ Regular 8-node 3-degree graphs



All 5 regular 8-node 3-degree graphs are depicted in Figure <ref>.
These graphs are not distinguishable by the Weisfeiler-Lehman test with dim=1 and consequently by many graph neural networks.










§ ATTENTION PLOTS









§ ACKNOWLEDGMENT

YM is funded by Research Council KU Leuven: C14/18/092 SymBioSys3; CELSA-HIDUCTION CELSA/17/032
Flemish Government:IWT: Exaptation, PhD grants
FWO 06260 (Iterative and multi-level methods for Bayesian multirelational factorization with features). This research received funding from the Flemish Government under the “Onderzoeksprogramma Artificiële Intelligentie (AI) Vlaanderen” program. EU: “MELLODDY” This project has received funding from the Innovative Medicines Initiative 2 Joint Undertaking under grant agreement No 831472. This Joint Undertaking receives support from the European Union’s Horizon 2020 research and innovation program and EFPIA. EDB is funded by a FWO-SB grant.































IEEEtran












Michael Shell
Biography text here.



John Doe
Biography text here.






Jane Doe
Biography text here.



















