MARS: A neurosymbolic approach for interpretable drug discovery

## Abstract

Abstract Neurosymbolic (NeSy) artificial intelligence describes the combination of logic or rule-based techniques with neural networks. Compared to neural approaches, NeSy methods often possess enhanced interpretability, which is particularly promising for biomedical applications like drug discovery. However, since interpretability is broadly defined, there are no clear guidelines for assessing the biological plausibility of model interpretations. To assess interpretability in the context of drug discovery, we devise a novel prediction task, called drug mechanism-of-action (MoA) deconvolution, with an associated, tailored knowledge graph (KG), MoA-net . We then develop the MoA Retrieval System (MARS) , a NeSy approach for drug discovery which leverages logical rules with learned rule weights. Using this interpretable feature alongside domain knowledge, we find that MARS and other NeSy approaches on KGs are susceptible to reasoning shortcuts, in which the prediction of true labels is driven by “degree-bias” rather than the domain-based rules. Subsequently, we demonstrate ways to identify and mitigate this. Thereafter, MARS achieves performance on par with current state-of-the-art models while producing model interpretations aligned with known MoAs.

## 1 Introduction

Drug discovery (DD), the search for novel drugs or chemical compounds to treat ailments, often involves the screening of thousands of small compounds (Lin et al., ). Many computational approaches have been developed to accelerate and streamline this screening process (Gottlieb et al., ; Gan et al., ). Specifically, hundreds of such approaches operate upon knowledge graphs (KGs), in which nodes representing drugs, proteins, or medical conditions are connected by edges, representing the relationships between them (Chen et al., ). Typically, DD is formulated on a KG as a link prediction task between drugs and the corresponding medical conditions (indications) to be treated (Schultz et al., ; Rivas-Barragan et al., ).

Moreover, it is important to understand each drug’s mechanism-of-action (MoA), the molecular processes by which it achieves its medicinal effect. For instance, one drug, cortisone acetate, upregulates the activity of the glucocorticoid (GC) receptor protein, which, in turn, downregulates the cyclooxygenase (COX) protein. Since the COX protein is involved in creating inflammation, its inhibition reduces inflammation (Gonzalez-Cavazos et al., ). Thus, MoAs typically involve chains of physical, molecular interactions induced by a drug. Uncovering these interactions informs researchers as to how each drug works and de-risks potential side effects (Palve et al., ; Green et al., ).

However, revealing MoAs alongside computational DD, a task we call MoA deconvolution, requires model interpretability: transparency into the processes or patterns which led to certain predictions (Molnar, ). Unfortunately, most state-of-the-art techniques on KGs rely on “black-box” models (Wu et al., ). Recently, neurosymbolic (NeSy) approaches, which combine logical rules with neural networks (DeLong et al., ) have been positioned as a promising avenue for MoA deconvolution. Specifically, NeSy methods tend to possess enhanced interpretability and better integrate structured domain knowledge, of which there is plenty in biomedicine (Acharya et al., ; Rigden & Fernández, ). However, since interpretability is broadly defined, there are no clear guidelines for assessing the plausibility of model interpretations, especially for this novel task (Molnar, ).

Thus, we propose MoA deconvolution as a prediction task for evaluating interpretable methods on KGs. We generate a tailored KG for benchmarking this task, MoA-net. To perform MoA deconvolution, we introduce a NeSy DD approach called the MoA Retrieval System (MARS). To create MARS, we draw inspiration from previous NeSy methods (Liu et al., ; Drancé et al., ). However, unlike its predecessors, MARS leverages logical rules with learned rule weights. Following training, rule weights reflect relative usefulness to MARS’ reasoning processes.

Alongside biomedical domain knowledge, MARS’ enhanced interpretability reveals a reasoning shortcut in which predictions are based on unintended semantics (Marconato et al., ). Essentially, predictions upon MoA-net are driven by “degree-bias”, an artifact of node degree variance (Zietz et al., ), rather than the rules representing domain knowledge. Therefore, to address this second challenge, we consider the desiderata from Marconato et al. () for making NeSy systems shortcut-aware: (1) calibration, high accuracy on concepts unaffected by reasoning shortcuts, (2) performance, high accuracy despite reasoning shortcuts being present, and (3) cost effectiveness achieved through simple mitigation strategies. Using these desiderata as guidelines, we make MARS shortcut-aware for more insightful predictions involving DD and MoA deconvolution. Ultimately, our study underscores the importance of evaluating the capabilities of NeSy models within applied domains: by evaluating model interpretations against specific domain knowledge, we can more easily identify and mitigate shortcuts.

Contributions. In summary, our work:

- 1.
devises a novel prediction task, MoA deconvolution, to assess the interpretability of NeSy and other KG-based approaches;
- 2.
provides a novel biomedical KG, MoA-net, on which MoA deconvolution can be benchmarked;
- 3.
presents the MoA retrieval system, MARS, a NeSy approach with enhanced interpretability through dynamic rule-weight updates;
- 4.
through MARS’ enhanced interpretability, identifies a key reasoning shortcut affecting some NeSy models;
- 5.
demonstrates ways to test and mitigate this reasoning shortcut, making MARS shortcut-aware and better aligned with domain knowledge.

## 2 Related Work and Background

Drug Discovery and MoA Deconvolution. Amongst previous open-source approaches for DD, many focus on drug repositioning (DR), as opposed to finding novel compounds (Rivas-Barragan et al., ; Urbina et al., ; Zheng et al., ). DR aims to discover new uses for drugs that have already been approved for distribution (Urbina et al., ). This prediction task allows researchers to leverage more available knowledge on approved drugs, including recorded clinical outcomes and side effects (Ratajczak et al., ; Liu et al., ). However, such relations are revealed during or after clinical trials, and are therefore rare or unknown for novel compounds. Furthermore, with such relations present, DR models oftentimes rely on associative patterns (e.g., drugs causing the same side effects also likely share indications) rather than mechanistic ones (e.g., involving physical, molecular interactions). Therefore, such models offer limited insight into drug MoAs.

NeSy AI for KGs. Several NeSy approaches for KGs involve logical rules in which the rule bodies reflect path-like patterns in the KG (e.g., $\texttt{treats}(A,D)\!\leftarrow\!\texttt{targets}(A,B)\land\texttt{regulates}
(B,C)\land\texttt{attenuates}(C,D)$) (Liu et al., ; Dash & Goncalves, ; Sen et al., ). This is ideal for investigating MoAs, which typically occur as consecutive chains of molecular interactions (Crino, ). Some of these approaches utilize deep reinforcement learning (RL), in which a neural network contributes toward the optimization of a reward function (Acharya et al., ). In these specific cases, the reward function is partially determined by whether an agent’s trajectory through the KG corresponds to a rule body, thereby satisfying a rule. Therefore, the agent is guided by the path-like rules, and predictions are expected to align, to some extent, with such rules. In this study, however, we also identify a major risk: the agent will ignore rules if it can exploit other semantics for reward optimization. This results in *reasoning shortcuts*.

Reasoning Shortcuts. Several NeSy approaches are designed to abide by rules and domain knowledge (Drancé et al., ; Dash & Goncalves, ), which might portray such approaches as more trustworthy than neural, black box ones (Gaur & Sheth, ). Recent studies, however, have found that NeSy approaches may suffer from reasoning shortcuts, in which a model predicts the correct outcome via unintended semantics (Marconato et al., ; Li et al., ). While reasoning shortcuts are not exclusive to NeSy methods (Jiang & Bansal, ; Li et al., ), they may be more easily overlooked when such approaches are portrayed as trustworthy.

## 3 MARS: a NeSy approach for MoA deconvolution

Here, we build the MoA Retrieval System (MARS)(^1^11https://github.com/laurendelong21/MARS) to perform MoA deconvolution. MARS improves upon a method called Policy-guided walks with logical rules (PoLo) (Liu et al., ) by introducing dynamic, learned rule weights. This differs from previous approaches, where weights are static and pre-computed (e.g., mined or literature-derived) (Liu et al., ; Drancé et al., ). As discussed further, these learned weights make MARS shortcut-aware.

MARS takes two major inputs. The first is a KG. A KG uses nodes to represent entities and edges between them to represent relationships. A KG triple comprises two nodes connected by an edge of some specific type, or relation. For example, $\texttt{interacts}(Protein,Protein)$ states that two $Protein$ nodes are connected via the interacts relation. Node degree describes the number of edges connected to a node.

Specifically, the input KG must contain triples involving some relation of interest. To understand each MoA as the biological response to drug administration, we aim to investigate relations between drugs and biological processes (BPs), such as signal transduction or inflammation (Consortium, ) (i.e., $\texttt{induces}(Drug,BP)$). As per our knowledge, this is a novel application in the KG field. Further details on our KG are introduced within section 3.4. The second input includes metapaths of the KG with corresponding weights.

### 3.1 Metapath-based rules

Metapaths are abstract representations of instantiated paths in a graph (Sun et al., ; Himmelstein et al., ; Noori et al., ). For example, given the following path, $P$, in our KG:

$$ $\textit{Cortisone acetate}\xrightarrow{upregulates}\textit{GC receptor} \xrightarrow{interacts}\textit{COX protein}\xrightarrow{participates}\textit{Inflammation}$ $$

the corresponding metapath, $\tilde{P}$, would be:

$$ $\textit{Drug}\xrightarrow{upregulates}\textit{Protein}\xrightarrow{interacts} \textit{Protein}\xrightarrow{participates}\textit{Biological Process}$ $$

Within this study, metapaths can be understood as a sequence of relations within the KG structure, making them inherently interpretable. In MARS, metapaths are used as the bodies of logical rules, in which the rule head, the left side of the implication arrow ($\Leftarrow$), constitutes a single, binary predicate representing the relation of interest between the first and last node types, e.g.,:

$\texttt{induces}(\textit{Drug},\textit{BP})\Leftarrow\texttt{upregulates}(
\textit{Drug},\textit{Protein}_{A})\land\\
\hphantom{\texttt{induces}(\textit{Drug},\textit{BP})}\texttt{interacts}(
\textit{Protein}_{A},\textit{Protein}_{B})\land\texttt{participates}(\textit{
Protein}_{B},\textit{BP})$

For each metapath-based rule, $M_{i}$, in the set of feasible rules, $\mathcal{M}=\{M_{1},M_{2},...,M_{m}\}$, we denote the rule weight by $w(M_{i})\in\mathbb{R}$, where $0\leq w(M_{i})\leq 1$. Such a weight indicates the relative usefulness of the metapath-based rule to the prediction task. In Section 3.3, we discuss how we initialize and compute these weights.

### 3.2 Overview of MARS.

Using a deep RL process, MARS trains an agent to take walks of length $L$ through the KG to connect pairs of nodes having the pre-defined relation of interest (e.g., induces); such edges are masked from the agent during training. Each walk generates a path, $P$, which can be understood as a series of $L$ transitions: $P:=(e_{c}\xrightarrow{r_{1}}e_{2}\xrightarrow{r_{2}}...\xrightarrow{r_{L}}e_{L
+1})$. The agent may also remain at its current node. Ultimately, the goal of the agent is episodic: to find paths in which the starting node, $e_{c}$, and the terminal node, $e_{L+1}$, make up one of the true input pairs. By training the agent to do so, it can identify node pairs with some desired relationship, thus generalizing beyond the training set to predict novel pairs. Through a process akin to a Markov Decision Process (Bellman, ), the agent makes decisions about its next move based on information about its current position and the next possible actions. Additionally, however, the history of the agent’s previous actions are encoded with an LSTM (Hochreiter & Schmidhuber, ; Sherstinsky, ), whose parameters are trained to optimize the reward function, $R(S_{L+1})$ (Eq. 1), which is evaluated each time the agent completes $L$ transitions from some starting node, reaching a state $S_{L+1}$:

$$ $R(S_{L+1})=\mathbbm{1}_{\{e_{L+1}=e_{d}\}}+\mathbbm{1}_{\{e_{L+1}=e_{d}\}} \lambda\sum_{i=1}^{m}w(M_{i})\mathbbm{1}_{\{\tilde{P}=M_{i}\}}$ (1) $$

$$ $\mathbbm{1}_{\{A\}}=\begin{cases}1&\text{if }A=\text{true}\\ 0&\text{if }A=\text{false}\end{cases}$ (2) $$

This reward function, adapted from Liu et al. (), quantifies how successful $P$ is according to two summands, where the hyperparameter $\lambda$ influences the balance between them. The first summand indicates whether the terminal node in the path, $e_{L+1}$, is one of the desired target (BP) nodes ($e_{d}$) that forms a true pair with the starting (drug) node, $e_{c}$. The second summand, contingent upon the first, indicates whether the corresponding metapath, $\tilde{P}$, matches any metapath-based rule, $M_{i}$. Therefore, the agent is not only encouraged to find connections between true pairs of nodes, but it is also guided toward paths which resemble known MoAs. Implementation details and hyperparameter selection are described in Appendices A.2 and A.4, respectively.

Of note, the second summand is proportional to some metapath-based rule weight, and MARS updates these weights during training. We accomplish this through a novel algorithm we call two-hop joint probability, or $P_{2H}$.

### 3.3 Weight updates via 2-hop joint probability

After MARS is executed, its learned rule weights reflect each rule’s relative usefulness in the prediction task. Additionally, during training, weight updates drive the agent toward more informative paths and bypass the assumptions that pre-assigned rule weights are correct and informative. This eliminates the need for pre-computed or literature-derived rule weights; therefore, we initialize all weights uniformly as 0.5, a medium level of importance.

The naive way to implement weight updates ($\text{MARS}_{\text{naive}}$) is to increase weights according to the frequency at which each metapath-based rule is satisfied. However, this would not account for partial metapath matches. As an alternative, we propose updates via two-hop joint probability ($P_{2H}$). This metric approximates the usefulness of metapath-based rules based on the observed frequencies of the two-hop fragments that they comprise. Since the metapaths constituting rule bodies typically contain three or four consecutive relations, each two-hop fragment is extracted as in the following example:

$M_{example}:=\texttt{induces}(A,E)\Leftarrow\texttt{upregulates}(A,B)\land
\texttt{interacts}(B,C)\land\\
\hphantom{M_{example}:=\texttt{induces}(A,E)\Leftarrow\texttt{upregul}}\texttt
{interacts}(C,D)\land\texttt{participates}(D,E)$

where two-hop fragments are pairs of binary predicates which share a variable:

- •
$Fragment1:\texttt{upregulates}(X,Y)\land\texttt{interacts}(Y,Z)$
- •
$Fragment2:\texttt{interacts}(X,Y)\land\texttt{interacts}(Y,Z)$
- •
$Fragment3:\texttt{interacts}(X,Y)\land\texttt{participates}(Y,Z)$

Here, the probability of each metapath-based rule is computed as the joint probability of its fragments. For example, the $P_{2H}$ metric for the above metapath-based rule would be computed as $P_{2H}(M_{example})=p(Fragment1)\times p(Fragment2)\times p(Fragment3)$. The probability of each fragment is determined by the ratio of its observed and expected frequencies (see Appendix A.3). We note that, while the fragments do not necessarily represent independent events, the way in which $P_{2H}$ is computed imposes this assumption. To compute the joint probability for dependent events, we require the conditional probabilities of the dependent events (Russell & Norvig, ). This can, in theory, be computed, but it becomes more complex with metapath length. For simplicity’s sake, we assume independence, and the $P_{2H}$ metric serves as an approximation for the empirical probabilities of metapath-based rules. Further formalization and pseudocode are in Appendix A.3. Ultimately, MARS with $P_{2H}$ updates ($\text{MARS}_{P_{2H}}$) uses all information from successful trajectories.

### 3.4 Datasets: MoA-net and its variants

We design our KG, MoA-net(^2^22https://github.com/laurendelong21/MoA-Net), specifically for MoA prediction. MoA-net consists of drugs, proteins, and BPs (Fig. 1-A). We assemble it using the causal relations between drugs and proteins from Custom KG (Rivas-Barragan et al., ) and OpenBioLink KG (Breit et al., ). The BP nodes come from molecular function annotations in UniProt (Consortium, ).

To predict drug-BP triples, which are unique to MoA-net, we make use of publicly available functional and biochemical assays in ChEMBL (v33), an open access database of bioactive compounds (Gaulton et al., ). Of the 1,622 drug-BP triples obtained, 48 also had known MoAs in DrugMechDB (Gonzalez-Cavazos et al., ), a manually-curated compendium of MoAs. Between the three node types, we define five unique edge types, or relations, shown in Fig. 1-B. We also include all inverse relations, running in the opposite direction of causality.

Using the hetnetpy package (Himmelstein et al., ), we extract all metapaths (see Section 3.1) from MoA-net which we considered to be valid MoAs: those comprising directed, mechanistic paths between drug and BP nodes (see Appendix A.1). Although other metapaths are possible, we exclude them from our set of metapath-based rules. Based on MoAs found in DrugMechDB, we limit metapaths to a maximum length of four relations (or hops).

Finally, we create variants of MoA-net. To investigate reasoning shortcuts, we use the Zietz et al. () implementation of XSwap (Hanhijärvi et al., ), which swaps edges in a KG without affecting the distribution of node degrees. We call the resultant KG MoA-net-permuted. Additionally, we implement an automatic trimming step, which reduces edges of each class to below a user-specified threshold by iteratively removing those between the highest-degree nodes. By setting the threshold to 10,000 (thereby reducing protein-protein interactions to $\sim 50\%$ of edges), our approach can work on a subgraph of the MoA-net, which we refer to as MoA-net-10k.

Figure: Figure 1: A) MoA-net schema. B) Distribution of node and relationship types of MoA-net. Note that the $y$-axis has a logarithmic scale.
Refer to caption: extracted/5924853/figures/Neurosymbolic_schema.png

### 3.5 Evaluation

We split the drug-BP triples within MoA-net into training (60%), validation (20%), and test (20%) sets. We evaluate the models using Hits@$k$, where $k\in\{1,3,10\}$ and mean reciprocal rank (MRR), optimizing for the latter. Hits@$k$ reports the proportion of times the correct results are in the top $k$ ranked entries, while MRR reports how highly ranked the first correct item is amongst ranked results (Chen et al., ). In addition to these standard metrics, we report the pruned metrics: these are computed on a subset of the predictions that utilized one of the pre-defined metapath-based rules (see Appendix A.1), excluding all predictions which did not satisfy a rule. Notably, pruned metrics help us assess the calibration desideratum, as introduced in Section 1, since rule-based predictions follow the expected model semantics.

We conduct an extensive benchmark of our method against nine different baseline KG embedding (KGE) models, two state-of-the-art NeSy methods, and one network measure. Based on previous biomedical KG benchmarks (Rivas-Barragan et al., ), the KGE models include RESCAL (Nickel et al., ), TransE (Bordes et al., ), DistMult (Yang et al., ), ER-MLP (Dong et al., ), ComplEx (Trouillon et al., ), HolE (Nickel et al., ), ConvE (Dettmers et al., ), RotatE (Sun et al., ) and MuRE (Balazevic et al., ). We also compare against PoLo (Liu et al., ) and its predecessor MINERVA (Das et al., ), which optimizes only the first summand of the reward function (i.e., $\lambda=0$ in Eq. 1). Additionally, we test prioritization of drug-BP triples based on degree-weighted path count (DWPC) using 0.4 as damping exponent (Himmelstein & Baranzini, ). We train and evaluate these models on the same data splits as MARS on MoA-net-10k.

Finally, $\text{MARS}_{P_{2H}}$ has two key interpretable features: firstly, all successful trajectories are recorded, serving as potential MoA predictions. This allows us to compare the predicted MoAs of 48 drug-BP pairs against their known MoAs (see Section 3.4). Secondly, learned rule weights serve as a proxy for the importance of each metapath-based rule. This helps determine whether agent trajectories are biased toward certain types of paths. Alongside the pruned metrics, these features help evaluate MARS’ alignment with domain knowledge.

## 4 Results

### 4.1 Associative patterns improve accuracy but offer limited practical use

In an initial set of experiments on MoA-net, we observed that pruned metrics were consistently lower than standard ones (Fig. 2-A), indicating that the metapath-based rules were not being utilized in most predictions. This can happen because rule-based rewards are contingent upon a true pair being found (Eq. 1). Additionally, amongst recorded trajectories, most did not follow our metapath-based rules; instead, most trajectories used the following metapath, involving inverse edges:

$\texttt{induces}(\textit{Drug}_{1},\textit{BP}_{2})\Leftarrow\texttt{induces}(
\textit{Drug}_{1},\textit{BP}_{1})\land\\
\hphantom{\texttt{induces}(\textit{Drug}_{1},\textit{BP})\Leftarrow}\texttt{
induces}(\textit{Drug}_{2},\textit{BP}_{1})\land\texttt{induces}(\textit{Drug}
_{2},\textit{BP}_{2})$

This associative pattern indicates that two drugs inducing a common BP also likely induce another BP. This type of pattern is also present in a previous NeSy RL study, Liu et al. (), in which the most used metapath was the following:

$\texttt{treats}(\textit{Drug}_{1},\textit{Disease})\Leftarrow\texttt{causes}(
\textit{Drug}_{1},\textit{Side Effect})\land\\
\hphantom{\texttt{treats}(\textit{Compo}_{1},\textit{Dis})\Leftarrow}\texttt{
causes}(\textit{Drug}_{2},\textit{Side Effect})\land\texttt{treats}(\textit{
Drug}_{2},\textit{Disease})$

However, when we reproduced the results from Liu et al. (), we achieved the same reported metrics even in the absence of the above, associative rule (See Appendix A.5). This suggests that, although the associative rule may serve as a plausible model explanation, it does not necessarily guide model training.

Previously, associative patterns have been regarded as informative in other KG prediction tasks such as DR and target prediction (Fernández-Torras et al., ; Ratajczak et al., ; Himmelstein et al., ; Renaux et al., ; Li et al., ; Muslu et al., ). However, associative patterns do not constitute biologically valid MoAs. As stated in Section 1, MoAs, like those in DrugMechDB, involve physical, molecular interactions. Therefore, an MoA path would not traverse an additional drug node. In other words, while such associations may be a valid proxy for other biomedical tasks, MoA deconvolution for novel compounds is inherently limited to information around its protein targets and related BPs.

Figure: Figure 2: Performance metrics for $\text{MARS}_{P_{2H}}$ compared to PoLo and $\text{MARS}_{\text{naive}}$ upon several variants of MoA-net. Each bar is the average and standard deviation across five independent training and testing iterations. From left to right: Little change between initial metrics upon MoA-net (A) in comparison to the standard MoA-net-permuted metrics (B) provides evidence that predictions are influenced by degree bias, resulting in a reasoning shortcut. Thereafter, inverse edges were removed to prohibit the reasoning shortcut, hindering performance (C). Performance was restored upon MoA-net-10k with the KG trimming step (D), with $\text{MARS}_{P_{2H}}$ showing the best standard and pruned metrics. Finally, $\text{MARS}_{P_{2H}}$ maintains high pruned metrics even when inverse edges (and reasoning shortcuts) are re-introduced (E).
Refer to caption: extracted/5924853/figures/hitsat1_barplot.png

### 4.2 P 2 ⁢ H subscript 𝑃 2 𝐻 P_{2H} italic_P start_POSTSUBSCRIPT 2 italic_H end_POSTSUBSCRIPT updates reveal reasoning shortcuts via degree bias

In addition to analyzing the agent trajectories, we used $P_{2H}$ to assess how informative each of the metapath-based rules is in making predictions. In particular, $\text{MARS}_{P_{2H}}$ weights showed that paths involving consecutive protein-protein interactions (PPIs) (i.e., $\texttt{interacts}(Protein,Protein)$), were consistently less important (Fig. 3). This indicated that the agent avoided exploring consecutive PPIs.

Figure: Figure 3: Metapath-based rule weights from $\text{MARS}_{P_{2H}}$ on MoA-net (Fig. 2-A). Each bar is the average and standard error across five independent training and testing iterations. Paths involving consecutive PPIs ($\texttt{interacts}(Protein,Protein)$), the most common relation type, have consistently lower weights.
Refer to caption: extracted/5924853/figures/confidences.png

Previous research on KGs has shown that node degree distribution, the number of adjacent edges for each KG node, can significantly bias predictions (Tang et al., ; Ju et al., ). Specifically, inspection bias, a type of degree bias, occurs when the KG is not uniformly inspected or sampled (Zietz et al., ). Since PPIs are the most common relation type in MoA-net (90% of edges) (Fig. 1-B), protein nodes have a higher degree distribution than other node types. We hypothesized, therefore, that the agent circumvents denser parts of the KG, creating an inspection bias. Although rule-based predictions merit a larger reward, the MARS agent exploits associative patterns for a more reliable reward. To confirm the existence of degree bias, we tested our approach upon MoA-net-permuted. As explained in Section 3.4, MoA-net-permuted is a variant of MoA-net in which edges are swapped while preserving node degree distribution. This tests the extent to which node degree drives predictions. Indeed, the lack of change amongst standard performance metrics suggested that node degree was largely responsible for predictions (Fig. 2-B). Put simply, the agent gets lost when exploring the PPIs, so it avoids them.

### 4.3 Identifying and mitigating degree bias improves performance

To temporarily prohibit the models from using associative patterns as in Section 4.1, we removed inverse edges from MoA-net and corresponding metapath-based rules. Consequently, the performance metrics were poor, (e.g., MRR consistently $<0.1$ (Fig. 2-C)). This confirmed that the models relied on associative patterns for predictions.

Next, we wanted to confirm that the agent was getting lost within the PPIs. As explained in Section 3.4, MoA-net-10k is a variant of MoA-net with fewer PPIs. We tested $\text{MARS}_{P_{2H}}$, $\text{MARS}_{\text{naive}}$, and PoLo with the same parameters upon on MoA-net-10k (Fig. 2-D). As before, we excluded inverse edges. Since we set trajectory length $L=4$, our approach automatically removed drug-BP triples from the validation/test sets that were no longer connected via directed paths of length $\leq 4$, resulting in 100 and 90 triples, respectively. Metrics were markedly improved for PoLo, $\text{MARS}_{\text{naive}}$, and particularly for $\text{MARS}_{P_{2H}}$, in comparison to the full MoA-net without inverse edges (Fig. 2-C). To ensure this improvement was not simply the result of a reduced test set, we also tested the approaches upon MoA-net with 100 sampled test triples, which showed no change (see Appendix A.6).

While removing inverse edges improved metrics, a shortcut-aware system should achieve high performance even with the shortcut present (Marconato et al., ). We addressed this next.

### 4.4 MARS P 2 ⁢ H subscript MARS subscript 𝑃 2 𝐻 \text{MARS}_{P_{2H}} MARS start_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT 2 italic_H end_POSTSUBSCRIPT end_POSTSUBSCRIPT retains performance amongst rule-based predictions

We re-introduced inverse edges to MoA-net-10k, thereby restoring the ability to use reasoning shortcuts. Thereafter, we tested each of $\text{MARS}_{P_{2H}}$, $\text{MARS}_{\text{naive}}$, and PoLo again (Fig. 2-E). While each approach was optimized for standard MRR, pruned metrics indicated how well positive predictions aligned with rules. In Fig. 2-E, we show that both MARS variants and PoLo achieved standard metrics on par with or better than MoA-net-10k without inverse edges (Fig. 2-D). However, $\text{MARS}_{P_{2H}}$ also achieved pruned metrics comparable to its standard metrics, showing improved calibration relative to PoLo and $\text{MARS}_{\text{naive}}$. Finally, as in section 4.2, we used XSwap on MoA-net-10k to assess the susceptibility of $\text{MARS}_{P_{2H}}$ to degree bias. Unlike in Section 4.2, we found no evidence for degree bias (see Appendix A.7).

### 4.5 External validation of MARS P 2 ⁢ H subscript MARS subscript 𝑃 2 𝐻 \text{MARS}_{P_{2H}} MARS start_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT 2 italic_H end_POSTSUBSCRIPT end_POSTSUBSCRIPT on MoA-net-10k

In comparison to baseline methods, $\text{MARS}_{P_{2H}}$’s metrics outperformed all but MINERVA’s, with which they were comparable (Table 1). However, since MINERVA does not, by design, utilize rules for guidance, it suffers the same reasoning shortcuts as PoLo and $\text{MARS}_{naive}$. In contrast to MINERVA, DWPC suffers the opposite limitation: predictions are based only on metapath-based rules. $\text{MARS}_{P_{2H}}$’s pruned metrics, which are directly comparable, also outperform DWPC. Finally, as mentioned in Section 3.4, several drug-BP pairs corresponding to known MoAs in DrugMechDB were included in the MoA-net test set. Of these, 33 pairs remained within MoA-net-10k’s test set, and $\text{MARS}_{P_{2H}}$ recovered the correct MoA for all of them. Thus, this comprehensive benchmark highlights MARS’ ability to achieve near state-of-the-art performance by effectively balancing domain-specific knowledge with the capacity to generalize beyond it.

**Table 1: Performance evaluations of MARS upon MoA-net-10k against baseline models. Metrics are presented as the (average, standard deviation) across five independent training/testing iterations. The best of each standard (top) and pruned (bottom) metric are in bold.**
| Model | metric type | Hits@1 | Hits@3 | Hits@10 | MRR |
| --- | --- | --- | --- | --- | --- |
| ComplEx | standard | 0.137 | 0.291 | 0.527 | 0.263 |
| ConvE | standard | 0.055 | 0.132 | 0.247 | 0.116 |
| DistMult | standard | 0.011 | 0.027 | 0.071 | 0.031 |
| ER-MLP | standard | 0.011 | 0.093 | 0.242 | 0.086 |
| HolE | standard | 0.038 | 0.077 | 0.203 | 0.096 |
| MuRE | standard | 0.071 | 0.181 | 0.253 | 0.177 |
| RESCAL | standard | 0.055 | 0.148 | 0.286 | 0.135 |
| RotatE | standard | 0.176 | 0.242 | 0.368 | 0.239 |
| TransE | standard | 0.154 | 0.297 | 0.538 | 0.274 |
| MINERVA | standard | (0.342, 0.016) | (0.516, 0.042) | (0.66, 0.066) | (0.45, 0.026) |
| PoLo | standard | (0.272, 0.041) | (0.462, 0.054) | (0.606, 0.061) | (0.387, 0.044) |
| $\text{MARS}_{\text{naive}}$ | standard | (0.33, 0.031) | (0.482, 0.066) | (0.664, 0.036) | (0.433, 0.036) |
| $\text{MARS}_{P_{2H}}$ | standard | (0.23, 0.007) | (0.492, 0.027) | (0.684, 0.03) | (0.395, 0.016) |
| Metapaths with DWPC | pruned | 0.370 | 0.560 | 0.780 | 0.508 |
| PoLo | pruned | (0.17, 0.049) | (0.228, 0.061) | (0.228, 0.061) | (0.198, 0.052) |
| $\text{MARS}_{\text{naive}}$ | pruned | (0.22, 0.049) | (0.238, 0.048) | (0.238, 0.048) | (0.229, 0.048) |
| $\text{MARS}_{P_{2H}}$ | pruned | (0.394, 0.026) | (0.644, 0.034) | (0.788, 0.018) | (0.535, 0.02) |

## 5 Discussion

NeSy approaches are sometimes portrayed as more trustworthy than their black-box counterparts, partially due to increased interpretability (Gaur & Sheth, ; DeLong et al., ). Within this study, we enhanced the interpretability of a NeSy RL approach to produce our model, $\text{MARS}_{P_{2H}}$. Specifically, through our novel algorithm, two-hop joint probabilities ($P_{2H}$), MARS learned weights corresponding to path-like rules, serving as a proxy for each rule’s importance. However, these insights revealed a new issue: NeSy RL approaches on KGs are susceptible to reasoning shortcuts. Specifically, in our study, predictions were driven by node degree bias. Ultimately, MARS’ enhanced interpretability called the trustworthiness of such approaches to question.

To address this, we considered Marconato et al. ()’s desiderata for a shortcut-aware NeSy system. Specifically, on MoA-net-10k, $\text{MARS}_{P_{2H}}$ showed both competitive performance as well as calibration in comparison to other models. Notably, however, measuring calibration is challenging in this domain. While rule-based predictions, measured through pruned metrics, follow the expected semantics for MoA deconvolution, we can not determine whether every other prediction follows unintended semantics. For example, in the classic MNIST addition task, popularly used to assess NeSy methods Manhaeve et al. (), a model is trained to determine the sum of two handwritten digits. In this toy example, the misclassification of a handwritten ‘2’ as ‘3’ and vice versa would still amount to the same sum. Thus, reasoning shortcuts can be objectively identified. On the contrary, while we provide evidence that predictions using associative patterns are largely affected by node degree bias, we can not determine whether such patterns always reflect a reasoning shortcut.

Finally, regarding cost effectiveness, $\text{MARS}_{P_{2H}}$ can be applied to any KG, serving as a generalizeable mitigation strategy. However, we also note that this was achieved upon MoA-net-10k, a trimmed version of MoA-net. While we automated this trimming step, such a strategy does not make use of all available information. To scale $\text{MARS}_{P_{2H}}$ to denser KGs and maintain its shortcut-aware status, several future directions could be explored. For instance, one could merge similar, high-degree nodes or rely upon domain knowledge, like the identification of promiscuous proteins (Copley, ), to make more informed choices about edge trimming or masking.

Our study highlights a key concern in which the behavior of some NeSy RL approaches could be attributed to node degree bias, rather than meaningful, domain-specific concepts. The interpretability of our approach, $\text{MARS}_{P_{2H}}$, allowed insight into this reasoning shortcut. Therefore, we question whether such shortcuts are identifiable amongst black-box approaches. Additionally, by testing a NeSy approach upon a novel applied task, MoA deconvolution, we could flag down patterns, like associative ones, which were plausible yet arguably less meaningful to biomedical researchers. Therefore, our study emphasizes the importance of testing interpretable models, like NeSy ones, in an applied domain. Finally, while our study honors the desiderata for shortcut-aware NeSy systems, we also examined the extent to which they were applicable to a biomedical domain.

## 6 Conclusions

We propose a novel prediction task for NeSy approaches on biomedical KGs: mechanism-of-action (MoA) deconvolution. In contrast to previous DD approaches, MoA deconvolution utilizes model interpretability to uncover the molecular mechanisms behind medicinal drugs. We also constructed a publicly available KG, MoA-net, for evaluating this task. To predict drug MoAs alongside indications, we designed the MoA Retrieval System (MARS). Relative to previous NeSy approaches, MARS has enhanced interpretability as it dynamically learns weights corresponding to logical rules. We showed that, with respect to the three desiderata for reasoning-aware NeSy systems, MARS has improved calibration and cost effectiveness compared to its predecessors, thereby enabling the identification and mitigation of a reasoning shortcut based on node degree bias.

#### Author Contributions

LND designed the methodology, wrote the software, and conducted the investigation. YG and LND did data curation and visualization. DDF and LND conceptualized the project. DDF, JDF, and PG supervised the project. All authors took part in writing the manuscript.

#### Acknowledgments

LND is funded by the University of Edinburgh Informatics Graduate School through the Global Informatics Scholarship. LND partially conducted this work during an internship at Enveda Biosciences Inc. JDF and PG are funded by the National Institute for Health Research (NIHR) Artificial Intelligence and Multimorbidity: Clustering in Individuals, Space and Clinical Context (AIM-CISC) grant NIHR202639. The views expressed are those of the author(s) and not necessarily those of the NIHR or the Department of Health and Social Care.

We thank these institutions for their support. We also thank Emile van Krieken and Guillermo Romero Moreno for their thoughtful feedback.

## Appendix A Appendix

### A.1 Selected Metapaths

**Table A1: Metapaths representing feasible MoAs. Drugs are represented with a $D$, proteins with a $P$, and biological processes with $BP$.**
| $\texttt{downregulates}(D,P)\rightarrow\texttt{participates}(P,BP)$ |
| --- |
| $\texttt{upregulates}(D,P)\rightarrow\texttt{participates}(P,BP)$ |
| $\texttt{downregulates}(D,P)\rightarrow\texttt{interacts}(P,P)\rightarrow<br>\texttt{participates}(P,BP)$ |
| $\texttt{upregulates}(D,P)\rightarrow\texttt{interacts}(P,P)\rightarrow\texttt{<br>participates}(P,BP)$ |
| $\texttt{downregulates}(D,P)\rightarrow\texttt{interacts}(P,P)\rightarrow<br>\texttt{interacts}(P,P)\rightarrow\texttt{participates}(P,BP)$ |
| $\texttt{upregulates}(D,P)\rightarrow\texttt{interacts}(P,P)\rightarrow\texttt{<br>interacts}(P,P)\rightarrow\texttt{participates}(P,BP)$ |

### A.2 Implementation

We implemented MARS using TensorFlow (version 2.10). The method is packaged in Python(^3^33https://github.com/laurendelong21/MARS). We used a grid search hyperparameter optimization (Feurer & Hutter, ); further details are within Appendix A.4.

The baseline KGE models have been trained using the PyKEEN framework (v1.10.1) (Ali et al., ). KGEMs were trained using PyKEEN’s hyperparameter optimization pipeline over 30 trials using as initial parameters the best configurations from (Rivas-Barragan et al., ). The evaluation in the hyperparameter optimization was conducted using Hits@10 for all the models on a link prediction task for the previously-described splits. Network algorithms were implemented in NetworkX (v3.1) (Hagberg et al., ) and metapaths were calculated using the hetnetpy Python package (Himmelstein et al., ; Himmelstein & Baranzini, ). Lastly, source code and data are openly available(^4^44https://github.com/laurendelong21/MoA-Net).

### A.3 P 2 ⁢ H subscript 𝑃 2 𝐻 P_{2H} italic_P start_POSTSUBSCRIPT 2 italic_H end_POSTSUBSCRIPT formalization

Here, we describe the weight updates in more detail. Both naive and $P_{2H}$ updates are normalized according to expected values. For naive updates, we record observed frequency, $O$, at which each metapath-based rule is satisfied. This is normalized by the batch-specific expected frequency, $E$, in which every rule has a uniform probability across the total number rule groundings in that batch. For instance, if there are 10 metapath-based rules, and the agent satisfied at least one rule 100 times in a given batch, $E$ would compute to 0.10, and each observed frequency $O_{M}$, would be normalized by $E$. If the agent found zero groundings, no weight updates would be made in that batch. Ultimately, this produces a metric, $\mu$ (Eq. 3) in which any value over one indicates usefulness (the agent used that rule more than others), and any value less than one indicates otherwise.

$$ $\mu_{M_{i}}=O_{M_{i}}/E_{M_{i}}$ (3) $$

Notably, the value of $\mu$ was bound to avoid division by zero and extreme values. To adjust weight updates relative to the batch size and rollouts, we defined the minimum and maximum bounds on $\mu$ as Eqs. 4 and 5, where $\rho$ is the total number of metapath-based rules:

$$ $\mu_{min}=\frac{\rho}{\text{batch size}\times\text{rollouts}}$ (4) $$

$$ $\mu_{max}=\rho\times\text{batch size}\times\text{rollouts}$ (5) $$

Using Eq. 6, denoted by $\Phi$, $\mu$ (Eq. 3) is used to update the weight of a rule, $w(M_{i})$. Eq. 6 is regularized by the hyperparameter $\alpha\in\mathbb{R}$, where $0\leq\alpha\leq 1$ to control how drastic the weight update is. An $\alpha$ closer to zero indicates subtle updates; this might be suitable if the user wants weights to act as a gauge on relative importance between rules. An $\alpha$ closer to one indicates drastic updates and could be suitable if the user has many rules and wants those which are less useful to be reduced to zero. If $\alpha=0$, no weight updates will be made. Therefore, $\alpha$ can be selected based on user needs or via hyperparameter optimization.

$$ $\Phi(\mu,w(M_{i}))=w(M_{i})\times 2\alpha(\frac{\mu-1}{\mu+1})$ (6) $$

The second, more complex method by which weight updates were done is based on a term we coined, two-hop joint probability, $P_{2H}$. Pseudocode for $P_{2H}$ can be found in Algorithm 1 below. This metric approximates the usefulness of metapath-based rules based on the observed frequencies of the two-hop fragments which they comprise. Here, the main differences are that (1) Eq. 3 is computed based on the observed and expected probabilities of two-hop fragments, rather than whole metapaths, and (2) $\rho$ within Equations 4 and 5 is the number of unique two-hop fragments possible.

Figure: Algorithm 1 $P_{2H}$ weight updates

### A.4 Hyperparameter selection

Here, we describe hyperparameter selection. Table A2 describes the hyperparameter search space for optimization, and Table A3 describes the hyperparameters which were fixed for every model.

**Table A2: Hyperparameter search space for grid search optimization (Feurer & Hutter, )**
| Hyperparameter | Description | Search space |
| --- | --- | --- |
| $\lambda$ (Lambda) | ratio at which the second summand is applied relative to the first summand in the reward function | {5, 8, 10} |
| $\alpha$ (alpha) | how dramatically weight updates should be made (if applicable) | {0.001, 0.01, 0.1} |
| learning rate | learning rate of the optimizer | {0.0001, 0.001, 0.01} |
| hidden size | size of hidden layers | {64 128 256} |
| batch size | size of sampled mini-batch for training | {128, 256} |
| rollouts | number of times each query (source-terminal node pair) is made or attempted during training | {50, 100} |
| $\gamma_{\text{baseline}}$ (gamma baseline) | discount factor for the baseline | {0.05, 0.5} |
| $\beta$ (beta) | entropy regularization factor | {0.025 0.05} |

**Table A3: Fixed hyperparameter settings**
| Hyperparameter | Description | Value |
| --- | --- | --- |
| embedding size | size of the relation and entity embeddings | 256 |
| LSTM layers | number of LTSM layers | 2 |
| test rollouts | number of times each query (source-terminal node pair) is made or attempted during testing | 50 |
| max branching | maximum number of outgoing edges per node shown to the agent in an episode | 150 |
| $\gamma$ (gamma) | discount factor | 1 |
| positive reward | reward for finding a true pair | 1 |
| negative reward | penalty for failing to find a true pair | 0 |

### A.5 PoLo metrics without associative rules

We ran PoLo on the Hetionet KG (Himmelstein et al., ) using the same parameters and data splits as reported by Liu et al. (). In contrast to Liu et al. (), we input all directed metapaths of length $L\leq 4$ as rule bodies (as in Appendix A.1). These metapaths served as the metapath-based rules for PoLo. Notably, these metapaths excluded the associative metapath mentioned in section 4.1:

$\textit{treats}(\textit{Drug}_{1},\textit{Disease})\Leftarrow\textit{causes}(
\textit{Drug}_{1},\textit{Side Effect})\land\\
\hphantom{\textit{treats}(\textit{Drug}_{1},\textit{Disease})\Leftarrow}
\textit{causes}(\textit{Drug}_{2},\textit{Side Effect})\land\textit{treats}(
\textit{Drug}_{2},\textit{Disease})$

Despite the most-used metapath-based rule being absent, PoLo achieved the same standard metrics as previously reported (Table A4).

**Table A4: Performance evaluations of PoLo upon Hetionet as reported in Liu et al. () (average across five independent training/testing iterations) and PoLo upon Hetionet without associative rules ((average, standard deviation) across four independent training/testing iterations.)**
| rule types | Hits@1 | Hits@3 | Hits@10 | MRR |
| --- | --- | --- | --- | --- |
| associative ((Liu et al., )) | 0.314 | 0.428 | 0.609 | 0.402 |
| mechanistic (this study) | (0.328, 0.046) | (0.465, 0.037) | (0.656, 0.044) | (0.431, 0.035) |

### A.6 Ablation study

Here, we tested the effects of reducing the test set size (n=100) on performance. The lack of change between Fig. 4-C and C (test=100) indicates that a reduction in test set size is not responsible for improvements observed in Fig. 4-D.

Figure: Figure 4: Performance evaluations upon MoA-net (no inverse edges) with a test set of 100 triples. Metrics are presented as the average, with error bars representing standard deviation across five independent training/testing iterations. The lack of change between C and C (test=100) indicates that a reduction in test set size is not responsible for improvements observed in D.
Refer to caption: extracted/5924853/figures/hitsat1_n100.png

### A.7 XSwap permutations: MARS P 2 ⁢ H subscript MARS subscript 𝑃 2 𝐻 \text{MARS}_{P_{2H}} MARS start_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT 2 italic_H end_POSTSUBSCRIPT end_POSTSUBSCRIPT on MoA-net-10k

Using the XSwap algorithm as in Section 4.2, we checked, once again, whether the prediction metrics achieved using $\text{MARS}_{P_{2H}}$ on MoA-net-10k were influenced by degree bias. This time, there was a stark decrease in performance metrics upon the permuted KG (Fig. 5). This showed that predictions made by $\text{MARS}_{P_{2H}}$ were due to factors beyond node degree bias.

Figure: Figure 5: Performance evaluations of $\text{MARS}_{P_{2H}}$ on MoA-net-10k as well as a permuted variant of MoA-net-10k via the XSwap algorithm. Metrics are presented as the average, with error bars representing standard deviation across five independent training/testing iterations. A drop in performance metrics (E (permuted)) indicates that node degree was not the main driver in predictions made in E.
Refer to caption: extracted/5924853/figures/hitsat1_xswap.png

## Peer Reviews

confidence: 5
contribution: 2
initial_rating: 5
presentation: 3
questions: - Could you include more training details e.g. hardware, training time, validation loss plots etc.?
- On what basis were the "feasible" MoA's chosen in Table A1?
- When comparing with PoLo, which rewards function(s) were used?
soundness: 3
strengths: - Using the $ P_{2H} $ update rule, the new method has been shown, via a series of ablation studies, to avoid reasoning shortcuts due to degree bias.
summary: *November 26: after first response, the "Presentation" score has been raised from 2 to 3*

This work proposes a new neurosymbolic (NeSy) model for drug discovery as well as an associated dataset, where the problem is cast as a prediction task for drug-biological process (BP) pairs. More specifically, biochemical knowledge is arranged into a knowledge graph (KG) where entities (vertices) such as drugs, proteins and BP are connected with each other via directed edges. By using a new update rule of edge weights, the new model was able to avoid reasoning shortcuts inherent in previous methods.
weaknesses: - The novelty of the manuscript appears to be limited:
    - The new dataset has fewer kinds of entities than PoLo, and in particular does not appears to include those relevant for therapeutics such as "disease" or "side effect".
    - Similarly, one would expect more kinds of entities beyond a generic "protein" (e.g. receptor, transcription factor, enzyme etc.).
    - Another issue is that all metapaths are chains, like $ \mathrm{interacts}(P_{k}, P_{k + 1}) $, whereas in real-world applications there are many instance of multi-protein complexes for biological functions.
    - The model reward function is almost the same as PoLo (Equation 4): $$ R(S_{L + 1}) = 1_{ e_{L + 1} = e_d  }  + b \lambda  \sum\limits_{i = 1}^{m} s(M_i) 1_{\tilde{P} = M_i} $$ In fact, that paper states that the hyparameter $ b $ can be "set to $ 1_{e_{L + 1} = e_d} $" (p. 382).
- The presentation of the manuscript is not as efficacious as one would hope:
    - Key concepts such as "deconvolution" or $ \land $ (presumably conjunction?) were not defined.
    - The 2-hop probability update ($ \S $A.4 and Alg. 1), which is the major novelty in the manuscript, should be in the main text.
    - The main predictive task was not explicitly formulated: from Figure 1A it appears to be link prediction.
    - Moreover, it would be useful to have the catalogue of different entities--entities interaction types (Figure 1B) in a table.
    - The background colour of Figure 1B makes it harder to read, plus the font sizes for the labels are a bit too small.
    - The citation for PoLo should use the published version (https://link.springer.com/chapter/10.1007/978-3-030-77385-4_22)

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: - What is the computation time of training the network? How susceptible is the model with inclusion of new data?
- What is the "confidence" in figure 3?
- What is the exact neural network architecture for MARS? What is a simple example of input and output pairs?
- How would we approach in finding drug discovery? Does the model provide any candidate paths or of any sort? If so, how would we approach justification of the result?
- What is the "enhanced interpretability"? Would there be an example of this?
soundness: 3
strengths: The paper is generally well-written and easy to follow. The objectives of the proposed framework clearly outlined and the model effectively handles several specific problems that other comparable methods hold. The proposed method shows competitive performances and provides a reasonable form of interpretability using knowledge graphs.
summary: The paper proposed a systematic framework for interpretable drug discovery. The model comprises of drug mechanism-of-action (MoA) deconvolution; MoA-net and its variants, a specifically designed knowledge graph based on public biomedical data; MARS, a Neural Symbolic approach with through dynamic rule-weight updates. The work provides enhanced interpretability and strong predictive performances compared to other baselines.
weaknesses: - The overview of MARS is unclear, and a figure delineating process would be helpful for reader's understanding. This should include the input-output pairs for the proposed model.
- It is difficult to interpret Figure 3, as the main body does not explain anything about the confidences.
- While the paper proposes the term "drug-discovery", the it unclear in how to find  "novel drugs or chemical compounds". Given that the model is trained on the full data (no test) it could provide some meaningful results.
- The computation cost in training the model could provide valuable insights.
- The term "enhanced interpretability" is difficult to understand. An example for the term (possibly from the data that the authors used) would help better understanding of the term, and the strength of the model.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 2
questions: - Testing the model on real-world datasets, such as drug-protein interactions or clinical outcome data, would strengthen the claims regarding its applicability in drug discovery.
- How does MARS handle drug discovery tasks beyond MoA deconvolution?
- How could biological interpretability be further improved? Considering biological plausibility, could the authors incorporate features like protein binding affinities or pharmacokinetic properties to add biological context to MARS’s interpretability?
- Given the computational cost of two-hop joint probability calculations, what specific modifications or future optimizations are planned to enhance scalability?
soundness: 3
strengths: - MARS applies NeSy methods to drug discovery, a field with a strong demand for interpretability. The study’s introduction of MoA deconvolution as a KG-based task and MoA-net as a novel KG for biomedical applications broadens the use of KGs in drug discovery.
- The authors conduct extensive testing of MARS on both synthetic and real data and provide detailed comparisons to existing NeSy methods, particularly regarding susceptibility to reasoning shortcuts.
- By enhancing interpretability in drug discovery, MARS could help researchers gain insights into drug mechanisms, potentially leading to safer and more effective treatments.
summary: This study presents the MARS (Mechanism-of-Action Retrieval System), a neurosymbolic (NeSy) approach designed for interpretable drug discovery. MARS aims to improve understanding of drug mechanisms of action (MoAs) using knowledge graphs (KGs) combined with logical rule-based inference. The model introduces a tailored KG called MoA-net, enabling MoA deconvolution. MARS enhances interpretability by assigning weights to logical rules and employing a two-hop joint probability (P2H) metric to improve model calibration and address potential reasoning shortcuts.
weaknesses: - The model is primarily tested on synthetic KG data and lacks real-world validation on clinical datasets or pharmacological records, which limits the practical assessment of its interpretability and generalization for drug discovery tasks.
- While the P2H metric is introduced to make MARS shortcut-aware, its reliance on rule weighting and complex computations might hinder scalability, especially in large, densely connected KGs that are typical in biomedical data.
- AMARS mainly addresses computational rather than biological significance.
- MARS is highly susceptible to node degree bias, leading to unintended reasoning shortcuts.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: * In Lines 235-239, why do the authors state, "Although other metapaths are possible, we exclude them from our set of metapath-based rules"? What's the reason for excluding possible metapaths?
* Following the previous question, what is the rationale for limiting the length of metapaths to 4? Could this lead to insufficient exploration of reasoning paths? If possible, could a comparative experiment on different lengths be conducted?
* In Lines 277-282, the comparison baselines seem to be the latest from 2020. Are there any recent works related to knowledge graphs?
* Why do a significant portion of the results in the experimental tables lack standard deviation? For example, the top 9 baseline models in Table 1.
* Please provide visualizations of the learned rule weights and corresponding analyses.

If the authors could adequately address my concerns, I will consider raising the score.
soundness: 2
strengths: * The integration of neurosymbolic approaches provides a better understanding of the evidence chain in drug discovery.
* The dynamically learned rule weights may offer an interpretable explanation for the potential impacts of reasoning shortcuts.
* The experimental results are relatively comprehensive and reasonable.
summary: The paper proposes MARS, a neurosymbolic model designed to improve interpretability in drug discovery through a task called Mechanism of Action (MoA) deconvolution, using a specialized knowledge graph named MoA-net. MARS intends to offer insights into the interactions among drug, protein, and biological process via learned rule weights.
weaknesses: * The paper is poorly structured and written, and quite difficult to follow.
* Baseline methods compared in the experiments appear to be outdated, with the latest only dating back to 2020.
* Although the authors claim that interpretable symbols can enhance the learning process, there is insufficient experimental support for this assertion, especially regarding the discovery of unseen reasoning paths.
* While this work is presented as a reasonable exploration of biologically meaningful evidence chain reasoning shortcuts, the authors only apply it to three types of entity relationships. Existing study [1] has researched more comprehensive biological pathway evidence chain mining using symbolic reasoning and reinforcement learning. A detailed comparison of the similarities and differences between the two approaches, along with a more substantial discussion of the advantages of your work, is needed.

**References**

[1] Sudhahar, S., Ozer, B., Chang, J. et al. An experimentally validated approach to automated biological evidence generation in drug discovery using knowledge graphs. Nat Commun 15, 5703 (2024).

Meta-review: This paper introduces a neurosymbolic system for drug discovery, which consists of drug mechanism-of-action (MoA) deconvolution; MoA-net and its variants, a specifically designed knowledge graph based on public biomedical data; MARS, a Neural Symbolic approach through dynamic rule-weight updates. The paper is about AI for science (biology and drug discovery), neurosymbolic approaches have the potential for interpretable drug discovery, and the paper is clearly written. However, reviewers raised concerns that the experiments conducted on a dataset don't sufficiently capture the biochemical complexity. For example, there is no differentiation between kinds of proteins (e.g. receptors, enzymes, transporters), drug-protein interactions should be more than just upregulates and downregulates, and drug acts by binding to multiple target proteins that are jointly responsible for a biological process. Thus, there lack of concrete case studies to demonstrate its real-world value in drug discovery.

Decision: Reject