# MPVNN: Mutated Pathway Visible Neural Network architecture for interpretable
prediction of cancer-specific survival risk

Gourab Ghosh Roy
School of Computer Science
University of Birmingham
Birmingham, B15 2TT, UK
School of Computing and Information Systems
University of Melbourne
Melbourne, 3052, Australia
<EMAIL_ADDRESS>
&Nicholas Geard
School of Computing and Information Systems
University of Melbourne
Melbourne, 3052, Australia
<EMAIL_ADDRESS>
&Karin Verspoor
School of Computing Technologies
RMIT University
Melbourne, 3000, Australia
<EMAIL_ADDRESS>
&Shan He
School of Computer Science
University of Birmingham
Birmingham, B15 2TT, UK
<EMAIL_ADDRESS>

Survival risk prediction using gene expression data is important in making
treatment decisions in cancer. Standard neural network (NN) survival analysis
models are black boxes with lack of interpretability. More interpretable
visible neural network (VNN) architectures are designed using biological
pathway knowledge. But they do not model how pathway structures can change for
particular cancer types. We propose a novel Mutated Pathway VNN or MPVNN
architecture, designed using prior signaling pathway knowledge and gene
mutation data-based edge randomization simulating signal flow disruption. As a
case study, we use the PI3K-Akt pathway and demonstrate overall improved
cancer-specific survival risk prediction results of MPVNN over standard non-NN
and other similar sized NN survival analysis methods. We show that trained
MPVNN architecture interpretation, which points to smaller sets of genes
connected by signal flow within the PI3K-Akt pathway that are important in
risk prediction for particular cancer types, is reliable.

##  1 Introduction

Cancer is a leading cause of death worldwide, and a substantial amount of
medical research is focused on survival analysis of cancer patients. The
suitability of a particular treatment method could be guided by the predicted
survival risk of the patient. The effectiveness of a treatment method could
also be measured using the predicted risk. Gene expression data has been
extensively used for cancer risk prediction, and different machine learning
methods have been applied for this survival analysis task by learning from
survival data [1]. Neural networks (NNs) are an effective category of machine
learning methods which have been used for gene expression-based cancer risk
prediction [2, 3].

One major challenge of some machine learning models including standard NNs is
lack of interpretability. Model interpretability refers to the degree to which
the model’s internal operations can be understood by a human [4], which here
denotes being understood in biological terms. This understanding would point
to how particular biological entities and relationships are used internally by
the model in mapping the input to the output. A more interpretable model would
inherently increase a user’s trust on the model, and hence be considered more
reliable for use in high-stakes clinical applications [5].

Visible NNs or VNNs are neural networks where biological meanings are attached
to intermediate neurons, and these VNNs provide increased model
interpretability compared to standard black box NNs [6]. Under a visible
architecture, neurons represent biological entities like genes, proteins,
pathways, cell subsystems, etc. and connections between the neurons represent
biological relationships. Here we refer to those NNs as VNNs where neurons in
each intermediate layer are associated with some explicit biological meaning.

An important source of biological knowledge for use in the design of VNNs is
that of biological pathways. VNNs whose design uses pathway knowledge with
neurons representing pathways include P-NET [7], GenNet [8] and DCell [9].
These VNNs are designed from hierarchical knowledge and have multiple
intermediate layers. Another VNN designed using pathway knowledge is
knowledge-primed NN or KPNN [10], which is a NN architecture analogous to the
structure of signaling pathways and also has multiple intermediate layers. A
VNN with one intermediate layer, designed using protein–protein (PPI) and
protein–DNA (PDI) interactions data, for predicting cell type from single cell
expression values has been explored [11]. However, none of these VNNs model
how a known biological pathway structure can change for a particular disease.

In this study, we propose a Mutated Pathway Visible Neural Network or MPVNN
architecture for cancer-specific survival risk prediction from gene expression
data. This is a simple VNN architecture, with one intermediate layer between
the input genotype and the output phenotype layers. Each neuron in the single
intermediate layer is assigned to represent perturbation at a gene. Here a
gene is considered to be perturbed if it is in the path of signal flow. A gene
perturbation is derived from the expression of the same gene and the
expression of the genes it shares edges with. Initially these edges are
signaling pathway edges from prior knowledge. However, the signal flow through
the known pathway edges can be disrupted in cancer. Simulating the change in
pathway structure for a particular cancer type, we use additional gene
mutation data for that cancer to replace a certain fraction of known pathway
edges with random gene connections. These new connections are aimed at
capturing signaling interactions not present in the used prior knowledge,
which can be important in survival risk prediction for that particular cancer
type. So the connections between the input layer and the intermediate layer
neurons are used to represent the signaling edges, either obtained from prior
knowledge or from the randomization using mutation data. Finally, from these
gene perturbations, the cancer-specific survival risk is predicted in the
MPVNN output.

Using prior knowledge of PI3K-Akt signaling pathway and mutation data in the
design, we assess the cancer-specific risk prediction performance of our
proposed MPVNN architecture on data for 10 cancer types from The Cancer Genome
Atlas or TCGA [12]. We compare its prediction performance against one standard
non-NN and other similar sized NN survival analysis methods. We show how
interpretation of trained MPVNN architecture points to signal flow involving
gene sets within the larger PI3K-Akt pathway which are important in cancer-
specific risk prediction. We also investigate how the interpreted signal flow
can be through some connections not present in the used PI3K-Akt pathway prior
knowledge. To assess the reliability of the MPVNN interpreted insights for two
cancer types, ovarian and liver cancer, we use validation from prior studies
in literature. All of this is aimed at evaluating the suitability of MPVNN for
use as a cancer-specific survival risk predictor in clinical applications.

Our main contributions in this paper are summarized as follows:

  1. 1.

A novel Mutated Pathway Visible Neural Network or MPVNN architecture designed
using signaling pathway knowledge and mutation data-based edge randomization
is proposed for cancer-specific survival risk prediction.

  2. 2.

Using the PI3K-Akt signaling pathway as a case study, we demonstrate the
overall improved risk prediction performance on real-world cancer datasets for
MPVNN over standard and other comparable NN survival analysis methods.

  3. 3.

We also show how MPVNN interpretation can provide insights on signal flow
important in risk prediction for particular cancer types, and assess the
reliability of some such insights using evidence from literature.

##  2 Methods

###  2.1 Problem

We address the survival analysis task as a regression problem of predicting
the survival risk. There are other approaches of addressing survival analysis.
For example, in a classification setting, the patients can be stratified into
risk categories [13]. However, such a stratification, based on a threshold
like the median survival time, might not be accurate in a patient population
with a mix of cancer stages corresponding to widely varying medians. In
another analysis approach, time-to-event distributions can be estimated [14].
Here our aim is to predict risk scores indicative of the survival times so
that ordering is correct. A clinician might be interested only in knowing if a
particular treatment causes an increase in survival time without the need to
know the exact value of the said time.

The survival data comprises two major components : the $i^{th}$ patient
survival time to event of interest $t_{i}$, and the censoring indicator
$l_{i}$ where a value of $1$ denotes that the event was observed and a value
of $0$ denotes censoring. In our risk prediction task, the objective is to
predict a relative risk score $f(x_{i})$ as output given the expression
profile $x_{i}\in\mathbb{R}^{N}$ consisting of expression values of $N$ genes
as input, where a larger survival time can be associated with a larger risk
score for cases that can be ordered.

Here we have used a standard survival analysis performance assessment metric
of the Concordance Index (CI) or c-index [15]. It is a measure of agreement
between the predicted survival risk and the observed survival. The CI is
defined as follows

$\text{CI}=\frac{1}{|\epsilon|}\sum_{(i,j)\in\epsilon}{1_{f(x_{i})<f(x_{j})}},\vspace*{-2pt}$
(1)

where $\epsilon=\\{(i,j)\hskip 2.84544pt|\hskip 2.84544ptl_{i}=1\hskip
2.84544pt\text{and}\hskip 2.84544ptt_{j}>t_{i}\\}$. The CI metric is a value
between $0$ and $1$, where a value of 0.5 denotes random prediction. The
objective here is to maximize the CI value. Since the CI itself is not
differentiable, we consider the following differentiable exponential lower
bound on the CI [16]:

$\frac{1}{|\epsilon|}\sum_{(i,j)\in\epsilon}{1-e^{f(x_{i})-f(x_{j})}}.\vspace*{-2pt}$
(2)

In our regression setting, we use the negative of this lower bound as the loss
function for minimization. As in [17], during optimization we ignore the
denominator of Eq. 2, and evaluate the loss over training batches.

The patient cancer survival outcome endpoint used in our experiments is the
disease-specific survival (DSS) endpoint [18]. A DSS event denotes death
specifically from the diagnosed cancer type, and the event time is measured
from the date of initial diagnosis until that of event. The censored time
denotes the time from the date of initial diagnosis until the date of death
due to another cause or the date of last contact. This DSS endpoint is
difficult to derive and is an approximation of the true DSS [18]. However, we
want to be able to predict the cancer-specific survival risk, which has a more
direct relationship with the cancer modeling. So we have selected DSS over a
more commonly used endpoint – the overall survival (OS), and used the TCGA
datasets for our experiments which provide this DSS data for multiple cancer
types.

###  2.2 Proposed Architecture

To solve this problem of risk prediction from gene expression profiles, we
propose a simple VNN architecture with one intermediate layer connecting the
input genotype layer and the output phenotype layer. The proposed MPVNN
architecture is presented in Fig. 11. The input and the intermediate layers
have $N$ neurons each, and the output layer has one. The connections between
the input and the intermediate layers is guided by the design algorithm in
Algorithm 1.

Figure 1: Proposed MPVNN architecture. Each intermediate layer neuron is
assigned to represent perturbation at a gene, which is derived from its own
expression and the expression of genes which are connected by either a pathway
or a randomly connected edge. For example, such an edge $G1\rightarrow G2$
between two genes leads to two NN connections, one from input layer $G1$
expression neuron to intermediate layer $G2$ perturbation neuron, and second
from $G2$ expression neuron to $G1$ perturbation neuron.

1:Input: List of all $N$ genes $N_{all}$; List of all $P$ pathways $P_{all}$,
Gene list $N_{p}$ & Edge list $E_{p}$ for every pathway $p\in P_{all}$;
Mutation data $M\in\\{0,1\\}^{K\times N}$;

2:Output: Connection matrix $W\in\\{0,1\\}^{N\times N}$ denoting connections
between input and intermediate layers

3:foreach gene $Gn\in N_{all}$ do

4: Number of mutated samples $M^{sum}_{Gn}\leftarrow\sum_{k=1}^{K}{M_{k,n}}$

5:end foreach

6:Most mutated genes $M^{mut}\leftarrow\\{Gm\hskip 2.84544pt\forall\hskip
2.84544ptM^{sum}_{Gm}\geq 0.99$ quantile of all non-zero $M^{sum}_{Gn}\\}$

7:foreach pathway $p\in P_{all}$ do

8: $frac_{p}\leftarrow\frac{|N_{p}|}{N}$, $thr_{p}\leftarrow\frac{|N_{p}\cap
M^{mut}|}{|M^{mut}|}$

9:end foreach

10:Connection matrix $W\leftarrow$ createW($N_{all}$, $N_{p}$, $E_{p}$,
$thr_{p}$, $frac_{p}$ for $p\in P_{all}$)

11:

12:function createW($N_{all}$, $N_{p}$, $E_{p}$, $thr_{p}$, $frac_{p}$ for
$p\in P_{all}$)

13: $W\leftarrow I_{N}$

14: $E_{new}\leftarrow\emptyset$

15: foreach pathway $p\in P_{all}$ do

16: foreach edge $a\rightarrow b\in E_{p}$ do

17: $rnum_{1}\leftarrow random()$

18: if $rnum_{1}<thr_{p}$ then

19: do

20: $rnum_{2}\leftarrow random()$

21: if $rnum_{2}\geq thr_{p}$ then

22: $a,b\leftarrow$ Select 2 unique genes randomly from $N_{p}$

23: else if $rnum_{2}<thr_{p}\frac{frac_{1\neq p}}{1-frac_{p}}$ then

24: if $random()\geq thr_{1\neq p}$ then

25: $a,b\leftarrow$ 2 unique genes randomly from $N_{1\neq p}$

26: else

27: $a,b\leftarrow$ 2 unique genes randomly from $N_{all}-N_{1\neq p}$

28: end if

29: ⋮

30: else if $rnum_{2}<thr_{p}\frac{frac_{1\neq p}+\cdots+frac_{P\neq
p}}{1-frac_{p}}$ then

31: if $random()\geq thr_{P\neq p}$ then

32: $a,b\leftarrow$ 2 unique genes randomly from $N_{P\neq p}$

33: else

34: $a,b\leftarrow$ 2 unique genes randomly from $N_{all}-N_{P\neq p}$

35: end if

36: else

37: $a,b\leftarrow$ 2 unique genes randomly from $N_{all}-\sum_{p\in
P_{all}}N_{p}$

38: end if

39: while $(a\rightarrow b\in\cup_{p\in P_{all}}E_{p})\lor(a\rightarrow b\in
E_{new})\lor(b\rightarrow a\in\cup_{p\in P_{all}}E_{p})\lor(b\rightarrow a\in
E_{new})$

40: Add $a\rightarrow b$ in $E_{new}$

41: end if

42: $W_{ab}=1$

43: $W_{ba}=1$

44: end foreach

45: end foreach

46: return $W$

47:end function

Algorithm 1 MPVNN design algorithm

In the architecture, each intermediate layer neuron representing the
perturbation at one gene, is connected to the input layer neuron representing
the expression of the same gene. When using just pathway knowledge and no
mutation data, each intermediate layer gene perturbation neuron is
additionally connected to other input layer neurons representing the
expression of genes which are its known neighbors in a signaling pathway. So a
gene perturbation is derived from its own expression and the expression of its
pathway neighbors, which are the genes with edges to or from the gene in
consideration. This VNN design is unchanged by autoregulatory edges, or by
cases where an edge is shared by two or more pathways or edges exist between
two genes in both directions within or across pathways. This is the PVNN
architecture, which is the version of the MPVNN architecture without the
mutation data-based edge randomization.

The MPVNN architecture is an extension of the PVNN architecture using another
type of data, non-silent gene mutation data. For a cancer type, we compute a
mutation fraction per signaling pathway, which is a rough estimate of how much
the signal flow in the pathway is disrupted. With the mutation fraction, we
randomly replace this same fraction of the known pathway edges with new edges
and connect the neurons accordingly. In this case, the intermediate layer
neurons can derive perturbations at a gene from its own expression, and
possibly from the expression of its neighbors as per known pathway knowledge,
and possibly from the expression of other genes that it is randomly connected
with, based on the mutation data, as per the design algorithm.

For the MPVNN architecture design, the input is the list of all $N$ genes and
$P$ pathways, the gene and the edge lists for every individual pathway and the
mutation data for the $N$ genes. We assume that over all the pathway edge
lists, no autoregulatory edge exists, an edge and its reversed edge together
do not exist, and an edge belongs to only one pathway. This keeps the number
of connections in PVNN and MPVNN same for fair performance comparison. The
design algorithm output is the matrix $W\in\\{0,1\\}^{N\times N}$ denoting
which connections exist between the N input layer gene expression neurons and
the N intermediate layer gene perturbation neurons. From the mutation data, we
compute the number of samples in which each gene $Gn$ is mutated –
$M^{sum}_{Gn}$. Then we obtain $M^{mut}$ consisting of most mutated genes,
each of whose above number is greater than or equal to the $0.99$ quantile of
all non-zero $M^{sum}_{Gn}$. For each pathway $p$ with $frac_{p}$ fraction of
total $N$ genes, we compute a mutation fraction $thr_{p}$ equal to the
fraction of genes from this pathway in $M^{mut}$. Here we assume that the
input gene expression profile does not consist only of genes all belonging to
one single pathway, or we would use the PVNN architecture instead. For each
known pathway edge, if a generated random number $rnum_{1}\in[0,1)$ is $\geq
thr_{p}$, we select the known edge, otherwise the edge is replaced and a new
edge is selected.

For replacing a known edge with a new edge, another random number
$rnum_{2}\in[0,1)$ is generated. If this number is $\geq thr_{p}$, we select
two pathway $p$ genes randomly for a new edge. Otherwise, based on the value
of the random number $rnum_{2}$, each remaining pathway $q\neq p$ is selected
for consideration $thr_{p}\frac{frac_{q}}{1-frac_{p}}$ fraction of times. If
the pathway $q$ is selected for consideration, based on whether another new
random number is $\geq thr_{q}$, two pathway $q$ genes are selected randomly
for a replacing edge, otherwise two genes are selected randomly from all genes
not in pathway $q$. This is motivated by the design goal that pathways that
have a high fraction of genes in most mutated genes and probably have greater
disruption of signal flow, should contribute less to the genes of replacing
edges. When no such pathway is selected dictated by the value of $rnum_{2}$,
two genes, not belonging to any pathway in the list of pathways used in the
design, are selected randomly for a replacing edge.

While obtaining a replacing edge in the MPVNN design, we only select two
unique genes, and ensure that the selected replacing edge or its reversed edge
do not belong either to the list of all known pathway edges or to the list of
all already selected replacing edges. This keeps the number of connections of
MPVNN same with that of PVNN, and also can help in finding important signaling
connections novel to what is present in the used prior knowledge. The
connection matrix $W$ is initially an identity matrix, and two entries where
the row and column interchangeably denote the two genes of every selected
edge, are made $1$.

Though MPVNN architecture is generalized to use prior knowledge and mutation
data for multiple pathways, here as a case study, we use one well-known cancer
related pathway–phosphatidylinositol 3’ -kinase (PI3K)-Akt signaling pathway
in its design. This pathway is a key regulator of processes involved in cell
growth, proliferation, survival and apoptosis. It has been observed to play
critical roles in various cancers [19]. So we have assessed the effectiveness
of our proposed architecture designed with the PI3K-Akt pathway obtained from
the KEGG Pathway database [20].

###  2.3 Interpretation

Compared to standard black box NN survival analysis methods, a major benefit
of our proposed VNN architecture is increased interpretability. Both the
trained VNN architectures – PVNN and MPVNN can be interpreted to obtain top
gene sets. These sets of genes within the larger signaling pathway are linked
by flow of signal that is important in the prediction of survival risk for a
particular cancer type, and are ranked in order of the importance associated
with a set. The benefit of MPVNN over PVNN is that from interpretation we can
obtain important signaling connections which are not present in the used prior
knowledge.

Here we describe a straightforward mean weight amplitude-based VNN
interpretation method, which fits with the signaling edge-based design of our
VNNs. We first obtain the top gene perturbations in the intermediate layer,
from the absolute weights connecting the intermediate layer gene perturbation
neurons to the output neuron, averaged over the VNN runs. From each top gene
perturbation, we find the top gene set connected by flow of signal associated
with the gene perturbations.

The first candidate gene is the one with the top gene perturbation, and then
the following process is repeated for every new gene added to the top gene
set. For every candidate gene, we evaluate which of the gene expression
neurons in the input layer, apart from the candidate gene itself and the last
gene added to the top gene set in consideration, has the highest mean absolute
weight connecting to the candidate gene perturbation neuron in the
intermediate layer. We check whether the highest weight is $\geq 0.85$
quantile of all the non-zero weights between the input layer and the
intermediate layer neurons. If this new gene already belongs to the top gene
set in consideration, we stop the process. Otherwise, we check whether the
weight connecting the expression neuron for the newly selected gene to its
perturbation neuron in the intermediate layer satisfies the weight threshold
above. We additionally check if the perturbation neuron for the newly selected
gene to the output layer phenotype neuron is $\geq 0.85$ quantile of all the
non-zero weights between the intermediate layer and the output layer neurons.
Then this new gene is added to the top gene set. This whole process is done
twice starting from the top gene perturbation to be able to possibly capture
the important signal flow into and out of the gene.

The thresholds used are applied to ensure that genes with low connection
weights and hence lower importance are not included in the top gene sets. We
used a value of $0.85$ quantile, which for the intermediate layer neurons,
roughly translates to that out of the 1440 gene perturbations, around 200 can
play some role of importance in risk prediction for a cancer type. The top
gene sets obtained from the interpretation would depend on the values of the
thresholds used.

##  3 Experiments

Our experimental TCGA data for 10 cancer types is obtained from the UCSC Xena
browser [21]. The data includes gene expression RNAseq data, binary gene-level
non-silent mutation data and disease-specific survival data. We consider
cancer samples, one per patient, which do not fall under the normal category
in TCGA. For the expression profile, apart from the genes in the PI3K-Akt
pathway, we also consider 1285 genes found to have systematic expression
change in cancer [22]. Finally our patient expression profile consists of 1440
genes, for which both expression and mutation data are available.

The cancer types used in our experiments are presented in detail in Table 1.
For each cancer type, the input-output data is split after a random shuffle
into training ($80\%$) and test ($20\%$) sets with stratification on the
censoring indicator/ event observed values. For the MPVNN architecture design,
the available mutation data is used as a whole. For every machine learning
method in our experiments, the input data features are standardized by mean
subtraction and scaling to unit variance, based on those metrics computed from
the training data.

Table 1: Our experimental datasets. Cancer |  Description |  Number of samples |  Events observed |  Censored
---|---|---|---|---
BLCA |  Bladder urothelial carcinoma |  390 |  119 |  271
BRCA |  Breast invasive carcinoma |  1070 |  83 |  987
COADREAD |  Colon and rectum adenocarcinoma |  354 |  41 |  313
GBM |  Glioblastoma multiforme |  147 |  113 |  34
HNSC |  Head and neck squamous cell carcinoma |  493 |  130 |  363
KIRC |  Kidney renal clear cell carcinoma |  522 |  109 |  413
LIHC |  Liver hepatocellular carcinoma |  361 |  79 |  282
LUNG |  Lung squamous cell carcinoma and adenocarcinoma |  914 |  199 |  715
OV |  Ovarian serous cystadenocarcinoma |  272 |  150 |  122
STAD |  Stomach adenocarcinoma |  386 |  95 |  291

We evaluated the performance of MPVNN against other comparable NN
architectures. First there is the PVNN architecture which is designed using
pathway knowledge and no mutation data-based randomization. In RaNN which is
another randomized version of PVNN, the connections between the hidden layer
and the input layer neurons are randomly shuffled. This represents a same
sized NN with the same number of connections but designed without using prior
knowledge or additional mutation data. We evaluated the performance of the
fully connected artificial neural network or ANN with the same number of
neurons. We have also compared the performance of these NN methods with a
standard survival analysis method – the semi-parametric Cox Proportional
Hazards (Cox-PH) model [23]. We drop the expression values of those genes
which give an initial warning of having very low variance in the Cox PH
regression fitter.

We used a 4-fold cross-validation on the training set to identify the optimal
hyperparameters. The optimal hyperparameter setting for a particular method
for a particular cancer type is selected from a certain number of random
searches on the list of given hyperparameter values. For this number, we used
around $10\%$ of the total number of possible combinations. So the number of
random searches was set to be 300 for NNs and 10 for Cox-PH model. The list of
values to choose from for NNs is given as follows :

  1. 1.

Activation function : tanh, sigmoid, ReLU

  2. 2.

Optimizer : Adam, SGD

  3. 3.

Learning Rate : 0.1, 0.01, 0.001

  4. 4.

Batch Size : 16, 32, 64, 128

  5. 5.

Epoch : 10, 25, 50, 100

  6. 6.

Regularizer : L1, L2

  7. 7.

Regularizer parameter : 0.0, 1e-5, 1e-3, 1e-2, 1e-1, 1.0.

The hyperparameters for the Cox-PH model to select from are penalizer values
in the set

$\\{0.0,0.0001,0.0005,0.001,0.005,0.01,0.05,0.1,0.5,1.0\\}$

and the L1 ratio values between $0$ and $1$ in intervals of $0.1$. For
randomized architectures MPVNN and RaNN, the optimal connection matrix is also
selected from the same 300 random searches.

The performance metric used is the CI value on the hold-out test set. We
present the mean and the standard deviation from 20 runs. The data and code
are available at https://github.com/gourabghoshroy/MPVNN.

##  4 Results

In Table 2 we present the results of our MPVNN architecture for cancer-
specific survival risk prediction, for each cancer type, and finally these
values macro-averaged over all cancers. We also show each method’s macro-
average weighted rank, where the weighted rank is calculated as the difference
between the maximum of mean CI metric values for all methods and the mean CI
metric value for a particular method in a cancer type. The lower the macro-
average weighted rank, the closer to the best mean performance a method’s mean
performance is on an average across cancer types.

Table 2: Cancer-specific survival risk prediction performance evaluation of MPVNN. The mean and standard deviation of the CI metric from 20 runs are shown for each cancer type, and then these values are macro-averaged over all cancer types. The macro-average weighted ranks are also presented. The mean best, that is the method with the highest mean value, individually and macro-averaged, and the one with the lowest macro-average weighted rank, are all marked in boldface. Cancer |  Cox-PH |  ANN |  RaNN |  PVNN |  MPVNN
---|---|---|---|---|---
BLCA |  0.6919$\pm$0.0000 |  0.7455$\pm$0.0107 |  0.7463$\pm$0.0046 |  0.7130$\pm$0.0089 |  0.7214$\pm$0.0029
BRCA |  0.4810$\pm$0.0142 |  0.6507$\pm$0.0174 |  0.5434$\pm$0.0082 |  0.6367$\pm$0.0151 |  0.6298$\pm$0.0040
COADREAD |  0.5868$\pm$0.0000 |  0.7731$\pm$0.0673 |  0.7658$\pm$0.0474 |  0.6085$\pm$0.0230 |  0.7804$\pm$0.0267
GBM |  0.5822$\pm$0.0041 |  0.5153$\pm$0.0521 |  0.5762$\pm$0.0421 |  0.4991$\pm$0.0168 |  0.5763$\pm$0.0115
HNSC |  0.4912$\pm$0.0006 |  0.5864$\pm$0.0092 |  0.5920$\pm$0.0091 |  0.6498$\pm$0.0025 |  0.5827$\pm$0.0081
KIRC |  0.7747$\pm$0.1159 |  0.8374$\pm$0.0257 |  0.8474$\pm$0.0063 |  0.7955$\pm$0.0049 |  0.8082$\pm$0.0063
LIHC |  0.5549$\pm$0.0042 |  0.7527$\pm$0.0292 |  0.7600$\pm$0.0060 |  0.7918$\pm$0.0052 |  0.7878$\pm$0.0050
LUNG |  0.6644$\pm$0.0000 |  0.5837$\pm$0.0124 |  0.5649$\pm$0.0035 |  0.6248$\pm$0.0031 |  0.6160$\pm$0.0057
OV |  0.5895$\pm$0.0305 |  0.5886$\pm$0.0334 |  0.5869$\pm$0.0130 |  0.6102$\pm$0.0068 |  0.6248$\pm$0.0050
STAD |  0.7429$\pm$0.0000 |  0.6207$\pm$0.0516 |  0.6699$\pm$0.0087 |  0.5502$\pm$0.0395 |  0.6800$\pm$0.0062
MACRO-AVERAGE |  0.6160$\pm$0.0170 |  0.6654$\pm$0.0309 |  0.6653$\pm$0.0149 |  0.6480$\pm$0.0126 |  0.6807$\pm$0.0081
MACRO-AVERAGE WEIGHTED RANK |  0.0921 |  0.0427 |  0.0428 |  0.0601 |  0.0273

Our results show that the overall mean cancer-specific risk prediction
performance of MPVNN is the best among the methods that were tested. Firstly,
it is much much better than the standard Cox-PH survival analysis model macro-
averaged across all cancer types. Second, its overall performance is better
than those of other comparable NN architectures – PVNN, RaNN and ANN. Compared
to the next best method, MPVNN on macro-average has a mean CI metric higher by
$0.015$. Even though MPVNN is not the best in all of the 10 cancer types, its
macro-averaged metrics demonstrate that the incorporation of signaling pathway
knowledge and gene mutation data-based edge randomization in the VNN design
can improve the overall prediction performance.

To demonstrate the interpretability of our proposed MPVNN architecture, we
have shown two top gene sets obtained from MPVNN interpretation in Fig. 22.
The top gene sets are given for 2 cancer types, ovarian and liver, in which
the VNNs outperform other methods (MPVNN performs the best and the second
best). We have selected the top gene set which has the first occurrence of any
signaling connection that is not present in the used PI3K-Akt pathway edge
list. For ovarian cancer, the top gene set consists of these genes :
GNB3–PPP2R1B–FGF2$\rightarrow$NGFR–PPP2R2B$\rightarrow$AKT1$\rightarrow$CHUK,
where the top gene perturbation is marked in bold and the novel signaling
interactions are marked by dashes. The shown top gene set for liver cancer
comprises ANGPT2$\rightarrow$FLT3–FLT4$\leftarrow$ANGPT2. An important novel
connection between genes FLT4 and FLT3 in the same gene group RTK is
interpreted from the MPVNN architecture, however these 2 genes are not
connected by any given PI3K-Akt signaling path.

Figure 2: Top gene set from MPVNN interpretation for Ovarian (OV) and Liver
(LIHC) cancers. Ovarian :
GNB3–PPP2R1B–FGF2$\rightarrow$NGFR–PPP2R2B$\rightarrow$AKT1$\rightarrow$CHUK,
Liver : ANGPT2$\rightarrow$FLT3–FLT4$\leftarrow$ANGPT2. The gene group in the
PI3K-Akt pathway diagram [20] for each gene in the top gene set is bordered by
a colored rectangle (Ovarian - Red, Liver - Blue). The novel connections which
are not present in the pathway knowledge are shown by colored curved lines.

##  5 Discussion

From our results, we demonstrate that the overall cancer-specific survival
risk prediction mean performance of our proposed MPVNN architecture is
substantially better than other survival analysis methods. However, the mean
performance of MPVNN is not the best in every single cancer type, given the
corresponding data. The Cox-PH model is found to be a better risk prediction
model than the NNs in some cancer types, though MPVNN is observed to be much
better on macro-averaged metrics. The high importance of the PI3K-Akt pathway
in risk prediction is highlighted for those particular cancer types where PVNN
performs the best. We also observe that for some cancers, the randomized RaNN
or the fully connected ANN gives the best mean CI metrics. Since these two
architectures lack interpretability, it is more difficult to understand in
biological terms how they work. When RaNN is the best, one possibility is that
there are other pathways involving the input genes, which are more effective
in risk prediction. Incorporating knowledge and mutation data for multiple
pathways in the architecture design might be helpful for performance. As
described previously, our VNN architectures can be applied to use multiple
pathways. For the cancer where ANN is the best, one possibility is that
pathways which can play significant roles in risk prediction are not well
represented in the input dataset genes, and a larger set of input genes, along
with multiple pathway knowledge and mutation data in the architecture design,
might improve the performance.

A major benefit of our MPVNN architecture is increased interpretability
compared to black box NN methods, making it more reliable for use in clinical
survival analysis. To validate some of the insights we obtained from MPVNN
interpretation, we looked for literature evidence. In the top gene set for
ovarian cancer, we observed important signaling interactions within the
PI3K-Akt pathway between parts of the growth factor, the chemokine, the Akt
and the NF$\kappa$B signaling pathways. Some of these interactions are novel
to what exists in the used prior knowledge. Signal flow involving these four
pathways together being important in the context of higher mortality in
ovarian cancer has been previously indicated [24, 25].

In the top gene set in liver cancer, we observed important pathway edges
ANGPT2$\rightarrow$FLT3 and ANGPT2$\rightarrow$FLT4 and a novel interaction
between RTK genes FLT3 and FLT4. Connection to ANGPT2 from a RTK gene via
calcineurin-NFAT has been studied in lung metastasis [26]. So such a signaling
path from any of the RTK genes FLT3 or FLT4 back to ANGPT2 can correspond to
the interpreted novel interaction. Both the FLT3 and FLT4 genes are among the
targets of the popular drug Sorafenib, which has been used effectively in the
treatment of advanced liver cancer for over a decade [27]. Interestingly, in
an extensive Sorafenib trial [28], ANGPT2 was found to be an independent
survival predictor in the overall and the drug treated population, but not
related with drug induced survival benefit. Compared to the baseline
concentration, ANGPT2 concentration after 12 weeks of treatment was observed
to increase in the placebo group, however it remained almost unchanged for
treatment with the drug. So it is possible that the signal flow in the MPVNN
interpreted top gene set is disrupted during treatment with Sorafenib by
targeting both FLT3 and FLT4 genes, as a result of which the ANGPT2
concentration remains almost constant. Later the signal flow involving these
genes probably again comes into effect.

Based on the above literature validation, we argue that MPVNN interpreted
insights are reliable, pointing to signal flow that have critical roles in
controlling cancer-specific survival risk. Further validation of these
insights would require experimental verification, which is outside the scope
of this study. These insights are flexible enough to include important
signaling connections that are not present in the used prior knowledge. MPVNN
interpreted insights can have correspondence with the mode of action of
existing drugs, and can provide directions for novel and more effective single
or combinatorial drug therapy.

Edge randomization based on additional mutation data, which models how a
pathway structure can change for a particular cancer type and replaces known
edges with new connections, is a novelty of our MPVNN architecture.
Heteroscedastic dropout with privileged or additional information is used in
[29], although this differs substantially from MPVNN. Dropout for a NN
architecture is used only during learning and not in prediction, whereas the
edges in the MPVNN architecture are randomized at the beginning, before
learning and subsequent prediction. Also in MPVNN, the use of the additional
mutation data in the supervised prediction setting is unsupervised, without
the requirement of being related to the input-output training data. This makes
the architecture more robust to bias that may arise in supervised machine
learning when some population groups are not well represented in the training
set, for instance due to lack of survival data here, but mutation data is
available for those groups and can be incorporated in the model.

As only the PI3K-Akt pathway used in this study is not sufficient in fully
demonstrating MPVNN predictive power, future work would look at conducting
larger scale experiments with larger gene expression profiles using knowledge
and data for multiple pathways in the MPVNN architecture design to further
evaluate cancer-specific survival risk prediction performance. In future, we
would also like to experimentally investigate the roles of the important
signal flow identified from MPVNN interpretation for different cancer types,
and importantly explore their relevance to cancer drug target identification.

## Funding

Gourab Ghosh Roy is supported by a Priestley Scholarship for joint study at
the University of Birmingham and the University of Melbourne.

## References

  * [1] Wei-Yi Cheng, Tai-Hsien Ou Yang, and Dimitris Anastassiou.  Development of a prognostic model for breast cancer survival in an open challenge environment.  Science translational medicine, 5(181):181ra50–181ra50, 2013.
  * [2] Jared L Katzman, Uri Shaham, Alexander Cloninger, Jonathan Bates, Tingting Jiang, and Yuval Kluger.  Deepsurv: personalized treatment recommender system using a cox proportional hazards deep neural network.  BMC medical research methodology, 18(1):1–12, 2018.
  * [3] Zhi Huang, Xiaohui Zhan, Shunian Xiang, Travis S Johnson, Bryan Helm, Christina Y Yu, Jie Zhang, Paul Salama, Maher Rizkalla, Zhi Han, et al.  Salmon: survival analysis learning with multi-omics neural networks on breast cancer.  Frontiers in genetics, 10:166, 2019.
  * [4] Or Biran and Courtenay Cotton.  Explanation and justification in machine learning: A survey.  In IJCAI-17 workshop on explainable AI (XAI), volume 8, pages 8–13, 2017.
  * [5] Cynthia Rudin.  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  Nature Machine Intelligence, 1(5):206–215, 2019.
  * [6] K Yu Michael, Jianzhu Ma, Jasmin Fisher, Jason F Kreisberg, Benjamin J Raphael, and Trey Ideker.  Visible machine learning for biomedicine.  Cell, 173(7):1562–1565, 2018.
  * [7] Haitham A Elmarakeby, Justin Hwang, Rand Arafeh, Jett Crowdis, Sydney Gang, David Liu, Saud H AlDubayan, Keyan Salari, Steven Kregel, Camden Richter, et al.  Biologically informed deep neural network for prostate cancer discovery.  Nature, 598(7880):348–352, 2021.
  * [8] Arno van Hilten, Steven A Kushner, Manfred Kayser, M Arfan Ikram, Hieab HH Adams, Caroline CW Klaver, Wiro J Niessen, and Gennady V Roshchupkin.  Gennet framework: interpretable neural networks for phenotype prediction.  bioRxiv, 2020.
  * [9] Jianzhu Ma, Michael Ku Yu, Samson Fong, Keiichiro Ono, Eric Sage, Barry Demchak, Roded Sharan, and Trey Ideker.  Using deep learning to model the hierarchical structure and function of a cell.  Nature methods, 15(4):290, 2018.
  * [10] Nikolaus Fortelny and Christoph Bock.  Knowledge-primed neural networks enable biologically interpretable deep learning on single-cell sequencing data.  Genome biology, 21(1):1–36, 2020.
  * [11] Chieh Lin, Siddhartha Jain, Hannah Kim, and Ziv Bar-Joseph.  Using neural networks for reducing the dimensions of single-cell rna-seq data.  Nucleic acids research, 45(17):e156–e156, 2017.
  * [12] Katarzyna Tomczak, Patrycja Czerwińska, and Maciej Wiznerowicz.  The cancer genome atlas (tcga): an immeasurable source of knowledge.  Contemporary oncology, 19(1A):A68, 2015.
  * [13] Yen-Chen Chen, Wan-Chi Ke, and Hung-Wen Chiu.  Risk classification of cancer survival using ann with gene expression data from multiple laboratories.  Computers in biology and medicine, 48:1–7, 2014.
  * [14] Paidamoyo Chapfuwa, Chenyang Tao, Chunyuan Li, Courtney Page, Benjamin Goldstein, Lawrence Carin Duke, and Ricardo Henao.  Adversarial time-to-event modeling.  In International Conference on Machine Learning, pages 735–744. PMLR, 2018.
  * [15] Frank E Harrell Jr, Kerry L Lee, and Daniel B Mark.  Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors.  Statistics in medicine, 15(4):361–387, 1996.
  * [16] Harald Steck, Balaji Krishnapuram, Cary Dehing-Oberije, Philippe Lambin, and Vikas C Raykar.  On ranking in survival analysis: Bounds on the concordance index.  In Advances in neural information processing systems, pages 1209–1216, 2008.
  * [17] Ellery Wulczyn, David F Steiner, Zhaoyang Xu, Apaar Sadhwani, Hongwu Wang, Isabelle Flament-Auvigne, Craig H Mermel, Po-Hsuan Cameron Chen, Yun Liu, and Martin C Stumpe.  Deep learning-based survival prediction for multiple cancer types using histopathology images.  PLoS One, 15(6):e0233678, 2020.
  * [18] Jianfang Liu, Tara Lichtenberg, Katherine A Hoadley, Laila M Poisson, Alexander J Lazar, Andrew D Cherniack, Albert J Kovatich, Christopher C Benz, Douglas A Levine, Adrian V Lee, et al.  An integrated tcga pan-cancer clinical data resource to drive high-quality survival outcome analytics.  Cell, 173(2):400–416, 2018.
  * [19] Ningni Jiang, Qijie Dai, Xiaorui Su, Jianjiang Fu, Xuancheng Feng, and Juan Peng.  Role of pi3k/akt pathway in cancer: the framework of malignant behavior.  Molecular biology reports, 47(6):4587–4629, 2020.
  * [20] Minoru Kanehisa and Susumu Goto.  Kegg: kyoto encyclopedia of genes and genomes.  Nucleic acids research, 28(1):27–30, 2000.
  * [21] Mary J Goldman, Brian Craft, Mim Hastie, Kristupas Repečka, Fran McDade, Akhil Kamath, Ayan Banerjee, Yunhai Luo, Dave Rogers, Angela N Brooks, et al.  Visualizing and interpreting cancer genomics data via the xena platform.  Nature biotechnology, 38(6):675–678, 2020.
  * [22] Aurora Torrente, Margus Lukk, Vincent Xue, Helen Parkinson, Johan Rung, and Alvis Brazma.  Identification of cancer related genes using a comprehensive map of human gene expression.  PloS one, 11(6):e0157484, 2016.
  * [23] David R Cox.  Regression models and life-tables.  Journal of the Royal Statistical Society: Series B (Methodological), 34(2):187–202, 1972.
  * [24] Yuan-Lin Dong, Syeda M Kabir, Eun-Sook Lee, and Deok-Soo Son.  Cxcr2-driven ovarian cancer progression involves upregulation of proinflammatory chemokines by potentiating nf-$\kappa$b activation via egfr-transactivated akt signaling.  PloS one, 8(12):e83789, 2013.
  * [25] Deok-Soo Son, Syeda M Kabir, Yuanlin Dong, Eunsook Lee, and Samuel E Adunyah.  Characteristics of chemokine signatures elicited by egf and tnf in ovarian cancer cells.  Journal of inflammation, 10(1):1–12, 2013.
  * [26] Takashi Minami, Shuying Jiang, Keri Schadler, Jun-ichi Suehiro, Tsuyoshi Osawa, Yuichi Oike, Mai Miura, Makoto Naito, Tatsuhiko Kodama, and Sandra Ryeom.  The calcineurin-nfat-angiopoietin-2 signaling axis in lung endothelium is critical for the establishment of lung metastases.  Cell reports, 4(4):709–723, 2013.
  * [27] Giorgia Marisi, Alessandro Cucchetti, Paola Ulivi, Matteo Canale, Giuseppe Cabibbo, Leonardo Solaini, Francesco G Foschi, Serena De Matteis, Giorgio Ercolani, Martina Valgiusti, et al.  Ten years of sorafenib in hepatocellular carcinoma: Are there any predictive and/or prognostic markers?  World journal of gastroenterology, 24(36):4152, 2018.
  * [28] Josep M Llovet, Carol EA Peña, Chetan D Lathia, Michael Shan, Gerold Meinhardt, and Jordi Bruix.  Plasma biomarkers as predictors of outcome in patients with advanced hepatocellular carcinoma.  Clinical Cancer Research, 18(8):2290–2300, 2012.
  * [29] John Lambert, Ozan Sener, and Silvio Savarese.  Deep learning under privileged information using heteroscedastic dropout.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 8886–8895, 2018.


# Motif-based Graph Self-Supervised Learning for Molecular Property Prediction

Zaixi Zhang1, Qi Liu1∗, Hao Wang1, Chengqiang Lu1, Chee-Kong Lee2
1: Anhui Province Key Lab of Big Data Analysis and Application,
School of Computer Science and Technology,
University of Science and Technology of China
2: Tencent America
{zaixi, wanghao3<EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS>Qi Liu and Chee-Kong Lee are corresponding authors.

###### Abstract

Predicting molecular properties with data-driven methods has drawn much
attention in recent years. Particularly, Graph Neural Networks (GNNs) have
demonstrated remarkable success in various molecular generation and prediction
tasks. In cases where labeled data is scarce, GNNs can be pre-trained on
unlabeled molecular data to first learn the general semantic and structural
information before being finetuned for specific tasks. However, most existing
self-supervised pre-training frameworks for GNNs only focus on node-level or
graph-level tasks. These approaches cannot capture the rich information in
subgraphs or graph motifs. For example, functional groups (frequently-occurred
subgraphs in molecular graphs) often carry indicative information about the
molecular properties. To bridge this gap, we propose Motif-based Graph Self-
supervised Learning (MGSSL) by introducing a novel self-supervised motif
generation framework for GNNs. First, for motif extraction from molecular
graphs, we design a molecule fragmentation method that leverages a
retrosynthesis-based algorithm BRICS and additional rules for controlling the
size of motif vocabulary. Second, we design a general motif-based generative
pre-training framework in which GNNs are asked to make topological and label
predictions. This generative framework can be implemented in two different
ways, i.e., breadth-first or depth-first. Finally, to take the multi-scale
information in molecular graphs into consideration, we introduce a multi-level
self-supervised pre-training. Extensive experiments on various downstream
benchmark tasks show that our methods outperform all state-of-the-art
baselines.

##  1 Introduction

The past decade has witnessed the remarkable success of deep learning in
natural language processing (NLP) [5], computer vision (CV) [13], and graph
analysis [19]. Inspired by these developments, researchers in the chemistry
domain have also tried to exploit deep learning methods for molecule-based
tasks such as retrosynthesis [47] and drug discovery [8]. To preserve the
internal structural information, molecules can be naturally modeled as graphs
where nodes represent atoms and edges denote the chemical bonds. Recently,
some works applied Graph Neural Network (GNN) and some of its variants for
molecular property prediction and obtained promising results [8, 25].

Though GNNs have achieved remarkable accuracy on molecular property
prediction, they are usually data-hungry, i.e. a large amount of labeled data
(i.e., molecules with known property data) is required for training [8, 12].
However, labeled molecules only occupy an extremely small portion of the
enormous chemical space since they can only be obtained from wet-lab
experiments or quantum chemistry calculations, which are time-consuming and
expensive. Moreover, directly training GNNs on small labeled molecule datasets
in a supervised fashion is prone to over-fitting and the trained GNNs can
hardly generalize to out-of-distribution data.

Similar issues have also been encountered in natural language processing and
computer vision. Recent advances in NLP and CV address them by self-supervised
learning (SSL) where a model is first pre-trained on a large unlabeled dataset
and then transferred to downstream tasks with limited labels [13, 2, 5]. For
example, the pre-trained BERT language model [5] is able to learn expressive
contextualized word representations through reconstructing the input text-next
sentence and masked language predictions so that it can significantly improve
the performance of downstream tasks. Inspired by these developments, various
self-supervised pre-training methods of GNNs have been proposed [14, 15, 52,
51, 32, 10, 30]. Based on how the self-supervised tasks are constructed, these
methods can be classified into two categories, contrastive methods and
predictive methods. Contrastive methods force views from the same graph (e.g.,
sampling nodes and edges from graphs) to become closer and push views from
different graphs apart. On the other hand, the predictive methods construct
prediction tasks by exploiting the intrinsic properties of data. For example,
Hu et.al. [14] designed node-level pre-training tasks such as predicting the
context of atoms and the attributes of masked atoms and bonds. [15] introduced
an attributed graph reconstruction task where the generative model predicts
the node attributes and edges to be generated at each step.

However, we argue that existing self-supervised learning tasks on GNNs are
sub-optimal since most of them fail to exploit the rich semantic information
from graph motifs. Graph motifs can be defined as significant subgraph
patterns that frequently occur [27]. Motifs usually contain semantic meanings
and are indicative of the characteristics of the whole graph. For example, the
hydroxide (–OH) functional group in small molecules typically implies higher
water solubility. Therefore, it is vital to design motif-level self-supervised
learning tasks which can benefit downstream tasks such as molecular property
prediction.

Designing motif-level self-supervised learning tasks brings unique challenges.
First, existing motif mining techniques could not be directly utilized to
derive expressive motifs for molecular graphs because they only rely on the
discrete count of subgraph structures and overlook the chemical validity [1,
18]. Second, most graph generation techniques generate graphs node-by-node
[23, 50], which are not suitable for our task of motif generation. Finally,
how to unify multi-level self-supervised pre-training tasks harmoniously
brings a new challenge. One naive solution to do pre-training tasks
sequentially may lead to catastrophic forgetting similar to continual learning
[7].

To tackle the aforementioned challenges, we propose Motif-based Graph Self-
Supervised Learning (MGSSL) and Multi-level self-supervised pre-training in
this paper. Firstly, MGSSL introduces a novel motif generation task that
empowers GNNs to capture the rich structural and semantic information from
graph motifs. To derive semantically meaningful motifs and construct motif
trees for molecular graphs, we leverage the BRICS algorithm [4] which is based
on retrosynthesis from the chemistry domain. Two additional fragmentation
rules are further introduced to reduce the redundancy of motif vocabulary.
Secondly, a general motif-based generative pre-training framework is designed
to generate molecular graphs motif-by-motif. The pre-trained model is required
to make topology and attribute predictions at each step and two specific
generation orders are implemented (breadth-first and depth-first).
Furthermore, to take the multi-scale regularities of molecules into
consideration, we introduce Multi-level self-supervised pre-training for GNNs
where the weights of different SSL tasks are adaptively adjusted by the Frank-
Wolfe algorithm [16]. Finally, by pre-training GNNs on the ZINC dataset with
our methods, the pre-trained GNNs outperforms all the state-of-the-art
baselines on various downstream benchmark tasks, demonstrating the
effectiveness of our design. The implementation is publicly available at
https://github.com/zaixizhang/MGSSL.

##  2 Preliminaries and Related Work

###  2.1 Molecular Property Prediction

Prediction of molecular properties is a central research topic in physics,
chemistry, and materials science [44]. Among the traditional methods, density
functional theory (DFT) is the most popular one and plays a vital role in
advancing the field [21]. However, DFT is very time-consuming and its
complexity could be approximated as $\mathcal{O}(N^{3})$, where $N$ denotes
the number of particles. To find more efficient solutions for molecular
property prediction, various machine learning methods have been leveraged such
as kernel ridge regression, random forest, and convolutional neural networks
[6, 26, 36]. To fully consider the internal spatial and distance information
of atoms in molecules, many works [8, 35, 48, 25] regard molecules as graphs
and explore the graph convolutional network for property prediction. To better
capture the interactions among atoms, [8] proposes a message passing framework
and [20, 48] extend this framework to model bond interactions. [25] builds a
hierarchical GNN to capture multilevel interactions. Furthermore, [32]
integrates GNNs into the Transformer-style architecture to deliver a more
expressive model.

###  2.2 Preliminaries of Graph Neural Networks

Recent years have witnessed the success of Graph Neural Networks (GNNs) for
modeling graph data [11, 19, 40, 46, 43, 42, 55]. Let $G=(V,E)$ denotes a
graph with node attributes $X_{v}$ for $v\in V$ and edge attributes $e_{uv}$
for $(u,v)\in E$. GNNs leverage the graph connectivity as well as node and
edge features to learn a representation vector (i.e., embedding) $h_{v}$ for
each node $v\in G$ and a vector $h_{G}$ for the entire graph $G$. Generally,
GNNs follows a message passing paradigm, in which representation of node $v$
is iteratively updated by aggregating the representations of $v$’s neighboring
nodes and edges. Specifically, there are two basic operators for GNNs:
AGGREGATE(·) and COMBINE(·). AGGREGATE(·) extracts the neighboring information
of node $v$; COMBINE(·) serves as the aggregation function of the neighborhood
information. After $k$ iterations of aggregation, $v$’s representation
captures the structural information within its $k$-hop network neighborhood.
Formally, the $k$-th layer of a GNN is:

$h_{v}^{(k)}=\textup{COMBINE}^{(k)}\left(h_{v}^{(k-1)},\textup{AGGREGATE}^{(k)}\left(\left\\{(h_{v}^{(k-1)},h_{u}^{(k-1)},e_{uv}):u\in\mathcal{N}(v)\right\\}\right)\right),$
(1)

where $h_{v}^{(k)}$ denotes the representation of node $v$ at the $k$-th
layer, $e_{uv}$ is the feature vector of edge between $u$ and $v$ and and
$\mathcal{N}(v)$ represents the neighborhood set of node $v$. $h_{v}^{(0)}$ is
initialized with $X_{v}$. Furthermore, to obtain the representation of the
entire graph $h_{G}$, READOUT functions are designed to pool node
representations at the final iteration $K$,

$h_{G}=\textup{READOUT}\left(\left\\{h_{v}^{(K)}|~{}v\in G\right\\}\right).$
(2)

READOUT is a permutation-invariant function, such as averaging, sum, max or
more sophisticated graph-level pooling functions [49, 53].

###  2.3 Self-supervised Learning of Graphs

Graph Self-supervised learning aims to learn the intermediate representations
of unlabeled graph data that are useful for unknown downstream tasks. [14, 15,
52, 51, 32, 10, 30, 39, 28, 41]. Traditional graph embedding methods [9, 29]
define different kinds of graph proximities, i.e., the vertex proximity, as
the self-supervised objective to learn vertex representations. Furthermore,
[39, 28, 41] proposed to use the mutual information maximization as the
optimization objective for GNNs. Recently, more self-supervised tasks for GNNs
have been proposed [14, 15, 52, 51, 32, 10, 54, 30]. Based on how the self-
supervised tasks are constructed, these models can be classified into two
categories, namely contrastive models and predictive models. Contrastive
models try to generate informative and diverse views from data instances and
perform node-to-context [14], node-to-graph [39] or motif-to-graph [54]
contrastive learning. On the other hand, predictive models are trained in a
supervised fashion, where the labels are generated based on certain properties
of the input graph data, i.e., node attributes [32], or by selecting certain
parts of the graph [15, 32].

However, few works try to leverage the information of graph motifs for graph
self-supervised learning. Grover [32] use traditional softwares to extract
motifs and treat them as classification labels. MICRO-Graph [54] exploits
graph motifs for motif-graph contrastive learning. Unfortunately, both methods
fail to consider the topology information of motifs.

##  3 Motif-based Graph Self-supervised Learning

In this section, we introduce the framework of motif-based graph self-
supervised learning (Figure 1). Generally, the framework consists of three
parts: chemistry-inspired molecule fragmentation, motif generation and multi-
level self-supervised pre-training for GNNs.

###  3.1 Chemistry-inspired Molecule Fragmentation

Given a dataset of molecules, the first step of our method is to decompose
molecules into several fragments/motifs. Based on these fragments, a molecule
graph can be converted into a motif tree structure where each node represents
a motif and the edges denotes the relative spatial relationships among motifs.
We choose to represent the structure of motifs as a tree because it
facilitates the motif generation task. Formally, given a molecule graph
$G=(V,E)$, a motif tree $\mathcal{T}(G)=(\mathcal{V},\mathcal{E},\mathcal{X})$
is a connected labeled tree whose node set is
$\mathcal{V}=\\{M_{1},...,M_{n}\\}$ and edge set is $\mathcal{E}$.
$\mathcal{X}$ refers to the induced motif vocabulary. Each motif
$M_{i}=(V_{i},E_{i})$ is a subgraph of $G$. There are many ways to fragment a
given graph while the designed molecule fragmentation method should achieve
the following goals: 1) In a motif tree $\mathcal{T}(G)$, the union of all
motifs $M_{i}$ should equals $G$. Formally, $\bigcup_{i}V_{i}=V$ and
$\bigcup_{i}E_{i}\bigcup\mathcal{E}=E$. 2) In a motif tree $\mathcal{T}(G)$,
the motifs should have no intersections. That is $M_{i}\cap M_{j}=\emptyset$.
3) The induced motifs should capture semantic meanings, e.g., similar to
meaningful functional groups in the chemistry domain. 4) The occurrence of
motifs in the dataset should be frequent enough so that the pre-trained GNNs
can learn semantic information of motifs that can be generalized to downstream
tasks. After the molecule fragmentation, the motif vocabulary of the molecule
dataset should have a moderate size.

In Figure 2, we show the overview of molecule fragmentation. Generally, there
are three procedures, BRICS fragmentation, further decomposition, and motif
tree construction. A motif vocabulary can be built via preprocessing the whole
molecule dataset following the molecule fragmentation procedures.

To fragment molecule graphs and construct motif trees, we firstly use the
Breaking of Retrosynthetically Interesting Chemical Substructures (BRICS)
algorithm [4] that leverages the domain knowledge from chemistry. BRICS
defines 16 rules and breaks strategic bonds in a molecule that match a set of
chemical reactions. "Dummy" atoms are attached to each end of the cleavage
sites, marking the location where two fragments can join together. BRICS
cleavage rules are designed to retain molecular components with valuable
structural and functional content, e.g. aromatic rings.

However, we find BRICS alone cannot generate desired motifs for molecule
graphs. This is because BRICS only breaks bonds based on a limited set of
chemical reactions and tends to generate several large fragments for a
molecule. Moreover, due to the combination explosion of graph structure, we
find BRICS produces many motifs that are variations of the same underlying
structure (e.g., Furan ring with different combinations of Halogen atoms). The
motif vocabulary is large (more than 100k unique fragments) while most of
these motifs appear less than 5 times in the whole dataset.

To tackle the problems above, we introduce a post-processing procedure to
BRICS. To alleviate the combination explosion, we define two rules operating
on the output fragments of BRICS: (1) Break the bond where one end atom is in
a ring while the other end not. (2) Select non-ring atoms with three or more
neighboring atoms as new motifs and break the neighboring bonds. The first
rule reduces the number of ring variants and the second rule break the side
chains. Experiments show that these rules effectively reduce the size of motif
vocabulary and improve the occurrence frequency of motifs.

.

Figure 1: Illustration of Motif-based Graph Self-supervised learning (MGSSL).
The multi-level pre-training consists of two layers, Atom layer and Motif
layer. In the Atom layer, we mask node/edge attributes and let GNNs predict
those attributes based on neighboring structures. In the Motif layer, we
construct motif trees and perform motif generative pre-training. In each step,
based on existing motifs and connections, topology and motif predictions are
made iteratively.

###  3.2 Motif Generation

Here, we present the framework of motif generation for the generative self-
supervised pre-training. The goal of the motif generation task is to let GNNs
learn the data distribution of graph motifs so that the pre-trained GNNs can
easily generalize to downstream tasks with several finetuning steps on graphs
from similar domains.

Given a molecule graph $G=(V,E)$ and a GNN model $f_{\theta}$, we first
convert the molecule graph to a motif tree
$\mathcal{T}(G)=(\mathcal{V},\mathcal{E},\mathcal{X})$. Then we can model the
likelihood over this motif tree by the GNN model as
$p(\mathcal{T}(G);\theta)$, representing how the motifs are labeled and
connected. Generally, our method aims to pretrain the GNN model $f_{\theta}$
via maximizing the likelihood of motif trees, i.e., $\theta^{*}={\rm
argmax}_{\theta}~{}p(\mathcal{T}(G);\theta)$. To model the likelihood of motif
trees, special predictions heads for topology and motif label predictions are
designed (as shown in the following sections) and optimized along with
$f_{\theta}$. After pre-training, only the GNN model $f_{\theta}$ is
transferred to downstream tasks.

We note that most existing works on graph generation [17, 23] follow the auto-
regressive manner to factorize the probability objective, i.e.,
$p(\mathcal{T}(G);\theta)$ in this work. For each molecular graph, they
decompose it into a sequence of generation steps. Similarly in this paper, we
interleave the addition of a new motif, and the addition of bonds to connect
the newly added motif to the existing partial motif tree. We denote a
permutation vector $\pi$ to determine the motif ordering, in which $i^{\pi}$
denotes the motif ID of $i$-th position in permutation $\pi$. Therefore, the
probability $p(\mathcal{T}(G);\theta)$ is equivalent to the expected
likelihood over all possible permutations, i.e.,

$p(\mathcal{T}(G);\theta)=\mathbb{E}_{\pi}\left[p_{\theta}(\mathcal{V}^{\pi},\mathcal{E}^{\pi})\right],$
(3)

where $\mathcal{V}^{\pi}$ denotes the permuted motif labels and
$\mathcal{E}^{\pi}$ denotes the edges among motifs.

Figure 2: Overview of molecule fragmentation. Generally, there are three
steps: 1) Firstly a molecule graph is cleaved based on BRICS. 2) Further
decomposition to reduce the redundancy of motifs 3) Construct motif trees from
molecule graphs. A motif vocabulary is built after preprocessing the whole
molecule dataset.

Our formalism permits a variety of orders. For simplicity, we assume that any
node ordering $\pi$ has an equal probability and we also omit the subscript
$\pi$ when illustrating the generative process for one permutation in the
following sections. Given a permutation order, the probability of generating
motif tree $\mathcal{T}(G)$ can be decomposed as follows:

${\rm
log}~{}p_{\theta}(\mathcal{V},\mathcal{E})=\sum_{i=1}^{|\mathcal{V}|}{\rm
log}~{}p_{\theta}(\mathcal{V}_{i},\mathcal{E}_{i}~{}|~{}\mathcal{V}_{<i},\mathcal{E}_{<i}).$
(4)

At each step $i$, we use motif attributes $\mathcal{V}_{<i}$ and structures
$\mathcal{E}_{<i}$ of all motifs generated before $i$ to generate a new motif
$\mathcal{V}_{i}$ and its connection with existing motifs $\mathcal{E}_{i}$.

Equation 4 describes the autoregressive generative process of motif trees.
Then the question is how to choose an efficient generation order and how to
model the conditional probability ${\rm
log}~{}p_{\theta}(\mathcal{V}_{i},\mathcal{E}_{i}~{}|~{}\mathcal{V}_{<i},\mathcal{E}_{<i})$.
In the following sections, we introduce two efficient generation orders,
breadth-first (BFS) and depth-first (DFS) orders and show the corresponding
auto-regressive generative models.

To generate a motif tree from scratch, we need to first choose the root of
motif tree. In our experiments, we simply choose the motif with the first atom
in the canonical order [34]. Then MGSSL generates motifs in DFS or BFS orders
(see Figure 3). In the DFS order, for every visited motif, MGSSL first makes a
topological prediction: whether this node has children to be generated. If a
new child motif node is generated, we predict its label and recurse this
process. MGSSL backtracks when there is no more children to generate. As for
the BFS order, MGSSL generates motif nodes layer-wise. For motifs nodes in the
$k$-th layer, MGSSL makes topological predictions and label predictions. If
all the children nodes of motifs in the $k$-th layer are generated, MGSSL will
move to the next layer. We also note that the motif node ordering in BFS and
DFS are not unique as the order within sibling nodes is ambiguous. In the
experiments, we pre-train in one order and leave this issue and other
potential generation orders for future exploration.

At each time step, a motif node receives information from other generated
motifs for making those predictions. The information is propagated through
message vectors $h_{ij}$ when motif trees are incrementally constructed.
Formally, let $\mathcal{\hat{E}}_{t}$ be the set of message at time $t$. The
model visits motif $i$ at time $t$ and $x_{i}$ denotes the embedding of motif
$i$, which can be obtained by pooling the atom embeddings in motif $i$. The
message ${\rm h}_{i,j}$ is updated through previous messages:

${\rm h}_{i,j}={\rm GRU}\left(x_{i},~{}\\{{\rm
h}_{k,i}\\}_{(k,i)\in\mathcal{\hat{E}}_{t},k\neq j}\right),$ (5)

where GRU is Gated Recurrent Unit [3] adapted for motif tree message passing:

$\displaystyle s_{i,j}$
$\displaystyle=\sum_{(k,i)\in\mathcal{\hat{E}}_{t},k\neq j}{\rm h}_{k,i}$ (6)
$\displaystyle z_{i,j}$ $\displaystyle=\sigma({\rm W}^{z}x_{i}+{\rm
U}^{z}s_{i,j}+b^{z})$ (7) $\displaystyle r_{k,i}$ $\displaystyle=\sigma({\rm
W}^{r}x_{i}+{\rm U}^{r}{\rm h}_{k,i}+b^{r})$ (8) $\displaystyle\tilde{\rm
h}_{i,j}$ $\displaystyle={\rm tanh}({\rm
W}x_{i}+U\sum_{k=\mathcal{N}(i)\backslash j}r_{k,i}\odot{\rm h}_{k,i})$ (9)
$\displaystyle{\rm h}_{i,j}$ $\displaystyle=(1-z_{ij})\odot
s_{ij}+z_{ij}\odot\tilde{\rm h}_{i,j}.$ (10)

Topology Prediction: When MGSSL visits motif $i$, it needs to make binary
predictions on whether it has children to be generated. We compute the
probability via a one hidden layer network followed by a sigmoid function
taking messages and motif embeddings into consideration:

$p_{t}=\sigma\left(U^{d}\cdot\tau(W_{1}^{d}x_{i}+W_{2}^{d}\sum_{(k,i)\in\mathcal{\hat{E}}_{t}}h_{k,i})\right),$
(11)

where $d$ is the dimension of the hidden layer.

.

Figure 3: Illustration of the motif generation orders. The first row is the
illustration of DFS order and the second row is the BFS order.

Motif Label Prediction: When motif $i$ generate a child motif $j$, we predict
the label of child $j$ with:

$q_{j}={\rm softmax}(U^{l}\tau(W^{l}h_{ij})),$ (12)

where $q_{j}$ is the distribution over the motif vocabulary $\mathcal{X}$ and
$l$ is the hidden layer dimension. Let $\hat{p}_{t}\in\\{0,1\\}$ and
$\hat{q}_{j}$ be the ground truth topological and motif label values, the
motif generation loss is the sum of cross-entropy losses of topological and
motif label predictions:

$\mathcal{L}_{motif}=\sum_{t}\mathcal{L}_{topo}(p_{t},\hat{p}_{t})+\sum_{j}\mathcal{L}_{pred}(q_{j},\hat{q}_{j}).$
(13)

In the optimization process, minimizing the above loss function corresponds to
maximizing the log likelihood in equation 4. Note that in the training
process, after topological and motif label prediction at each step, we replace
them with their ground truth so that the MGSSL makes predictions based on
correct histories.

###  3.3 Multi-level Self-supervised Pre-training

To capture the multi-scale information in molecules, MGSSL is designed to be a
hierarchical framework including Atom-level and Motif-level tasks (Figure 1).
For Atom-level pre-training, we leverage attribute masking to let GNNs firstly
learn the regularities of the node/edge attributes. In attribute masking,
randomly sampled node and bond attributes (e.g., atom numbers, bond types) are
replaced with special masked indicators. Then we apply GNNs to obtain the
corresponding node/edge embeddings (edge embeddings can be obtained as a
combination of node embeddings of the edge’s end nodes). Finally, a fully
connected layer on top of the embeddings predicts the node/edge attributes.
The cross-entropy prediction losses are denoted as $\mathcal{L}_{atom}$ and
$\mathcal{L}_{bond}$ respectively.

To avoid catastrophic forgetting in sequential pre-training, we unify multi-
level tasks and aim to minimize the hybrid loss in pre-training:

$\mathcal{L}_{ssl}=\lambda_{1}\mathcal{L}_{motif}+\lambda_{2}\mathcal{L}_{atom}+\lambda_{3}\mathcal{L}_{bond},$
(14)

where $\lambda_{i}$ are weights of losses. However, it is time-consuming to do
a grid search to determine the optimal weights. Here, we adapts the MGDA-UB
algorithm [37] from multi-task learning to efficiently solve the optimization
problem (Equation 14). Since MGDA-UB calculates the weights $\lambda_{i}$ by
Frank-Wolfe algorithm [16] at each training step, we do not have to give
weights explicitly. The pseudo codes of the training process is included in
the Appendix.

##  4 Experimental Results

###  4.1 Experimental Settings

Datasets and Dataset Splittings. In this paper, we mainly focus on the
molecular property prediction tasks, where large-scale unlabeled molecules are
abundant whereas downstream labeled data is scarce. Specifically, we use 250k
unlabeled molecules sampled from the ZINC15 database [38] for self-supervised
pre-training tasks. As for the downstream finetune tasks, we consider 8 binary
classification benchmark datasets contained in MoleculeNet [45]. The detailed
dataset statistics are summarized in the Appendix. We use the open-source
package RDKit [22] to preprocess the SMILES strings from various datasets. To
mimic the real-world use case, we split the downstream dataset by scaffold-
split [14, 31], which splits the molecules according to their structures. We
apply 3 independent runs on random data splitting and report the means and
standard deviations.

Baselines. We comprehensively evaluate the performance of MGSSL against five
state-of-the-art self-supervised pre-training methods for GNNs:

  * •

Deep Graph Infomax [41] maximizes the mutual information between the
representations of the whole graphs and the representations of its sampled
subgraphs.

  * •

Attribute masking [14] masks node/edge features and let GNNs predict these
attributes.

  * •

GCC [30] designs the pretraining task as discriminating ego-networks sampled
from a certain node ego-networks sampled from other nodes.

  * •

Grover [32] predicts the contextual properties based on atom embeddings to
encode contextual information into node embeddings.

  * •

GPT-GNN [15] is a generative pretraining task which predicts masked edges and
node attributes.

In experiments, we also consider GNN without pre-training (direct supervised
finetuning) and MGSSL with different generation orders (BFS and DFS).

Model Configuration. In the following experiments, we select a five-layer
Graph Isomorphism Networks (GINs) as the backbone architecture, which is one
of the state-of-the-art GNN methods. Mean pooling is used as the Readout
function for GIN. In MGSSL, sum pooling is used to get the embedding of graph
motifs. Atom number and chirality tag are input as node features and bond type
and direction are regarded as edge features. In the process of pre-training,
GNNs are pre-trained for 100 epochs with Adam optimizer and learning rate
0.001. In the finetuning stage, we train for 100 epochs and report the testing
score with the best cross-validation performance. The hidden dimension is set
to 300 and the batch size is set to 32 for pre-training and finetuning. The
split for train/validation/test sets is $80\%:10\%:10\%$. All experiments are
conducted on Tesla V100 GPUs.

Table 1: Test ROC-AUC (%) performance on molecular property prediction benchmarks using different pre-training strategies with GIN. The rightmost column averages the mean of test performance across the 8 datasets. The best result for each dataset are bolded. SSL methods | muv | clintox | sider | hiv | tox21 | bace | toxcast | bbbp | Avg.
---|---|---|---|---|---|---|---|---|---
No pretrain |  71.7$\pm 2.3$ |  58.2$\pm 2.8$ |  57.2$\pm 0.7$ |  75.4$\pm 1.5$ |  74.3$\pm 0.5$ |  70.0$\pm 2.5$ |  63.3$\pm 1.5$ |  65.5$\pm 1.8$ | 67.0
Infomax |  75.1$\pm 2.8$ |  73.0$\pm 3.2$ |  58.2$\pm 0.5$ |  76.5$\pm 1.6$ |  75.2$\pm 0.3$ |  75.6$\pm 1.0$ |  62.8$\pm 0.6$ |  68.1$\pm 1.3$ | 70.6
Attribute masking |  74.7$\pm 1.9$ |  77.5$\pm 3.1$ |  59.6$\pm 0.7$ |  77.9$\pm 1.2$ | 77.2$\mathbf{\pm 0.4}$ |  78.3$\pm 1.1$ |  63.3$\pm 0.8$ |  65.6$\pm 0.9$ | 71.8
GCC |  74.1$\pm 1.4$ |  73.2$\pm 2.6$ |  58.0$\pm 0.9$ |  75.5$\pm 0.8$ |  76.6$\pm 0.5$ |  75.0$\pm 1.5$ |  63.5$\pm 0.4$ |  66.9$\pm 0.7$ | 70.4
GPT-GNN |  75.0$\pm 2.5$ |  74.9$\pm 2.7$ |  59.3$\pm 0.8$ |  77.0$\pm 1.7$ |  76.1$\pm 0.4$ |  78.5$\pm 0.9$ |  63.1$\pm 0.5$ |  67.5$\pm 1.3$ | 71.4
Grover |  75.8$\pm 1.7$ |  76.9$\pm 1.9$ |  60.7$\pm 0.5$ |  77.8$\pm 1.4$ |  76.3$\pm 0.6$ |  79.5$\pm 1.1$ |  63.4$\pm 0.6$ |  68.0$\pm 1.5$ | 72.3
MGSSL (DFS) |  78.1$\pm 1.8$ |  79.7$\pm 2.2$ |  60.5$\pm 0.7$ | 79.5$\mathbf{\pm 1.1}$ |  76.4$\pm 0.4$ | 79.7$\mathbf{\pm 0.8}$ |  63.8$\pm 0.3$ | 70.5$\mathbf{\pm 1.1}$ | 73.5
MGSSL (BFS) | 78.7$\mathbf{\pm 1.5}$ | 80.7$\mathbf{\pm 2.1}$ | 61.8$\mathbf{\pm 0.8}$ |  78.8$\pm 1.2$ |  76.5$\pm 0.3$ |  79.1$\pm 0.9$ | 64.1$\mathbf{\pm 0.7}$ |  69.7$\pm 0.9$ | 73.7

(a) clintox (b) sider (c) hiv (d) bace

Figure 4: Training and testing curves of different pre-training strategies on
GINs. Solid and dashed lines indicate training and testing curves
respectively.

###  4.2 Results and Analysis

Results on Downstream Tasks.  In Table 1, we show the testing performance on
downstream molecular prediction benchmarks using different self-supervised
pre-training strategies with GIN. We have the following observations: 1)
Generally, GNN models can benefit from various self-supervised pre-training
tasks. The average prediction ROC-AUC of all the pre-trained models are better
than GNN with direct finetuning. 2) MGSSL methods achieve the best performance
on 7 out of 8 benchmarks, demonstrating the effectiveness of motif-based self-
supervised pre-training. 3) We also show MGSSL with two motif generation
orders in Table 1. Both methods show significant performance improvements on
downstream tasks and BFS has a small edge over DFS on average. This may be
explained by the fact that MGSSL is required to generate motifs layer-wise in
BFS orders, which helps GNNs learn more structural information of motifs. In
the following experiments, we use BFS order as the default setting.

In Figure 4, we further show the training and testing curves of MGSSL (BFS and
DFS) and the selected baselines. Due to the page limits, we select 4 benchmark
datasets here. Beyond predictive performance improvement, our pre-trained GNNs
have faster training and testing convergence than baseline methods. Since the
pre-training is a one-time-effort, once pre-trained with MGSSL, the pre-
trained GNNs can be used for various downstream tasks with minimal finetuning
overhead.

Table 2: Compare pre-training gains with different GNN architectures, averaged ROC-AUC ($\%$) on 8 benchmark datasets Model | GCN | GIN | RGCN | DAGNN | GraphSAGE
---|---|---|---|---|---
No pretrain | 68.8 | 67.0 | 68.3 | 67.1 | 68.3
MGSSL (BFS) | 72.7 | 73.7 | 73.0 | 72.3 | 73.4
Relative gain | 5.7$\%$ | 10.0$\%$ | 6.9$\%$ | 7.7 $\%$ | 7.5$\%$

Influence of the Base GNN. In Table 2 we show that MGSSL is agnostic to the
GNN architectures by trying five popular GNN models including GIN [46], GCN
[19], RGCN[33], GraphSAGE[11] and DAGNN [24]. We report the average ROC-AUC
and the relative gains on 8 benchmarks. We observe that all these GNN
architectures can benefit from motif-based pre-training tasks. Moreover, GIN
achieves the largest relative gain and the best performance after pre-
training.

(a) MGSSL (BFS) (b) MGSSL (DFS)

Figure 5: Influence of the size of motif vocabulary

Influence of Molecule Fragmentation.  Here we show proper molecule
fragmentation methods are vital for the motif-based pre-training. Given
different molecule fragmentation methods, different motif vocabulary are
generated with varying sizes. Other than the fragmentation method introduced
in this paper, we also try other fragmentation schemes with different
granularities [17, 4]. BRICS alone [4] tends to generate motifs with large
numbers of atoms. Due to the combinatorial explosion, its generated motif
vocabulary has a size over 100k while more than 90$\%$ motifs have frequencies
less than 5. On the other hand, JT-VAE [17] fragments molecules into rings and
bonds and has a motif vocabulary size of less than 800. Our methods generate
around 12k distinct motifs. By combining different fragmentation strategies,
we are able to fragment molecules with intermediate granularities. In Figure
5, we show the influence of the size of motif vocabulary on 5 benchmark
datasets. We can observe that the pre-trained models achieve the optimal
performance with the motif vocabulary generated by our method. This may be
explained by the following reasons: 1) When the motif segmentation is too
coarse and the motif vocabulary is too large, the generated motif trees have
fewer nodes. It is harder for GNNs to capture the structural information of
motifs. Moreover, the generated motifs have low occurrence frequencies, which
prevents GNNs from learning the general semantic information of motifs that
can be generalized to downstream tasks. 2) When the motif segmentation is too
fine, many generated motifs are single atoms or bonds, which inhibits GNNs
from learning higher level semantic information through motif generation
tasks.

Methods | Avg. ROC-AUC
---|---
w/o atom-level | 73.0
Sequential pre-training | 73.4
Multi-level | 73.7
Table 3: Ablation studies on multi-level self-supervised pre-training.

Ablation Studies on Multi-level Self-supervised Pre-training. We perform
ablation studies to show the effectiveness of the multi-level pre-training. In
Table 3, w/o atom-level denotes pre-training GNNs with Motif-level tasks only
and the sequential pre-training denotes performing the Motif-level tasks after
the Atom-level tasks. As observed from Table 3, the multi-level pre-training
has larger average ROC-AUC than the two variants. We can have the following
interesting insights: 1) The Atom-level pre-training tasks enables GNNs to
first capture the atom-level information, which can benefit higher level,
i.e., motif-level tasks. 2) Since our multi-level pre-training unifies multi-
scale pre-training tasks and adaptively assigns the weights for hierarchical
tasks, it can achieve better performance than sequential pre-training.

##  5 Conclusion and Future Works

In this paper, we proposed Motif-based Graph Self-supervised Learning (MGSSL),
which pre-trains GNNs with a novel motif generation task. Through pre-
training, MGSSL empowers GNNs to capture the rich semantic and structural
information in graph motifs. First, a retrosynthesis-based algorithm with two
additional rules are leveraged to fragment molecule graphs and derive semantic
meaningful motifs. Second, a motif generative pre-training framework is
designed and two specific generation orders are considered (BFS and DFS). At
each step, the pre-trained GNN is required to make topology and motif label
predictions. Furthermore, we designed a multi-level pre-training to unify
hierarchical self-supervised tasks. Finally, we conducted extensive
experiments to show that MGSSL overperforms all the state-of-the-art baselines
on various downstream benchmark tasks. Interesting future work includes 1)
Designing more self-supervised pre-training tasks based on graph motifs. 2)
Exploring motif-based pre-training in other domains other than molecules.

## Acknowledgments and Disclosure of Funding

We thank the anonymous reviewers for valuable feedback. This research was
supported by grants from the National Natural Science Foundation of China
(Grants No. 61922073 and U20A20229) and 2021 Tencent Rhino-Bird Research Elite
Training Program.

## References

  * [1] Jin Chen, Wynne Hsu, Mong Li Lee, and See-Kiong Ng.  Nemofinder: Dissecting genome-wide protein-protein interactions with meso-scale network motifs.  In SIGKDD, pages 106–115, 2006.
  * [2] Ting Chen, Simon Kornblith, Mohammad Norouzi, and Geoffrey Hinton.  A simple framework for contrastive learning of visual representations.  In ICML, pages 1597–1607. PMLR, 2020.
  * [3] Junyoung Chung, Caglar Gulcehre, KyungHyun Cho, and Yoshua Bengio.  Empirical evaluation of gated recurrent neural networks on sequence modeling.  arXiv preprint arXiv:1412.3555, 2014.
  * [4] Jörg Degen, Christof Wegscheid-Gerlach, Andrea Zaliani, and Matthias Rarey.  On the art of compiling and using’drug-like’chemical fragment spaces.  ChemMedChem: Chemistry Enabling Drug Discovery, 3(10):1503–1507, 2008.
  * [5] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  Bert: Pre-training of deep bidirectional transformers for language understanding.  ACL, 2018.
  * [6] Felix A Faber, Luke Hutchison, Bing Huang, Justin Gilmer, Samuel S Schoenholz, George E Dahl, Oriol Vinyals, Steven Kearnes, Patrick F Riley, and O Anatole Von Lilienfeld.  Prediction errors of molecular machine learning models lower than hybrid dft error.  Journal of chemical theory and computation, 13(11):5255–5264, 2017\.
  * [7] Robert M French.  Catastrophic forgetting in connectionist networks.  Trends in cognitive sciences, 3(4):128–135, 1999.
  * [8] Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl.  Neural message passing for quantum chemistry.  In ICML, pages 1263–1272. PMLR, 2017.
  * [9] Aditya Grover and Jure Leskovec.  node2vec: Scalable feature learning for networks.  In SIGKDD, pages 855–864, 2016.
  * [10] Hakim Hafidi, Mounir Ghogho, Philippe Ciblat, and Ananthram Swami.  Graphcl: Contrastive self-supervised learning of graph representations.  arXiv preprint arXiv:2007.08025, 2020.
  * [11] William L Hamilton, Rex Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  NeurIPS, 2017.
  * [12] Zhongkai Hao, Chengqiang Lu, Zhenya Huang, Hao Wang, Zheyuan Hu, Qi Liu, Enhong Chen, and Cheekong Lee.  Asgn: An active semi-supervised graph neural network for molecular property prediction.  KDD ’20, 2020.
  * [13] Kaiming He, Haoqi Fan, Yuxin Wu, Saining Xie, and Ross Girshick.  Momentum contrast for unsupervised visual representation learning.  In CVPR, pages 9729–9738, 2020.
  * [14] Weihua Hu*, Bowen Liu*, Joseph Gomes, Marinka Zitnik, Percy Liang, Vijay Pande, and Jure Leskovec.  Strategies for pre-training graph neural networks.  In ICLR, 2020.
  * [15] Ziniu Hu, Yuxiao Dong, Kuansan Wang, Kai-Wei Chang, and Yizhou Sun.  Gpt-gnn: Generative pre-training of graph neural networks.  In SIGKDD, pages 1857–1867, 2020.
  * [16] Martin Jaggi.  Revisiting frank-wolfe: Projection-free sparse convex optimization.  In ICML, pages 427–435. PMLR, 2013.
  * [17] Wengong Jin, Regina Barzilay, and Tommi Jaakkola.  Junction tree variational autoencoder for molecular graph generation.  In ICML, pages 2323–2332. PMLR, 2018.
  * [18] Nadav Kashtan, Shalev Itzkovitz, Ron Milo, and Uri Alon.  Efficient sampling algorithm for estimating subgraph concentrations and detecting network motifs.  Bioinformatics, 20(11):1746–1758, 2004.
  * [19] Thomas N Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  ICLR, 2017.
  * [20] Johannes Klicpera, Janek Groß, and Stephan Günnemann.  Directional message passing for molecular graphs.  In ICLR, 2020.
  * [21] Walter Kohn and Lu Jeu Sham.  Self-consistent equations including exchange and correlation effects.  Physical review, 140(4A):A1133, 1965.
  * [22] Greg Landrum.  Rdkit documentation.  Release, 1(1-79):4, 2013.
  * [23] Renjie Liao, Yujia Li, Yang Song, Shenlong Wang, Charlie Nash, William L Hamilton, David Duvenaud, Raquel Urtasun, and Richard S Zemel.  Efficient graph generation with graph recurrent attention networks.  NeurIPS, 2019.
  * [24] Meng Liu, Hongyang Gao, and Shuiwang Ji.  Towards deeper graph neural networks.  In SIGKDD, pages 338–348, 2020.
  * [25] Chengqiang Lu, Qi Liu, Chao Wang, Zhenya Huang, Peize Lin, and Lixin He.  Molecular property prediction: A multilevel quantum interactions modeling perspective.  In AAAI, volume 33, pages 1052–1060, 2019.
  * [26] James L McDonagh, Arnaldo F Silva, Mark A Vincent, and Paul LA Popelier.  Machine learning of dynamic electron correlation energies from topological atoms.  Journal of chemical theory and computation, 14(1):216–224, 2017\.
  * [27] Ron Milo, Shai Shen-Orr, Shalev Itzkovitz, Nadav Kashtan, Dmitri Chklovskii, and Uri Alon.  Network motifs: simple building blocks of complex networks.  Science, 298(5594):824–827, 2002.
  * [28] Zhen Peng, Wenbing Huang, Minnan Luo, Qinghua Zheng, Yu Rong, Tingyang Xu, and Junzhou Huang.  Graph representation learning via graphical mutual information maximization.  In Proceedings of The Web Conference 2020, pages 259–270, 2020\.
  * [29] Bryan Perozzi, Rami Al-Rfou, and Steven Skiena.  Deepwalk: Online learning of social representations.  In SIGKDD, pages 701–710, 2014.
  * [30] Jiezhong Qiu, Qibin Chen, Yuxiao Dong, Jing Zhang, Hongxia Yang, Ming Ding, Kuansan Wang, and Jie Tang.  Gcc: Graph contrastive coding for graph neural network pre-training.  In SIGKDD, pages 1150–1160, 2020.
  * [31] Bharath Ramsundar, Peter Eastman, Patrick Walters, and Vijay Pande.  Deep learning for the life sciences: applying deep learning to genomics, microscopy, drug discovery, and more.  " O’Reilly Media, Inc.", 2019.
  * [32] Yu Rong, Yatao Bian, Tingyang Xu, Weiyang Xie, Ying Wei, Wenbing Huang, and Junzhou Huang.  Self-supervised graph transformer on large-scale molecular data.  NeurIPS, 33, 2020.
  * [33] Michael Schlichtkrull, Thomas N Kipf, Peter Bloem, Rianne Van Den Berg, Ivan Titov, and Max Welling.  Modeling relational data with graph convolutional networks.  In European semantic web conference, pages 593–607. Springer, 2018\.
  * [34] Nadine Schneider, Roger A Sayle, and Gregory A Landrum.  Get your atoms in order an open-source implementation of a novel and robust molecular canonicalization algorithm.  Journal of chemical information and modeling, 55(10):2111–2120, 2015.
  * [35] Kristof T Schütt, Farhad Arbabzadah, Stefan Chmiela, Klaus R Müller, and Alexandre Tkatchenko.  Quantum-chemical insights from deep tensor neural networks.  Nature communications, 8(1):1–8, 2017.
  * [36] Kristof T Schütt, Pieter-Jan Kindermans, Huziel E Sauceda, Stefan Chmiela, Alexandre Tkatchenko, and Klaus-Robert Müller.  Schnet: A continuous-filter convolutional neural network for modeling quantum interactions.  NeurIPS, 2017.
  * [37] Ozan Sener and Vladlen Koltun.  Multi-task learning as multi-objective optimization.  NeurIPS, 2018.
  * [38] Teague Sterling and John J Irwin.  Zinc 15–ligand discovery for everyone.  Journal of chemical information and modeling, 55(11):2324–2337, 2015.
  * [39] Fan-Yun Sun, Jordan Hoffman, Vikas Verma, and Jian Tang.  Infograph: Unsupervised and semi-supervised graph-level representation learning via mutual information maximization.  In ICLR, 2019.
  * [40] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.  Graph attention networks.  arXiv preprint arXiv:1710.10903, 2017.
  * [41] Petar Velickovic, William Fedus, William L Hamilton, Pietro Liò, Yoshua Bengio, and R Devon Hjelm.  Deep graph infomax.  In ICLR, 2019.
  * [42] Hao Wang, Defu Lian, Hanghang Tong, Qi Liu, Zhenya Huang, and Enhong Chen.  Hypersorec: Exploiting hyperbolic user and item representations with multiple aspects for social-aware recommendation.  ACM Trans. Inf. Syst., 2021.
  * [43] Hao Wang, Tong Xu, Qi Liu, Defu Lian, Enhong Chen, Dongfang Du, Han Wu, and Wen Su.  Mcne: An end-to-end framework for learning multiple conditional network representations of social network.  KDD ’19, 2019.
  * [44] Junmei Wang and Tingjun Hou.  Application of molecular dynamics simulations in molecular property prediction ii: diffusion coefficient.  Journal of computational chemistry, 2011.
  * [45] Zhenqin Wu, Bharath Ramsundar, Evan N Feinberg, Joseph Gomes, Caleb Geniesse, Aneesh S Pappu, Karl Leswing, and Vijay Pande.  Moleculenet: a benchmark for molecular machine learning.  Chemical science, 9(2):513–530, 2018.
  * [46] Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?  In ICLR, 2019.
  * [47] Chaochao Yan, Qianggang Ding, Peilin Zhao, Shuangjia Zheng, JINYU YANG, Yang Yu, and Junzhou Huang.  Retroxpert: Decompose retrosynthesis prediction like a chemist.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, NeurIPS, volume 33, pages 11248–11258. Curran Associates, Inc., 2020.
  * [48] Kevin Yang, Kyle Swanson, Wengong Jin, Connor Coley, Philipp Eiden, Hua Gao, Angel Guzman-Perez, Timothy Hopper, Brian Kelley, Miriam Mathea, et al.  Analyzing learned molecular representations for property prediction.  Journal of chemical information and modeling, 59(8):3370–3388, 2019\.
  * [49] Zhitao Ying, Jiaxuan You, Christopher Morris, Xiang Ren, Will Hamilton, and Jure Leskovec.  Hierarchical graph representation learning with differentiable pooling.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, NeurIPS, volume 31, 2018.
  * [50] Jiaxuan You, Rex Ying, Xiang Ren, William Hamilton, and Jure Leskovec.  Graphrnn: Generating realistic graphs with deep auto-regressive models.  In ICML, pages 5708–5717. PMLR, 2018.
  * [51] Yuning You, Tianlong Chen, Yongduo Sui, Ting Chen, Zhangyang Wang, and Yang Shen.  Graph contrastive learning with augmentations.  NeurIPS, 33, 2020.
  * [52] Yuning You, Tianlong Chen, Zhangyang Wang, and Yang Shen.  When does self-supervision help graph convolutional networks?  In ICML. PMLR, 2020.
  * [53] Muhan Zhang, Zhicheng Cui, Marion Neumann, and Yixin Chen.  An end-to-end deep learning architecture for graph classification.  In AAAI, volume 32, 2018.
  * [54] Shichang Zhang, Ziniu Hu, Arjun Subramonian, and Yizhou Sun.  Motif-driven contrastive learning of graph representations.  arXiv preprint arXiv:2012.12533, 2020.
  * [55] Zaixi Zhang, Qi Liu, Zhenya Huang, Hao Wang, Chengqiang Lu, Chuanren Liu, and Enhong Chen.  Graphmi: Extracting private graph data from graph neural networks.  In Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI-21.