==== Front Bioinformatics Bioinformatics bioinformatics Bioinformatics 1367-4803 1367-4811 Oxford University Press 37387164 10.1093/bioinformatics/btad207 btad207 Systems Biology and Networks AcademicSubjects/SCI01060 ArkDTA: attention regularization guided by non-covalent interactions for explainable drug–target binding affinity prediction https://orcid.org/0000-0002-6458-7723 Gim Mogan Department of Computer Science and Engineering, Korea University, Seoul 02841, Republic of Korea https://orcid.org/0000-0001-9548-7146 Choe Junseok Department of Computer Science and Engineering, Korea University, Seoul 02841, Republic of Korea https://orcid.org/0000-0002-3231-478X Baek Seungheun Department of Computer Science and Engineering, Korea University, Seoul 02841, Republic of Korea https://orcid.org/0000-0002-0613-120X Park Jueon Department of Computer Science and Engineering, Korea University, Seoul 02841, Republic of Korea https://orcid.org/0000-0003-3778-6177 Lee Chaeeun Department of Computer Science and Engineering, Korea University, Seoul 02841, Republic of Korea https://orcid.org/0009-0003-4048-4089 Ju Minjae LG CNS, AI Research Center, Seoul 07795, Republic of Korea Lee Sumin LG AI Research, Seoul 07795, Republic of Korea https://orcid.org/0000-0001-6798-9106 Kang Jaewoo Department of Computer Science and Engineering, Korea University, Seoul 02841, Republic of Korea AIGEN Sciences, Seoul 04778, Republic of Korea Corresponding author. Department of Computer Science, Korea University, Seoul, 02841, Republic of Korea. E-mail: kangj@korea.ac.kr (J.K.) This work was done while the author was a graduate student at Korea University Computer Science Department. 6 2023 30 6 2023 30 6 2023 39 Suppl 1 ISMB/ECCB 2023 Proceedings i448i457 © The Author(s) 2023. Published by Oxford University Press. 2023 https://creativecommons.org/licenses/by/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract Motivation Protein–ligand binding affinity prediction is a central task in drug design and development. Cross-modal attention mechanism has recently become a core component of many deep learning models due to its potential to improve model explainability. Non-covalent interactions (NCIs), one of the most critical domain knowledge in binding affinity prediction task, should be incorporated into protein–ligand attention mechanism for more explainable deep drug–target interaction models. We propose ArkDTA, a novel deep neural architecture for explainable binding affinity prediction guided by NCIs. Results Experimental results show that ArkDTA achieves predictive performance comparable to current state-of-the-art models while significantly improving model explainability. Qualitative investigation into our novel attention mechanism reveals that ArkDTA can identify potential regions for NCIs between candidate drug compounds and target proteins, as well as guiding internal operations of the model in a more interpretable and domain-aware manner. Availability ArkDTA is available at https://github.com/dmis-lab/ArkDTA Contact kangj@korea.ac.kr National Research Foundation of Korea 10.13039/501100003725 NRF-2023R1A2C3004176 Ministry of Health & Welfare, Republic of Korea HR20C0021(3) Ministry of Science and ICT 10.13039/501100014188 Institute for Information & communications Technology Planning & Evaluation ==== Body pmc1 Introduction Identification of drug–target interactions (DTIs) is a central task in drug design and development. Due to the costly and labor-intensive nature of traditional drug development process based on in vivo and in vitro experiments, deep learning models for protein–ligand binding affinity prediction have gained recognition (Chen et al. 2018). However, limited model explainability remains an obstacle to the adoption of such models by domain experts (Preuer et al. 2019). With the unique ethical and regulatory requirements, there is a growing demand for interpretable deep models in the field of biomedicine. In recent works, attention-based methods were studied to address the issue of explainability (Liang et al. 2021). Critical domain knowledge should be integrated to ensure that the model’s implicit assumptions are compatible with expert opinions (Dash et al. 2022). In DTI, one such key concept is that of protein–ligand non-covalent interactions (NCIs). NCIs are essential for understanding how proteins and ligands interact and form complexes with each other, which affects the mechanism of action for drug compounds (Tang et al. 2017; Chen et al. 2019; Anighoro 2020; Aljoundi et al. 2020). Most drug compounds are small organic molecules that act as ligands and interact with proteins to carry out their functions. The majority of drugs deliver their effects by forming noncovalent bonds with their biological targets. NCIs induce conformational changes in target proteins which influences the overall binding affinity. This is crucial for the stabilization of the protein–ligand complex in its final form (Davis and Phipps 2017; Aljoundi et al. 2020). Despite being highlighted as a fundamental concept in protein–ligand affinity prediction task, few studies have addressed the importance of protein–ligand NCIs. While MONN (Li et al. 2020) explicitly utilized NCIs in its auxiliary task, the resulting pairwise interaction matrix between all protein residues and all ligand atoms is limited in its capacity to differentiate active and inactive binding sites. On the other hand, AttentionDTA sought to distinguish active and inactive residues using parameterized weights in their attention mechanism without explicitly using NCI labels. However, its attention mechanism is based on the convoluted features from each of its respective modality-wise encoder modules. We addressed the limitations of both previous works by utilizing NCI labels in the attention mechanism to identify active protein residues (Zhao et al. 2022a). We present ArkDTA, an explainable deep DTI prediction model with NCI-aware attention regularization. Taking as input a set of protein residues and a set of ligand substructures, our novel regularization method modulates the distribution of cross-modal attention weights from the residues to chemical substructures in a manner that allows a distinction between active and inactive residues. Our modified cross-modal attention module appends a pseudo-substructure embedding to the set of key chemical substructures and focuses the attention on the pseudo embedding where the query protein residue is inactive. Examining the final attention weights yields qualitative insights into the model’s internal operations. Experimental results on three benchmark datasets reveal that ArkDTA achieves predictive performance comparable to the current state-of-the-art models while significantly improving model explainability. Qualitative investigation into the attention maps demonstrates our model’s ability to identify NCI-forming regions in seen and unseen protein–ligand complexes as well as highlight chemical substructures commonly used as pharmaceutical agents. 2 Materials and methods 2.1 Dataset Three different benchmark datasets were used in this study which are PDBbind version 2020 (PDBbind), Davis et al. (Davis), Metz et al. (Metz) to conduct experiments on ArkDTA and baseline models (Davis et al. 2011; Metz et al. 2011; Liu et al. 2017). The ith data instance Xi in each of these datasets consists of a protein–ligand pair with its binding affinity score, expressed in one of the following measurement types: inhibition constant (Ki), dissociation constant (Kd) and inhibitory concentration 50 (IC50). For the purpose of this study, the PDBbind was sub-divided into two subsets according to the binding affinity measurement type. The KIKD subset consists of all protein–ligand instances whose binding affinity scores are expressed as Ki or Kd value, and the remaining instances whose affinity scores are expressed as IC50 value were combined to form the IC50 subset. The Davis and Metz dataset contains protein–ligand pairs with only KIKD-based affinity scores. We applied several data curation methods to our constructed datasets. For each dataset, protein–ligand pairs whose number of amino acids in the protein sequence exceeds 1000 or the exact affinity value is unavailable (e.g. expressed as inequality “>50000 M”) were excluded. We then normalized binding affinity scores in each dataset into values in unit “M” and subsequently transformed them into log space for consistent comparison (Öztürk et al. 2018). Table 1 shows the total number of proteins, ligands, and curated data instances in two measurement types (KIKD, IC50) for each dataset. Overall, each ith data instance Xi in the binding affinity datasets is defined as the following, where pi, ci, yi∈R are input protein, ligand, and annotated binding affinity value, respectively. Table 1 shows the statistics of each dataset. (1) Xi=(pi,ci,yi) Table 1. Statistics of the DTI datasets used in experiments. Dataset PDBbind Davis et al. Metz et al. No. of proteins 10 162 311 121 No. of ligands 13 015 68 240 No. of KIKD instances 9327 21 331 13 669 No. of IC50 instances 6 593 0 0 To regularize ArkDTA’s attention mechanism, we further augmented the preprocessed PDBbind dataset with NCI labels. We used Protein–Ligand Interaction Profiler (PLIP) to extract the NCI labels from each binding complex structure contained in the original PDBbind dataset (Adasme et al. 2021). The NCI labels for each protein–ligand pair are represented as a m×n 2D binary matrix to indicate the presence of any type of NCIs (e.g. hydrogen bonding, salt bridges) where m and n are the numbers of amino acid residues and atoms in ligand, respectively. Since the attention mechanism in ArkDTA is based on cross-modal interactions between the protein residues and chemical substructures, we converted the binary matrix into a m-dimensional binary vector k→ where each residue in protein is labeled as 1 if it has NCI with at least one atom in its ligand partner. We use k→ as ground truth residue-wise NCI labels for attention regularization in ArkDTA. The augmented ith data instance X’i in the PDBbind dataset is defined as the following, where k→i∈[0,1]n is the converted binary vector indicating the presence of each residue’s NCI with the input ligand ci. (2) X’i=(pi,ci,yi,k→i) Despite having the least data instances, our preprocessed PDBbind dataset is the primary dataset of this study since it contains the residue-wise NCI labels. We randomly partitioned the PDBbind dataset into 5-folds where 5% of the training instances in each fold were used for validation. This split method yields an ensemble of five models trained on different folds of this dataset. The Davis and Metz dataset used for fine-tuning purposes was partitioned into training and test instances (8:2) where 5% of the training instances were also used for validation. The purpose of these two datasets is to fine-tune each of the five models previously trained on the PDBbind dataset. 2.2 Model architecture 2.2.1 Overview Our ArkDTA model consists of ‘Protein Encoder Module’, ‘Ligand Encoder Module’, ‘Protein-Ligand Integration Module’, and ‘Affinity Prediction Module’. The first two modules encode input data into protein residue-wise and chemical substructure-wise representations, respectively. The ‘Protein-Ligand Integration Module’ refines the residue-wise representations based on their attention weights given the substructure-wise representations and subsequently aggregates them into a single binding complex representation. Finally, the ‘Affinity Prediction Module’ takes the binding complex representation as input and predicts the binding affinity value as output. The formal definition for ArkDTA is the following, where the input are protein p and ligand (compound) c, and output is the predicted binding affinity value y^. Figure 1 shows the overall model architecture of ArkDTA. (3) y^=ArkDTA(p,c) 2.2.2 Protein Encoder Module The ‘Protein Encoder Module’ takes a protein p as input and encodes it into a set of m d-dimensional residue embeddings as output. The initial representation for input protein is its FASTA sequence. While such 1D-based representations may have limitations in representing proteins, large-scale language models have been introduced to alleviate these issues (Rives et al. 2019). These models have shown promising results in protein structure and function prediction tasks. We imported a pre-trained protein language model called Evolutionary Scale Model (ESM) and its tokenizer to obtain residue embeddings (Rives et al. 2019; Lin et al. 2022). The tokenizer converts the input protein p into a sequence of tokens subsequently fed to the ESM model. Finally, the ESM model converts the input tokens to a set of protein residue embeddings R∈Rm×d. The model version of its pretrained weights is ESM-2 (8M) where its number of layers is 6. 2.2.3 Ligand Encoder Module The ‘Ligand Encoder Module’ takes a ligand c as input and encodes it into a set of n d-dimensional substructure embeddings as output. The initial representation for input ligand is its SMILES string. The SMILES string is first converted into a Morgan fingerprint represented as a 1024-dimensional bit vector f→∈[0,1]1024. Each bit position in the vector indicates the presence of its corresponding chemical substructure. The ‘Ligand Encoder Module’ leverages this information by gathering the positional indices of that vector where its bit is 1 and uses a lookup table to obtain a set of trainable d-dimensional chemical substructure embeddings S={s1,s2…sn|si∈S} where n is the number of chemical substructures extracted from f→. S is a set of 1024 trainable chemical substructures where each of them corresponds to its bit position in the Morgan fingerprint. 2.2.4 Protein–Ligand Integration Module The ‘Protein-Ligand Integration Module’ consists of a Multihead Attention Block (MAB) and a Pooling Layer. The MAB refines an input set of protein residues R from ‘Protein Encoder Module’ based on our novel attention mechanism with another input set of chemical substructures S from ‘Ligand Encoder Module’. The Pooling Layer subsequently aggregates the refined residue-wise embeddings into one single binding complex embedding. The MAB’s operations reflect the conformational transitions proteins undergo when bound to a ligand, while the Pooling Layer’s output corresponds to the final protein–ligand complex that determines the binding affinity value. The MAB in ‘Protein-Ligand Integration Module’ employs multihead attention mechanism and produces a set of ‘refined’ residues given R and S as inputs (Vaswani et al. 2017; Lee et al. 2019). Following the definitions made by previous works, the MAB takes R, S, and S as queries, keys, and values, respectively. Let A∈Rm×n be the calculated attention weights between the query and key linear projections of R and S, respectively. For each ith residue (i∈{1,2,…,m}), the attention weights are distributed across all n corresponding chemical substructures where ∑j=1nAi,j=1. However, as most residues do not form NCIs with the incoming ligand, it may be undesirable to utilize all calculated attention weights. Our modified version of MAB first appends a trainable universal d-dimensional pseudo-substructure embedding p→∈Rd to current set of chemical substructure embeddings S. The main purpose is to regularize the attention between the query protein residues and key-value chemical substructures based on their NCIs. Specifically, we devised a strategy that makes the attention weights from non-binding query residues (i.e. residues having no NCIs with ligand) skewed toward the key pseudo-substructure embedding. For binding query residues (i.e. residues having NCIs with ligand), the attention weights are prevented from being skewed toward the pseudo-substructure but distributed to actual chemical substructure embeddings in an unsupervised fashion. We denote this modification as Attention Regularization based on NCIs in MAB (ARK-MAB). The ARK-MAB that takes R and S as input is mathematically expressed as follows, where R∗∈Rm×d is a set of ‘refined’ residue embeddings, ‘LayerNorm’ is layer-wise normalization method (Ba et al. 2016) and RFF is row-wise feedforward layer. MultiAttn is k-headed attention layer where X is linearly projected to query vectors and Y is linearly projected to both key and value vectors. For the calculation of attention weights, we employed Additive Attention originally proposed by Bahdanau et al. (2014) and used four attention heads. (4) S+=S∪{p→} (5) R∗=ARKMAB(S+,R) (6) ARKMAB(S+,R)=LayerNorm(H+RFF2(R)) (7) H=LayerNorm(R+RFF1(MultiAttn(R,S+,S+))) We adopted Pooling by Multihead Attention (PMA) from the Set Transformer framework (Lee et al. 2019) for the Pooling Layer. The m refined residue embeddings R∗∈Rm×d are aggregated based on a set of u trainable seed vectors U∈Ru×d into a set of u aggregated residue embeddings Ra∈Ru×d. Following the explanation in the Set Transformer paper, the PMA layer is built based on the MAB that takes U, R∗, and R∗ as queries, keys, and values, respectively. Subsequently, the aggregated residues are concatenated vector-wise and reduced to a d-dimensional binding complex embedding via a simple linear layer. The order of vector-wise concatenation is determined by the fixed order of seed vectors U. The PMA layer that takes the refined residues R∗ as input is mathematically expressed as follows, where C∈R1×d is the binding complex embedding built from vector-wise concatenation (⊕) of the aggregated residues {ria∈Ra|i=1,2,…,u}. Linear is linear layer without nonlinear activation that reduces the binding complex embedding’s expanded dimension to d where the weights and bias are Wlinear∈Rdu˙×d, blinear∈Rd, respectively. Figure 2 shows the detailed description of ARK-MAB and PMA. (8) Ra=PMA(R∗) (9) PMA(X)=MAB(U,X) (10) C=Linear([r1a⊕⋯⊕rua]) Figure 1. Overview of ArkDTA. A refers to the Protein Encoder Module, B refers to the Ligand Encoder Module, C refers to the Protein–Ligand Integration Module and D refers to the Affinity Prediction Module. y^ is the predicted binding affinity value. Figure 2. Detailed description of ARK-MAB (a) and PMA (b). (a) The input queries is a set of m residues while the key values is a set of n chemical substructures appended with a pseudo-embedding p→. The output of this sub-module is a set of m refined residues R∗∈Rm×d. (b) The input queries is a set of u seed vectors while the key values is a set of m refined residues. The output of this sub-module is a set of u aggregated residues concatenated to each other. Note that the order of vector-wise concatenation is determined by the fixed order of seed vectors. 2.2.5 Affinity Prediction Module The ‘Affinity Prediction Module’ that takes the binding complex embedding C as input for predicting the binding affinity score y^ is mathematically expressed as follows, where y^∈R. MLPscore is a two-layered MLP where the intermediate layers use Dropout and ReLU as nonlinear activation. The weights, bias in the linear layers from top to bottom are Wscore1∈Rd ×d, bscore1∈Rd, Wscore2∈Rd×1, bscore2∈R1, respectively. (11) y^=MLPscore(C) 2.2.6 Training and optimization The loss objective for training ArkDTA consists of two terms which are the main and auxiliary loss objective. The main loss objective is based on root mean squared error (RMSE) between the binding affinity predictions and each of their corresponding values. The auxiliary loss objective was specially designed to impose regularization on the attention mechanism utilized in the ‘Protein-Ligand Integration Module’ using binary cross entropy as its criterion. The batch-wise main loss objective for binding affinity prediction is mathematically expressed as follows, where Y^ is a b-sized batch of predicted binding affinities, Y is a b-sized batch of ground truth binding affinities, and MSE is mean squared loss criterion for binding affinity prediction. (12) L1(Y^,Y)=MSE(Y^,Y) (13) Y^=[y^1⋯y^b] (14) Y=[y1⋯yb] For attention regularization described in the ‘Protein-Ligand Integration Module’, let A¯+∈Rm×(n+1) be the calculated attention weight matrix averaged head-wise, given the set of m protein residue embeddings R∗∈Rm×d as queries and set of n chemical substructure embeddings S∈R(n+1)×d appended with pseudo-embedding p→ as keys. For each residue in A¯+, the summation of n attention weights corresponding to n chemical substructures is equivalent to the NCI score deemed as the predicted class probability having NCI with the ligand compound. On the contrary, the attention weight corresponding to the pseudo-substructure is deemed as the predicted class probability having no such interactions. Figure 3 illustrates how the attention mechanism in the ARK-MAB works. If the ith residue does not have any NCI with the ligand, the ARK-MAB is guided to generate attention weights A¯+∈Rm×(n+1) where ∑j=1nA¯i,j