1 Background↩︎

Drug discovery (DD), the search for novel drugs to treat ailments, often involves screening thousands of small chemical compounds [1]. Many computational approaches have been developed to accelerate and streamline this screening process [2], [3]. 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 [4]. Typically, DD is formulated on a KG as a link prediction task between drugs and the corresponding medical conditions (indications) to be treated [5], [6].

In addition to predicting drug indications, it is also important to understand a drug’s mechanism-of-action (MoA), the molecular processes by which it achieves its medicinal effect. As depicted in Fig. 1, MoAs typically involve chains or paths of physical, molecular interactions induced by a drug [7]. Revealing a drug’s MoA, which we call MoA deconvolution, informs researchers as to how each drug works, de-risks potential side effects [8], [9], and helps to predict clinical trial success [10]. Unfortunately, traditional approaches for MoA deconvolution involve additional laboratory assays, which do not scale well to the magnitudes at which DD screening takes place [1], [11]. Our approach, called the MoA Retrieval System (MARS) [12], could serve as a unique, AI-based alternative.

Specifically, MARS is a neurosymbolic (NeSy) approach. NeSy approaches combine logical rules with neural networks [14], so they tend to possess enhanced interpretability compared to state-of-the-art techniques on KGs [15]. Model interpretability involves transparency into the processes or patterns which led to certain predictions [16]. Therefore, NeSy approaches have already been positioned as a promising avenue for MoA deconvolution [17], [18]. However, interpretability is broadly defined [16], which poses an additional challenge: there are no clear guidelines for assessing the plausibility of model interpretations, especially for this novel task.

Although some previous studies present explainable or interpretable pipelines [17], [19], [20], the corresponding explanations leverage associative patterns: two nodes with mutual connections are likely to share other connections [21]. For example, such methods utilize associations regarding a drug’s pharmacological class [19], side effects [17], or known indications [20]. These associative patterns, which are discovered during or after clinical trials, might be useful for repurposing clinically approved drugs [5], [6]. However, associative patterns are rare or absent for novel compounds which have not yet undergone clinical trials. Furthermore, such patterns can not represent the MoA of a drug; instead, an MoA involves mechanistic patterns. For example, Figure 1 shows the MoA of cortisone acetate, which involves physical, molecular interactions [13]. Therefore, we define computational MoA deconvolution as the prediction of mechanistic paths between drugs and their biological effects, such as biological processes (BPs). BPs are a class of the Gene Ontology (GO) [22] representing the outcomes of molecular interactions.

To benchmark this novel task, we generate a tailored KG, called MoA-net [23], from real-world, experimental data. Unlike existing KGs that are designed for DR [24], [25], the entities and relations in MoA-net exclusively comprise molecular interactions which could be captured in a typical laboratory setting. To perform MoA deconvolution, MARS is used upon MoA-net to predict the links between nodes representing drugs and the BPs they cause.

Several interpretable features of MARS provide insight into the drug’s MoA, represented as a path between the two nodes. However, MARS’ interpretable features also reveal a reasoning shortcut. A reasoning shortcut describes a situation in which predictions are based on unintended semantics [26], [27]. For example, take the classic MNIST addition task, popularly used to assess NeSy methods [28]. To accomplish the MNIST addition task, a model is trained to determine the sum of two handwritten digits. In some cases, models were found to predict the correct sum, despite having misclassified the digits [29]. In this study, we found that predictions upon MoA-net are often driven by degree-bias, rather than the rules representing domain knowledge. Degree-bias occurs when a KG’s node degree distribution, describing the number of edges connected to each node, varies greatly [30]. Consequently, this affects how accessible, discoverable, and influential certain nodes might be to various statistical and AI-based methodologies [31], [32].

To address this reasoning shortcut, we consider the desiderata from Marconato et al. [33] 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.

1.1 Related Work↩︎

1.1.1 Interpretable methods for drug discovery↩︎

Several previous methods for interpretable, computational DD lay the groundwork for MARS. Specifically, several of such methods use metapaths. Metapaths are abstract representations of paths in a KG [24], [34], and they are particularly useful for capturing long-range dependencies between nodes which are several hops away from one another. Project Rephetio [24], in particular, was foundational for the use of metapaths in DD. The authors of Project Rephetio prioritized drug-disease associations by evaluating the prevalence of specific metapath patterns across a large, heterogeneous KG. To do so, they used degree-weighted path counts (DWPC) [35]. In short, DWPCs enumerate all possible paths between source and target nodes, taking node degree into account. From this information, they predicted new drug-disease relations, and the metapaths used to make such predictions served as model interpretations. Another study by Kawichai et al. [36] used a similar approach, but they extracted metapaths involving relevant GO entities, like BPs (see Section 1). Their approach could effectively predict novel drug-disease associations to inform research on drug repurposing. Becuase metapaths were shown to be both useful and interpretable tools for this task, other studies applied them in different ways. For example, Ratajczak et al. [19] leveraged metapaths from Project Rephetio, but they used them, instead, to filter the KG down to a more informative subgraph for drug repurposing. Another method, called Bioteque [20], obtained representative embeddings of KG paths following each metapath. Like the aforementioned studies, they also found that their embeddings were useful for predicting novel drug-disease associations. Although our article discusses the use of mechanistic versus associative patterns used by KG-based methods, a variety of methods exist which capture different concepts. For example, the GFlowNet [37] uses paths in a graph to represent potential synthesis routes, involving chemical reactions, between reagents and a resultant drug. For a wider scope of interpretable computational DD, we refer the interested reader to broader reviews on this topic [38], [39].

1.1.2 NeSy AI for MoA deconvolution↩︎

Several NeSy approaches involve logical rules reflecting path-like patterns in biomedical KGs for DD-related link prediction tasks. For example, an approach by Tian et al. [40] used a combination of path patterns (similar to that in Section 2.3) and graph neural networks (GNNs) to predict drug side effects in an interpretable manner. Similarly, to predict drug indications, Sudhahar et al. [41] investigated evidence chains, paths explaining associations between drugs and diseases. However, these explanations were derived separately from indication predictions, via an additional rule-mining model [42], much like that of Renaux et al. [43], who used rule mining methods to unveil disease-causing gene interactions. Other approaches [17], [18], [44] accomplished biomedical link prediction for a variety of DD-related tasks through deep reinforcement learning (RL), in which a neural network contributes toward the optimization of a reward function [45]. One prominent example of such an approach is PoLo [17], which is based on MINERVA [18]. MINERVA involves a deep RL agent which traverses paths in a KG, receiving a reward when the first and last nodes of its path are known to have some relation of interest. PoLo expanded upon this by including an additional reward when the agent’s trajectories utilized a set of predefined, informative metapaths. Critically, they applied a static, literature-derived weight to each of the metapaths to quantify how important each one should be to the agent. Consequently, after training, predictions were expected to utilize and align, to some extent, with the most important metapaths. Another study by Drancé et al. [44] expanded on this idea by deriving the metapaths and their corresponding weights through rule-mining. Although weights remained static, they were mined directly from the KG, rather than sourced from an external study. In contrast to these studies, MARS learns weights for each metapath, so the resultant weights are a more direct measure of relative importance toward predictions. Additionally, while all these studies found that metapaths representing associative patterns were most frequently used, our study focuses upon paths involving mechanistic, molecular relations. In this study, we also identify a major risk: the approach may neglect to utilize rules in favor of other semantics for reward optimization. This results in reasoning shortcuts.

1.1.3 Trustworthy NeSy methods and reasoning shortcuts↩︎

Like the deep RL approaches described above, many other NeSy approaches are designed to abide by rules and domain knowledge [44], [46]. For example, a study called Biomedical KG refinement with Embedding and Rules, or BioGRER [47], used domain knowledge on the COVID-19 disease to construct relevant rules about it. By alternating between rules and a representative KG embeddings, BioGRER predicted novel drug indications and other biologically relevant information on COVID-19. Another NeSy approach, called Walking RDF and OWL [48], presented a pipeline for predicting drug indications in which ontological rules helped to predict new relations and augment the KG. Subsequently, a GNN was used to refine the augmented KG. Unlike the aforementioned approaches which utilize metapaths, methods which utilize GNNs or KG embeddings often suffer from local receptive fields [49], experiencing peak performances at two layers [50], [51]. Nonetheless, they behave similarly in that they use rules to guide or inform a neural network as well as to explain biomedical link predictions. Thus, the use of rules and domain knowledge in NeSy AI often portray such approaches as more trustworthy than neural, black box ones [52]. Recent studies, however, have found that NeSy approaches may suffer from reasoning shortcuts, in which a model predicts the correct outcome via unintended semantics [26], [27]. While reasoning shortcuts are not exclusive to NeSy methods [53], [54], they may be more easily overlooked when such approaches are portrayed as trustworthy.

2 Methods↩︎

Here, we present the MoA Retrieval System (MARS)¹ [12] to perform MoA deconvolution. MARS improves PoLo [17] (Section 1.1) by introducing dynamic, learned rule weights. This differs from previous approaches, where weights are static and pre-computed (e.g., mined or literature-derived) [17], [44]. As discussed further, these learned weights also make MARS shortcut-aware. To test MARS, we designed a novel KG, called MoA-net, comprising data that is particularly relevant to predicting drug MoAs for new compounds, as opposed to approved ones.

2.1 Datasets: MoA-net and its variants↩︎

We designed our KG, MoA-net² [23], specifically for MoA deconvolution. MoA-net consists of triples, two nodes connected by an edge, representing relations between drugs, proteins, and BPs (Table 1). In this study, we depict triples in MoA-net as binary predicates. For example, \(\texttt{interacts}(\textsl{Protein}, \textsl{Protein})\) states that two \(\textsl{Protein}\) nodes are connected via the interacts relation. Such protein-protein interactions (PPIs) are derived from several real-world datasets comprising experimental data, including Custom KG [55] and OpenBioLink KG [25]. In addition to PPIs, MoA-net contains causal relations between drugs and proteins, such as upregulates and downregulates, also from Custom KG and OpenBioLink KG. The BP nodes as well as the causal relations between proteins and BPs come from experimentally-derived and expert-curated GO annotations available in the UniProt database [56] (release 2023.05), a gold-standard collection of protein data. Such GO annotations are entered into UniProt by several groups of designated researchers, developers, and curators [57], and they denote which GO terms, including BPs, pertain to each protein in UniProt.

While other biomedical KGs exist [24], [25], MoA-net is unique from other KGs because it comprises drug-BP triples. These are derived from publicly available functional and biochemical assays in ChEMBL (v33), an open access database of bioactive compounds [58]. BPs are the fundamental processes for the regulation, defense, and communication of cells, including processes like signal transduction or inflammation [22]. Thus, MoA-net is particularly useful for DD within the context of novel compounds because it comprises entities and relations which could be discovered in a laboratory setting before clinical trials. By representing a drug’s therapeutic effect as a BP, framing DD within a pre-clinical context is possible.

Between the three node types, five unique edge types, or relations, are defined and shown in Table 1. Notably, the inverse relations of all edge types in Table 1, which run in the opposite direction of causality, are also included. Of the 1,622 drug-BP triples in MoA-net, 48 also had known MoAs in a database called DrugMechDB [13], a manually-curated compendium of known drug MoAs. DrugMechDB contains evidence-derived, mechanistic paths, or MoAs, between drugs and their respective molecular and physiological effects. Specifically, like MoA-net and other widely used biomedical KGs [24], [25], paths in DrugMechDB involve biomedical entities sourced from gold-standard biomedical databases.

Finally, we created two variants of MoA-net. To investigate reasoning shortcuts, we used the Zietz et al. [30] implementation of XSwap [59], which swaps edges in a KG without affecting the distribution of node degrees. We called the resultant KG MoA-net-permuted. Additionally, MARS, described in Section 2.2, includes an optional, automatic KG trimming step to sparsify the KG. In essence, this step reduces edges of each class to below a user-specified threshold by iteratively removing those between the highest-degree nodes. To create a sparse variant of the MoA-net, we set this threshold to 10,000 PPIs. This specific value was selected because it reduced the PPIs to \(\sim 50\%\) of the total edges, striking a balance between mitigating the skewness of the node degree distribution and preserving the global connectivity of the KG. We refer to this trimmed subgraph as MoA-net-10k (see Section 3.3).

2.2 MARS: a NeSy approach for MoA deconvolution↩︎

As shown in Fig. 2, MARS takes two major inputs. The first involves a KG. In this case, we test MARS on MoA-net and its variants (see Section 2.1). Specifically, the input KG into MARS must contain triples involving some relation of interest. As shown in Figure 1, some information about drug candidates, such their indications, is discovered during or after clinical trials, so this information is typically unavailable for novel compounds. Therefore, to understand each MoA as the biological response to drug administration, we investigate MoAs as the relations between drugs and BPs (i.e., \(\texttt{induces}(\textsl{Drug}, \textsl{BP})\)). Specifically, we formulate this as a link prediction task. In other words, we task MARS with determining whether edges of type \(\texttt{induces}\) exist between \(\textsl{Drug}\) and \(\textsl{BP}\) nodes in MoA-net. Thereafter, new predictions regarding the \(\texttt{induces}\) relation serve as potential therapeutic outcomes for the chemical compound represented by the \(\textsl{Drug}\) node. As discussed in Section 2.1, this is a novel application in the KG field.

The second input into MARS, as depicted in Figure 2, includes metapath-based rules with corresponding weights.

2.3 Metapath-based rules↩︎

As mentioned in Section 1.1, metapaths are abstract representations of instantiated paths in a graph [24], [34], [60]. For example, given the following path, \(P\), in MoA-net:

\[\begin{align} \textsl{Cortisone acetate} &\xrightarrow{upregulates} \textsl{GC receptor} \notag \\ & \xrightarrow{interacts} \textsl{COX protein} \xrightarrow{participates} \textsl{Inflammation} & \label{cortisone95acetate} \end{align}\tag{1}\]

the corresponding metapath, \(\tilde{P}\), would be: \[\begin{align} \textsl{Drug} \xrightarrow{upregulates} \textsl{Protein} \xrightarrow{interacts} \textsl{Protein} \xrightarrow{participates} \textsl{Biological Process (BP)}\label{cortisone95acetate95moa} \end{align}\tag{2}\]

Within this study, metapaths can be understood as a sequence of triples within the KG structure, making them inherently interpretable. In MARS, metapaths are used as the bodies of logical rules, in which triples are connected by logical conjunctions (\(\land\)). Conjunctions indicate that, if all triples in the rule body are true, then the rule head is evaluated as true. In other words, the rule is satisfied. The rule head, the left side of the implication arrow (\(\Leftarrow\)), is a single triple representing the relation of interest between the first and last node types of the metapath, e.g.,: \[\begin{align} \texttt{induces}(\textsl{Drug}, \textsl{BP}) &\Leftarrow \texttt{upregulates}(\textsl{Drug}, \textsl{Protein}_A) \land \notag \\ &\hphantom{\Leftarrow} \texttt{interacts}(\textsl{Protein}_A, \textsl{Protein}_B) \land \texttt{participates}(\textsl{Protein}_B, \textsl{BP}) & \label{example95rule} \end{align}\tag{3}\]

For each metapath-based rule, \(M_i\), in a set of rules, \(\mathcal{M} = \{M_1, M_2, ..., M_m\}\), the rule weight is denoted 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 2.5, we discuss how we initialize and compute these weights.

Using the hetnetpy package [61], we extracted all valid MoA patterns as metapaths from MoA-net (see Section 2.1). We defined valid MoA patterns as metapaths comprising directed, mechanistic paths between drug and BP nodes (see Table 3 in Appendix 6). Based on MoAs found in DrugMechDB, metapath length was limited to a maximum of four relations (or hops). We excluded metapaths depicting associative patterns, such as those leveraging information about shared BP targets, from the set of metapath-based rules. Specifically, associative patterns were identified based on paths in which a node was receiving directed, incoming edges from nodes of the same type. Finally, each of these metapaths served as the body of a logical rule with \(\texttt{induces}(\textsl{Drug}, \textsl{BP})\) as the rule head, as in formula 3 .

2.4 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. Here, that relation is \(\texttt{induces}(\textsl{Drug}, \textsl{BP})\), which are masked from the agent during training. Each walk generates a path, \(P\), such as the one in the previous section 2.3. This path, \(P\), can also be understood as a series of \(L\) transitions: \(P := (e_{Drug} \xrightarrow{r_1} e_2 \xrightarrow{r_2} ... \xrightarrow{r_L} e_{L+1})\). With respect to classic RL terminology, each node can be considered a state, and each relation or transition between nodes can be considered an action. 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_{Drug}\) (the drug node), and the terminal node, \(e_{L+1}\), have the induces relation. By training the agent to do so, it can identify node pairs with the desired relationship, thus generalizing beyond the training set to predict novel pairs in a holdout, test set. In other words, while true positive predictions in the test set serve as validation, false positives are positioned as potentially novel \(\texttt{induces}(\textsl{Drug}, \textsl{BP})\) predictions.

Similarly to a Markov Decision Process [62], the agent moves based on its current position and its next possible actions, with no information about its target destination. Additionally, at each state or node in the agent’s trajectory, the history of the agent’s previous actions are encoded with a two-layer LSTM (see Appendix 8). Given the current node, \(e_c\), its history would comprise entity and relation embeddings (\(\boldsymbol{\textit{e}}\) and r, respectively) of each entity and relation traversed to that point:

\[\boldsymbol{\textit{h}}_c = \text{LSTM}(\boldsymbol{\textit{h}}_{c-1}, \boldsymbol{\textit{r}}_{c-1}, \boldsymbol{\textit{e}}_{c-1}) \label{lstm}\tag{4}\]

In this case, entity and relation embeddings can be understood as type-specific, learnable parameters of the LSTM. Using Eq. 4 , the agent’s next action can be modeled as a probability distribution, \(\boldsymbol{\textit{a}}_c\), where \(\boldsymbol{\textit{W}}_1\) and \(\boldsymbol{\textit{W}}_2\) are learnable weight matrices, and \(\boldsymbol{\textit{A}}_c\) denotes the action space (adjacent nodes) from the current node:

\[\boldsymbol{\textit{a}}_c = softmax(\boldsymbol{\textit{A}}_c(\textit{\boldsymbol{W}}_2 \text{ReLU}(\boldsymbol{\textit{W}}_1 [\textit{\boldsymbol{h}}_c; \boldsymbol{\textit{e}}_c]))) \label{action95dist}\tag{5}\]

Ultimately, MARS’ policy, is to draw actions from a categorical probability distribution parametrized by Eq. 5 . The parameters are trained to optimize the reward function, \(R(S_{L+1})\) (Eq. 6 ), which is evaluated each time the agent completes \(L\) transitions from some starting node:

\[R(S_{L+1}) = \mathbb{1}_{\{e_{L+1} = e_{BP}\}} + \mathbb{1}_{\{e_{L+1} = e_{BP}\}}\lambda\sum_{i=1}^{m} w(M_i) \mathbb{1}_{\{\tilde{P} = M_i\}}\label{eq:reward}\tag{6}\]

\[\mathbb{1}_{\{A\}} = \begin{cases} 1 & \text{if } A = \text{true} \\ 0 & \text{if } A = \text{false} \end{cases}\]

In short, the reward function, originally from Liu et al. [17], quantifies how successful \(P\) is according to two rewards (Figure 2). The first, base reward indicates whether the terminal node in the path, \(e_{L+1}\), is one of the desired target (BP) nodes (\(e_{BP}\)) that forms a known triple or true pair with the starting (drug) node, \(e_{Drug}\). Put simply: “given an \(\texttt{induces}(\textsl{Drug}, \textsl{BP})\) triple, did the agent make a successful traversal between the drug and BP nodes?"

The second, supplementary reward, contingent upon the first, indicates whether the corresponding metapath, \(\tilde{P}\), matches any metapath-based rule, \(M_i\). The second reward is also proportional to the rule’s corresponding weight, and MARS updates these weights during training. These updates are accomplished through a novel algorithm called two-hop joint probability, or \(P_{2H}\) (Section 2.5). 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. Thus, MARS has two key interpretable features: (1) paths between nodes which serve as potential MoA predictions, and (2) learned rule weights which serve as a proxy for the importance of each rule.

2.5 MARS dynamically updates rule weights↩︎

After MARS is executed, its learned rule weights reflect each rule’s relative usefulness in the prediction task. Rule weights are updated after every batch, which includes a subset of the training triples. This means that, during training, the weight updates drive the agent toward more informative paths and bypass the assumptions that pre-assigned rule weights are correct. This eliminates the need for pre-computed or literature-derived rule weights; thus, we initialize all weights uniformly as 0.5, a medium level of importance. We tested two different approaches for updating rule weights.

2.5.1 Naive weight updates↩︎

The first method for updating rule weights is called naive updates, and it lays the groundwork for the \(P_{2H}\) updates, covered in the next section. 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, \(M_i\) is satisfied. In other words, frequency serves as a heuristic for how often the agent makes successful predictions with a given metapath. Formally, this involves recording the observed frequency, \(O_{M_i}\), at which each metapath-based rule is satisfied in each batch. The observed frequency is then normalized by the batch-specific expected frequency, \(E\), which assumes that every rule has a uniform probability across the total number of occurrences in that batch. Therefore, \(E\) is the same for every metapath-based rule. Ultimately, this produces a metric, \(\mu_{M_i}\) (Eq. 7 ) in which \(\mu_{M_i} > 1\) indicates usefulness (the agent used that rule more than others), and \(\mu_{M_i} < 1\) indicates otherwise.

\[\mu_{M_i}= O_{M_i} / E\label{eq:normalization}\tag{7}\]

Then, using Eq. 8 , denoted by \(\Phi\), and Eq. 9 , \(\mu_{M_i}\) is used to update the weight of a rule, \(w(M_i)\).

\[\Phi(\mu_{M_i}, w(M_i)) = w(M_i) \times 2\alpha(\frac{\mu_{M_i}-1}{\mu_{M_i}+1})\label{eq:update}\tag{8}\]

\[w(M_i)' = w(M_i) + \Phi(\mu_{M_i}, w(M_i))\label{eq:new95weight}\tag{9}\]

Several measures are taken to allow weight updates to happen gradually. Firstly, Eq. 8 is regularized by the hyperparameter \(\alpha \in \mathbb{R}\), where \(0 \leq \alpha \leq 1\), to control how subtle or drastic the weight update is, respectively. If \(\alpha = 0\), no weight updates are made. We selected \(\alpha\) via hyperparameter optimization (Appendix ¿tbl:hyperparameters?). Secondly, given that \(\mu_{M_i} \geq 0\), the range of Eq. 8 is \([-w(M_i)2\alpha, w(M_i)2\alpha]\). Therefore, if \(\alpha=1\), a weight can only, at most, triple during a single update. Thirdly, weight updates are executed after each batch, but batches only capture a fraction of possible trajectories through the KG. To increase the variety of metapaths seen, MARS’ agent performs multiple traversals, or rollouts, for each triple in a batch. Despite this, however, certain metapaths may still go under- or over- sampled in a single batch, resulting in drastic weight fluctuations.

To avoid extreme values that may result from sampling biases, \(\mu_{M_i}\) is further restricted as described in Eq. 12 by lower and upper bounds in Eqs. 10 and 11 . Specifically, \({\text{batch size} \times \text{rollouts}}\) describes the number of opportunities for which the agent could traverse a metapath, so the bounds are designed relative to this. For naive updates, \(\rho\) is the total number of metapath-based rules.

\[\mu_{min}=\frac{\rho}{\text{batch size} \times \text{rollouts}} \label{eq:lowerbound}\tag{10}\]

\[\mu_{max}=\rho \times \text{batch size} \times \text{rollouts} \label{eq:upperbound}\tag{11}\]

\[\mu_{M_i} = \text{min}(\mu_{max}, \text{max}(\mu_{min}, \mu_{M_i})) \label{eq:bound95mu}\tag{12}\]

To understand the intuition behind Eqs. 12  and 10 , recall that \(\mu_{M_i}\) is a ratio to approximate the usefulness of a rule. If \(\mu_{M_i}<1\), the corresponding weight, \(w(M_i)\), will decrease, and if \(\mu_{M_i}>1\), it will increase. If \(\rho\) is large, it will be more difficult for the agent to sample all metapath-based rules. Therefore, to reduce the impact of under-sampling on a weight update, Eq. 10 defines a larger lower bound in such cases. If \(\rho\) is small, the agent is likely to encounter each metapath-based rule more frequently within a single batch. In this scenario, over-sampling is a risk, so Eq. 11 scales the upper bound constrains a positive weight update by the relative size of the rule space. Finally, if the agent finds zero occurrences of any metapath-based rule, no weight updates are made at all.

2.5.2 \(P_{2H}\) weight updates↩︎

The second, more complex method to update weights, which is also a novel component of this work, is coined two-hop joint probability (\(P_{2H}\)). In essence, this metric approximates the usefulness of metapath-based rules based on full and partial matches. Pseudocode for \(P_{2H}\) can be found in Algorithm 1 below. In summary, each two-hop fragment is extracted from every path traversed by the agent (lines 2-8).

Two-hop fragments are defined as the metapaths involving two consecutive relations in a given path or trajectory. For example, the two-hop fragments of the path expressed in Eq. 1 would be:

• \(\texttt{upregulates}(\textsl{Drug}, \textsl{Protein}) \land \texttt{interacts}(\textsl{Protein}, \textsl{Protein})\), as well as

• \(\texttt{interacts}(\textsl{Protein}, \textsl{Protein}) \land \texttt{participates}(\textsl{Protein}, \textsl{BP})\)

After all the two-hop fragments are extracted from the agent’s trajectories in a given batch, weight updates can be computed similarly to naive updates (Section 2.5.1). Here, \(E\) denotes the expected frequency of each fragment, computed as the reciprocal of the number unique fragments encountered (line 9). The observed frequencies are also adjusted in a similar manner, describing the number of times each two-hop fragment, \(f_j \in \mathcal{F}\), is encountered (\(O_{f_j}\)) (lines 10-12). With these alterations, Eq. 7 is used to compute \(\mu_{f_j}\).

To compute weight updates, the two-hop fragments are also determined for each of the metapath-based rules (line 14). For example, take the following metapath-based rule, \(M_{example}\): \[\begin{align} M_{example} &:= \texttt{induces}(\textsl{A}, \textsl{E}) \Leftarrow \texttt{upregulates}(\textsl{A}, \textsl{B}) \land \texttt{interacts}(\textsl{B}, \textsl{C}) \land \notag \\ &\hphantom{:= \texttt{induces}(\textsl{A}, \textsl{E}) \Leftarrow \texttt{upregulates}} \texttt{interacts}(\textsl{C}, \textsl{D}) \land \texttt{participates}(\textsl{D}, \textsl{E}) & \label{example95metapath} \end{align}\tag{13}\] where two-hop fragments are pairs of binary predicates which share a variable:

• \(\textsl{Fragment 1}: \texttt{upregulates}(\textsl{X}, \textsl{Y}) \land \texttt{interacts}(\textsl{Y}, \textsl{Z})\)

• \(\textsl{Fragment 2}: \texttt{interacts}(\textsl{X}, \textsl{Y}) \land \texttt{interacts}(\textsl{Y}, \textsl{Z})\)

• \(\textsl{Fragment 3}: \texttt{interacts}(\textsl{X}, \textsl{Y}) \land \texttt{participates}(\textsl{Y}, \textsl{Z})\)

Here, the probability of each metapath-based rule is computed as the joint probability of its fragments (line 15 of Algorithm 1). For example, the \(P_{2H}\) metric for the metapath-based rule shown in formula 13 would be computed as \(P_{2H}(M_{example}) = p(\textsl{Fragment 1}) \times p(\textsl{Fragment 2}) \times p(\textsl{Fragment 3})\). Notably, there are two caveats. Firstly, to account for partial metapath matches, the definition of conjunction here is relaxed, allowing truth to be evaluated on the fragment level. Secondly, metapath fragments do not necessarily represent independent events, but the way in which the \(P_{2H}\) metric is computed suggests otherwise. However, to avoid complex computation involving conditional probabilities for dependent events [63], independence is assumed, and the \(P_{2H}\) metric serves as an approximation for the empirical probabilities of metapath-based rules. From this \(P_{2H}\) metric, weights can be updated using Eqa. 8  and 9 from Section 2.5.1 (line 16 of Algorithm 1). Notably, for \(P_{2H}\) updates, \(\rho\) in Eqs. 11 and 10 is the total number of unique two-hop fragments possible. Ultimately, MARS with \(P_{2H}\) updates (\(\text{MARS}_{P_{2H}}\)) uses all information from successful trajectories.

We describe implementation details and hyperparameter selection in Appendices 7 and 8, respectively.

2.6 Evaluation↩︎

As mentioned in Section 2.2, we test MARS upon MoA-net, and performance is assessed based on link prediction over drug-BP triples. Therefore, we split the drug-BP triples within MoA-net into training (60%), validation (20%), and test (20%) sets. Models are evaluated 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 [4]. In addition to these standard metrics, we report pruned metrics: these are computed solely on a subset of the predictions that utilized one of the pre-defined metapath-based rules (see Table 3), 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 MARS against five different baseline KGE models, two state-of-the-art NeSy RL methods, and one network measure. Based on previous biomedical KG benchmarks [6], the KGE baseline models include ComplEx [64], RotatE [65], MuRE [66], CompGCN [67], and PairRE [68]. The two RL baseline models include PoLo [17] and its predecessor MINERVA [18] (Section 1.1), which is not guided by rules. We describe hyperparameter selection for baseline models in Appendix 8. Lastly, the network measure baseline includes prioritization of drug-BP triples based on DWPC (see Section 1.1). We train and evaluate all baseline 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 for a comparison between MARS’ predicted MoAs of 48 drug-BP pairs against their known MoAs (see Section 2.1). Secondly, learned rule weights serve as a proxy for the importance of each metapath-based rule. This helps to determine whether agent trajectories are biased toward certain types of paths. Alongside the pruned metrics, these features help us evaluate MARS’ alignment with domain knowledge.

3 Results↩︎

Within the next sections, we present the results of \(\text{MARS}_{P_{2H}}\) and compare them against relevant baseline models. Specifically, we report standard and pruned metrics (see Section 2.6) to evaluate performance as well as Marconato et al.’s [33] shortcut-aware desiderata, introduced in Section 1. In summary, we initially found that both variants of MARS, as well as PoLo, make predictions which, primarily, avoid using metapath-based rules. Instead, all three methods utilized paths following associative patterns. After several subsequent experiments, we verified that this reasoning shortcut was likely due to node degree bias in MoA-net. Thereafter, we found that \(\text{MARS}_{P_{2H}}\) can mitigate this reasoning shortcut on a trimmed version of MoA-net.

3.1 Associative patterns improve accuracy but offer limited practical use↩︎

After an initial set of experiments on MoA-net, the results were already different than expected. Specifically, the pruned metrics were consistently lower than standard ones (Figure 4-\(\alpha\)), 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. 6 ). Additionally, amongst recorded trajectories, most did not follow the metapath-based rules; instead, most trajectories used the following associative pattern, involving inverse edges: \[\begin{align} \texttt{induces}(\textsl{Drug}_1, \textsl{BP}_2) &\Leftarrow \texttt{induces}(\textsl{Drug}_1, \textsl{BP}_1) \land \texttt{induces}(\textsl{Drug}_2, \textsl{BP}_1) \notag\\ &\hphantom{\Leftarrow \texttt{induces}(\textsl{Drug}_1, \textsl{BP}_1)} \land \texttt{induces}(\textsl{Drug}_2, \textsl{BP}_2) & \label{mars95assoc95pattern} \end{align}\tag{14}\] This associative pattern indicates that two drugs inducing a common BP also likely induce another BP. This type of pattern was also present in the original PoLo study by Liu et al. [17], in which the most used pattern was the following: \[\begin{align} \texttt{treats}(\textsl{Drug}_1, \textsl{Disease}) &\Leftarrow \texttt{causes}(\textsl{Drug}_1, \textsl{Side Effect}) \land \texttt{causes}(\textsl{Drug}_2, \textsl{Side Effect}) \notag\\ &\hphantom{\Leftarrow \texttt{causes}(\textsl{Drug}_1, \textsl{Side Effect})} \land \texttt{treats}(\textsl{Drug}_2, \textsl{Disease}) & \label{polo95assoc95pattern} \end{align}\tag{15}\]

In Liu et al., the idea that pattern 15 was commonly used by the agent had served as validation that the agent was, indeed, using the metapath-based rules. This is because pattern 15 was one of the metapath-based rules which could be satisfied to achieve the supplementary reward (see Section 2.4). Additionally, in its original study, PoLo had obtained high standard and pruned metrics, thereby reinforcing this notion. However, high performance metrics do not guarantee biologically meaningful results. When the results from Liu et al. [17] are reproduced, the same standard metrics can be achieved – even when pattern 15 is not used to guide the agent (see Appendix 9). Additionally, pattern 15 is still the most commonly used pattern, even when it is not used as a metapath-based rule. This suggests that, although this associative rule may serve as a plausible model explanation, it does not necessarily guide model training. From this, it was clear that neither PoLo’s nor MARS’ agents were using the supplementary reward in Eq. 6 . In the next section, MARS’ interpretable features helped to understand why this was happening.

3.2 \(P_{2H}\) updates reveal reasoning shortcuts via degree bias↩︎

In addition to analyzing the agent trajectories, we used the rule weights with \(P_{2H}\) updates to assess how informative each of the metapath-based rules was in making predictions. In particular, we found that \(\text{MARS}_{P_{2H}}\) weights for paths involving consecutive PPIs (see Section 2.1) were consistently less important (Figure 5). This indicated that the agent avoided exploring consecutive PPIs.

Previous research on KGs has shown that node degree distribution, the number of adjacent edges for each KG node, can significantly bias predictions [31], [32]. Specifically, inspection bias, a type of degree bias, occurs when the KG is not uniformly inspected or sampled [30]. Since PPIs are the most common relation type in MoA-net (90% of edges) (Table 1), protein nodes have a higher degree distribution than other node types. From this evidence, we hypothesized that the agent was circumventing 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 MARS upon MoA-net-permuted. As explained in Section 2.1, 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 (Figure 4-\(\beta\)). Put simply, the agent gets lost when exploring the PPIs, so it avoids them.

3.3 Identifying and mitigating degree bias improves performance↩︎

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

Next, to confirm that the agent was getting lost within the PPIs, we employed MARS’ automatic trimming step. As explained in Section 2.1, MoA-net-10k is a variant of MoA-net with fewer PPIs. We tested each of \(\text{MARS}_{P_{2H}}\), \(\text{MARS}_{\text{naive}}\), and PoLo with the same parameters upon on MoA-net-10k (Figure 4-\(\gamma\)). As before, inverse edges were excluded. Since trajectory length was set to \(L=4\), MARS’ trimming approach also 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 (Figure 4-\(\alpha\)). To ensure this improvement was not simply the result of a reduced test set, we tested each of the three approaches upon MoA-net with 100 sampled test triples, which showed no change (see Appendix 10).

While removing inverse edges improved metrics, a shortcut-aware system should achieve high performance even with the shortcut present [33]. We addressed this next.

3.4 \(\text{MARS}_{P_{2H}}\) retains performance amongst rule-based predictions↩︎

To restore MARS’ option to use reasoning shortcuts, we re-introduced inverse edges to MoA-net-10k. Thereafter, we tested each of \(\text{MARS}_{P_{2H}}\), \(\text{MARS}_{\text{naive}}\), and PoLo again (Figure 4-\(\epsilon\)). While each approach was optimized for standard MRR, pruned metrics indicated how well positive predictions aligned with rules. In Figure 4-\(\epsilon\), one can see that both MARS variants and PoLo achieved standard metrics on par with or better than MoA-net-10k without inverse edges (Figure 4-\(\gamma\)). 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 3.2, XSwap was used on MoA-net-10k to assess the susceptibility of \(\text{MARS}_{P_{2H}}\) to degree bias. Unlike in Section 3.2, there was now no evidence for degree bias (see Appendix 11).

3.5 External validation of \(\text{MARS}_{P_{2H}}\) on MoA-net-10k↩︎

In comparison to baseline methods, MARS’ standard metrics were on par with that of \(\text{MARS}_{naive}\) and PoLo, and sometimes outperformed by that of MINERVA (Table 2). However, \(\text{MARS}_{P_{2H}}\)’s pruned metrics outperformed both standard and pruned metrics for all baseline models. Since \(\text{MARS}_{P_{2H}}\)’s pruned metric intervals overlapped with reported intervals for MINERVA, PoLo, and \(\text{MARS}_{naive}\), significance was verified based on Mann-Whitney U tests [69] (\(\alpha=0.05\), one-sided) with Holm-Bonferroni correction for multiple testing [70], and \(p\)-values are reported in Appendix 12.

Since MINERVA does not, by design, utilize rules for guidance, it suffers the same reasoning shortcuts as PoLo and \(\text{MARS}_{naive}\). This is evident by its low pruned metrics reported in Table 2. 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 2.1, 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. Based on the sequence of drugs, proteins, and BPs in an MoA, \(\text{MARS}_{P_{2H}}\) recovered the correct MoA for all 33 pairs. Notably, \(\text{MARS}_{P_{2H}}\) correctly predicted MoAs which consisted of three and four hops, including that of norethisterone (Eq. 16 ) to treat heavy menstrual bleeding, and dexamethasone acetate (Eq. 17 ) to treat a cancer known as Multiple myeloma [13]. In contrast, MINERVA failed to connect the drugs and BPs for these longer-range MoAs. 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.

\[\begin{align} \textsl{Norethisterone} &\xrightarrow{upregulates} \textsl{Progesterone receptor} \notag \\ &\xrightarrow{interacts} \textsl{Gonadotropin releasing hormone} \notag \\ &\xrightarrow{participates} \textsl{Luteinizing hormone secretion} & \label{norethisterone} \end{align}\tag{16}\]

\[\begin{align} \textsl{Dexamethasone Acetate} &\xrightarrow{upregulates} \textsl{GC receptor} \xrightarrow{interacts} \textsl{Annexin A1} \notag \\ &\xrightarrow{interacts} \textsl{Phospholipase} \xrightarrow{participates} \textsl{Leukotriene biosynthesis} & \label{dex95acetate} \end{align}\tag{17}\]

4 Discussion↩︎

NeSy approaches are sometimes portrayed as more trustworthy than their black-box counterparts, partially due to increased interpretability [14], [52]. Here, we presented a NeSy RL approach, \(\text{MARS}_{P_{2H}}\), which promotes interpretability by deconvoluting drug MoAs. Specifically, through our novel algorithm, two-hop joint probabilities (\(P_{2H}\)), MARS learned weights corresponding to rules representing MoA patterns; each weight served 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’ interpretability called the trustworthiness of such approaches to question.

To address this, we considered Marconato et al.’s [33] 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 introduced in Section 1, 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 generalizable mitigation strategy.

4.1 Limitations and prospective directions↩︎

MARS also has several limitations. Firstly, we note that the cost effectiveness achieved by \(\text{MARS}_{P_{2H}}\) was limited to 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 [71], to make more informed choices about edge trimming or masking. Notably, a recent study [72] developed a NeSy RL method, REx, which attenuates structural KG biases by incorporating information content (IC), which encodes information about node degree and type, directly into the reward function. The integration of IC into MARS’ reward function could serve as an alternative mitigation strategy for node degree bias. Importantly, it could, potentially, reduce the need for KG pruning beforehand, thereby improving the generalizability of MARS to KGs which are denser or have imbalanced node and edge type distributions.

Additionally, as mentioned in Section 2.5.2, the \(P_{2H}\) updates assume independence between fragments. While this assumption was made to reduce computational complexity, it does not necessarily reflect the reality of the biology, in which such molecular interactions are not typically independent. Therefore, rewarding partial metapath matches without such an assumption is a key future direction. In addition to addressing methodological limitations, prospective studies could explore biological improvements. For example, as the DrugMechDB [13] continues to curate MoAs, future studies could utilize more ground-truth MoAs than the 33 used here. This could also, potentially, allow the incorporation of more complex MoAs. Alternatively, other future directions might involve the inclusion of binding or expression values, or more specific protein subclasses, which could help to better understand the patterns found here. For instance, the weights reported in Figure 5 show that paths involving a \(\texttt{downregulates}\) relation were often weighted higher than paths of equal length with an \(\texttt{upregulates}\) relation. Including such information could help one understand whether this occurs due to a structural bias in the KG or a potential biological explanation.

In summary, 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.

5 Conclusion↩︎

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. However, since shortcut mitigation was achieved upon a pruned, task-specific KG, a key future direction includes the generalization of MARS upon larger, denser KGs.

Declarations↩︎

5.1 Ethics approval and consent to participate↩︎

Not applicable.

5.2 Consent for publication↩︎

Not applicable.

5.3 Availability of data and materials↩︎

The MoA-net KG used in this study is publicly available in the MoA-net repository: https://github.com/laurendelong21/MoA-Net[23].

MARS code is publicly available in the MARS repository: https://github.com/laurendelong21/MARS[12].

5.4 Competing interests↩︎

YG and DDF were employees of Enveda during the course of this work and have real or potential interest in the company. LND was a temporary employee of Lyzeum Ltd but began employment after the course of this work. PG and LND are current employees of Cancer Research UK (CRUK), but CRUK is not affiliated or associated with this work. All other authors declare no competing interests.

5.5 Funding↩︎

During completion of this study LND was funded by the University of Edinburgh Informatics Graduate School through the Global Informatics Scholarship as well as the Alan Turing Institute Enrichment Scheme. LND partially conducted this work during an internship at Enveda. JDF and PG were 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.

5.6 Authors’ 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.

5.7 Acknowledgements↩︎

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

References↩︎

---

1. https://github.com/laurendelong21/MARS↩︎

2. https://github.com/laurendelong21/MoA-Net↩︎

3. https://github.com/laurendelong21/MARS↩︎