# Explainable Biomedical Recommendations via Reinforcement Learning Reasoning
on Knowledge Graphs

Gavin Edwards1<EMAIL_ADDRESS>Sebastian Nilsson2
<EMAIL_ADDRESS>Benedek Rozemberczki1
<EMAIL_ADDRESS>Eliseo Papa1<EMAIL_ADDRESS>1Biological Insight Knowledge Graph (BIKG)  Research D&A  R&D IT  AstraZeneca
Cambridge  UK  2Biological Insight Knowledge Graph (BIKG)  Research D&A  R&D
IT  AstraZeneca  Gothenburg  Sweden

###### Abstract

For Artificial Intelligence to have a greater impact in biology and medicine,
it is crucial that recommendations are both accurate and transparent. In other
domains, a neurosymbolic approach of multi-hop reasoning on knowledge graphs
has been shown to produce transparent explanations. However, there is a lack
of research applying it to complex biomedical datasets and problems. In this
paper, the approach is explored for drug discovery to draw solid conclusions
on its applicability. For the first time, we systematically apply it to
multiple biomedical datasets and recommendation tasks with fair benchmark
comparisons. The approach is found to outperform the best baselines by 21.7%
on average whilst producing novel, biologically relevant explanations.

##  1 Introduction

Drug development costs are ten times what the industry spent per year in the
1980s and continously growing. The average new drug is now expected to take
ten or more years to develop and cost over $2 billion (Austin and Hayford,
2021). To improve on this trend, there is increasing interest in enhancing
human decision making with Artificial Intelligence (AI) recommendation
systems.

AI recommendation systems commonly use black-box models that solely focus on
optimising performance for tasks such as product recommendation and
personalised content. These are high-frequency, low-risk recommendations that
receive fast feedback to improve and build trust. However, the application of
recommendations in the drug discovery space suffers from significant
limitations (Gogleva et al., 2021): decisions are infrequent, high-risk and
can have a large impact.

Figure 1: Transparent reasoning for why Everolimus could treat Unipolar
Depression.

AI may pave the way for innovative new therapies, but pursuing the
recommendations means devoting time and resources. First, experts are needed
to label and rank recommendation lists, and then laboratory scientists must
undertake costly and often time-consuming experiments to provide further
validation. Finally, stakeholders need to persevere when provided
recommendations are not successful, which is certain to happen, given the
notoriously high failure rates of new drugs. The whole process demands robust,
sustained trust from humans. Therefore, a human-centric approach is essential.
We argue that, in tandem with accurate recommendations, providing transparent
reasoning, as shown in Figure 1, is a key part to building the trust required.

Drug discovery recommendation tasks include, but are not limited to, drug
repositioning, disease target identification and personalised medicine. Some
recent successes have been obtained in the application of recommendation
systems on top of networks and knowledge graphs of biological knowledge
(Barabási et al., 2011). For instance, machine learning models trained on
graphs have been used to find novel drug candidates for COVID-19 (Zhou et al.,
2020; Gysi et al., 2021), identify efficacious drug combination therapies
(Cheng et al., 2019) and discover novel drivers of drug resistance (Gogleva et
al., 2021).

A Knowledge Graph (KG) is a data representation that links diverse types of
data together into a single, unified model. It enables complex and nuanced
relationships to be captured. Due to the increasing breadth and depth of
biomedical data, biomedical KGs have become a popular way to integrate the
complex and diverse data into a single unified representation.

Biomedical KGs can contain and integrate data ranging from disease networks
and protein interactions to healthcare records and scientific knowledge
extracted from literature. The intrinsic patterns, uncertainty and bias within
biomedical data presents many unique challenges for reasoning (Bonner et al.,
2021). Biomedical KGs are also structurally different from standard KGs
(Figure 4). Because of this, the common assumptions and methods for graph
reasoning often fail to generalise (Liu et al., 2021).

###  1.1 Main contributions:

This paper makes the following novel contributions:

  1. 1.

Applying multi-hop neural-driven recommendation to complex biomedical KGs 10x
bigger than previously used.

  2. 2.

Benchmarking multiple algorithms on representative biomedical datasets and
tasks with standardised evaluations.

  3. 3.

Validating if biologically relevant explanations are produced with biomedical
expert end users.

  4. 4.

Finding that multi-hop reasoning has the potential to generate explanations
and boost the performance of black-box methods.

##  2 Related Work

###  2.1 Knowledge Graphs

Figure 2: Basic graph (a) vs knowledge graph (b)

At the fundamental level, graphs contain information about a set of entities
(also referred to as nodes or mathematically as vertex) and the set of
relationships (also known as edges) between them. A Knowledge Graph (KG) is a
heterogeneous, multi-relation graph that can contain directed and undirected
edges. This means it has many types of entities and relations, can have
multiple edges of the same type between entities, and the edges may indicate
direction (Figure 2).

The core component of a graph is a triple. A single triple defines how two
entities are linked. It consists of a head entity $h$, relation $r$ and target
entity $t$. A triple is therefore defined by $(h,r,t)$. An example of a triple
would be $(London,\text{ }capital\text{ }of,\text{ }United\text{ }Kingdom)$.
When multiple triples are combined, it creates a graph.

A common graph reasoning task is link prediction. It is used to discover new
links within incomplete KGs. Many classical AI problems such as
recommendations, reasoning and question answering can be rephrased as link
prediction problems. The task is to find new links by identifying entities
that correctly complete a query. For example, predicting drugs to treat
COVID-19 could be defined as finding tail entities that correctly complete the
triple $(COVID\text{-}19,treated\text{ }by,?)$.

In industry, KGs have been successfully deployed in many applications such as
search engines, social graphs, fraud detection and product recommendation. KGs
commonly used for benchmarking within research include FB15K-237 (Toutanova et
al., 2015), NELL-995 (Xiong et al., 2017) and WN18RR (Dettmers et al., 2018).

However, the lack of large and reliable benchmark graphs for research purposes
is well known (Hu et al., 2020). In particular, standard, representative and
fair benchmark graphs are still limited for specialist problems such as
biomedical recommendation. Initiatives such as the Open Graph Benchmark (Hu et
al., 2020) are good attempts at solving this, but still not representative of
large biomedical KGs used in industry. OGB’s biomedical KG is relatively small
and has only 5 node types (appendix:dataset_details). It also differs
structurally when compared to the low degree graphs with NLP edges extracted
from literature (Figure 4).

###  2.2 Biomedical Knowledge Graphs

Biomedical KGs are a specialist type of KG that focus on the biomedical
domain. They include information between entities such as genes, diseases,
pathways and compounds. The data is usually collected from a wide variety of
data sources (Bonner et al., 2021). An example ontology can be seen in Figure
3. The biomedical data is large, noisy, incomplete and contains contradictory
observations. Structurally, the graphs exhibit long-range dependencies, many
high-degree hub entities, a higher density of links, long tails of weakly
connected entities and higher heterogeneity than standard benchmark graphs
(Figure 4).

Biomedical KGs are becoming increasingly popular within industry and as open-
source initiatives. Examples include, BioKG (Walsh et al., 2020), DRKG
(Ioannidis et al., 2020), Clinical KG (Santos et al., 2020), Hetionet
(Himmelstein et al., 2017), OpenBioLink (Breit et al., 2020), OGBL-BIOKG (Hu
et al., 2020) and PharmKG (Zheng et al., 2020b).

Figure 3: Ontology of Hetionet biomedical knowledge graph (Himmelstein et al.,
2017). Source: https://het.io/about/

[] $10^{4}$$10^{5}$$10^{6}$$10^{7}$$10^{5}$$10^{6}$$10^{7}$$10^{8}$Number of
entitiesNumber of triplesGraph Size$10^{1}$$10^{2}$$10^{1}$$10^{2}$Number of
entity typesNumber of edge typesGraph Heterogeneity4 []
$0$$50$$100$$150$WN18RRNELL-995FB15K-237ClinicalKGOpenBioLink-HQBIKG
Hetionet+BioKGOGBL-BIOKGDRKGHetionetBIKG HetionetPharmKG-8kDegree

Figure 4: Comparison of biological and generic KG datasets. (a) Graph size and
heterogeneity. Triangles ($\Delta$) represent benchmark KGs, circles ($\circ$)
note biomedical KGs. (b) Skewness of degree distribution. Box plots represent
the quartiles of the degree distributions and the mean degree is denoted by
the diamonds ($\diamondsuit$). Biomedical KGs have a highly skewed degree
distribution caused by a long-tail of poorly connected entities and a large
number of highly connected hub-like entities.

###  2.3 Existing Reasoning Approaches

Existing approaches to reasoning on KGs originate from different theoretical
perspectives with varying levels of explainability. The approaches can be
categorised into symbolic, neural and neurosymbolic reasoning approaches
(Zhang et al., 2021).

#### Symbolic

Symbolism assumes symbols are the fundamental unit of human intelligence. It
treats cognition purely as a series of inferences upon symbolic
representations (Haugeland, 1989). Symbolic reasoning approaches deduce
general logic rules to use when predicting new links.

Importantly, as they follow clear symbolic rules, the reasoning is fully
transparent to the user. AnyBURL (Meilicke et al., 2019) and pLogicNet (Qu and
Tang, 2019) both involve the learning of rules. Recent updates to AnyBURL
added a more neuro-symbolic style. They were shown to yield impressive but not
superior performance on the biomedical KG Hetionet (Liu et al., 2021).

Due to the strict matching and discrete logic operations, symbolic methods are
thought to be intolerant to the ambiguity and noise of biomedical data.
However, the inflexible rules can restrict the expressiveness. Additionally,
it was hypothesised that large amounts of high-degree entities makes the
learning and applying of logical rules more difficult (Liu et al., 2021).

#### Neural

In contrast, connectionism is subsymbolic and imitates the neurons within
brains to build models (Rosenblatt, 1958; Rumelhart et al., 1986). Instead of
symbolic representations, neural reasoning approaches map entities and edges
in KGs to low-dimensional vector representations known as embeddings. This
approach is the focus of a large body of recent research (Wu et al., 2020).

Despite achieving state-of-the-art results on benchmark tasks, a fundamental
problem with the neural approach is its non-transparent nature. Predictions
are single-hop within vector space and the contributing features driving the
prediction are hidden to humans.

Traditional embeddings methods are trained to minimise the reconstruction
error in the immediate first-order neighbourhood while discarding higher-order
proximities. There is also a homophily assumption in proximity preserving
embeddings. However, most explanatory metapaths for drug repurposing were
found to have length two or more (Liu et al., 2021) and biomedical graphs
display heterophily.

More recently, Graph Neural Networks (GNNs) learn node embeddings by
aggregating incoming messages from neighbouring nodes. In theory, these
methods are capable of modeling longer-term dependencies. However, in the
presence of many high-degree nodes, like in biomedical KGs, the size of
neighbourhoods leads to many uninformative signals and causes over-smoothing
(Kipf and Welling, 2016).

GNNs are also limited by poor expressiveness of long-tail entities with small
neighbourhoods (Guo et al., 2019). This is common in biomedical KGs,
demonstrated by the average degree’s 25th percentile of just 3 neighbours and
positive skews in Figure 4.

#### Neurosymbolic

The neural-symbolic approach attempts to synergise the neural and symbolic
approaches. It is being increasingly referred to by the colloquial term
”neurosymbolic”. In graph reasoning, this leverages a symbolic graph
representation to reduce the search space, cleanse uninformative signals and
provide interpretability, whilst using neural approaches to increase
expressiveness and adapt the reasoning to ambiguity.

These attributes are highly relevant for reasoning on biomedical KGs. Multi-
hop reasoning with neural networks, a path based approach to KG reasoning,
falls into the subcategory of neural-driven symbolic reasoning.

###  2.4 Reinforcement Learning Reasoning

Formulating KG reasoning as a neural-driven multi-hop problem naturally lends
itself to Reinforcement Learning (RL). RL is known to be effective in using
neural approaches for solving sequential decision-making problems. Alternative
neurosymbolic approaches such as NTP (Rocktäschel and Riedel, 2017) and Neural
LP (Yang et al., 2017) fail to scale to large knowledge graphs (Liu et al.,
2021).

RL methods have achieved superhuman levels of performance in sequential
decision making, famously for arcade games (Mnih et al., 2013), board games
(Silver et al., 2017), online multiplayer games (Berner et al., 2019) and
imperfect-information games (Brown et al., 2020). These successes demonstrated
RL’s ability to operate on ambiguous data, understand complex environments and
infer high-level causal relationships. We believe these to be very relevant
attributes for the application to biomedical KGs and drug discovery
recommendations in general.

The approach of using RL for recommendation systems is increasing in
popularity in industry and academia (Afsar et al., 2021; Chen et al., 2019).
RL has already been applied to the tasks of question answering (Das et al.,
2018), fact-checking (Xiong et al., 2017) and recommendations (Liu et al.,
2021).

Most of the approaches use path based RL reasoning, which sequentially
traverses triples in a graph to make predictions. Because of this, it provides
transparent reasoning for every prediction. An issue with this approach can be
convergence and generating relevant reasoning paths. To help, demonstrations
of reasoning can be provided by metapaths which acts as a form of weak
supervision (Zhao et al., 2020; Liu et al., 2021). RL reasoning applied to the
standard benchmarks graphs used in research has not shown promising
performance and thus has not been the focus of much research. A small sample
of path-based models and their characteristics are shown in Table 1.

Table 1: Sample of RL reasoning algorithms Algorithm | Reference | Approach
---|---|---
DeepPath | Xiong et al. (2017) | RL path finding for fact checking
MINERVA | Das et al. (2018) | Policy-based RL for QA
M-Walk | Shen et al. (2018) | Monte Carlo Tree Search when sampling
MultiHop | Lin et al. (2018) | Soft reward and action dropout for MINERVA
PGPR | Xian et al. (2019) | Policy-based RL walks
RuleGuider | Lei et al. (2020) | Discovers rules to use for RL walks
R2D2 | Hildebrandt et al. (2020) | MINERVA with debate dynamics
PoLo | Liu et al. (2021) | MINERVA with metapaths

MINERVA is a path based RL algorithm for the similar task of question
answering. Two extensions of MINERVA, R2D2 and PoLo, applied the method to a
biomedical KG (Hildebrandt et al., 2020; Liu et al., 2021). PoLo highlighted
the potential for multi-hop reasoning to be especially suited to biomedical
KGs.

Whilst promising in their approach, these works used Hetionet which is a
relatively small KG that lacks the full complexity, size and noise found in
industrial biomedical KGs. Additionally, Hetionet has not been updated since
2017. PoLo’s results came from testing a limited number (151) of Compound-
treats-Disease links and the evaluation strategy was not standardised as
pruning was not tested across all models. The ambiguity of how metrics are
calculated also creates the possibility for unfair comparisons (Sun et al.,
2020; Berrendorf et al., 2021). To overcome the limitations of these studies,
we extend their work to larger biomedical KGs and offer a set of standardised
evaluations.

##  3 Preliminaries

The algorithms MINERVA and PoLo are used for path-based RL reasoning to
directly build on the claims in Liu et al. (2021) that the approach is suited
to biomedical KGs. Both algorithms use the same definitions of the knowledge
graph environment and task, which are defined below.

###  3.1 Knowledge Graph

A triple links two entities together. It is defined by $(h,r,t)$ where $h$ is
the head entity, $r$ is the relation, and $t$ is the tail entity. An example
of a triple would be $(Ibuprofen,treats,Headache)$. E represents the full set
of entities and R the set full of relations in a graph. Every entity has a
unique type T. The KG is defined as a set of $(h,r,t)$ triples
$\set{KG}\subset\set{E}\times\set{R}\times\set{E}$. To enable RL agents to
travel bidirectionally across triples, the inverse triples are added to the
$KG$. The inverse relations are denoted as $r^{-1}$. For example for
$(e_{1},r,e_{2})$ the triple $(e_{2},r^{-1},e_{1})$ is also added. All triples
in the KG are treated as true facts.

###  3.2 Metapaths

A metapath is a logical reasoning rule defining a general pathway between two
entities. For example, a desired metapath for predicting
$(Compound,treats,Disease)$ triples might be:

$(Compound\overset{binds}{\rightarrow}Gene\overset{associates}{\rightarrow}Disease)$

###  3.3 Recommendation

The task of recommendation is posed as a link prediction problem via
sequential graph traversal. Given a query $(e_{1},r,?)$, such as
$(COVID\text{-}19,treated\text{ }by,?)$, the task is to sequentially traverse
triples in the KG to find the correct answer entities. At each step $s$ in the
graph, the set of possible actions is the outgoing triples of the current
entity $e_{s}$. Each action is represented by the triples relation and tail
entity $(r,t)$. Vector space embeddings are used to represent the relations
and entities for input to the models. A ’NO_OP’ action is augmented to each
node, which effectively allows for no action to be taken and staying at the
current node. A maximum number of steps is set to prevent the model from
traversing indefinitely.

##  4 Experiments

We ran multiple experiments to establish if the path based RL reasoning
approach is suited to biomedical recommendation tasks. We applied it to three
different KGs on two separate tasks to draw concrete conclusions on its
suitability. Ablations are performed to understand the approach in more depth.

###  4.1 Tasks

The models are evaluated on two common biomedical recommendation tasks: _drug
repurposing_ and _drug-target interaction_. In the first task, the goal is to
find $(Compound,treats,Disease)$ triples. Formally the query is defined as
$(Compound,treats,?)$. For the second task, the goal is predicting
$(Compound,binds,Gene)$ triples, with the formal query being defined as
$(Compound,binds,?)$.

###  4.2 Datasets

Three datasets are used. Hetionet as a small open-source baseline graph and
two projections of Biological Insights Knowledge Graph (BIKG) (Geleta et al.,
2021). BIKG is an internal biomedical KG actively used by organization name
anonymised for peer review for recommendations. The two projections are ‘BIKG
Hetionet’ and ‘BIKG Hetionet+’. The specific projections of BIKG are chosen as
they represent the two separate structural types of biomedical KGs, which can
be seen in Figure 4. Importantly, these specific datasets were chosen as they
share the same core ontology as Hetionet (Figure 3), which allows for a unique
systematic analysis.

#### Hetionet

Hetionet (Himmelstein et al., 2017) is an open-source dataset developed in
2017 as part of a study looking at drug repurposing using a KG. It is the most
commonly used biomedical graph for benchmarking. Hetionet combines data from
29 public databases and contains over 47k entities of 11 types, including
genes, compounds and diseases. These are linked by over 2.2M edges of 24
different relation types. The full ontology is shown in Figure 3.

#### BIKG Hetionet

BIKG (Geleta et al., 2021) unifies data from different public, private and
internal sources into a single, comprehensive knowledge graph. BIKG Hetionet
is a projection of BIKG using the same Hetionet schema in Figure 3, except
with almost five times more edges and three times more nodes. It represents
the type of biomedical graph with a higher density and average degree (Figure
4). Usually this is because they do not include edges mined from literature.
Exact dataset statistics can be found in appendix:experiment_details.

#### BIKG Hetionet+

Like other biomedical KGs, a large part of the full BIKG is from Natural
Language Processing (NLP) pipelines which extract triples from literature.
These are inherently more noisy and ambiguous edges. BIKG Hetionet+ extends on
BIKG Hetionet by adding ‘has link’ NLP edges and 22 other new edge types. This
nearly doubles the number of edge types. One extra node type ‘Mechanism of
Action’ is also added. These additions result in a much noisier and more
complex KG, representing the sparser biomedical graphs that contain literature
edges and a lower average degree caused by a longer tail of poorly connected
entities (Figure 4). Exact statistics can be found in
appendix:experiment_details.

###  4.3 Algorithms

For the RL algorithm, MINERVA (Das et al., 2018) is selected. PoLo (Liu et
al., 2021), which extends MINERVA with extra rewards for following pre-defined
metapaths, is also used on the drug repurposing task. Using PoLo tests the
effect of weak supervision from injecting prior domain knowledge. The
metapaths are re-used from the original paper.

The following baseline algorithms were chosen to cover the different
theoretical approaches and the most competitive models from the initial Liu et
al. (2021) results on Hetionet – TransE (Wang et al., 2014), DistMult (Yang et
al., 2014), R-GCN (Schlichtkrull et al., 2018) and AnyBURL (JUNO version)
(Meilicke et al., 2019).

###  4.4 Evaluation Strategy

Ambiguity exists in the evaluation of link prediction methods. Metrics are not
uniformly calculated across individual algorithm implementations and common
software libraries. This leads to unfair comparisons in reported results (Sun
et al., 2020; Berrendorf et al., 2021). Therefore, we implement a independent
and standardised approach to fairly evaluate predictions. To the best of our
knowledge, this is the first time it has been done for biomedical KG
recommendations across symbolic, neural and neuro-symbolic methods.

All metrics are filtered (Bordes et al., 2013), deduplicated and evaluated for
tail-sided predictions against all nodes in the graph. Pruning is also applied
which filters predictions to only nodes of the correct target type, as this
would be done before giving recommendation lists to scientists. Standard
metrics for recommendation tasks are used, HITS@1, HITS@3, HITS@10, and mean
reciprocal rank (MRR).

An initial hyperparameter search is applied for every dataset, task and model
to find sensible hyperparameters (appendix:experiment_details). The best
hyperparameters are used for training 5-fold cross-validation with seeded
randomly split triples. Finally, the mean of the five runs is reported
alongside standard errors.

##  5 Results

The results are shown in Table 2. The interpretable RL methods, MINERVA and
PoLo, performed the best. Averaged across all datasets and tasks, MINERVA
outperformed the best baselines by 41.8% for HITS@1, 24.3% for HITS@3, 7.3%
for HITS@10 and 13.5% for MRR. Generating a average performance gain of 21.7%
across metrics.

When metapaths were used in PoLo, it improved MINERVA’s average performance by
a further 4.4%. We found applying pruning to RL results was especially
important; leading to significant performance boosts, as it effectively
eliminates incomplete and failed walks from predictions (appendix:pre-
pruned_results for pre-pruned).

The RL approaches superior performance notably decreased as the graphs
complexity increased. Falling from 41.7% performance gain in BIKG Hetionet to
just 6.8% gain in BIKG Hetionet+. The improvement from metapaths also degraded
from 5.4% on BIKG Hetionet to 2.7% across metrics for BIKG Hetionet+ as the
metapaths became less relevant for the expanded ontology.

The baseline neural and symbolic methods varied in performance, whereas
neurosymbolic methods (MINERVA and PoLo) were able to perform consistently
across datasets and tasks. The embedding models (TransE & DistMult) struggled
with larger biomedical KGs. However, TransE performed notably well for drug-
target interaction. R-GCN performed relatively consistently but achieved a
middle-of-the-pack performance. Finally, AnyBURL was competitive on larger
graphs with the higher ambiguity. The symbolic rules managed to cut through
the noise and perform well. However, in other graphs AnyBURL’s rigid rules
limited its expressiveness.

In summary, we found that neural methods were better suited to smaller,
cleaner biomedical KGs and symbolic methods were suited to larger, noisy
biomedical KGs. However the neurosymbolic methods were able to perform
consistently well across both.

###  5.1 Ablations

#### Impact of Multi-Hop Reasoning

An ablation was ran to isolate the impact of single-hop vector space reasoning
vs multi-hop symbolic reasoning. It was found that using multi-hop symbolic
reasoning had better performance (Table 3) This ablation boosted the
performance of black-box TransE recommendations whilst also providing the much
needed explanations. It highlights the potential for multi-hop reasoning to
work synergistically with graph representations usually used for single-hop
vector space predictions. Multi-hop reasoning algorithms such as MINERVA could
be a way to add explainability on top of black-box approaches like GNNs and
traditional embeddings without sacrificing performance. However, further
research is required into the extent that this is true.

#### Impact of Neural-Driven Reasoning

A second ablation was run to see how the expressiveness of neural driven
reasoning impacts performance. With superior symbolic rules as guidance, the
extra expressiveness of neural-driven reasoning could not outperform the
purely symbolic rules (Table 3). This shows how symbolic rules are very
effective when there are high levels of noise and complexity in the KG. Even
when being guided by the symbolic rules, the dynamic capability of the neural
driven approach did not provide any benefits.

Table 2: Experiment results. green indicates the best result, blue indicates second best result. The interpretable, path based RL methods PoLo and MINERVA performed best. Task | Dataset | Model | HITS@1 | HITS@3 | HITS@10 | MRR
---|---|---|---|---|---|---
Drug Repurposing | Hetionet | TransE | $.193\pm.027$ | $.364\pm.031$ | $.608\pm.034$ | $.312\pm.022$
DistMult | $.028\pm.002$ | $.138\pm.022$ | $.373\pm.039$ | $.129\pm.002$
RGCN | $.047\pm.028$ | $.194\pm.018$ | $.492\pm.043$ | $.221\pm.040$
AnyBURL | $.116\pm.013$ | $.288\pm.024$ | $.572\pm.036$ | $.258\pm.012$
MINERVA | $\mathit{.351\pm.040}$ | $\mathit{.559\pm.058}$ | $\mathit{.790\pm.057}$ | $\mathit{.463\pm.041}$
PoLo | $\mathbf{.378\pm.023}$ | $\mathbf{.575\pm.025}$ | $\mathbf{.845\pm.009}$ | $\mathbf{.479\pm.015}$
BIKG Hetionet | TransE | $.016\pm.003$ | $.034\pm.005$ | $.088\pm.001$ | $.045\pm.004$
DistMult | $.004\pm.002$ | $.013\pm.002$ | $.036\pm.007$ | $.024\pm.003$
RGCN | $.056\pm.005$ | $.108\pm.008$ | $.239\pm.008$ | $.117\pm.005$
AnyBURL | $.000\pm.000$ | $.015\pm.001$ | $.046\pm.004$ | $.023\pm.001$
MINERVA | $\mathit{.121\pm.002}$ | $\mathit{.194\pm.010}$ | $\mathit{.249\pm.013}$ | $\mathit{.160\pm.004}$
PoLo | $\mathbf{.124\pm.005}$ | $\mathbf{.204\pm.010}$ | $\mathbf{.273\pm.016}$ | $\mathbf{.167\pm.008}$
BIKG Hetionet+ | TransE | $.017\pm.002$ | $.052\pm.004$ | $.119\pm.011$ | $.057\pm.002$
DistMult | $.004\pm.002$ | $.014\pm.004$ | $.039\pm.001$ | $.026\pm.004$
RGCN | $.039\pm.008$ | $.104\pm.016$ | $\mathbf{.222\pm.015}$ | $\mathbf{.115\pm.005}$
AnyBURL | $.000\pm.000$ | $.056\pm.006$ | $\mathit{.205\pm.014}$ | $.072\pm.002$
MINERVA | $\mathbf{.076\pm.008}$ | $\mathit{.124\pm.011}$ | $.166\pm.013$ | $\mathit{.113\pm.010}$
PoLo | $\mathit{.074\pm.006}$ | $\mathbf{.133\pm.010}$ | $.173\pm.012$ | $\mathbf{.115\pm.007}$
Drug-Target Interaction | Hetionet | TransE | $\mathbf{.287\pm.008}$ | $\mathbf{.546\pm.013}$ | $\mathbf{.918\pm.007}$ | $\mathbf{.331\pm.005}$
DistMult | $.039\pm.004$ | $.128\pm.007$ | $.367\pm.012$ | $.115\pm.005$
RGCN | $.060\pm.025$ | $.177\pm.057$ | $.376\pm.102$ | $.126\pm.031$
AnyBURL | $.169\pm.043$ | $.323\pm.081$ | $.493\pm.124$ | $.210\pm.051$
MINERVA | $.186\pm.015$ | $\mathit{.420\pm.028}$ | $\mathit{.807\pm.039}$ | $\mathit{.296\pm.019}$
BIKG Hetionet | TransE | $.193\pm.001$ | $.376\pm.012$ | $\mathit{.674\pm.012}$ | $.239\pm.004$
DistMult | $.044\pm.007$ | $.111\pm.014$ | $.269\pm.014$ | $.098\pm.007$
RGCN | $.173\pm.016$ | $.338\pm.013$ | $.539\pm.021$ | $.235\pm.008$
AnyBURL | $\mathit{.215\pm.005}$ | $\mathit{.408\pm.004}$ | $.623\pm.009$ | $\mathit{.264\pm.004}$
MINERVA | $\mathbf{.235\pm.005}$ | $\mathbf{.516\pm.011}$ | $\mathbf{.983\pm.015}$ | $\mathbf{.305\pm.004}$
BIKG Hetionet+ | TransE | $.054\pm.004$ | $.098\pm.007$ | $.171\pm.011$ | $.008\pm.004$
DistMult | $.035\pm.006$ | $.084\pm.011$ | $.206\pm.029$ | $.084\pm.001$
RGCN | $\mathit{.188\pm.008}$ | $.323\pm.019$ | $.539\pm.010$ | $.241\pm.006$
AnyBURL | $\mathbf{.215\pm.003}$ | $\mathbf{.408\pm.008}$ | $\mathit{.622\pm.005}$ | $\mathbf{.263\pm.003}$
MINERVA | $.181\pm.006$ | $\mathit{.368\pm.008}$ | $\mathbf{.628\pm.019}$ | $\mathit{.243\pm.004}$
Table 3: Multi-Hop Ablation. Retraining MINERVA using the same embeddings from TransE where they initially performed best (Drug-Target Interaction on Hetionet) improved performance. Neural-Driven Ablation. Where AnyBURLs symbolic rules performed best (Drug-Target Interaction on BIKG Hetionet+) the same rules were reused as metapaths in PoLo but it failed to improve performance. Ablation | Model | HITS@1 | HITS@3 | HITS@10 | MRR
---|---|---|---|---|---
Multi-Hop | TransE | $\mathit{.287\pm.008}$ | $\mathit{.546\pm.013}$ | $\mathit{.918\pm.007}$ | $\mathit{.331\pm.005}$
MINERVA | $.186\pm.015$ | $.420\pm.028$ | $.807\pm.039$ | $.296\pm.019$
MINERVA* | $\mathbf{.327\pm.014}$ | $\mathbf{.622\pm.013}$ | $\mathbf{.938\pm.018}$ | $\mathbf{.350\pm.006}$
Neural-Driven | AnyBURL | $\mathbf{.215\pm.003}$ | $\mathbf{.408\pm.008}$ | $\mathit{.622\pm.005}$ | $\mathbf{.263\pm.003}$
MINERVA | $.181\pm.006$ | $\mathit{.368\pm.008}$ | $\mathbf{.628\pm.019}$ | $.243\pm.004$
PoLo* | $\mathit{.191\pm.009}$ | $.366\pm.012$ | $.596\pm.019$ | $\mathit{.250\pm.010}$

##  6 Explanations

Although the results are promising, Lv et al. (2021) claim the reasoning
produced by multi-hop models can have gaps in interpretability. If the
explanations are unclear to the end user, it will be detrimental to trust. To
address this aspect, we run inference on the full BIKG Hetionet and BIKG
Hetionet+ graphs, using the trained MINERVA models. We then assess the
resulting high-level metapaths and individual reasoning paths for biological
relevance.

###  6.1 Reasoning Metapaths

We found that the high-level metapaths displayed sensible biological logic.
The reasoning ranges from clear interactions between genes, mechanisms of
action and biological pathways to more vague similarity-based reasons such as
compounds that resemble each other (Figure 5 & 6). However, inspection of the
metapaths also revealed that models can get fixated on less successful
policies, almost exclusively focusing on mechanisms of action entities in BIKG
Hetionet+ and using more cyclic paths, which ultimately led to lower
performance. The top 10 metapaths and their frequency can be found in
appendix:explanations.

Figure 5: MINERVA’s most common metapaths for drug repurposing Figure 6:
MINERVA’s most common metapaths for drug-target interaction

###  6.2 Case Studies

Although the high-level metapaths appear to be relevant, the low-level
reasoning could still be biologically flawed. Therefore, together with subject
matter experts we investigated a sample of the individual reasoning paths to
check the validity.

Novel predictions were made for Lixisenatide, Dulaflutide and similar
compounds to treat non-alcoholic fatty liver disease (NAFLD). It was able to
correctly reason that they act as a receptor agonist (RA) for the gene
Glucagon-like peptide-1 receptor (GLP1R). There is literature supporting the
role of GLP-1 RA’s as a treatment for NAFLD (Seghieri et al., 2018; Lv et al.,
2020), as well as clinical trials and case studies for NAFLD treatment (Gluud
et al., 2014; Seko et al., 2017). The next highest prediction for NAFLD was
for Secukinumab, which follows a similar reasoning but instead via
Interleukin-17A inhibition. There is a clinical trial currently recruiting to
study this on ClinicalTrials.gov (NCT04237116).

Everolimus, usually used to treat cancer, was unconventionally predicted to
treat unipolar depression (Figure 1). However there is case study evidence to
support this theory (Mir et al., 2018). The reasoning was due to links with
Cytochrome P450, which has known links to depression via the mTor pathway. It
also links to Gene PLXNA3, which is important for neurological development and
has studies linking it to depression (Wray et al., 2007). This is a more
speculative, but still biologically relevant, reasoning that would require
additional validation work from drug discovery scientists.

Finally, a prediction was made of Simvastatin binding to gene UGT1A1. This
example highlights how the predictions and explanations are susceptible to
errors in the KG data. The reasoning used the edges of Simvastatin binding to
enzymatic genes CYP2D6, SLCO1B1 and CYP3A5, which are involved in drug
metabolic pathways.UGT1A1 was also shown to be a participant in these
pathways, and therefore it reasoned that Simvastatin was likely to bind to
UGT1A1. All these enzymatic genes are indeed involved in Simvastatin
metabolism (Iwuchukwu et al., 2014). However, they do not actually ’bind’ to
these genes, the edges are misclassified and should be ’affects’. The benefit
of the transparent explanations is that a scientist can easily identify this
fault and dismiss the recommendation. It also means data issues can easily be
identified and fixed, which creates higher quality KG datasets. Whereas for
black-box recommendations the fault in the reasoning would be invisible to end
users and likely go unnoticed.

These case studies are typical of the types of biological reasoning that can
be found throughout the predictions, with some novel and interesting
explanations. There is also a proportion of predictions with more ambiguity,
linking entities together via similar genes, compounds and pharmacological
classes. We provide more examples in appendix:explanations.

##  7 Implications

The superior performance we observed and the accompanying transparent
explanations make the RL approach well suited to biomedical recommendations.
The general approach can be extended to produce even greater explanations.
RL’s reward function enables users to have more control over the reasoning
used and it’s possible to provide reasons a recommendation might not be true
(Hildebrandt et al., 2020).

Additionally, the multi-hop approach allows for an interactive, human-in-the-
loop reasoning process. Humans and AI can collaborate by combining knowledge
to collaboratively walk across the knowledge graph to make predictions. All
these attributes mean the approach can create better informed, more trusted
biomedical recommendations without compromising on accuracy.

##  8 Conclusions

This paper highlights that biomedical KG recommendation is a unique challenge
that cannot be approached in the same way as other recommendation tasks. The
reality of the problem demands recommendations to be human-centric and it is
shown that biomedical KGs unique structure differentiates them as a distinct
problem. Different theoretical types of algorithms were evaluated across
multiple representative datasets and tasks. It was found that multi-hop
neural-driven symbolic reasoning, using RL, outperformed other approaches.
Additionally, using this approach we were able to produce novel and relevant
biological explanations. An ablation also highlighted the potential for this
approach to create explanations for black-box models whilst boosting
performance. Therefore, we believe the approach is well suited to biomedical
recommendations and encourage more research in this direction.

## Acknowledgments

The authors would like to thank Fatima Lugtu, Greet De Baets and Piotr
Grabowski for their assistance in evaluating the recommendations and
explanations. Thanks also goes to the SCP (Scientific Compute Platform) team
at AstraZeneca.

## References

  * Abadi et al. (2016) Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et al.  TensorFlow: A System for Large-Scale Machine Learning.  In _12th $\\{$USENIX$\\}$ Symposium on Operating Systems Design and Implementation ($\\{$OSDI$\\}$ 16)_, pages 265–283, 2016.
  * Afsar et al. (2021) M Mehdi Afsar, Trafford Crump, and Behrouz Far.  Reinforcement learning based recommender systems: A survey.  _arXiv preprint arXiv:2101.06286_ , 2021.
  * Akiba et al. (2019) Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama.  Optuna: A next-generation hyperparameter optimization framework.  In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data mining_, pages 2623–2631, 2019.
  * Austin and Hayford (2021) David Austin and Tamara Hayford.  Research and development in the pharmaceutical industry, Apr 2021.  URL https://www.cbo.gov/publication/57025.
  * Barabási et al. (2011) Albert-László Barabási, Natali Gulbahce, and Joseph Loscalzo.  Network Medicine: a Network-based Approach to Human Disease.  _Nature Reviews Genetics_ , 12(1):56–68, 2011\.
  * Berner et al. (2019) Christopher Berner, Greg Brockman, Brooke Chan, Vicki Cheung, Przemysław Debiak, Christy Dennison, David Farhi, Quirin Fischer, Shariq Hashme, Chris Hesse, et al.  Dota 2 with large scale deep reinforcement learning.  _arXiv preprint arXiv:1912.06680_ , 2019.
  * Berrendorf et al. (2021) Max Berrendorf, Evgeniy Faerman, Laurent Vermue, and Volker Tresp.  On the Ambiguity of Rank-Based Evaluation of Entity Alignment or Link Prediction Methods, 2021.
  * Bonner et al. (2021) Stephen Bonner, Ian P Barrett, Cheng Ye, Rowan Swiers, Ola Engkvist, Andreas Bender, Charles Tapley Hoyt, and William Hamilton.  A review of biomedical datasets relating to drug discovery: A knowledge graph perspective, 2021.
  * Bordes et al. (2013) Antoine Bordes, Nicolas Usunier, Alberto Garcia-Duran, Jason Weston, and Oksana Yakhnenko.  Translating Embeddings for Modeling Multi-relational Data.  In C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems_ , volume 26. Curran Associates, Inc., 2013.
  * Breit et al. (2020) Anna Breit, Simon Ott, Asan Agibetov, and Matthias Samwald.  OpenBioLink: A Benchmarking Framework for Large-Scale Biomedical Link Prediction.  _Bioinformatics_ , 36(13):4097–4098, 2020.
  * Brown et al. (2020) Noam Brown, Anton Bakhtin, Adam Lerer, and Qucheng Gong.  Combining deep reinforcement learning and search for imperfect-information games.  _arXiv preprint arXiv:2007.13544_ , 2020.
  * Chen et al. (2019) Minmin Chen, Alex Beutel, Paul Covington, Sagar Jain, Francois Belletti, and Ed H Chi.  Top-k off-policy correction for a reinforce recommender system.  In _Proceedings of the Twelfth ACM International Conference on Web Search and Data Mining_ , pages 456–464, 2019.
  * Cheng et al. (2019) Feixiong Cheng, István A Kovács, and Albert-László Barabási.  Network-based Prediction of Drug Combinations.  _Nature Communications_ , 10(1):1–11, 2019\.
  * Das et al. (2018) Rajarshi Das, Shehzaad Dhuliawala, Manzil Zaheer, Luke Vilnis, Ishan Durugkar, Akshay Krishnamurthy, Alex Smola, and Andrew McCallum.  Go for a Walk and Arrive at the Answer: Reasoning Over Paths in Knowledge Bases Using Reinforcement Learning.  In _International Conference on Learning Representations_ , 2018.  URL https://openreview.net/forum?id=Syg-YfWCW.
  * Dettmers et al. (2018) Tim Dettmers, Pasquale Minervini, Pontus Stenetorp, and Sebastian Riedel.  Convolutional 2D Knowledge Graph Embeddings.  In _Thirty-Second AAAI Conference on Artificial Intelligence_ , 2018\.
  * Fey and Lenssen (2019) Matthias Fey and Jan Eric Lenssen.  Fast graph representation learning with PyTorch Geometric.  _arXiv preprint arXiv:1903.02428_ , 2019.
  * Geleta et al. (2021) David Geleta, Andriy Nikolov, Gavin Edwards, Anna Gogleva, Richard Jackson, Erik Jansson, Andrej Lamov, Sebastian Nilsson, Marina Pettersson, Vladimir Poroshin, Benedek Rozemberczki, Timothy Scrivener, Michael Ughetto, and Eliseo Papa.  Biological insights knowledge graph: an integrated knowledge graph to support drug development, 2021.  URL https://doi.org/10.1101/2021.10.28.466262.
  * Gluud et al. (2014) Lise L Gluud, Filip K Knop, and Tina Vilsbøll.  Effects of lixisenatide on elevated liver transaminases: systematic review with individual patient data meta-analysis of randomised controlled trials on patients with type 2 diabetes.  _BMJ open_ , 4(12):e005325, 2014.
  * Gogleva et al. (2021) Anna Gogleva, Dimitris Polychronopoulos, Matthias Pfeifer, Vladimir Poroshin, Michaël Ughetto, Ben Sidders, Jonathan Dry, Miika Ahdesmäki, Ultan McDermott, Eliseo Papa, et al.  Knowledge graph-based recommendation framework identifies novel drivers of resistance in egfr mutant non-small cell lung cancer.  _bioRxiv_ , 2021.
  * Guo et al. (2019) Lingbing Guo, Zequn Sun, and Wei Hu.  Learning to exploit long-term relational dependencies in knowledge graphs.  In _International Conference on Machine Learning_ , pages 2505–2514. PMLR, 2019.
  * Gysi et al. (2021) Deisy Morselli Gysi, Ítalo Do Valle, Marinka Zitnik, Asher Ameli, Xiao Gan, Onur Varol, Susan Dina Ghiassian, JJ Patten, Robert A Davey, Joseph Loscalzo, et al.  Network Medicine Framework for Identifying Drug-Repurposing Opportunities for COVID-19.  _Proceedings of the National Academy of Sciences_ , 118(19), 2021.
  * Haugeland (1989) John Haugeland.  _Artificial intelligence: The very idea_.  MIT press, 1989.
  * Hildebrandt et al. (2020) Marcel Hildebrandt, Jorge Andres Quintero Serna, Yunpu Ma, Martin Ringsquandl, Mitchell Joblin, and Volker Tresp.  Reasoning on knowledge graphs with debate dynamics.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pages 4123–4131, 2020.
  * Himmelstein et al. (2017) Daniel Scott Himmelstein, Antoine Lizee, Christine Hessler, Leo Brueggeman, Sabrina L Chen, Dexter Hadley, Ari Green, Pouya Khankhanian, and Sergio E Baranzini.  Systematic Integration of Biomedical Knowledge Prioritizes Drugs for Repurposing.  _Elife_ , 6:e26726, 2017.
  * Hu et al. (2020) Weihua Hu, Matthias Fey, Marinka Zitnik, Yuxiao Dong, Hongyu Ren, Bowen Liu, Michele Catasta, and Jure Leskovec.  Open Graph Benchmark: Datasets for Machine Learning on Graphs.  _arXiv preprint arXiv:2005.00687_ , 2020.
  * Ioannidis et al. (2020) Vassilis N Ioannidis, Xiang Song, Saurav Manchanda, Mufei Li, Xiaoqin Pan, Da Zheng, Xia Ning, Xiangxiang Zeng, and George Karypis.  DRKG - Drug Repurposing Knowledge Graph for COVID-19, 2020.
  * Iwuchukwu et al. (2014) Otito F Iwuchukwu, QiPing Feng, Wei-Qi Wei, Lan Jiang, Min Jiang, Hua Xu, Joshua C Denny, Russell A Wilke, Ronald M Krauss, Dan M Roden, et al.  Genetic variation in the ugt1a locus is associated with simvastatin efficacy in a clinical practice setting.  _Pharmacogenomics_ , 15(14):1739–1747, 2014.
  * Kingma and Ba (2017) Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization, 2017.
  * Kipf and Welling (2016) Thomas N Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  _arXiv preprint arXiv:1609.02907_ , 2016.
  * Lei et al. (2020) Deren Lei, Gangrong Jiang, Xiaotao Gu, Kexuan Sun, Yuning Mao, and Xiang Ren.  Learning Collaborative Agents with Rule Guidance for Knowledge Graph Reasoning.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP)_ , pages 8541–8547, 2020.
  * Lin et al. (2018) Xi Victoria Lin, Richard Socher, and Caiming Xiong.  Multi-hop knowledge graph reasoning with reward shaping.  _arXiv preprint arXiv:1808.10568_ , 2018.
  * Liu et al. (2021) Yushan Liu, Marcel Hildebrandt, Mitchell Joblin, Martin Ringsquandl, Rime Raissouni, and Volker Tresp.  Neural Multi-Hop Reasoning With Logical Rules on Biomedical Knowledge Graphs.  In _European Semantic Web Conference_ , pages 375–391. Springer, 2021\.
  * Lv et al. (2020) Xiaodan Lv, Yongqiang Dong, Lingling Hu, Feiyu Lu, Changyu Zhou, and Shaoyou Qin.  Glucagon-like peptide-1 receptor agonists (glp-1 ras) for the management of nonalcoholic fatty liver disease (nafld): A systematic review.  _Endocrinology, diabetes & metabolism_, 3(3):e00163, 2020.
  * Lv et al. (2021) Xin Lv, Yixin Cao, Lei Hou, Juanzi Li, Zhiyuan Liu, Yichi Zhang, and Zelin Dai.  Is multi-hop reasoning really explainable? towards benchmarking reasoning interpretability.  _arXiv preprint arXiv:2104.06751_ , 2021.
  * Meilicke et al. (2019) Christian Meilicke, Melisachew Wudage Chekol, Daniel Ruffinelli, and Heiner Stuckenschmidt.  Anytime Bottom-Up Rule Learning for Knowledge Graph Completion.  In _IJCAI_ , pages 3137–3143, 2019.
  * Mir et al. (2018) Olivier Mir, Alexandre Salvador, Sarah Dauchy, Stanislas Ropert, Cédric Lemogne, and Raphaël Gaillard.  Everolimus induced mood changes in breast cancer patients: a case-control study.  _Investigational new drugs_ , 36(3):503–508, 2018\.
  * Mnih et al. (2013) Volodymyr Mnih, Koray Kavukcuoglu, David Silver, Alex Graves, Ioannis Antonoglou, Daan Wierstra, and Martin Riedmiller.  Playing atari with deep reinforcement learning.  _arXiv preprint arXiv:1312.5602_ , 2013.
  * Nair and Hinton (2010) Vinod Nair and Geoffrey E Hinton.  Rectified linear units improve restricted boltzmann machines.  In _Icml_ , 2010.
  * Paszke et al. (2017) Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer.  Automatic differentiation in pytorch.  _NIPS_ , 2017.
  * Qu and Tang (2019) Meng Qu and Jian Tang.  Probabilistic logic neural networks for reasoning.  _arXiv preprint arXiv:1906.08495_ , 2019.
  * Rocktäschel and Riedel (2017) Tim Rocktäschel and Sebastian Riedel.  End-to-end differentiable proving.  _arXiv preprint arXiv:1705.11040_ , 2017.
  * Rosenblatt (1958) Frank Rosenblatt.  The perceptron: a probabilistic model for information storage and organization in the brain.  _Psychological review_ , 65(6):386, 1958.
  * Rumelhart et al. (1986) David E Rumelhart, Geoffrey E Hinton, and Ronald J Williams.  Learning representations by back-propagating errors.  _nature_ , 323(6088):533–536, 1986.
  * Santos et al. (2020) Alberto Santos, Ana Rita Colaço, Annelaura Bach Nielsen, Lili Niu, Philipp Emanuel Geyer, Fabian Coscia, Nicolai Jacob Wewer Albrechtsen, Filip Mundt, Lars Juhl Jensen, and Matthias Mann.  Clinical Knowledge Graph Integrates Proteomics Data Into Clinical Decision-Making.  _bioRxiv_ , 2020.
  * Schlichtkrull et al. (2018) 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.
  * Seghieri et al. (2018) Marta Seghieri, Alexander S Christensen, Andreas Andersen, Anna Solini, Filip K Knop, and Tina Vilsbøll.  Future perspectives on glp-1 receptor agonists and glp-1/glucagon receptor co-agonists in the treatment of nafld.  _Frontiers in endocrinology_ , 9:649, 2018.
  * Seko et al. (2017) Yuya Seko, Yoshio Sumida, Saiyu Tanaka, Kojiroh Mori, Hiroyoshi Taketani, Hiroshi Ishiba, Tasuku Hara, Akira Okajima, Atsushi Umemura, Taichiro Nishikawa, et al.  Effect of 12-week dulaglutide therapy in japanese patients with biopsy-proven non-alcoholic fatty liver disease and type 2 diabetes mellitus.  _Hepatology Research_ , 47(11):1206–1211, 2017\.
  * Shen et al. (2018) Yelong Shen, Jianshu Chen, Po-Sen Huang, Yuqing Guo, and Jianfeng Gao.  M-Walk: Learning to Walk Over Graphs Using Monte Carlo Tree Search.  In _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , pages 6787–6798, 2018.
  * Silver et al. (2017) David Silver, Julian Schrittwieser, Karen Simonyan, Ioannis Antonoglou, Aja Huang, Arthur Guez, Thomas Hubert, Lucas Baker, Matthew Lai, Adrian Bolton, et al.  Mastering the game of go without human knowledge.  _nature_ , 550(7676):354–359, 2017.
  * Srivastava et al. (2014) Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov.  Dropout: a simple way to prevent neural networks from overfitting.  _The journal of machine learning research_ , 15(1):1929–1958, 2014.
  * Sun et al. (2020) Zhiqing Sun, Shikhar Vashishth, Soumya Sanyal, Partha Talukdar, and Yiming Yang.  A Re-evaluation of Knowledge Graph Completion Methods, 2020.
  * Toutanova et al. (2015) Kristina Toutanova, Danqi Chen, Patrick Pantel, Hoifung Poon, Pallavi Choudhury, and Michael Gamon.  Representing Text for Joint Embedding of Text and Knowledge Bases.  In _Proceedings of the 2015 Conference on Empirical Methods in Natural Language Processing_ , pages 1499–1509, 2015.
  * Walsh et al. (2020) Brian Walsh, Sameh K Mohamed, and Vít Nováček.  BioKG: A Knowledge Graph for Relational Learning On Biological Data.  In _Proceedings of the 29th ACM International Conference on Information & Knowledge Management_, pages 3173–3180, 2020.
  * Wang et al. (2014) Zhen Wang, Jianwen Zhang, Jianlin Feng, and Zheng Chen.  Knowledge Graph Embedding by Translating on Hyperplanes.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 28, 2014.
  * Wray et al. (2007) Naomi R Wray, Michael R James, Steven P Mah, Matthew Nelson, Gavin Andrews, Patrick F Sullivan, Grant W Montgomery, Andrew J Birley, Andreas Braun, and Nicholas G Martin.  Anxiety and comorbid measures associated with plxna2.  _Archives of general psychiatry_ , 64(3):318–326, 2007.
  * Wu et al. (2020) Zonghan Wu, Shirui Pan, Fengwen Chen, Guodong Long, Chengqi Zhang, and S Yu Philip.  A comprehensive survey on graph neural networks.  _IEEE transactions on neural networks and learning systems_ , 32(1):4–24, 2020.
  * Xian et al. (2019) Yikun Xian, Zuohui Fu, Shan Muthukrishnan, Gerard De Melo, and Yongfeng Zhang.  Reinforcement Knowledge Graph Reasoning for Explainable Recommendation.  In _Proceedings of the 42nd International ACM SIGIR Conference on Research and Development in Information Retrieval_ , pages 285–294, 2019.
  * Xiong et al. (2017) Wenhan Xiong, Thien Hoang, and William Yang Wang.  DeepPath: A Reinforcement Learning Method for Knowledge Graph Reasoning.  In _Proceedings of the 2017 Conference on Empirical Methods in Natural Language Processing_ , pages 564–573, 2017.
  * Yang et al. (2014) Bishan Yang, Wen-tau Yih, Xiaodong He, Jianfeng Gao, and Li Deng.  Embedding entities and relations for learning and inference in knowledge bases.  _arXiv preprint arXiv:1412.6575_ , 2014.
  * Yang et al. (2017) Fan Yang, Zhilin Yang, and William W Cohen.  Differentiable learning of logical rules for knowledge base reasoning.  _arXiv preprint arXiv:1702.08367_ , 2017.
  * Zhang et al. (2021) Jing Zhang, Bo Chen, Lingxi Zhang, Xirui Ke, and Haipeng Ding.  Neural, Symbolic and Neural-symbolic Reasoning on Knowledge Graphs.  _AI Open_ , 2:14–35, 2021.
  * Zhao et al. (2020) Kangzhi Zhao, Xiting Wang, Yuren Zhang, Li Zhao, Zheng Liu, Chunxiao Xing, and Xing Xie.  Leveraging Demonstrations for Reinforcement Recommendation Reasoning Over Knowledge Graphs.  In _Proceedings of the 43rd International ACM SIGIR Conference on Research and Development in Information Retrieval_ , pages 239–248, 2020.
  * Zheng et al. (2020a) Da Zheng, Xiang Song, Chao Ma, Zeyuan Tan, Zihao Ye, Jin Dong, Hao Xiong, Zheng Zhang, and George Karypis.  DGL-KE: Training Knowledge Graph Embeddings at Scale.  In _Proceedings of the 43rd International ACM SIGIR Conference on Research and Development in Information Retrieval_ , pages 739–748, 2020a.
  * Zheng et al. (2020b) Shuangjia Zheng, Jiahua Rao, Ying Song, Jixian Zhang, Xianglu Xiao, Evandro Fei Fang, Yuedong Yang, and Zhangming Niu.  PharmKG: A Dedicated Knowledge Graph Benchmark for Biomedical Data Mining.  _Briefings in Bioinformatics_ , 2020b.
  * Zhou et al. (2020) Yadi Zhou, Yuan Hou, Jiayu Shen, Yin Huang, William Martin, and Feixiong Cheng.  Network-based Drug Repurposing for Novel Coronavirus 2019-nCoV/SARS-CoV-2.  _Cell Discovery_ , 6(1):1–18, 2020.

##  Appendix A Experimental Details

###  A.1 Model Implementations

Our experiments with PoLo and MINERVA utilized the TensorFlow (Abadi et al.,
2016) codebase released by the authors of PoLo. AnyBURL (Meilicke et al.,
2019) uses the latest version (AnyBURL-JUNO) from the authors. The DistMult
and TransE (L2 norm) embeddings were learned with the DGL-KE (Zheng et al.,
2020a) on top of the PyTorch (Paszke et al., 2017) automatic differentiation
backend library. The hyperparameters of these embedding models were tuned on
the validation set with the open-source Optuna framework (Akiba et al., 2019).
Finally, the R-GCN baselines were generated with the PyTorch Geometric library
(Fey and Lenssen, 2019) using Adam optimiser Kingma and Ba (2017), Dropout
Srivastava et al. (2014) and ReLU Nair and Hinton (2010).

###  A.2 Hyperparameters

To get the best results for each model, the Table X and Y show hyperparameters
that are searched. For MINERVA, $\lambda$ is set to 0 which recovers MINERVA
from PoLo.

Table 4: Hyperparameter search space for MINERVA and PoLo. For MINERVA $\lambda$ is set to 0. Hyperparameter | Values
---|---
Embedding size | {128, 256}
Hidden layer size | {256,512}
Learning rate | {0.0001,0.001}
$\lambda$ | {0.1,1}
$\beta$ | {0.01,0.1}
Table 5: Hyperparameter search space for knowledge graph embedding techniques TransE and DistMult. Hyperparameter | Values
---|---
Embedding size | {128, 256, …, 1024}
Learning rate | [0.001,0.2]
Negative samples | {1,…,500}
Max step | {12000, 22000, …, 200000}
Table 6: Hyperparameter search space for R-GCN. Hyperparameter | Values
---|---
Embedding size | $\\{32,64,128\\}$
Learning rate | $\\{0.01,0.05\\}$
Negative samples | $\\{10\\}$
Max step | $\\{0.01\\}$

#### AnyBURL

We use the most recent (JUNO) version of AnyBURL, which is more of a
neurosymbolic approach as it implements Reinforcement Learning to learn rules.
The maximum length of rules is set to 3 and the rules are learned for a total
of 500 seconds. The AnyBURL algorithm does not need hyperparameter tuning to
achieve good results, therefore we kept the default values of remaining
hyperparameters.

###  A.3 Dataset Statistics

Table 7: Statistics of each tasks dataset sizes Task | Dataset | Relation | Nodes | Edges | Train | Valid | Test
---|---|---|---|---|---|---|---
Drug Repurposing | Hetionet | CtD | 47,031 | 2,250,197 | 483 | 121 | 151
BIKG Hetionet | CtD | 131,668 | 11,289,713 | 2,062 | 516 | 664
BIKG Hetionet+ | CtD | 556,219 | 25,762,310 | 2,114 | 529 | 661
Drug-Target Interaction | Hetionet | CbG | 47,031 | 2,250,197 | 7,404 | 1,852 | 2,315
BIKG Hetionet | CbG | 131,668 | 11,289,713 | 7,351 | 1,838 | 2,298
BIKG Hetionet+ | CbG | 556,219 | 25,762,310 | 7,351 | 1,838 | 2,298
Table 8: Statistics of biomedical and benchmark knowledge graphs Dataset | Domain | Node Types | Nodes | Edge Types | Edges | Mean Degree
---|---|---|---|---|---|---
ClinicalKG | Biological | 35 | 19,251,579 | 57 | 208,177,953 | 21.63
BIKG | Biological | 27 | 10,948,027 | 52 | 83,849,252 | 17.87
BIKG Hetionet+ | Biological | 12 | 556,219 | 47 | 25,762,310 | 92.63
BIKG Hetionet | Biological | 11 | 131,668 | 24 | 11,289,713 | 171.49
DRKG | Biological | 13 | 97,238 | 107 | 5,874,258 | 120.81
OGBL-BIOKG | Biological | 5 | 45,085 | 51 | 4,762,677 | 108.52
OpenBioLink-HQ | Biological | 7 | 184,667 | 30 | 4,778,683 | 51.80
Hetionet | Biological | 11 | 47,031 | 24 | 2,250,197 | 95.83
BioKG | Biological | 10 | 105,524 | 17 | 2,067,998 | 39.19
PharmKG-8k | Biological | 3 | 7,262 | 29 | 479,902 | 132.16
FB15K-237 | General | 571 | 14,505 | 237 | 310,116 | 42.50
WN18RR | Lexical | 4 | 40,945 | 11 | 90,003 | 4.53
NELL-995 | General | 269 | 75,492 | 200 | 154,213 | 4.09

##  Appendix B Explanations

The following section contains metapath statistics and explanations generated
by MINVERA on the BIKG Hetionet and BIKG Hetionet+ datasets.

###  B.1 Metapaths Statistics

Table 9: MINERVA’s top 10 frequent reasoning metapaths for drug repurposing (pruned) Dataset | Metapath | %
---|---|---
BIKG Hetionet | $Compound\xrightarrow{\text{treats}}Disease\xrightarrow{\text{associates}^{-1}}Gene\xrightarrow{\text{associates}}Disease$ | 71.7
$Compound\xrightarrow{\text{includes}^{-1}}Pharmacologic\ Class\xrightarrow{\text{includes}}Compound\xrightarrow{\text{treats}}Disease$ | 5.1
$Compound\xrightarrow{\text{treats}}Disease\xrightarrow{\text{upregulates}}Gene\xrightarrow{\text{associates}}Disease$ | 4.3
$Compound\xrightarrow{\text{treats}}Disease\xrightarrow{\text{treats}^{-1}}Compound\xrightarrow{\text{treats}}Disease$ | 2.7
$Compound\xrightarrow{\text{palliates}}Disease\xrightarrow{\text{upregulates}}Gene\xrightarrow{\text{associates}}Disease$ | 2.1
$Compound\xrightarrow{\text{resembles}}Compound\xrightarrow{\text{resembles}}Compound\xrightarrow{\text{treats}}Disease$ | 1.5
$Compound\xrightarrow{\text{binds}}Gene\xrightarrow{\text{associates}}Disease$ | 1.2
$Compound\xrightarrow{\text{resembles}}Compound\xrightarrow{\text{resembles}^{-1}}Compound\xrightarrow{\text{treats}}Disease$ | 1.2
$Compound\xrightarrow{\text{palliates}}Disease\xrightarrow{\text{associates}^{-1}}Gene\xrightarrow{\text{associates}}Disease$ | 1.2
$Compound\xrightarrow{\text{treats}}Disease\xrightarrow{\text{localizes}}Anatomy\xrightarrow{\text{localizes}^{-1}}Disease^{\prime}$ | 0.9
BIKG Hetionet+ | $Compound\xrightarrow{\text{has mechanism}}Mechanism\ of\ Action\xrightarrow{\text{has mechanism}^{-1}}Compound\xrightarrow{\text{has link}}Disease$ | 51.5
$Compound\xrightarrow{\text{has mechanism}}Mechanism\ of\ Action\xrightarrow{\text{has assay target}}Gene\xrightarrow{\text{associates}}Disease$ | 34.2
$Compound\xrightarrow{\text{has link}}Disease\xrightarrow{\text{has link}^{-1}}Disease\xrightarrow{\text{has link}}Disease$ | 9.2
$Compound\xrightarrow{\text{has mechanism}}Mechanism\ of\ Action\xrightarrow{\text{has assay target}}Gene\xrightarrow{\text{has link}}Disease$ | 1.6
$Compound\xrightarrow{\text{treats}}Disease\xrightarrow{\text{has link}^{-1}}Compound\xrightarrow{\text{has link}}Disease$ | 1.3
$Compound\xrightarrow{\text{treats}}Disease\xrightarrow{\text{associates}^{-1}}Gene\xrightarrow{\text{associates}}Disease$ | 0.4
$Compound\xrightarrow{\text{resembles}}Compound\xrightarrow{\text{binds}}Gene\xrightarrow{\text{associates}}Disease$ | 0.3
$Compound\xrightarrow{\text{has link}}Disease\xrightarrow{\text{involved in}^{-1}}Biological\ Process\xrightarrow{\text{involved in}}Disease$ | 0.2
$Compound\xrightarrow{\text{resembles}^{-1}}Compound\xrightarrow{\text{binds}}Gene\xrightarrow{\text{associates}}Disease$ | 0.2
$Compound\xrightarrow{\text{treats}}Disease\xrightarrow{\text{involved in}^{-1}}Biological\ Process\xrightarrow{\text{involved in}}Disease$ | 0.1
Table 10: MINERVA’s top 10 frequent reasoning metapaths for drug target interaction (pruned) Dataset | Metapath | %
---|---|---
BIKG Hetionet | $Compound\xrightarrow{\text{binds}}Gene\xrightarrow{\text{participates}}Pathway\xrightarrow{\text{participates}^{-1}}Gene$ | 14.2
$Compound\xrightarrow{\text{resembles}^{-1}}Compound\xrightarrow{\text{resembles}}Compound\xrightarrow{\text{binds}}Gene$ | 13.8
$Compound\xrightarrow{\text{includes}^{-1}}Pharmacologic\ Class\xrightarrow{\text{includes}}Compound\xrightarrow{\text{binds}}Gene$ | 13.2
$Compound\xrightarrow{\text{resembles}}Compound\xrightarrow{\text{resembles}}Compound\xrightarrow{\text{binds}}Gene$ | 12.4
$Compound\xrightarrow{\text{resembles}}Compound\xrightarrow{\text{resembles}^{-1}}Compound\xrightarrow{\text{binds}}Gene$ | 12.2
$Compound\xrightarrow{\text{binds}}Gene\xrightarrow{\text{participates}}Molecular\ Function\xrightarrow{\text{participates}^{-1}}Gene$ | 4.8
$Compound\xrightarrow{\text{resembles}}Compound\xrightarrow{\text{binds}}Gene$ | 3.5
$Compound\xrightarrow{\text{resembles}}Compound\xrightarrow{\text{binds}}Gene$ | 3.3
$Compound\xrightarrow{\text{binds}}Gene\xrightarrow{\text{participates}}Biological\ Process\xrightarrow{\text{participates}}Gene$ | 3.2
$Compound\xrightarrow{\text{resembles}^{-1}}Compound\xrightarrow{\text{resembles}^{-1}}Compound\xrightarrow{\text{binds}}Gene$ | 3.2
BIKG Hetionet+ | $Compound\xrightarrow{\text{has mechanism}}Mechanism\ of\ Action\xrightarrow{\text{has assay target}}Gene\xrightarrow{\text{has link}}Gene$ | 47.6
$Compound\xrightarrow{\text{has mechanism}}Mechanism\ of\ Action\xrightarrow{\text{has assay target}}Gene$ | 18.8
$Compound\xrightarrow{\text{has mechanism}}Mechanism\ of\ Action\xrightarrow{\text{has assay target}}Gene\xrightarrow{\text{interacts}}Gene$ | 10.1
$Compound\xrightarrow{\text{binds}}Gene\xrightarrow{\text{has link}}Gene\xrightarrow{\text{has link}}Gene$ | 5.9
$Compound\xrightarrow{\text{has link}}Gene\xrightarrow{\text{has link}}Gene\xrightarrow{\text{has link}}Gene$ | 3.0
$Compound\xrightarrow{\text{binds}}Gene\xrightarrow{\text{has link}}Gene$ | 2.2
$Compound\xrightarrow{\text{associates}}Pathway\xrightarrow{\text{involved in}}Disease\xrightarrow{\text{associates}^{-1}}Gene$ | 1.4
$Compound\xrightarrow{\text{resembles}^{-1}}Compound\xrightarrow{\text{binds}}Gene\xrightarrow{\text{has link}}Gene$ | 1.4
$Compound\xrightarrow{\text{upregulates}}Gene\xrightarrow{\text{has link}}Gene\xrightarrow{\text{has link}}Gene$ | 1.0
$Compound\xrightarrow{\text{has mechanism}}Mechanism\ of\ Action\xrightarrow{\text{has assay target}}Gene\xrightarrow{\text{has link}}Gene$ | 0.8

###  B.2 Explanations

Figure 7: Explainable reasoning for recommending Benralizumab in the treatment
of Asthma. Figure 8: Explainable reasoning for recommending Inotuzumab
ozogamicin in the treatment of B-cell Chronic Lymphocytic Leukemia. Figure 9:
Explainable reasoning for recommending Naproxen Sodium in the treatment of
rheumatic diseases. Figure 10: Explainable reasoning for recommending
Quinapril as an ACE inhibitor.

##  Appendix C Pre-pruned Results

Below are the results before pruning of expected target node types was
applied. This effectively keeps failed walks and predictions in the results
sets.

Table 11: Drug repurposing results (pre-pruning) Dataset | Model | HITS@1 | HITS@3 | HITS@10 | MRR
---|---|---|---|---|---
Hetionet | TransE | $.193\pm.027$ | $.364\pm.031$ | $.608\pm.034$ | $.312\pm.022$
DistMult | $.020\pm.014$ | $.086\pm.023$ | $.270\pm.022$ | $.096\pm.015$
RGCN | $.053\pm.027$ | $.200\pm.016$ | $.454\pm.036$ | $.217\pm.041$
AnyBURL | $.116\pm.013$ | $.288\pm.024$ | $.572\pm.036$ | $.258\pm.012$
MINERVA | $\mathit{.281\pm.027}$ | $\mathit{.462\pm.055}$ | $\mathit{.701\pm.055}$ | $\mathit{.391\pm.034}$
PoLo | $\mathbf{.306\pm.017}$ | $\mathbf{.481\pm.029}$ | $\mathbf{.742\pm.013}$ | $\mathbf{.406\pm.015}$
BIKG Hetionet | TransE | $.016\pm.003$ | $.034\pm.005$ | $.088\pm.001$ | $.045\pm.004$
DistMult | $.012\pm.002$ | $.021\pm.003$ | $.046\pm.006$ | $.029\pm.003$
RGCN | $.056\pm.005$ | $.108\pm.008$ | $\mathit{.229\pm.011}$ | $.115\pm.004$
AnyBURL | $.000\pm.000$ | $.056\pm.005$ | $.205\pm.014$ | $.072\pm.002$
MINERVA | $\mathit{.070\pm.007}$ | $\mathit{.137\pm.011}$ | $.204\pm.020$ | $\mathit{.116\pm.009}$
PoLo | $\mathbf{.089\pm.005}$ | $\mathbf{.162\pm.007}$ | $\mathbf{.239\pm.008}$ | $\mathbf{.142\pm.005}$
BIKG Hetionet+ | TransE | $.017\pm.002$ | $.052\pm.004$ | $.119\pm.011$ | $.057\pm.002$
DistMult | $.004\pm.002$ | $.013\pm.004$ | $.036\pm.009$ | $.025\pm.004$
RGCN | $.039\pm.008$ | $.091\pm.014$ | $\mathit{.191\pm.012}$ | $\mathbf{.109\pm.060}$
AnyBURL | $.000\pm.000$ | $.065\pm.006$ | $\mathbf{.205\pm.008}$ | $.070\pm.009$
MINERVA | $\mathbf{.069\pm.008}$ | $\mathit{.115\pm.010}$ | $.151\pm.011$ | $\mathit{.105\pm.009}$
PoLo | $\mathit{.068\pm.007}$ | $\mathbf{.117\pm.011}$ | $.157\pm.011$ | $\mathit{.105\pm.008}$
Table 12: Drug repurposing results (pruned) Dataset | Model | HITS@1 | HITS@3 | HITS@10 | MRR
---|---|---|---|---|---
Hetionet | TransE | $.193\pm.027$ | $.364\pm.031$ | $.608\pm.034$ | $.312\pm.022$
DistMult | $.028\pm.002$ | $.138\pm.022$ | $.373\pm.039$ | $.129\pm.002$
RGCN | $.047\pm.028$ | $.194\pm.018$ | $.492\pm.043$ | $.221\pm.040$
AnyBURL | $.116\pm.013$ | $.288\pm.024$ | $.572\pm.036$ | $.258\pm.012$
MINERVA | $\mathit{.351\pm.040}$ | $\mathit{.559\pm.058}$ | $\mathit{.790\pm.057}$ | $\mathit{.463\pm.041}$
PoLo | $\mathbf{.378\pm.023}$ | $\mathbf{.575\pm.025}$ | $\mathbf{.845\pm.009}$ | $\mathbf{.479\pm.015}$
BIKG Hetionet | TransE | $.016\pm.003$ | $.034\pm.005$ | $.088\pm.001$ | $.045\pm.004$
DistMult | $.004\pm.002$ | $.013\pm.002$ | $.036\pm.007$ | $.024\pm.003$
RGCN | $.056\pm.005$ | $.108\pm.008$ | $.239\pm.008$ | $.117\pm.005$
AnyBURL | $.000\pm.000$ | $.015\pm.001$ | $.046\pm.004$ | $.023\pm.001$
MINERVA | $\mathit{.121\pm.002}$ | $\mathit{.194\pm.010}$ | $\mathit{.249\pm.013}$ | $\mathit{.160\pm.004}$
PoLo | $\mathbf{.124\pm.005}$ | $\mathbf{.204\pm.010}$ | $\mathbf{.273\pm.016}$ | $\mathbf{.167\pm.008}$
BIKG Hetionet+ | TransE | $.017\pm.002$ | $.052\pm.004$ | $.119\pm.011$ | $.057\pm.002$
DistMult | $.004\pm.002$ | $.014\pm.004$ | $.039\pm.001$ | $.026\pm.004$
RGCN | $.039\pm.008$ | $.104\pm.016$ | $\mathbf{.222\pm.015}$ | $\mathbf{.115\pm.005}$
AnyBURL | $.000\pm.000$ | $.056\pm.006$ | $\mathit{.205\pm.014}$ | $.072\pm.002$
MINERVA | $\mathbf{.076\pm.008}$ | $\mathit{.124\pm.011}$ | $.166\pm.013$ | $\mathit{.113\pm.010}$
PoLo | $\mathit{.074\pm.006}$ | $\mathbf{.133\pm.010}$ | $.173\pm.012$ | $\mathbf{.115\pm.007}$
Table 13: Drug-target interaction results (pre-pruning) Dataset | Model | HITS@1 | HITS@3 | HITS@10 | MRR
---|---|---|---|---|---
Hetionet | TransE | $\mathbf{.287\pm.008}$ | $\mathbf{.546\pm.013}$ | $\mathbf{.918\pm.007}$ | $\mathbf{.331\pm.005}$
DistMult | $.037\pm.004$ | $.120\pm.006$ | $.355\pm.010$ | $.111\pm.004$
RGCN | $.041\pm.017$ | $.123\pm.035$ | $.320\pm.087$ | $.104\pm.024$
AnyBURL | $.068\pm.018$ | $.128\pm.033$ | $.195\pm.049$ | $.210\pm.051$
MINERVA | $\mathit{.159\pm.015}$ | $\mathit{.373\pm.037}$ | $\mathit{.756\pm.045}$ | $\mathit{.269\pm.021}$
BIKG Hetionet | TransE | $\mathbf{.193\pm.010}$ | $\mathbf{.376\pm.012}$ | $\mathbf{.674\pm.012}$ | $\mathit{.239\pm.004}$
DistMult | $.043\pm.007$ | $.106\pm.015$ | $.258\pm.016$ | $.095\pm.007$
RGCN | $.057\pm.025$ | $.142\pm.017$ | $.270\pm.028$ | $.123\pm.012$
AnyBURL | $.088\pm.003$ | $.167\pm.004$ | $.254\pm.007$ | $\mathbf{.264\pm.004}$
MINERVA | $\mathit{.104\pm.011}$ | $\mathit{.224\pm.014}$ | $\mathit{.589\pm.013}$ | $.172\pm.008$
BIKG Hetionet+ | TransE | $.054\pm.004$ | $.098\pm.007$ | $.171\pm.011$ | $.080\pm.004$
DistMult | $.035\pm.006$ | $.083\pm.011$ | $.202\pm.028$ | $.083\pm.001$
RGCN | $\mathit{.058\pm.025}$ | $\mathit{.142\pm.017}$ | $\mathbf{.270\pm.022}$ | $\mathit{.123\pm.012}$
AnyBURL | $\mathbf{.086\pm.002}$ | $\mathbf{.162\pm.003}$ | $.251\pm.007$ | $\mathbf{.263\pm.003}$
MINERVA | $.026\pm.002$ | $.079\pm.006$ | $\mathit{.255\pm.027}$ | $.085\pm.004$
Table 14: Drug-target interaction results (pruned) Dataset | Model | HITS@1 | HITS@3 | HITS@10 | MRR
---|---|---|---|---|---
Hetionet | TransE | $\mathbf{.287\pm.008}$ | $\mathbf{.546\pm.013}$ | $\mathbf{.918\pm.007}$ | $\mathbf{.331\pm.005}$
DistMult | $.039\pm.004$ | $.128\pm.007$ | $.367\pm.012$ | $.115\pm.005$
RGCN | $.060\pm.025$ | $.177\pm.057$ | $.376\pm.102$ | $.126\pm.031$
AnyBURL | $.169\pm.043$ | $.323\pm.081$ | $.493\pm.124$ | $.210\pm.051$
MINERVA | $\mathit{.186\pm.015}$ | $\mathit{.420\pm.028}$ | $\mathit{.807\pm.039}$ | $\mathit{.296\pm.019}$
BIKG Hetionet | TransE | $.193\pm.001$ | $.376\pm.012$ | $\mathit{.674\pm.012}$ | $.239\pm.004$
DistMult | $.044\pm.007$ | $.111\pm.014$ | $.269\pm.014$ | $.098\pm.007$
RGCN | $.173\pm.016$ | $.338\pm.013$ | $.539\pm.021$ | $.235\pm.008$
AnyBURL | $\mathit{.215\pm.005}$ | $\mathit{.408\pm.004}$ | $.623\pm.009$ | $\mathit{.264\pm.004}$
MINERVA | $\mathbf{.235\pm.005}$ | $\mathbf{.516\pm.011}$ | $\mathbf{.983\pm.015}$ | $\mathbf{.305\pm.004}$
BIKG Hetionet+ | TransE | $.054\pm.004$ | $.098\pm.007$ | $.171\pm.011$ | $.008\pm.004$
DistMult | $.035\pm.006$ | $.084\pm.011$ | $.206\pm.029$ | $.084\pm.001$
RGCN | $\mathit{.188\pm.008}$ | $.323\pm.019$ | $.539\pm.010$ | $.241\pm.006$
AnyBURL | $\mathbf{.215\pm.003}$ | $\mathbf{.408\pm.008}$ | $\mathit{.622\pm.005}$ | $\mathbf{.263\pm.003}$
MINERVA | $.181\pm.006$ | $\mathit{.368\pm.008}$ | $\mathbf{.628\pm.019}$ | $\mathit{.243\pm.004}$


# An asymptotic approximation for the Riemann zeta function revisited

R. B. Paris
Division of Computing and Mathematics,
Abertay University, Dundee DD1 1HG, UK

###### Abstract

We revisit a representation for the Riemann zeta function $\zeta(s)$ expressed
in terms of normalised incomplete gamma functions given by the author and S.
Cang in Methods Appl. Anal. 4 (1997) 449–470. Use of the uniform asymptotics
of the incomplete gamma function produces an asymptotic-like expansion for
$\zeta(s)$ on the critical line $s=1/2+it$ as $t\to+\infty$. The main term
involves the original Dirichlet series smoothed by a complementary error
function of appropriate argument together with a series of correction terms.
It is the aim here to present these correction terms in a more user-friendly
format by expressing then in inverse powers of $\omega$, where $\omega^{2}=\pi
s/(2i)$, multiplied by coefficients involving trigonometric functions of
argument $\omega$.

Mathematics subject classification (2010): 11M06, 33B20, 34E05, 41A60

Keywords: Riemann zeta function, critical line, asymptotic expansion,
incomplete gamma function, modified complementary error function

1\. Introduction

The computation of the Riemann zeta function on the critical line
$s=\mbox{${\textstyle\frac{1}{2}}$}+it$, $t\geq 0$ is normally carried out
using the real function
$Z(t)=e^{i\vartheta(t)}\zeta(\mbox{${\textstyle\frac{1}{2}}$}+it)$, where the
phase angle

$\vartheta(t)=\arg\,\Gamma(\mbox{${\textstyle\frac{1}{4}}$}+\mbox{${\textstyle\frac{1}{2}}$}it)-\mbox{${\textstyle\frac{1}{2}}$}t\log\,\pi.$

In [3], a representation of $Z(t)$ on the critical line was given in the form

$Z(t)=2\Re
e^{i\vartheta(t)}\biggl{\\{}\sum_{n=1}^{\infty}n^{-s}Q(\mbox{${\textstyle\frac{1}{2}}$}s,\pi
n^{2}i)-\frac{\pi^{s/2}e^{\pi
is/4}}{s\Gamma(\mbox{${\textstyle\frac{1}{2}}$}s)}\biggr{\\}},$ (1.1)

wherein $s=\mbox{${\textstyle\frac{1}{2}}$}+it$, $t\geq 0$ and
$Q(a,z)=\Gamma(a,z)/\Gamma(a)$ is the normalised (upper) incomplete gamma
function. The behaviour of $Q(\mbox{${\textstyle\frac{1}{2}}$}s,\pi n^{2}i)$
for large $t$ changes abruptly in the neighbourhood of its transition point
$\mbox{${\textstyle\frac{1}{2}}$}s=\pi n^{2}i$; that is, when the summation
index $n$ roughly equals the Riemann-Siegel cut-off value $N_{t}$ given by

$N_{t}=[(t/2\pi)^{1/2}],$ (1.2)

where square brackets denote the integer part. Then for large $t$, the
function $|Q(\mbox{${\textstyle\frac{1}{2}}$}s,\pi n^{2}i)|\simeq 1$ when
$n\,\raisebox{-3.44444pt}{\mbox{$\stackrel{{\scriptstyle\textstyle<}}{{\sim}}$}}\,N_{t}$
and decays (algebraically) to zero when
$n\,\raisebox{-3.44444pt}{\mbox{$\stackrel{{\scriptstyle\textstyle>}}{{\sim}}$}}\,N_{t}$.
As a consequence, the absolutely convergent sum111The sum in (1.1) converges
absolutely since $\Gamma(a,z)\sim z^{a-1}e^{-z}$ as $|z|\to\infty$ in
$|\arg\,z|<3\pi/4$, and so late terms behave like $(-)^{n}/n^{2}$. in (1.1)
represents the original Dirichlet series “smoothed” by the incomplete gamma
function.

The uniform asymptotic expansion of $Q(a,z)$ valid for $a\to\infty$ when
$|z|\in[0,\infty)$ in the domains $|\arg\,z|<\pi$ and $|\arg\,(z/a)|<2\pi$ was
employed in [3] to derive an asymptotic approximation for $Z(t)$ valid as
$t\to\infty$. The resulting expansion involved the Dirichlet sum smoothed by a
complementary error function of appropriate argument together with asymptotic
correction terms in inverse powers of $\omega$, where $\omega^{2}=\pi s/(2i)$,
decorated by functions $S_{k}(\omega)$ involving $\csc\omega$ and its
derivatives. In [2] an attempt was made to establish the asymptotic nature of
this expansion; this was not successful, however, due to the use of an
insufficiently precise error bound for the incomplete gamma function.

In [3] the quantities $S_{k}(\omega)$ were treated as computable functions.
The purpose of the present paper is to revisit the expansion (1.1) with the
aim of presenting the coefficients in the correction terms in a more user-
friendly format. This is achieved using the expansion of the functions
$S_{k}(\omega)$ as polynomials in $1/\omega$ of degree $k-1$.

2\. The asymptotic approximation for $Z(t)$

To make this paper self-contained, we summarise in this section the main steps
in the derivation of the asymptotic approximation derived from (1.1) and
presented in [3]. We employ the uniform asymptotic expansion of the normalised
incomplete gamma function valid as $a\to\infty$ in $|\arg\,a|<\pi$ in the form
[1, p. 181]

$Q(a,z)=\frac{1}{2}\mbox{erfc}\,(\eta\sqrt{a/2})+\frac{e^{-\frac{1}{2}a\eta^{2}}}{\sqrt{2\pi
a}}\biggl{\\{}\sum_{k=0}^{m-1}c_{k}(\eta)a^{-k}+a^{-m}R_{m}(a,\eta)\biggr{\\}}$
(2.1)

for $m=1,2,\ldots$, where

$\mbox{${\textstyle\frac{1}{2}}$}\eta^{2}=\lambda-1-\log\,\lambda,\quad\lambda=\frac{z}{a},\quad\mu=\lambda-1$
(2.2)

and the quantity $R_{m}(a,\eta)$ is the remainder term in the expansion
truncated after $m$ terms. The choice of the square root branch for
$\eta(\lambda)$ is made such that $\eta(\lambda)$ and $\lambda-1$ have the
same sign when $\lambda>0$; we then have $\eta\simeq\lambda-1$ when
$\lambda\simeq 1$. The coefficients $c_{k}(\eta)$ are given by

$c_{k}(\eta)=\frac{(-)^{k}W_{k}(\mu)}{\mu^{2k+1}}-\frac{(-2)^{k}(\mbox{${\textstyle\frac{1}{2}}$})_{k}}{\eta^{2k+1}},\qquad
W_{k}(\mu):=\sum_{r=0}^{2k}\alpha_{r,k}\mu^{r},$ (2.3)

where $W_{k}(\mu)$ is a polynomial in $\mu$ of degree $2k$ and the
coefficients $\alpha_{r,k}$ are discussed in [3, Appendix]; see also [4]. The
first few coefficients are given by

$\displaystyle\alpha_{00}$ $\displaystyle=$ $\displaystyle 1$
$\displaystyle\alpha_{01}$ $\displaystyle=$ $\displaystyle 1,\ \
\alpha_{11}=1,\ \ \alpha_{21}=\mbox{${\textstyle\frac{1}{12}}$}$
$\displaystyle\alpha_{02}$ $\displaystyle=$ $\displaystyle 3,\ \
\alpha_{12}=5,\ \ \alpha_{22}=\mbox{${\textstyle\frac{25}{12}}$},\ \
\alpha_{32}=\mbox{${\textstyle\frac{1}{12}}$},\ \
\alpha_{42}=\mbox{${\textstyle\frac{1}{288}}$}$ $\displaystyle\alpha_{03}$
$\displaystyle=$ $\displaystyle 15,\ \alpha_{13}=35,\ \
\alpha_{23}=\mbox{${\textstyle\frac{105}{4}}$},\ \
\alpha_{33}=\mbox{${\textstyle\frac{77}{12}}$},\ \
\alpha_{43}=\mbox{${\textstyle\frac{49}{288}}$},\ \
\alpha_{53}=\mbox{${\textstyle\frac{1}{288}}$},\ \
\alpha_{63}=-\mbox{${\textstyle\frac{139}{51840}}$}.$

The modified complementary error function is introduced by

$\mbox{erfc}\,(z;m)=\mbox{erfc}\,z-\frac{e^{-z^{2}}}{\sqrt{\pi}}\sum_{r=0}^{m-1}\frac{(-)^{r}(\mbox{${\textstyle\frac{1}{2}}$})_{r}}{z^{2r+1}}\qquad(z\neq
0;\ m=1,2,\ldots\,),$ (2.4)

which corresponds to the deletion from $\mbox{erfc}\,z$ of the first $m$ terms
of its asymptotic expansion for $|z|\to\infty$ in
$|\arg\,z|<\mbox{${\textstyle\frac{3}{4}}$}\pi$. Then the expansion (2.1) can
rewritten in the form

$Q(a,z)=\frac{1}{2}\mbox{erfc}\,(\eta\sqrt{a/2};m)+\frac{e^{-\frac{1}{2}a\eta^{2}}}{\sqrt{2\pi
a}}\biggl{\\{}\sum_{k=0}^{m-1}\frac{(-)^{k}W_{k}(\mu)}{\mu^{2k+1}}a^{-k}+a^{-m}R_{m}(a,\eta)\biggr{\\}}.$
(2.5)

We now substitute the expansion (2.5) into (1.1) when we identify the
parameters $a$, $\lambda\equiv\lambda_{n}$ and $\mu\equiv\mu_{n}$ (and
correspondingly $\eta\equiv\eta_{n}$) with

$a=\mbox{${\textstyle\frac{1}{2}}$}s=\mbox{${\textstyle\frac{1}{4}}$}+\mbox{${\textstyle\frac{1}{2}}$}it,\quad\lambda_{n}=2\pi
n^{2}i/s,\quad\mu_{n}=(\pi n/\omega)^{2}-1,\quad\omega^{2}=\pi s/(2i).$ (2.6)

We also make use of the well-known expansion

$\frac{1}{\Gamma(\mbox{${\textstyle\frac{1}{2}}$}s)}=(2\pi)^{-1/2}(\mbox{${\textstyle\frac{1}{2}}$}s)^{(1-s)/2}e^{s/2}\biggl{\\{}\sum_{r=0}^{m-1}\gamma_{r}(\mbox{${\textstyle\frac{1}{2}}$}s)^{-r}+(\mbox{${\textstyle\frac{1}{2}}$}s)^{-m}H_{m}(s)\biggr{\\}}\qquad(t\to\infty),$

where the $\gamma_{r}$ are the Stirling coefficients with $\gamma_{1}=1$,
$\gamma_{2}=-\mbox{${\textstyle\frac{1}{12}}$}$,
$\gamma_{2}=\mbox{${\textstyle\frac{1}{288}}$}$,
$\gamma_{3}=\mbox{${\textstyle\frac{139}{51840}}$},\ldots\ ,$ and $H_{m}(s)$
is a remainder term. Since

$e^{-\frac{1}{2}a\eta_{n}^{2}}=(-)^{n}n^{s}e^{\pi is/4}(2\pi e/s)^{s/2},$

we then find that

$Z(t)=\Re
e^{i\vartheta(t)}\biggl{\\{}\sum_{n=1}^{\infty}n^{-s}\mbox{erfc}\,(\mbox{${\textstyle\frac{1}{2}}$}\eta_{n}\sqrt{s};m)-\frac{e^{\pi
is/4}}{\sqrt{\pi s}}\,\biggl{(}\frac{2\pi
e}{s}\biggr{)}^{s/2}\biggl{(}\sum_{r=0}^{m-1}\frac{(-)^{r}A_{r}(s)}{(\mbox{${\textstyle\frac{1}{2}}$}s)^{r}}+(\mbox{${\textstyle\frac{1}{2}}$}s)^{-m}{\hat{R}}_{m}\biggr{)}\biggr{\\}},$
(2.7)

where ${\hat{R}}_{m}$ is the remainder that contains contributions from the
remainder term in $Q(a,z)$ and that in the expansion of
$1/\Gamma(\mbox{${\textstyle\frac{1}{2}}$}s)$. We do not consider this term
here; an attempt was made in [2] to bound ${\hat{R}}_{m}$ and demonstrate that
(2.7) is an asymptotic expansion, but the bounds employed on the incomplete
gamma function were not sharp enough to achieve this.

The coefficients $A_{r}(s)$ are given by

$A_{r}(s)=(-)^{r}\gamma_{r}-2\sum_{n=1}^{\infty}\frac{(-)^{n}W_{r}(\mu_{n})}{\mu_{n}^{2r+1}}=\sum_{n=-\infty}^{\infty}\frac{(-)^{n-1}W_{r}(\mu_{n})}{\mu_{n}^{2r+1}}$

where $W_{r}(\mu_{0})=W_{r}(-1)=(-)^{r}\gamma_{r}$; see [3, (A.1)]. From the
definitions of $W_{k}(\mu_{n})$ and $\mu_{n}$ in (2.3) and (2.6), we then find
that

$A_{r}(s)=\sum_{k=0}^{2r}(-)^{k}\alpha_{r,k}\sum_{n=-\infty}^{\infty}(-)^{n}\biggl{(}\frac{\omega^{2}}{\omega^{2}-(\pi
n)^{2}}\biggr{)}^{\\!2r+1-k}.$

The inner sum may be expressed in terms of the functions $S_{k}(\omega)$
defined by

$S_{k}(\omega):=2^{k-1}\sum_{n=-\infty}^{\infty}\frac{(-)^{n}\omega^{k}}{(\omega^{2}-(\pi
n)^{2})^{k}}\qquad k=1,2,\ldots$ (2.8)

so that

$A_{r}(s)=2^{-2r}\omega^{2r+1}B_{r}(\omega),\qquad
B_{r}(\omega)=\sum_{k=0}^{2r}(-)^{k}\alpha_{r,k}\,\frac{S_{2r+1-k}(\omega)}{(\omega/2)^{k}}.$
(2.9)

After some routine algebra we finally obtain the expansion on the critical
line for large $t$ in the form [3]:

###### Theorem 1

$\\!\\!\\!.$ As $t\to+\infty$ on the critical line
$s=\mbox{${\textstyle\frac{1}{2}}$}+it$, we have the expansion

$Z(t)=\Re
e^{i\vartheta(t)}\biggl{\\{}\sum_{n=1}^{\infty}n^{-s}\mbox{erfc}\,(\mbox{${\textstyle\frac{1}{2}}$}\eta_{n}\sqrt{s};m)-\frac{1}{\sqrt{2i}}\biggl{(}\frac{\pi
e^{1/2}}{\omega}\biggr{)}^{\\!s}\biggl{(}T_{m}(\omega)+\omega^{-2m+1}{\hat{R}}_{m}\biggr{)}\biggr{\\}},$
(2.10)

where $\omega^{2}=\pi s/(2i)$ and

$T_{m}(\omega):=\sum_{r=0}^{m-1}(\pi i/4)^{r}B_{r}(\omega)$ (2.11)

with the coefficients $B_{r}(\omega)$ defined in (2.9).

In [3] the functions $S_{k}(\omega)$ appearing in the coefficients
$B_{r}(\omega)$ were treated as computable functions. In the next section we
express $B_{r}(\omega)$, and hence $T_{m}(\omega)$, as a finite series in
inverse powers of $\omega$ with coefficients containing certain trigonometric
functions.

3\. The coefficients in the expansion

In this section we examine the correction term $T_{m}(\omega)$ defined in
(2.11) and express it as a finite series in descending powers of $\omega$. To
achieve this we first observe that the functions $S_{k}(\omega)$ defined in
(2.8) satisfy

$S_{1}(\omega)=\csc\,\omega,\qquad
S_{k+1}(\omega)=\frac{1}{\omega}S_{k}(\omega)-\frac{1}{k}\,\frac{dS_{k}(\omega)}{d\omega}\quad(k\geq
1).$

From this we find with the help of Mathematica that

$S_{k+1}(\omega)=\frac{\csc\omega}{k!}\sum_{r=0}^{k}\frac{b_{r,k}}{\omega^{k}},$
(3.1)

where the coefficients $b_{r,k}$ for $0\leq r,k\leq 6$ are

$\displaystyle b_{00}$ $\displaystyle=$ $\displaystyle 1;$ $\displaystyle
b_{01}$ $\displaystyle=$ $\displaystyle\cot\,\omega,\ \ b_{11}=1;$
$\displaystyle b_{02}$ $\displaystyle=$
$\displaystyle\cot^{2}\omega+\csc^{2}\omega,\ \ b_{12}=3\cot\omega,\ \
b_{22}=3;$ $\displaystyle b_{03}$ $\displaystyle=$
$\displaystyle\cot^{3}\omega+5\cot\omega\csc^{2}\omega,\ \
b_{13}=6(\cot^{2}\omega+\csc^{2}\omega),\ \ b_{23}=15\cot\omega,\ \
b_{33}=15;$ $\displaystyle b_{04}$ $\displaystyle=$
$\displaystyle\cot^{4}\omega+18\cot^{2}\omega\csc^{2}\omega+5\csc^{4}\omega,\
\ b_{14}=10(\cot^{3}\omega+5\cot\omega\csc^{2}\omega),$ $\displaystyle b_{24}$
$\displaystyle=$ $\displaystyle 45(\cot^{2}\omega+\csc^{2}\omega),\ \
b_{34}=105\cot\omega,\ \ b_{44}=105;$ $\displaystyle b_{05}$ $\displaystyle=$
$\displaystyle\cot^{5}\omega+58\cot^{3}\omega\csc^{2}\omega+61\cot\omega\csc^{4}\omega,\
\ b_{15}=15(\cot^{4}\omega+18\cot^{2}\omega\csc^{2}\omega+5\csc^{4}\omega),$
$\displaystyle b_{25}$ $\displaystyle=$ $\displaystyle
105(\cot^{3}\omega+5\cot\omega\csc^{2}\omega),\ \
b_{35}=420(\cot^{2}\omega+\csc^{2}\omega),\ \ b_{45}=945\cot\omega,\ $
$\displaystyle\ b_{55}$ $\displaystyle=$ $\displaystyle 945;$ $\displaystyle
b_{06}$ $\displaystyle=$
$\displaystyle\cot^{6}\omega+179\cot^{4}\omega\csc^{2}\omega+479\cot^{2}\omega\csc^{4}\omega+61\csc^{6}\omega,$
$\displaystyle b_{16}$ $\displaystyle=$ $\displaystyle
21(\cot^{5}\omega+58\cot^{3}\omega\csc^{2}\omega+61\cot\omega\csc^{4}\omega),\
$ $\displaystyle b_{26}$ $\displaystyle=$ $\displaystyle
210(\cot^{4}\omega+18\cot^{2}\omega\csc^{2}\omega+5\csc^{4}\omega),\ \
b_{36}=1260(\cot^{3}\omega+5\cot\omega\csc^{2}\omega),\ $ $\displaystyle
b_{46}$ $\displaystyle=$ $\displaystyle 4725(\cot^{2}\omega+\csc^{2}\omega),\
\ b_{56}=10395\cot\omega,\ \ b_{66}=10395.$

Then, from (2.9) and (3.1), we find $B_{0}(\omega)=\csc\omega$ and

$B_{1}(\omega)=\alpha_{10}S_{3}(\omega)-\frac{2}{\omega}\alpha_{11}S_{2}(\omega)+\frac{4}{\omega^{2}}\alpha_{12}S_{1}(\omega)$
$=\csc\omega\biggl{(}\frac{\alpha_{10}}{2}\sum_{r=0}^{2}\frac{b_{r,2}}{\omega^{r}}-2\alpha_{11}\sum_{r=0}^{1}\frac{b_{r,1}}{\omega^{r+1}}+\frac{4}{\omega^{2}}\alpha_{12}b_{00}\biggr{)}=\csc\omega\sum_{k=0}^{2}\frac{d_{1,k}}{\omega^{k}},$

where

$d_{10}=\mbox{${\textstyle\frac{1}{2}}$}\alpha_{10}b_{02}=\mbox{${\textstyle\frac{1}{2}}$}(\cot^{2}\omega+\csc^{2}\omega),\qquad
d_{11}=\mbox{${\textstyle\frac{1}{2}}$}\alpha_{10}b_{12}-2\alpha_{11}b_{01}=-\mbox{${\textstyle\frac{1}{2}}$}\cot\omega,$
$d_{13}=\mbox{${\textstyle\frac{1}{2}}$}\alpha_{10}b_{22}-2\alpha_{11}b_{11}+4\alpha_{12}b_{00}=-\mbox{${\textstyle\frac{1}{6}}$}.$

In this manner we obtain after some effort

$B_{r}(\omega)=\csc\omega\sum_{k=0}^{2r}\frac{d_{r,k}}{\omega^{k}}\qquad(r\geq
1),$

where

$d_{r,k}=\sum_{p=0}^{k}\frac{(-2)^{k-p}\alpha_{r,k-m}}{(p+2r-k)!}\,b_{m,p^{*}}\qquad(p^{*}=2r\\!-\\!k\\!+\\!p\\!+\\!1).$

It then follows upon reversing the order of summation that

$T_{m}(\omega)=\csc\omega\sum_{r=0}^{m-1}(\pi
i/4)^{r}\sum_{k=0}^{2r}\frac{d_{r,k}}{\omega^{k}}=\csc\omega\sum_{k=0}^{2m-2}\frac{1}{\omega^{k}}\sum_{r=\lceil
k/2\rceil}^{m-1}(\pi i/4)^{r}d_{r,k}.$

Finally, we obtain the correction term in the following form:

###### Theorem 2

$\\!\\!\\!.$ The correction term $T_{m}(\omega)$ appearing in the expansion
(2.10) has the form

$T_{m}(\omega)=\csc\omega\sum_{k=0}^{2m-2}\frac{D_{k}(m)}{\omega^{k}},$ (3.2)

where the coefficients $D_{k}(m)\equiv D_{k}(m;\omega)$ are given by

$D_{k}(m)=\sum_{r=\lceil k/2\rceil}^{m-1}(\pi i/4)^{r}d_{r,k}.$ (3.3)

The coefficients $D_{k}(\omega)$ are easily seen to satisfy

$D_{k}(m+1)=D_{k}(m)+(\pi i/4)^{m}d_{m,k}\qquad(0\leq k\leq 2m-4)$

and involve the trigonometric functions $\csc\omega$ and $\cot\omega$. We now
present the explicit representation222The O$(\epsilon^{3})$ contribution in
$D_{2}(4)$ is zero. of the coefficients $D_{k}(m)$ for $m=1,2,3$ and 4, where
for brevity we have set $\epsilon=\pi i/4$.

$\displaystyle m=1:\qquad$ $\displaystyle D_{0}(1)=1;$ $\displaystyle
m=2:\qquad$ $\displaystyle
D_{0}(2)=D_{0}(1)+\mbox{${\textstyle\frac{1}{2}}$}\epsilon(\cot^{2}\omega+\csc^{2}\omega)$
$\displaystyle D_{1}(2)=-\mbox{${\textstyle\frac{1}{2}}$}\epsilon\cot\omega,\
\ \ D_{2}(2)=-\mbox{${\textstyle\frac{1}{6}}$}\epsilon;$ $\displaystyle
m=3:\qquad$ $\displaystyle
D_{0}(3)=D_{0}(2)+\mbox{${\textstyle\frac{1}{8}}$}\epsilon^{2}(\cot^{4}\omega+18\cot^{2}\omega\csc^{2}\omega+5\csc^{4}\omega)$
$\displaystyle
D_{1}(3)=D_{1}(2)-\mbox{${\textstyle\frac{5}{12}}$}\epsilon^{2}(\cot^{3}\omega+5\cot\omega\csc^{2}\omega)$
$\displaystyle
D_{2}(3)=D_{2}(2)-\mbox{${\textstyle\frac{5}{24}}$}\epsilon^{2}(\cot^{2}\omega+\csc^{2}\omega)$
$\displaystyle
D_{3}(3)=-\mbox{${\textstyle\frac{1}{24}}$}\epsilon^{2}\cot\omega,\ \ \
D_{4}(3)=\mbox{${\textstyle\frac{1}{72}}$}\epsilon^{2};$ $\displaystyle
m=4:\qquad$ $\displaystyle
D_{0}(4)=D_{0}(3)+\mbox{${\textstyle\frac{1}{48}}$}\epsilon^{3}(\cot^{6}\omega+179\cot^{4}\omega\csc^{2}\omega+479\cot^{2}\omega\csc^{4}\omega+61\csc^{6}\omega)$
$\displaystyle
D_{1}(4)=D_{1}(3)-\mbox{${\textstyle\frac{7}{48}}$}\epsilon^{3}(\cot^{5}\omega+58\cot^{3}\omega\csc^{2}\omega+61\cot\omega\csc^{4}\omega)$
$\displaystyle D_{2}(4)=D_{2}(3)$ $\displaystyle
D_{3}(4)=D_{3}(3)+\mbox{${\textstyle\frac{7}{36}}$}\epsilon^{3}(\cot^{3}\omega+5\cot\omega\csc^{2}\omega)$
$\displaystyle
D_{4}(4)=D_{4}(3)+\mbox{${\textstyle\frac{49}{144}}$}\epsilon^{3}(\cot^{2}\omega+\csc^{2}\omega)$
$\displaystyle
D_{5}(4)=\mbox{${\textstyle\frac{47}{144}}$}\epsilon^{3}\cot\omega,\ \ \
D_{6}(4)=\mbox{${\textstyle\frac{1003}{6480}}$}\epsilon^{3}.$

4\. Numerical results

In this section we describe computations using the expansion (2.10) with the
correction term $T_{m}(\omega)$ given in (3.2) and (3.3). The terms in the
main sum in (2.10) decay rapidly beyond $n\simeq N_{t}$. In the neighbourhood
of the transition point of $Q(\mbox{${\textstyle\frac{1}{2}}$}s,\pi n^{2}i)$
given by $\lambda_{n}=1$, we have

$\eta_{n}\simeq\lambda_{n}-1=\frac{\pi^{2}}{\omega^{2}}\biggl{(}n^{2}-\frac{\omega^{2}}{\pi^{2}}\biggr{)}\simeq\frac{2\pi
i}{s}\biggl{(}n^{2}-\frac{t}{2\pi}\biggr{)},$

so that the argument of the modified complementary error function in (2.10)
for large $t$ when $n\simeq N_{t}$ is

$\mbox{${\textstyle\frac{1}{2}}$}\eta_{n}\sqrt{s}\simeq\pi t^{-1/2}e^{\pi
i/4}(n^{2}-N_{t}^{2})\simeq 2\pi t^{-1/2}e^{\pi i/4}N_{t}(n-N_{t})$
$\simeq\sqrt{2\pi i}\,(n-N_{t}).$

From the asymptotic behaviour

$\mbox{erfc}\,(z;m)\sim\frac{e^{-z^{2}}}{\sqrt{\pi}}\,\frac{(-)^{m}(\mbox{${\textstyle\frac{1}{2}}$})_{m}}{z^{2m+1}}\qquad(|z|\to\infty,\
\ |\arg\,z|<\mbox{${\textstyle\frac{3}{4}}$}\pi),$

it is then found that when $n=N_{t}+K$, where integer $K\ll N_{t}$, the
magnitude of the terms in the main sum is given approximately by

$|n^{-s}\mbox{erfc}\,(\mbox{${\textstyle\frac{1}{2}}$}\eta_{n}\sqrt{s};m)|=(4\pi
K^{2})^{-m-1/2}O((t/2\pi)^{-1/4}).$

The decay of the terms is therefore controlled by $K^{-2m-1/2}$ together with
a scaling factor depending weakly on $t$ like $t^{-1/4}$. Thus if $K=10$, for
example, the magnitude of the terms in the main sum is approximately of order
$10^{-14}t^{-1/4}$ when $m=4$.

Table 1: Details of the computations when $t=2600$, $N_{t}=20$,
$Z(t)=-0.63210232$ for different values of $m$ and truncation index $K$ in the
main sum. $m=2$ Correction term=$+0.17012\ 33165\ 89694\ 33$
---
$K$ | Main Sum | $|Z_{\mbox{\footnotesize{approx}}}-Z(t)|$
10 | $-0.46197\ 90063\ 18404\ 31$ | $2.269\times 10^{-08}$
20 | $-0.46197\ 90041\ 44898\ 25$ | $9.502\times 10^{-11}$
30 | $-0.46197\ 90040\ 66916\ 27$ | $1.704\times 10^{-11}$
$m=3$ Correction term=$+0.08577\ 29470\ 58489\ 99$
$K$ | Main Sum | $|Z_{\mbox{\footnotesize{approx}}}-Z(t)|$
10 | $-0.54632\ 93736\ 04227\ 97$ | $2.315\times 10^{-11}$
20 | $-0.54632\ 93735\ 81347\ 73$ | $2.693\times 10^{-13}$
30 | $-0.54632\ 93735\ 81099\ 88$ | $2.146\times 10^{-14}$
$m=4$ Correction term=$-0.53258\ 74830\ 13204\ 32$
$K$ | Main Sum | $|Z_{\mbox{\footnotesize{approx}}}-Z(t)|$
10 | $-1.16468\ 98036\ 52714\ 66$ | $5.792\times 10^{-14}$
20 | $-1.16468\ 98036\ 52772\ 40$ | $1.874\times 10^{-16}$
30 | $-1.16468\ 98036\ 52772\ 58$ | $6.799\times 10^{-18}$
Table 2: Details of the computations when $t=2\times 10^{5}$, $N_{t}=178$,
$Z(t)=-3.51142011$ for different values of $m$ and truncation index $K$ in the
main sum. $m=2$ Correction term=$+0.04418\ 05095\ 82215\ 20$
---
$K$ | Main Sum | $|Z_{\mbox{\footnotesize{approx}}}-Z(t)|$
10 | $-3.46723\ 96042\ 56646\ 04$ | $5.980\times 10^{-10}$
20 | $-3.46723\ 96036\ 78939\ 00$ | $2.032\times 10^{-11}$
30 | $-3.46723\ 96036\ 61437\ 02$ | $2.815\times 10^{-12}$
$m=3$ Correction term=$+0.03215\ 07834\ 37290\ 07$
$K$ | Main Sum | $|Z_{\mbox{\footnotesize{approx}}}-Z(t)|$
10 | $-3.47926\ 93298\ 09067\ 88$ | $5.521\times 10^{-12}$
20 | $-3.47926\ 93298\ 03596\ 26$ | $4.958\times 10^{-14}$
30 | $-3.47926\ 93298\ 03549\ 81$ | $3.130\times 10^{-15}$
$m=4$ Correction term=$-0.03047\ 11685\ 29774\ 29$
$K$ | Main Sum | $|Z_{\mbox{\footnotesize{approx}}}-Z(t)|$
10 | $-3.54189\ 12817\ 70598\ 54$ | $1.251\times 10^{-14}$
20 | $-3.54189\ 12817\ 70611\ 01$ | $2.921\times 10^{-17}$
30 | $-3.54189\ 12817\ 70611\ 04$ | $8.377\times 10^{-19}$

In the computations we define the truncated main sum as

$\Re
e^{i\vartheta(t)}\sum_{k=1}^{N_{t}+K}n^{-s}\mbox{erfc}\,(\mbox{${\textstyle\frac{1}{2}}$}\eta_{n}\sqrt{s};m)$
$=\Re
e^{i\vartheta(t)}\biggl{\\{}\sum_{k=1}^{N_{t}}n^{-s}\biggl{(}2-\mbox{erfc}\,(\mbox{${\textstyle\frac{1}{2}}$}\eta_{n}\sqrt{s};m)\biggr{)}+\sum_{k=N_{t}+1}^{N_{t}+K}n^{-s}\mbox{erfc}\,(\mbox{${\textstyle\frac{1}{2}}$}\eta_{n}\sqrt{s};m)\biggr{\\}},$
(4.1)

where333We note that the factor 2 in the sum over $1\leq n\leq N_{t}$ in (4.1)
yields the standard Riemann-Siegel sum
$2\sum_{n=1}^{N_{t}}n^{-1/2}\cos(\vartheta(t)-t\log\,n)$. $N_{t}$ is the
Riemann-Siegel cut-off value defined in (1.2) and we have made use of the
result $\mbox{erfc}\,(-z;m)=2-\mbox{erfc}\,(z;m)$. The correction term is
given by

$\Re\frac{e^{i\vartheta(t)}}{\sqrt{2i}}\,\biggl{(}\frac{\pi
e^{1/2}}{\omega}\biggr{)}^{\\!s}T_{m}(\omega).$

The difference between these two contributions then yields the value
$Z_{\mbox{\footnotesize{approx}}}$. An example of the results is displayed in
Tables 1 and 2, where the value $Z(t)$ was obtained by Mathematica using the
command RiemannSiegelZ[t].. The values in the tables confirm that the accuracy
increases as both $m$ and the truncation index $K$ increase.

5\. Concluding remarks

We have revisited an expansion derived in [3] for $Z(t)$ on the critical line
$s=\mbox{${\textstyle\frac{1}{2}}$}+it$ as $t\to+\infty$ in which the main sum
is the original Dirichlet series smoothed by a complementary error function.
The correction term in this expansion has been expressed as a series in
inverse powers of $\omega$, where $\omega^{2}=\pi s/(2i)$, multiplied by
coefficients involving $\csc\omega$ and its derivatives. Numerical results are
presented to illustrate the accuracy achievable with this expansion.

A difficulty arises in the use of the expansion (2.10) when $(t/2\pi)^{1/2}$
is close to an integer i.e., at a discontinuity in $N_{t}$. This arises
because

$\frac{\omega}{\pi}=\biggl{(}\frac{t}{2\pi}\biggr{)}^{\\!1/2}\biggl{(}1-\frac{i}{2t}\biggr{)}^{\\!1/2},$

so that $\csc\omega$ and $\cot\omega$ become large, but never singular as
$\omega$ always has a small imaginary part (that decreases with increasing
$t$). This results in the term in the main sum and correction term
corresponding to $\eta_{n}\simeq 0$ becoming large. A means of overcoming this
problem, by removing and combining these terms, is given in [3, Section 5].

## References

  * [1] F.W.J. Olver, D.W. Lozier, R.F. Boisvert and C.W. Clark (eds.), NIST Handbook of Mathematical Functions, Cambridge University Press, Cambridge, 2010.
  * [2] R.B. Paris, New asymptotic formulas for the Riemann zeta function on the critical line, Special Functions, Proceedings of the International Workshop 21–25 June 1999, (eds. C. Dunkl, M. Ismail and R. Wong) 247–261, World Scientific, Singapore, 2000.
  * [3] R.B. Paris and S. Cang, An asymptotic representation for $\zeta(\mbox{${\textstyle\frac{1}{2}}$}+it)$, Methods Applic. Anal. 4 (1997) 449–470.
  * [4] N.M. Temme, The asymptotic expansions of the incomplete gamma functions, SIAM J. Math. Anal. 10 (1979) 757–766.