PAMol: Pocket-Aware Drug Design Method with Hypergraph Representation of Protein Pocket Structure and Feature Fusion

1 Introduction

Drug design aims to efficiently generate molecules that have both significant potential for clinical application and precise treatment of disease [Wong et al., 2024]. It relies on in-depth analysis of structures and biochemical properties of existing drugs or target proteins. Traditional drug design is a complex process with high costs, long cycles and high risks. It costs about 2.5 billion dollars to design a new drug, and the development process can take up to 10 to 15 years [Bano et al. ]. The chemical space for drugs ranges between approximately 1023 and 1060 molecules [Medina-Franco and L´opez-L´opez, 2024]. There will be more than 1015 kinds of diverse and new compounds that can be synthesized [Sadybekov and Katritch, 2023]. In such a large, discrete and disorganised chemistry space, it is a very difficult task to find molecules that interact with disease targets and conform to specific physicochemical properties. The application of deep learning in the field of drug design has received increasing attention [Zhang and Chen, 2022]. Compared with traditional methods, deep learning can learn molecular and protein features from massive data, accelerating the drug discovery process. Currently, drug design methods are usually divided into ligand-based and structure-based methods.

∗Corresponding author

†Corresponding author

Ligand-based methods are based on the fact that compounds with the same physicochemical properties or structures should have the same activity or similar targets [Fenglei et al., 2021]. [Wang et al., 2021] proposed a generation model that satisfies multiple constraints by combining the knowledge distillation, conditional Transformer and reinforcement learning. [Iwata et al., 2023] combined variation graph autoencoder and Monte Carlo Tree Search to capture structural features of molecules. [Mao et al., 2023] proposed a novel data-driven self-supervised pre-trained model to generate molecules, which extends the SMILES molecule generation space to optimize the generated molecules from a chemical semantic perspective. These methods have limitations in prediction accuracy and reliability due to ignoring the structural information of the target protein.

Structure-based methods are currently dominated by protein pocket-based drug design [Zhang et al., 2024], which relies on known structures of protein pockets. Currently, the structure of protein pockets is mainly represented in the form of graphs. [Peng et al., 2022] used the graph neural network (GNN) to capture the spatial relations of binding pockets, and generated molecules that satisfy geometric and chemical constraints. [Guan et al., 2023] represented protein pockets as sets of atomic points in 3D space, and used GNN to generate target-aware molecules in continuous space. [Zhang et al., 2023] proposed a fragment-based generation framework that encodes contextual information and used GNN to generate molecules. [Zhang and Liu, 2023] captured the interactions between sub-pockets and molecular motifs by learning sub-pocket prototypes, and constructed a global interaction graph to generate molecules. [Lin et al., 2024] represented pocket amino acids and molecular functional groups as fragments to generate molecules. [Qian et al., 2024] introduced a scoring function for binding affinity to generate molecules that bind with high affinity to specific targets. [Huang et al., 2024b] incorporated protein-ligand interaction priors to generate molecules with high affinity. [Huang et al., 2024a] proposed a pocket-based molecular diffusion model that incorporates protein pocket information to generate drug-like molecules. The limitations of structure-based methods are mainly in two aspects: (1) 3D structures of target proteins are often difficult to be obtained directly by experimental methods, which are demanding and time-consuming in terms of computational resources. (2) Even though some proteins have been experimentally resolved, their critical pocket structure information may not be fully annotated or compiled into protein databases, which restricts their applications.

These methods have made great progress in generating molecules, but there are still challenges. Ligand-based methods are limited by the available chemical space, making it difficult to generate molecules with novel structures. Structurebased methods require in-depth knowledge of protein structures, but obtaining 3D structures is expensive, and still has much to be explored. Although the protein structure prediction can be performed using AlphaFold [Jumper et al., 2021], the accuracy of the prediction cannot be fully guaranteed. In addition, most methods neglect the high-order information in the structure of protein pockets, which leads to an incomplete understanding of the properties of protein pockets. This limitation may restrict the model's ability to provide a comprehensive understanding of complex biological systems.

To address these issues, we propose a pocket-aware drug design framework (PAMol) to generate molecules, which constructs the hypergraph of protein pockets to represent the spatial structure. It can capture high-order relations and neighborhood information of protein pockets. This framework also fuses multi-modal embeddings from proteins and molecules, including the structure and sequence of protein pockets, fingerprint features and physicochemical properties of molecules. A multi-level cross fusion module integrates the structure and sequence of protein pockets to obtain fused features, which contain high-order structural information. The fused features serve as constraints for the conditional molecule generation module, helping to improve the quality of generated molecules for specific targets. In addition, the fused features of protein pockets and molecules provide more comprehensive information for the supervised discriminator, enabling it to optimize the quality of the generated molecules. The contribution of this work can be concluded as follows:

• To capture the high-order structural information in protein pockets, we proposed a pocket-aware drug design method with hypergraph representation of protein pocket structure. The proposed method helps to improve the performance of targeted drug design.
• We fused different modal embeddings from structure and sequence of proteins, fingerprint and physicochemical properties of molecules. We also developed a conditional molecule generation module that incorporates an unsupervised discriminator and a supervised discriminator. It uses the fused features of pockets that contain high-order structural information, as constraints to guide and optimize the process of molecule generation.
• We demonstrated the effectiveness of PAMol on CrossDocked dataset. PAMol outperforms related state-of-

the-art methods in terms of vina score, QED, Lipinski, diversity, and time, showing the feasibility for targeted drug design.

2 Methods

Figure 1 shows the framework of PAMol model. First, for a given protein pocket, the spacial structure and sequence are represented by HGNN and ProteinBERT, respectively, as shown in Figure 1 (a). Hypergraph structure and sequence features of protein pockets are fused by multi-level cross fusion module, as shown in Figure 1 (b). For a given molecule, fingerprint features and physicochemical properties are represented separately. Figure 1 (c) illustrates the process of obtaining and fusing these molecular features. Finally, the Conditional Molecule Generation module (Figure 1 (d)) uses fused features of pockets as conditions to guide the molecule generation.

2.1 Hypergraph Representation of Protein Pocket Structure

The spatial structure of proteins is the basis of their functions. We construct a hypergraph that contains the structural hyperedges of multiple amino acids. It can reflect the spatial relation between amino acids in protein pocket and thus represent the higher-order structural information of protein pocket more accurately.

The hypergraph of a protein pocket can be defined as G = (V,E,W), where V = {v1, v2, . . . , vn} denotes the set of nodes, with each node representing an amino acid in the protein pocket. E = {e1, e2, . . . , em} denotes the spatial structure hyperedge set. Suppose the coordinates of the central carbon atoms of the i-th and j-th amino acid are (xi, yi, zi) and (xj , yj , zj), respectively. The distance between them can be calculated by:

where j ∈ {1, 2, . . . , n} and i ̸= j. If Dij < 5˚ A, the j-th amino acid is added to the hyperedge ei.

In a hypergraph G, each hyperedge ei ∈ E is assigned a weight w(ei) that indicates the importance of its connectivity relations within the entire hypergraph. These weights are organized into a diagonal matrix W, defined as follows:

where diag(W) denotes the diagonal of matrix. Each diagonal element w(ei) corresponds to the weight of the i-th hyperedge ei in the hyperedge set E. This represents the individual importance of each hyperedge in the hypergraph.

To specifically describe the relation between nodes and hyperedges, the hypergraph of a protein pocket can be further represented as an association matrix Hp ∈ {0, 1}|V |×|E| , which is defined as:

where Hp(v, e) = 1 indicates that node v is a member of hyperedge e. Conversely, Hp(v, e) = 0 indicates that node v does not belong to hyperedge e.

To learn the high-order structural information in protein pockets, we apply hypergraph neural network [Feng et al., 2019] to encode the protein pocket hypergraph through its unique hypergraph convolution layer. The hypergraph of protein pocket is first mapped to a feature matrix Xh ∈ R|V |∗|E| . In the hypergraph convolution process, the representation of each node is updated based on the features of its connected hyperedges and neighboring nodes. The hypergraph convolution operation is defined as:

where λv ∈ R|V |∗|V | and λe ∈ R|E|∗|E| denote the diagonal matrices of the node degree and hyperedge degree, respectively. Hp is the association matrix. W is the weight matrix of hyperedges for protein pockets. Θ(l) is the trainable parameter matrix at layer l, which is used to capture the complex structural and attribute relations of protein pockets. P (l) h is the feature matrix of protein pocket nodes at layer l with P (0) h = Xh. This initial feature matrix Xh passes through two convolutional layers and one average pooling layer to obtain a final representation Ph ∈ R|V |∗ds , where ds =768.

By multi-layer aggregation and propagating mechanisms, hypergraph neural network is able to learn high-order relations within the protein pocket, capturing both local and global structural information.

2.2 Embedding Representation of Protein Pocket Sequence

The protein pocket sequences can map discrete amino acid sequences into a low-dimensional continuous vector space, generating embedding vectors. These embedding vectors are used as feature representations of pocket sequences and can capture key information in the sequence, such as the relative positions between amino acids, thus improving the performance of the model in drug design tasks. A pre-trained ProteinBERT model [Rao et al., 2019] is used to obtain the sequence features. First, the protein pocket sequences are encoded by IUPAC Tokenizer to get the token sequences.

The resulting token sequences are passed through a 12-layer Transformer with a hidden layer size of 512 units and 8 attention heads, to obtain the final representation Ps ∈ Rnp∗ds , where np is the number of amino acids in the protein pocket and ds is the embedding dimension of 768.

2.3 Multi-Level Cross Fusion

In this work, we fuse features of different modalities and different scales, improving the model's capability to represent features. Figure 1 (b) shows the process of feature fusion across the structure and sequence of protein pockets based on Multi-level Cross Fusion (MCF).

Before cross fusion operation, a fully connected layer maps the feature vectors into a unified embedding space. Specifically, Ph is the feature of protein pocket hypergraph structures, and Ps is the feature of protein pocket sequences. Then, the feature vectors are represented as h and s, corresponding to Ph and Ps, respectively.

where w is trainable weight, and b is bias, respectively.

Multi-Scale Feature Fusion. This block first constructs a cross matrix by the cross-product operation [Chen et al., 2021]. Let the structure representation and sequence representation of the protein pocket after transformation be denoted as vectors hi = [hi1, hi2, . . . , hid] and si = [si1, si2, . . . , sid], respectively, where d=768. hi and si denote the i-th row in vectors h and s. The cross matrix Ci ∈ Rd×d represents the interaction between hi and si, which is defined as:

Then, it extracts local and global features from the cross matrix at different scales for comprehensive understanding. The CNNmodel incorporates a pooling layer to capture localized interactive patterns, denoted as feature Flocal:

The flatten operation on the cross matrix Ci allows learning global features.Fglobal = Linear(flatten(Ci))

Flocal and Fglobal are passed through a MLP to obtain the final multi-scale fused feature representation Flg of protein pocket.

Scalar-Based Multi-Feature Fusion. First, the feature interaction between hi and si obtained from structure and sequence representations of the protein pocket is encoded by an element-wise product operation. Then, the element-wise vector is passed through MLP to obtain the scalar fusion feature Fsc, which is defined as:

Flg and Fsc are concatenated to obtain the final representation Fp with the embedding dimension of 512, which fuses the hypergraph structure features and sequence features of the protein pocket.Fp = MLP (Concat(Flg , Fsc)) (9)

2.4 Embedding Representation of Molecule

This module can extract fingerprint features and physicochemical properties of molecules [Kotsias et al., 2020], which is useful for optimizing the performance of generated molecules. The embedding vectors of structural features are obtained by Morgan fingerprint, which are generated by considering the topology structure of molecules. In the generation process, each atom and its neighboring atoms are iteratively considered until a predetermined radius is reached. This iterative process captures the connection patterns and distances between atoms, thereby reflecting the structure of molecules. Based on this structural data, we obtain a one-dimensional vector, denoted as Ms. The physicochemical properties obtained by RDKit [Landrum and others, 2013], including Octanol-Water Partition coefficient (LogP), Topological Polarity Surface Area (TPSA), Molecular Weight (MW), Number of Hydrogen Bond Acceptors (HBA) and Number of Hydrogen Bond Donors (HBD). The embedding representation of physicochemical properties is denoted as a one-dimensional vector Mp. The fingerprint features Ms and physicochemical features Mp are concatenated to obtain Fm with the embedding dimension of 512.

2.5 Conditional Molecule Generation (CMG) Module

The Conditional Molecule Generation (CMG) module uses the fusion features Fp of protein pockets that contain highorder structural information as generative conditions, to facilitate the discovery of potential drug molecules against specific targets. CMG module includes two discriminators, one generator, and one decoder, as shown in Figure 1 (d).

First, the 512-dimensional noise vector z sampled from a normal distribution [Chen et al., 2023] is concatenated with the fused features Fp of the protein pocket, serving as the input to the generator. It can enhance the expressive power of the generator, enabling it to produce molecules that are both diverse and conform to specific biological properties. Through the processing of a 3-layer neural network, the latent feature vectors Mf of the molecules are obtained.

In the process of molecule generation, to more accurately optimize the fit between the generated molecules and the protein pockets, the real molecule features Fm and the protein pocket features Fp are concatenated as fused feature Fmp:

We also fuse the latent feature vector Mf of generated molecules and protein pocket features Fp, obtaining the fusion feature Ffp:

Then, an unsupervised discriminator captures global features of real molecules to refine the overall quality of the generated ones. The features Fm of real molecules and the latent feature vectors Mf of generated molecules are passed through an unsupervised discriminator, to compute the unsupervised loss. This loss measures the difference or similarity between the features of generated molecules and real molecules. The quality of generated molecules is optimized in the back propagation, which motivates the generator to produce molecules that are closer to the real data distribution.

Meanwhile, a supervised discriminator, by focusing on protein pocket features including the high-order structural information, learns and optimizes relevant features of generated molecules to ensure compatibility with those pockets. The fusion features Fmp and Ffp are served as the inputs to the supervised discriminator. During the back propagation process, the discriminator learns the relevant features between the generated molecules and protein pockets. By calculating the supervised loss, it performs back propagation and updates its internal neural network parameters based on the loss gradients. This optimization process aims to improve the quality of the generated molecules for targeting the protein pockets.

An integrated loss function combines the discriminatory capabilities of the two discriminators, enabling the generative model to integrate the information of molecules and protein pockets during the training process. The loss is defined as:

where Ereal represents the expectation over the real data distribution. Efake represents the expectation over the generated data distribution. gp is a gradient penalty. λgp represents the weight coefficient of the gradient penalty. This loss function effectively balances the learning process of two discriminators, thus improving the overall performance of the model.

3 Experiments and Results

3.1 Dataset and Preprocessing

We used the same dataset CrossDocked with [Luo et al., 2021]. This dataset removes the protein-ligand pairs with a binding pose RMSD of less than 1 ˚A, leading to a total of 183,468 pairs. To avoid overlap between training and test sets, [Luo et al., 2021] first clustered the data based on protein sequence similarities. Then, 100,000 protein-ligand pairs were randomly selected from the clustered data for training. For the test set, 100 proteins were randomly selected from the remaining clusters, ensuring no overlap with the training set. In this work, to construct the structural hypergraph of protein pockets, we parsed the PDB files containing information about protein pockets. The files that could not be parsed were removed from both the training and test sets. Finally, we obtained 53,268 protein-ligand pairs for training and 47 protein-ligand pairs for testing.

3.2 Implementation Details

PAMol has been performed based on Python 3.8, Tensorflow and Keras. The hardware setup consisted of an NVIDIA GeForce RTX 3090 with CUDA and cuDNN. We set the training process to run for 2000 epochs, with a batch size of 64, and a low learning rate of 0.0001 to ensure smooth and stable convergence. Adam optimizer was used to optimize the training process of PAMol, ensuring an efficient adjustment of learning rates and a well-behaved convergence.

3.3 EvaluationMetrics

We evaluated the performance of PAMol, with common metrics [Luo et al., 2021; Polykovskiy et al., 2020] including: (1) Vina Score, it estimates the binding affinity between the ligand and target protein, which is a crucial measure to evaluate how well the generated molecule fits into the target protein pocket. (2) High Affinity, it represents the percentage of molecules whose Vina Score is higher than that of the ground truth molecule in the test set. (3) QED, it evaluates the drug-likeness of a molecule by combining multiple desirable molecular properties. (4) SA (Synthetic Accessibility), it measures the synthetic difficulty of the molecule. (5) LogP, it indicates the octanol-water partition coefficient, which should be between -0.4 and 5.6 [Ghose et al., 1999] for a good drug candidate. (6) Lipinski, it measures how well the molecule complies with Lipinski's five rules. (7) Diversity, it quantifies the average pairwise Tanimoto dissimilarity of the generated molecules for each target pocket. (8) Time, it represents the average time required to generate 100 samples for each pocket across all targets.

3.4 Ablation Study

The ablation studies were performed with different feature combinations, including the molecule features (MolF), the pocket sequence features (SeqF), the pocket structure features (StruF), and the pocket fused features (CrossF), to investigate their impacts. As shown in Table 1, MolF+CrossF achieves the best vina score of -7.646 among all models, indicating that it can generate molecules with higher binding affinity. MolF+CrossF performs second-best on the high affinity (0.840), with only a 0.001 gap from MolF+SeqF. In addition, MolF+CrossF obtains the highest score on both QED (0.778) and diversity (0.823), and has a logP value of 3.149 within the acceptable range, which indicates that it can improve the drug-likeness and diversity of generated molecules. SA of MolF+CrossF is higher than that of MolF+SeqF, but lower than that of two other combinations. MolF+CrossF and MolF+SeqF both score 5.000, complying with Lipinski's Rule of Five. Despite MolF+CrossF slightly longer runtime compared to the second best model, it is a promising model due to its superior performance on several key metrics.

To evaluate the global view of the chemical structures of generated molecules, Figure 2 shows the proportion of different ring sizes among molecules generated by PAMol. We notice that PAMol tends to generate molecules containing relatively more stable rings (5-atom ring and 6-atom ring), and few unstable rings. This is consistent with the regular principles of drug design, indicating that PAMol has high effectiveness in generating molecules with drug potential. In addition, we present the distribution of the nine common covalent bonds in the dataset and the generated molecules, including C-C, C=C, C-O, C=O, C-N, C-F, C=N, C-Cl and C-S. As shown in Figure 3, for all nine covalent bonds, the bond distributions of molecules generated by MolF+CrossF (PAMol) are closer to those of dataset (CrossDocked2020).

3.5 Performance Comparison and Analysis

Table 2 shows the comparison results of different models. PAMol outperforms other models on general metrics except vina score and SA. PAMol can generate the molecules with high affinity (0.840), which has improved by 6.06% compared to the second-best method. PAMol improves by 26.92% over the second-best method in QED, indicating significantly enhanced drug-likeness of generated molecules.

Figure 4: Docking results of the generated molecules with targets KRAS and EGFR.

The logP value (3.149) of PAMol within the acceptable range (-0.4 to 5.6) indicates that the generated molecules are more potential as drug candidates. PAMol scores 5.000 on the Lipinski criteria, indicating that generated molecules meet all the conditions of the Rule of Five. PAMol improves diversity by at least 11.22%, showing its ability to generate diverse and novel molecular structures. As in Table 2, vina score and SA of PAMol are not optimal among all models but higher than those in the test set, indicated that PAMol has certain potential. In addition, PAMol achieves a significant reduction in time, thereby enhancing efficiency and saving time costs in molecule generation, and showing competitiveness.

3.6 Case Study

To further validate the capability of PAMol model to generate molecules, the case experiments were conducted. We selected two key targets in pancreatic cancer, KRAS (UniProt ID: P01116, PDB ID: 8onv) and EGFR (UniProt ID: P00533, PDB ID: 8a27). Figure 4 shows the docking results of the generated molecule Mol K1 with target KRAS and the generated molecule Mol E1 with target EGFR. Light cyan indicates carbon, sulfur, hydrogen, and fluorine atoms. Black indicates nitrogen atoms, and yellow indicates oxygen atoms. It can be seen that the generated molecules bind to the target at distances between 2.4 ˚ A and 3.1 ˚ A, with the promising vina scores. It indicates that there are the strong interactions between the molecules and the targets. In addition, it is clear that the shapes of the generated molecules are ideally suited to the shapes of the active pockets.

To evaluate physicochemical properties of targeted molecules generated by PAMol, we used ADMET [Fu et al., 2024] to obtain radar charts of some basic properties, including MW, nRig, fChar, nHet, MRing, nRing, nRot, TPSA, nHD, nHA, logD, logS and logP, as shown in Figure 5. The green line represents the physicochemical property scores of the generated molecule, the blue outline indicates the upper limit, and the red outline indicates the lower limit. It can be seen that four targeted molecules (Mol K1, Mol K2, Mol E1 and Mol E2) meet most of the physicochemical property standards.

4 Conclusion

This paper proposes a pocket-aware drug design framework, namely PAMol, which captures the high-order structural information of protein pockets to generate molecules for specific targets. We constructed the hypergraph to represent the intricate spatial structure of protein pockets, aiming to capture high-order relations among residues and detailed neighborhood information that reflects the local environment within the pocket. We also fused the cross-modal embeddings from protein pockets and molecules to guide and optimize the process of molecule generation. In addition, we designed a Conditional Molecule Generation (CMG) module that focuses on the features of protein pockets including the high-order structural information. It learns the latent features and distribution patterns of molecules through an unsupervised discriminator, and uses a supervised discriminator to optimize relevant features of generated molecules to ensure compatibility with specific pockets. Experiments show that PAMol can efficiently generate molecules for specific targets. In the future, we will consider incorporating the molecule hypergraph and further optimize the quality of generated molecules.

Acknowledgments

The authors thank the members of Machine Learning and Artificial Intelligence Laboratory, School of Computer Science and Technology, Wuhan University of Science and Technology, for their helpful discussion within seminars. This work was supported by Hubei Province Natural Science Foundation of China (No.2024AFB865) and National Natural Science Foundation of China (No.61972299, 62372342).