Method for predicting protein-ligand interaction using protein-ligand interaction prediction model

Abstract

This method for predicting a protein-ligand interaction executed by an information processing system including a processor and a memory comprises the steps of: predicting a protein-ligand interaction probability by a GNN-based first model, wherein a protein graph and a ligand graph in which a protein pocket and a ligand of a target protein-ligand pair are respectively structured as individual graphs are input to the first model to output the protein-ligand interaction probability; predicting protein-ligand binding free energy by a GNN-based second or third model wherein a protein-ligand graph in which the structure of a target protein-ligand complex is structured is inputted to the second model or third model to output a protein-ligand binding free energy prediction value; and generating a protein-ligand interaction prediction value by multiplying the protein-ligand interaction probability by the protein-ligand binding free energy prediction value.

Description

Protein-ligand interaction prediction method using a protein-ligand interaction prediction model

[0001] The present invention relates to a method for predicting protein-ligand interactions, and more specifically, to a method for predicting protein-ligand interactions using a GNN-based model, a learning method for the model, and a system therefor.

[0002] The cost of drug discovery is steadily increasing, while the return on research and development (R&D) costs is declining. Predicting the binding affinity of protein-ligand complexes is crucial for drug design and discovery. Typically, drug screening involves identifying drug candidates that strongly bind to a target protein. However, experimentally measuring protein-ligand binding affinity is a time-consuming and labor-intensive process, posing a major obstacle to the drug discovery process. If the affinity between a target protein and a specific ligand could be predicted quickly and accurately, the efficiency of virtual screening could be significantly improved. Therefore, various computational methods for binding affinity prediction are being studied to optimize the drug discovery process.

[0003] In general, computational methods for predicting the interaction probability or binding affinity of a ligand for a target protein can be divided into three categories: 1) physics-based, 2) empirical, and 3) structural data-based methods.

[0004] Physics-based methods primarily rely on force-field models to faithfully calculate binding degrees of freedom. Their greatest advantage is that they can predict binding free energies for arbitrary small-molecule ligands to proteins with an average error of approximately 1 kcal / mol using state-of-the-art force-field models. However, rigorous binding free energy calculations require significant computational resources. This computational burden hinders the use of molecular dynamics (MD) and also hinders the calculation of free energies for molecules such as drugs, which require high-throughput screening.

[0005] Empirical scoring functions have been widely used in many protein-ligand docking programs and virtual screening processes. They approximate protein-ligand interactions using equations composed of various physics-based conditions, such as van der Waals interactions, solvation free energy, and electrostatic interactions. The parameters of these physics-based conditions are typically fitted to experimental data to reproduce measured binding affinity values. Due to their computational simplicity and close approximation to physics-based interactions, empirical scoring functions are still actively under development.

[0006] Recently, deep learning methods have enabled more accurate, data-driven predictions in various scientific fields. Many deep learning-based methods have been proposed for predicting protein-ligand binding affinity. In particular, a structure-based machine learning scoring function has been proposed and shows a high correlation with experimental binding affinity values.

[0007] For example, the literature (Accurate Protein-Ligand Binding Affinity Prediction Using an Ensemble of 3D-Convolutional Neural Networks by Kwon, Y.: Shin, W.-H.; Ko, J.; Lee, J. Int. J. Mol. Sci. 2020, 21, 8424.) proposed a structure-based binding affinity prediction model using 3D-CNN. The literature (Meli R, Morris GM and Biggin PC (2022) Scoring Functions for Protein-Ligand Binding Affinity Prediction Using Structure-based Deep Learning: A Review.) proposed a 3D-CNN-based model to predict binding affinity. The literature (Chao Shen, Ye Hu, Zhe Wang, Xujun Zhang, Jinping Pang, Gaoang Wang, Haiyang Zhong, Lei Xu, Dongsheng Cao, Tingjun Hou Beware of the generic machine learning-based scoring functions in structure-based virtual screening, Briefings in Bioinformatics, Volume 22, Issue 3, May 2021, bbaa070,) evaluates the effectiveness of several machine learning-based scoring functions.

[0008] However, conventional models focus on predicting the binding affinity between proteins and ligands themselves, and the limitations of the training data sets for these binding prediction models have resulted in many false positives. For virtual screening, a binding prediction model that overcomes these issues and enables more accurate and faster calculations is needed.

[0009] Accurate binding free energy also requires good initial poses. Many proteins have multiple apo conformations, and small conformational changes at the binding site can significantly affect relative binding energy calculations. Autodock-GPU enables faster protein-ligand docking, and libraries are growing. Protein-ligand docking generates multiple binding poses, and performing binding free energy calculations for all of the protein-ligand docking results is expensive. Therefore, a more cost-effective and accurate method for predicting protein-ligand interactions is needed.

[0010] According to one aspect of the present invention, the present invention aims to provide a learning method for a protein-ligand interaction prediction model with improved accuracy by solving the above problem and a prediction method using the model.

[0011] According to another aspect of the present invention, it is an object to provide an efficient protein-ligand interaction prediction model with improved accuracy and a prediction method using the same.

[0012] To solve the above problem, a protein-ligand interaction prediction method executed by an information processing device according to one aspect of the present invention is provided.

[0013] The above protein-ligand interaction prediction method is a protein-ligand interaction prediction method executed by an information processing system including a processor and a memory, comprising: a step of predicting a protein-ligand interaction probability by a first model based on a GNN, wherein a protein graph and a ligand graph, each of which structures a protein pocket and a ligand of a target protein-ligand pair as individual graphs, are input into the first model to output a protein-ligand interaction probability; a step of predicting a protein-ligand binding free energy by a second or third model based on a GNN, wherein a protein-ligand graph, each of which structures a structure of a target protein-ligand complex as a graph, is input into the second or third model to output a protein-ligand binding free energy prediction value; and a step of generating a protein-ligand interaction prediction value by multiplying the protein-ligand interaction probability by the protein-ligand binding free energy prediction value. It is preferable that the above first, second and third models are neural network models based on GNN.

[0014] According to another aspect of the present invention, an ensemble score is generated by adding a physical scoring term multiplied by a parameter to the protein-ligand interaction prediction value, and it is preferable that the physical scoring term is an output value of the AutoDock GPU model.

[0015] The input data of the first, second and third models can each be described as an undirected graph G=(V, E), where V represents a group of nodes, one node represents one atom of a protein or ligand, and E is a group of edges representing the characteristics of interactions between nodes.

[0016] The protein graph and ligand graph, which are input values ​​of the first model, have a 73-dimensional vector as a node feature and a 12-dimensional vector as an edge feature, and the protein-ligand graph, which are input values ​​for the second model and the third model, have a 73-dimensional vector as a node feature and a 24-dimensional vector as an edge feature, and the first model uses the protein binding pocket and the ligand separately as input data without considering the binding pose of the target protein-ligand pair in order to predict the interaction probability of the target protein-ligand pair.

[0017] It is desirable for the second and third models to use protein-ligand complexes with binding poses as input data to predict binding free energy.

[0018] The above first model (Akscore2_NonDock model) inputs the protein binding pocket and ligand separately without considering the binding poses of the protein and ligand to predict the protein-ligand interaction probability.

[0019] The first model for predicting protein-ligand interaction probability is trained using the structural data of each protein (pocket) and ligand of the target protein-ligand pair as input values, and the actual interaction of the protein-ligand pair as output values.

[0020] The second model takes structural data of a protein-ligand complex as input, the corresponding binding free energy and RMSD as output, and includes a graph neural network.

[0021] The above third model takes the structural data of a protein-ligand complex as input, the corresponding binding free energy as output, includes a graph neural network, and applies the predicted energy and RMSD to the loss function to perform parameter optimization.

[0022] The step of predicting the protein-ligand interaction probability by the first model is:

[0023] A step of constructing a protein graph and a ligand graph by structuring the protein pockets and ligands of a target protein and ligand pair into individual graphs;

[0024] A step of inputting the protein graph and ligand graph into an embedding layer and outputting an embedding vector;

[0025] A step of independently inputting the embedding vector of the protein graph and the embedding vector of the ligand graph into the protein GNN and the ligand GNN, respectively;

[0026] A step of pooling the output values ​​of the protein GNN and the output values ​​of the ligand GNN and then fusing them together;

[0027] A step of inputting the output value of the above fusion step into a multi-layer perceptron and outputting the output value of the binary classification; and

[0028] It includes a step of predicting the protein-ligand interaction probability by applying a sigmoid function to the output value of the above binary classification.

[0029] The first model includes a protein GNN into which a protein graph is input, a ligand GNN into which a ligand graph is input, a pooling layer for pooling the output values ​​of the protein GNN and the ligand GNN, a concatenation layer for concatenating the output values ​​of the pooling layers, and a multilayer perceptron for inputting the output values ​​of the concatenation layer and outputting the output values ​​of a binary classification. It is preferable that the protein GNN and the ligand GNN each include one or more GAT layers, and the output values ​​of the binary classification are input to a binary cross entropy (BCE) loss function to optimize the parameters of the neural network model.

[0030] The training data set of the first model includes a native data set and an augmented data set, wherein the native data set includes pairs of protein pockets and ligands of existing native protein-ligand complexes having actual binding affinities. The augmented data set includes a cross-decoy data set and a random decoy data set, wherein the augmented data set includes pairs of protein pockets and ligands generated by randomly selecting ligands other than the native ligands of the native protein-ligand complex and redocking them with the protein pockets.

[0031] The step of predicting the protein-ligand binding free energy by the second model is:

[0032] A step of constructing a protein-ligand complex graph by structuring the structure of a protein-ligand complex of a target protein-ligand pair into a graph;

[0033] A step of inputting the protein-ligand complex graph into an embedding layer and outputting a protein-ligand embedding vector;

[0034] A step of inputting the above protein-ligand embedding vector into a protein-ligand GNN;

[0035] A step of inputting the output value of the protein-ligand GNN to a pooling layer and pooling; and

[0036] It includes a step of inputting the output value of the global mean pooling layer into an energy multi-layer perceptron and an RMSD multi-layer perceptron, respectively, to output the energy and RMSD, respectively.

[0037] The second model includes a protein-ligand GNN, a pooling layer for pooling the output values ​​of the protein-ligand GNN, an energy multilayer perceptron that inputs the output values ​​of the pooling layer and outputs energy, and an RMSD multilayer perceptron that inputs the output values ​​of the pooling layer and outputs RMSD, and the protein-ligand GNN is composed of a 5-layer GATv2.

[0038] The step of predicting protein-ligand binding free energy using the third model is:

[0039] A step of constructing a protein-ligand graph by structuring the structure of a protein-ligand complex of a target protein-ligand pair into a graph;

[0040] Step of inputting the protein-ligand graph into the protein-ligand GNN;

[0041] A step of inputting the output of the protein-ligand GNN to a pooling layer and pooling; and

[0042] It includes a step of inputting the output of the above pooling layer into a multi-layer perceptron to output binding free energy.

[0043] The third model includes a protein-ligand GNN, a pooling layer for pooling the output values ​​of the protein-ligand GNN, and a multilayer perceptron into which the output values ​​of the pooling layer are input, and the parameters of the neural network are optimized through regression analysis.

[0044] According to one aspect of the present invention, a more accurate combination prediction is possible by solving the false positive problem of the prediction model.

[0045] According to another aspect of the present invention, the accuracy of prediction of binding of a protein-ligand pair is improved and computational resources can be utilized efficiently.

[0046] The training data set for training the second model includes a native data set of protein-ligand complexes and an augmented data set, and the augmented data set includes a structural decoy data set, a Decoys_Cross data set, and a Decoys_Random data set.

[0047] Figure 1 is a schematic diagram of a prediction system according to one embodiment of the present invention.

[0048] FIG. 2 is a graph showing the RMSD distribution in a training data set according to one embodiment of the present invention.

[0049] Figure 3 is a schematic diagram showing the overall architecture of a network according to one embodiment of the present invention.

[0050] FIG. 4 is a schematic diagram showing the structure of the Akscore2_NonDock model according to one embodiment of the present invention.

[0051] FIG. 5 is a schematic diagram showing the structure of the Akscore2_DockS model according to one embodiment of the present invention.

[0052] FIG. 6 is a schematic diagram showing the structure of the Akscore2_DockC model according to one embodiment of the present invention.

[0053] FIG. 7 is a diagram illustrating an ensemble strategy according to one embodiment of the present invention.

[0054] FIG. 8 is a diagram illustrating an ensemble strategy according to another embodiment of the present invention.

[0055] The present invention will be more clearly understood by reference to the embodiments described below in detail, taken in conjunction with the drawings attached to this specification. However, the present invention is not limited to the embodiments disclosed below, but may be implemented in various different forms. These embodiments are provided solely to ensure that the disclosure of the present invention is complete and to fully inform those skilled in the art of the scope of the invention. The present invention is defined solely by the scope of the claims.

[0056] Meanwhile, the terms used in this specification are intended to describe embodiments and are not intended to limit the present invention. In this specification, singular forms also include plural forms, unless specifically stated otherwise.

[0057] As used herein, the term "comprises" does not exclude the presence or addition of one or more other components or steps in addition to the components or steps mentioned. While terms such as "first," "second," etc. may be used to describe various components, the components should not be limited by the terms. The terms are used solely to distinguish one component from another.

[0058] Hereinafter, various embodiments of the present invention will be described by way of example with reference to the drawings.

[0059] Figure 1 is a block diagram illustrating an example of a prediction system configuration according to an embodiment, conceptually illustrating the relevant parts of the embodiment. Each component may be implemented on a single information processing device (computer) and executed on a single device, but is not limited thereto. The prediction system may be implemented on multiple information processing devices connected to each other via a communications network, and may also include operations performed on each information processing device.

[0060] The external server (20) can be connected to the prediction system (10) via a network and can provide information such as protein structure information, ligand structure information, protein-ligand complex structure information, binding affinity information according to protein-ligand complex binding sites, genetic information, intermolecular interaction information, and / or protein structure similarity information. For example, the external server (20) can be a database for prediction processing of the prediction system (10) or a server providing the same.

[0061] The prediction system (10) may include a control unit (11), a communication unit (12), an input / output interface unit (13), and a memory unit (14). The control unit (11) is a component that controls the entire prediction system (10), and may include, for example, a processing unit such as a CPU or a GPU. The control unit (11) may train models to be described later using information stored in the memory unit (14), and may also perform prediction value calculation for a new input using the learned model. Specifically, the control unit (11) may control a model that calculates a prediction of protein-ligand binding affinity. To this end, the control unit (11) may include a control program such as an operating system (OS), a program that defines various processing orders, and an internal memory for storing data. In addition, the control unit (11) may perform information processing to execute various processes by these programs, etc.

[0062] In addition, the communication unit (12) may include an interface that can be connected to a communication device such as a router connected to a communication line, etc., and may control communication between the prediction system (10) and an external server (20). The input / output interface unit (13) may be an interface connected to the input unit (15) and / or the display unit (16).

[0063] The user can communicate with the prediction system (10) through the input / output interface unit (13). For example, the display unit (16) may be a display means (e.g., a display, monitor, touch panel, etc. composed of a liquid crystal display or organic EL, etc.) that displays a display screen of an application, etc. In addition, the input unit (15) may be, for example, a key input unit, a touch panel, a control pad (e.g., a touch pad, a game pad, etc.), a mouse, a keyboard, a microphone, etc.

[0064] In addition, the memory unit (14) may be a device that stores various databases or tables, etc. For example, the memory unit may include information such as protein structure information, ligand structure information, protein-ligand complex structure information, binding affinity information according to protein-ligand complex binding sites, genetic information, intermolecular interaction information, and / or protein structure similarity information.

[0065] According to an embodiment of the present invention, three graph neural network (GNN)-based models were developed to predict interactions, binding affinity, or binding free energy.

[0066] The first model (Akscore2_NonDock model) inputs the structure of the protein binding pocket and the ligand of the protein-ligand pair separately without considering their binding poses to predict the protein-ligand interaction probability.

[0067] The second model (Akscore2_DockS model) and the third model (Akscore2_DockC model) input a protein-ligand complex with binding poses to predict binding free energy (kcal / mol) and RMSD (Root Mean Square Deviation, unit Å), or their combination, and further estimate binding affinity.

[0068] Table 1 lists the neural network training and test datasets. The training data is the general set (excluding CASF-2016) from the PDBbind v2020 dataset, along with the augmented datasets Decoys_Dock, Decoys_Cross, and Decoys_Random. The CASF2016 dataset and the LIT-PCBA dataset were used as benchmark test data.

[0069] Training DataProtein-Ligand Complex (Native)Decoys_Dock(Created)Decoys_Cross(Created)Decoys_Random(Created)PDBbind refined set(Without CASF2016)17,225900,9101,720,9581,721,583Benchmark Test DataProteinActiveInactiveCASF2016 coreset285CASF2016 docking data2855,10016,865CASF2016 screening data5727,8391,561,424LIT-PCBA15501,500139,750,100

[0070]

[0071] 1. Neural network training data sets

[0072] For neural network training, we used the PDBbind v2020 dataset, which contains protein-ligand complex structures with experimentally determined binding affinity data (Ki / Kd). Specifically, the general dataset of PDBbind v2020 was used as training data, but the core set of PDBbind 2016, which is a redundant sample, was excluded from the training data. In addition, samples that did not dock when using Autodock-GPU (Accelerating AutoDock4 with GPUs and Gradient-Based Local Search, J. Chem. Theory Comput. 2021, 10.1021 / acs.jctc.0c01006) were also excluded.

[0073] We used 4961 samples as the training data set (selected from the data set in Table 1).

[0074] To reduce computational cost, the protein was converted into a protein binding pocket, where the protein binding pocket is defined as the residues within 5.0 Å around the crystallized ligand.

[0075] Complexes whose protein pockets were not recognized by the RDKit toolkit (RDKit: Open-source cheminformatics. https: / / www.rdkit.org) were excluded. The RDKit toolkit is one of the major software development toolkits in the field of cheminformatics and is available under the BSD-3-Clause License. Details can be found at https: / / www.rdkit.org.

[0076] To ensure efficient training, we generated three different types of data sets: a native set, a conformational decoy set, and a cross-docked decoy set.

[0077] A subset of the native set and a structural decoy set were generated by redocking native ligands, which are real ligands, into the native binding pocket. To generate decoy structures, Autodock-GPU (Accelerating AutoDock4 with GPUs and Gradient-Based Local Search, J. Chem. Theory Comput. 2021, 10.1021 / acs.jctc.0c01006) was used to generate decoy structures by generating 50 poses around the native binding pose of the ligand using a grid box of 22.5Å × 22.5Å × 22.5Å, resulting in 900,910 complex poses (Decoys_Dock). However, for some targets, a sufficient number of poses were not generated, or the RMSD values ​​of all docked poses were less than 2Å. In these cases, an additional 900 poses were generated to maximize diversity. These poses were clustered using Autodock-GPU, with a maximum of 100 poses per target. Poses with RMSD values ​​less than 2 Å were classified as the native set, while poses with higher values ​​were classified as the structural decoy set (Decoys_Dock).

[0078] Figure 2 shows the RMSD distribution of docked poses in the training dataset. Figure 2(a) shows the RMSD distribution for each protein in the Decoys_Dock pose generated by Autodock GPU in the training dataset, where the vertical axis represents proteins and the horizontal axis represents RMSD. Figure 2(b) shows the RMSD distribution of the Decoys_Dock pose generated by Autodock GPU in the training dataset, where the vertical axis represents percentage and the horizontal axis represents RMSD.

[0079] We randomly selected 100 ligands (excluding the sample's own ligands) from the remaining complexes in the PDBbind dataset and redocked them, generating 1,720,958 complexes to create a cross-docking set (Decoy-Cross). Among the generated cross-docking complexes, 70,000 complex structures were randomly selected and used as the final cross-docking decoy set.

[0080] Additionally, to cover a wider chemical space, 100 ligands were randomly selected from the Enamine Discovery Diversity Set (DDS-50) (Enamine, https: / / enamine.net / ) and docked, generating 1,721,583 complexes. This is the Decoys_Random set.

[0081] 2. Model training method

[0082] Figure 3 illustrates the overall learning model structure according to one embodiment of the present invention. In the present invention, three independent neural network models were developed: Akscore2_NonDock, Akscore2_DockS, and Akscore2_DockC.

[0083] The Akscore2_NonDock model inputs protein binding pockets and ligands separately and outputs the interaction probabilities between them.

[0084] The Akscore2_DockS model takes a protein-ligand complex as input and outputs the energy (kcal / mol) and RMSD (Å), respectively, where RMSD corresponds to the native pose.

[0085] The Akscore2_DockC model also takes a protein-ligand complex as input, but outputs a single energy value (kcal / mol) with an RMSD penalty. The method of the present invention is described in detail below.

[0086] (1) Graph Construction

[0087] The inputs of all three models can be described by an undirected graph G = (V, E), where V represents a group of nodes, each representing an atom of a protein (protein pocket) or a ligand. E represents a group of edges that connect the nodes through interactions. Edges are constructed to describe the interactions between a protein and a ligand in a protein-ligand complex. There are two types of edges: covalent and noncovalent. A covalent edge is formed when two atoms in a protein or ligand share one or more electron pairs. A noncovalent bond is formed when the distance between a protein atom and a ligand atom is less than 8 Å.

[0088] The protein-ligand complex graph, which is the input data of the Akscore2_DockS and Akscore2_DockC models, has a 73-dimensional vector as a node feature and a 24-dimensional vector as an edge feature.

[0089] The graph for individual proteins and ligands, which is the input to Akscore2_NonDock, also has 73-dimensional vectors as node features, but 12-dimensional vectors as edge features, because it only has covalent bond type edges.

[0090] Initial node and edge features according to one embodiment of the present invention are described in Tables 2 and 3, respectively. Table 2 describes node features, their sizes, and a corresponding list for each node feature, while Table 3 discloses edge features, their sizes, and a corresponding list for each type. A detailed node and edge feature construction approach is described below.

[0091] Feature (size)listAtom type (20)C, H, O, N, S, P, F, Cl, Br, I, B, Se, Fe, Ru, Mn, Co, Ni, Cu, Zn, other (one-hot)Degree of atom (7)0, 1, 2, 3, 4, 5, 6 (one-hot)Num Hs (5)0, 1, 2, 3, 4 (one-hot)Implicit Valence (6)0, 1, 2, 3, 4, 5 (one-hot)Explicit Valence (6)0, 1, 2, 3, 4, 5 (one-hot)Implicit Hs (5)0, 1, 2, 3, 4 (one-hot)Explicit Hs (2)0, 1 (one-hot)Hybridization (5)sp, sp2, sp3, sp3d, sp3d2 (one-hot)Radical (3)0,1,2 (one-hot)Formal charge (7)-2, -1, 0, 1, 2, 3, 4 (one-hot)Aromatic, Ring, Hydrophobic, H Donor, H Accepter, Acidic, Basic (7)Binary

[0092] Edge typeFeature (size)listCovalent bondBond type (4)Single, double, triple, aromatic (one-hot)Stereo type (6)Any, Cis, Trans, E, Z, none (one-hot)Ring (1), Conjugating (1)BinaryNon-covalent bondPL bond type (5)Hydrophobic, PL Hydrogen, LP Hydrogen, PL ionic, LP ionic (binary)Distance (4)<2, <4, <6, <8 (binary)Atom-Atom connection (3)Protein-Protein, Ligand-Ligand, Protein-Ligand (one-hot)

[0093] Ligand molecules and proteins are represented as graphs with atoms and two valence bonds. The properties of nodes (atoms) are represented by 73 node features (Table 2), and the properties of bonds are represented by 24 edge features (Table 3). First, the edge connection type is used to indicate whether the bond is within the ligand, within the protein, or between the ligand and the protein. The 20 atom types (C, H, O, N, S, P, F, Cl, Br, I, B, Se, Fe, Ru, Mn, Co, Ni, Cu, Zn, and others) are represented as one-hot vectors. Various other properties of atoms were obtained using RDKit tools and represented in one-hot and binary. There are five features (Hydrophobic, H Donor, H Acceptor, Acidic, Basic) that could not be obtained with RDKit tools. Among these, Hydrophobic is represented separately as 0 if the charge of each atom calculated by Autodock Tools (MGLTools) is less than 0.2 and 1 if it is greater. H Donor and H Acceptor are represented as 0 and 1 depending on the number of H that N, O, S, and P atoms have, and Acidic and Basic are represented as 0 and 1 for carboxyl and amino groups.

[0094] The types of bonds used are covalent bonds between the ligand and the protein, and non-covalent bonds that interact between the ligand and protein atoms. Edge features are expressed as one-hot vectors according to bond type, stereotype, ring, conjunction, and distance. Covalent bonds are expressed as one-hot vectors of four types (single, double, triple, and aromatic) and five stereogenic directions (Any, Cis, Trans, E, Z, and None). Ring and conjugation were then checked. A non-covalent bond is formed between a ligand atom and a protein atom when the distance between the two atoms is within 8 Å. There are two types of non-covalent bonds: PL (protein-ligand) bond type and distance bond type. PL bond types include hydrophobic, hydrogen bond, and ionic bond. If both bonding atoms are hydrophobic, it is a hydrophobic bond. If an atom has one H acceptor and one H donor, it is a hydrogen bond. Similarly, if each atom has an acidic and basic property, it is an ionic bond. Hydrogen bonds and ionic bonds are distinguished from ligand-to-protein or protein-to-ligand. The distance between a ligand and a protein atom is expressed in four binary numbers (<2, <4, <6, <8).

[0095] (2) Reinforcement by embedding layer

[0096] After the initial graph is constructed, the initial node features are enhanced by an embedding layer before being fed into a graph neural network (GNN). Each structured graph is input to an embedding layer, which outputs an embedding vector.

[0097] The embedding layer, which is applied to all three models above, is a simple linear layer that projects the 73-dimensional node feature vector into a 256-dimensional node embedding.

[0098] (3) Model training

[0099] The enhanced node embedding (embedding vector) is fed into the GNN along with the edge features to update the node representation with more information.

[0100] All GNNs in the above three models are built by stacking multiple GATv2 (Graph Attention Network version 2) layers, implemented in the PyTorch Geometric library (PyTorch Geometric, https: / / pytorch-geometric.readthedocs.io / en / latest / ), with default settings. For more information on GATv2, see How Attentive are Graph Attention Networks? (https: / / arxiv.org / abs / 2105.14491).

[0101] As shown in Figures 3 and 4, for training the Akscore2_NonDock model, the enhanced protein binding pocket graph and ligand graph are input to the protein GNN (GNN-Protein) and ligand GNN (GNN-Ligand), respectively.

[0102] For training the Akscore2_DockS model, the enriched protein-ligand complex graph (embedding vector) with non-covalent bonds is input to the protein-ligand GNN (GNN-Protein Ligand Complex) (see Figures 3 and 5).

[0103] For training the Akscore2_DockC model, the enriched protein-ligand complex graph (embedding vector) with non-covalent bonds is input to the protein-ligand GNN (GNN-Protein Ligand Complex) (see Figures 3 and 6).

[0104] Protein GNN (GNN-Protein) and Protein-Ligand GNN (GNN-Protein Ligand Complex) each have the same network structure consisting of five layers of Gatv2, while Ligand GNN (GNN-Ligand) consists of three layers of Gatv2 because it has far fewer atoms and bonds than the protein binding pocket. Relu activation layers and dropout layers are sequentially added following each Gatv2 layer.

[0105] After updating the node embeddings of the graph by GNN, a global mean pooling layer is applied to average all node embeddings of the graph and obtain a single multidimensional vector, the graph-level embedding.

[0106] (4) Model output and loss function

[0107] The above process is described as being applied identically or similarly to the three models, but the output value and loss function of each model are different depending on the model and are described individually.

[0108]

[0109] Below, the training and inference processes for each model are described in detail.

[0110] Model 1: Akscore2_NonDock model

[0111] Referring to Fig. 4, the structure of the Akscore2_NonDock model, which inputs the protein binding pocket graph and the ligand graph separately without considering the binding pose and outputs the interaction probability between them, is described in detail.

[0112] As explained above, in the initial graph construction step, the input protein forget graph and ligand graph are constructed individually from the protein forget structure and the ligand structure, respectively, and each graph has a 73-dimensional vector as a node feature and a 12-dimensional vector as an edge feature.

[0113] After initial graph construction, node features are enhanced (from 73 dimensions to 256 dimensions) by an embedding layer before being fed into a graph neural network (GNN). At this stage, the structured protein and ligand graphs are each input to the embedding layer, which outputs enhanced embeddings.

[0114] The enhanced protein pocket graph and ligand graph, which are the outputs of the embedding layer, are input to the protein GNN (GNN-Protein) and ligand GNN (GNN-Ligand), respectively. As described above, the protein GNN consists of five layers of Gatv2 layers, Relu activation layers, and dropout layers, and the ligand GNN consists of three layers of Gatv2 layers, Relu activation layers, and dropout layers.

[0115] After the node embeddings of the above-mentioned enhanced protein and ligand graphs are updated by the protein GNN (GNN-Protein) and ligand GNN (GNN-Ligand), respectively, the output values ​​are applied (input) to a global mean pooling layer to average all node embeddings of the graph and obtain a single multidimensional vector, the graph-level embedding.

[0116] The protein binding pocket graph embedding and ligand graph embedding output from the above global average pooling layer are concatenated to form a single combined graph embedding.

[0117] The above single coupled graph embedding is input to a multi-layer perceptron, which outputs a binary classification value (active or inactive). The multi-layer perceptron consists of three linear layers with a ReLu activation layer and a dropout layer between the linear layers (linear layer-ReLu-Dropout-linear layer-ReLu-Dropout-linear layer). The loss function used for binary classification here is the binary cross entropy (BCE) loss function.

[0118] The Akscore2-NonDock model is designed to measure interactions between protein binding pockets and ligands, i.e., activation or inactivity. Therefore, binary cross entropy (BCE) was used as the loss function when training this model.

[0119] In Equation 1 below, L ND represents the loss value of this model, and P pred represents the output value of this model, and P true represents the ground truth. P of the native protein-ligand complex in the PDBind data set true The value is 1 and the P of the native protein-ligand complexes of Decoys_Cross and Decoys_Random true The value is 0. Decoys_Dock is not used when training this model because it is designed not to consider binding poses.

[0120] L ND = BCE (P pred , P true ) (1)

[0121] During the training phase, native protein-ligand complexes, Decoys_Cross and Decoys_Random, were fed to the model in a 1:2:2 ratio in one batch.

[0122] To summarize, the above training process includes the following steps:

[0123] 1. Step 1: Graphically construct protein pockets and ligands of target protein-ligand pairs in the training data set independently.

[0124] 2. Step of inputting the structured protein graph and ligand graph into the embedding layer and outputting the embedding vector respectively.

[0125] 3. Step of independently inputting the embedding vector of the protein graph and the embedding vector of the ligand graph into the protein GNN and the ligand GNN, respectively.

[0126] 4. Step of performing global mean pooling on the output values ​​of protein GNN and ligand GNN respectively.

[0127] 5. Fusion step that concatenates the global average pooling outputs of the protein GNN and ligand GNN.

[0128] 6. A step of inputting the output of the above fusion step into a multi-layer perceptron to output the binary classification result, which is either active or inactive.

[0129] 7. Step 7: Optimizing the parameters of the first model using binary cross entropy (BCE) as a loss function.

[0130] Steps 1 to 7 above are repeated to update the first model.

[0131] In the inference step, the sigmoid function outputs P predIt is applied to output a value in the range of 0 to 1, with a larger value indicating a higher interaction probability between the protein binding pocket and the ligand. Therefore, this model outputs the interaction probability (ranging from 0 to 1) between the protein binding pocket and the ligand.

[0132] To train this model, a binary classification model was created using both a data set with known binding affinity and a native set (negative set) that does not actually bind, as training data, so the probability of binding is predicted with high accuracy.

[0133]

[0134] Model 2: Akscore2_DockS model

[0135] The Akscore2_DockS model is designed to predict protein-ligand interactions, and takes a protein-ligand complex with a binding pose as input data and independently outputs binding free energy (kcal / mol) and RMSD (Å).

[0136] Hereinafter, the training method of the Akscore2_DockS model is described with reference to FIG. 5, but the description of the same configuration as the Akscore2-NonDock model is omitted.

[0137] The training method of the Akscore2_DockS model includes the following steps.

[0138] 1. A step of constructing a protein-ligand complex graph by structuring the protein-ligand complexes of the training data set (input: protein-ligand complex, output: binding free energy and RMSD) into a graph.

[0139] 2. A step of inputting the above protein-ligand complex graph into the embedding layer to output a protein-ligand embedding vector.

[0140] 3. Step 3: Input the protein-ligand embedding vector into the protein-ligand GNN.

[0141] 4. A step of inputting the output value of the above protein-ligand GNN into a global mean pooling layer and performing global mean pooling.

[0142] 5. Step of inputting the output values ​​of the global mean pooling layer into the energy multi-layer perceptron and the RMSD multi-layer perceptron to output the energy and RMSD, respectively.

[0143] 6. Step for calculating the loss function using the energy, which is the output value of the energy multilayer perceptron, and the RMSD, which is the output value of the RMSD multilayer perceptron.

[0144] 7. A step of optimizing the parameters of the model by executing regression analysis to minimize the above loss function.

[0145] The Akscore2_DockS model is trained using the above steps, and the model parameters are optimized by performing regression analysis using a loss function. Furthermore, steps 1 through 7 are repeated as additional training data is added.

[0146] The Akscore2_DockS model has two multilayer perceptrons (neural networks) to individually predict the energy and root mean square deviation (RMSD) of protein-ligand complexes. Here, RMSD represents the relative distance of the ligand to its native pose, if such a pose exists. In other words, RMSD is a value that indicates the extent to which the 3D structure of the protein-ligand complex input for training or inference deviates from the actual structure of the protein-ligand complex. A value of 0 indicates that the structure is identical to that of a naturally occurring protein-ligand, while a larger value indicates a greater difference from the actual structure.

[0147] The task of Akscore2_DockS is modeled as a regression problem predicting real-valued quantities.

[0148] Native protein-ligand complexes, Decoys_Cross, and Decoys_Random datasets were used to predict binding free energies (kcal / mol) in this model.

[0149] When the input data is a native protein-ligand complex dataset, its loss L DS_E_native is calculated by the mean squared error (MSE) loss function (Equation 2), and e pred is the binding free energy predicted by the model, and e true is the observed binding free energy provided by the PDBbind data set.

[0150] In general, since random molecular ligands are not true binders for a specific protein, we treat all protein-ligand pairs in Decoys_Cross and Decoys_Random as inactive and have energies higher than the threshold of -5 kcal / mol. Therefore, the loss L for these inputs (Decoys_Cross, Decoys_Random) DS_E_cross_random is calculated by Equation 3, which is similar to PIGNet (PIGNet: a physics-informed deep learning model toward generalized drug-target interaction predictions). This loss function is e pred If e is lower than the threshold -5, the difference from the threshold is considered a loss, and conversely, e pred If is above the threshold, it is considered to provide a correct prediction, so the loss is explained as 0.

[0151] To predict RMSD (Å), we used the native protein-ligand complex and Decoys_Dock datasets. When the input is a native protein-ligand complex, the model should try to output 0, since the ground truth RMSD of these samples is always 0. Therefore, the loss for the native protein-ligand complex, L DS_R_native is calculated using MSE as written in Equation 4, where r pred is the predicted RMSD value by this model. For the samples in Decoys_Dock, native ligands with various different binding poses were redocked into the protein binding pocket. The loss of Decoys_Dock, L DS_R_dock It is calculated using the mean squared error (MSE) loss function as described in Equation 5. Here, r true is the RMSD value calculated compared to the native pose.

[0152] Final total loss L of the Akscore2_DockS model DS_total L listed in Equations 2 to 5 DS_E_native , L DS_E_cross_random , L DS_R_native , L DS_R_dock It is the sum of the four losses (Equation 6). L DS_E_cross_random Note that we added weights by multiplying the factor by 5.

[0153]

[0154] During the training phase, native protein-ligand complexes, Decoys_Dock, Decoys_Cross, and Decoys_Random, were fed to the model in batches with a ratio of 1:4:2:2. As can be seen in the loss function above, when calculating the energy loss, the output samples of Decoys_Dock were ignored because they lack accurate ground truth energy values. Similarly, when calculating the loss for RMSD, the output samples of Decoys_Cross and Decoys_Random were ignored because they lack true binding poses for calculating RMSD.

[0155] The step of generating a protein-ligand binding free energy prediction value using the second model includes the steps of structuring the structure of the protein-ligand complex into a graph to construct a protein-ligand graph; inputting the protein-ligand graph into a protein-ligand GNN; inputting the output of the protein-ligand GNN into a pooling layer and pooling it; and inputting the output of the pooling layer into a multi-layer perceptron and outputting the binding free energy. At this time, a step of inputting the protein-ligand graph into an embedding layer and outputting an embedding vector may be added before the step of inputting it into the protein-ligand GNN. In this case, the embedding vector is input into the protein-ligand GNN.

[0156] During the inference phase, the predicted RMSD value is treated as a penalty, with a high RMSD indicating a poor binding pose. Consequently, the two separately predicted values ​​(energy, RMSD) must be added to obtain the final energy for model evaluation.

[0157] In this regard, the inventors thought about whether it would be possible to directly predict energy values ​​with RMSD penalty rather than predicting energy values ​​and RMSD separately, and thus additionally designed the Akscore2_DockC model.

[0158] Third model: Akscore2_DockC model

[0159] The Akscore2_DockC model can be viewed as the Akscore2_DockS model with the RMSD prediction multilayer perceptron removed and only the energy regression branch.

[0160] The Akscore2_DockC model predicts protein-ligand interactions by taking as input a protein-ligand complex with a binding pose and outputting the binding free energy (kcal / mol). However, this output value is the binding free energy with the RMSD penalized.

[0161] For example, suppose that the ligand that actually binds to a specific protein pocket (correct) is in the correct pose and has an RMSD of 0 and a binding free energy of -12. If we randomly change the pose to have an RMSD of 4 Å to create augmented data, the correct output that Akscore2_DockC should predict is the penalized binding free energy of -12 + 4 = -8. On the other hand, the Akscore2_DockS model predicts values ​​of -12 and 4 in this case, respectively.

[0162] Referring to Fig. 6, the training method of the Akscore2_DockC model includes the following steps.

[0163] 1. Step of constructing a protein-ligand complex graph by structuring the protein-ligand complexes of the training data set into a graph.

[0164] 2. A step of inputting the above protein-ligand complex graph into the embedding layer to output a protein-ligand embedding vector.

[0165] 3. Step 3: Input the protein-ligand embedding vector into the protein-ligand GNN.

[0166] 4. A step of inputting the output value of the above protein-ligand GNN into a global mean pooling layer and performing global mean pooling.

[0167] 5. Step of inputting the output value of the global mean pooling layer to a multi-layer perceptron.

[0168] 6. A step for calculating the loss function using the energy and RMSD, which are the output values ​​of the multilayer perceptron (here, RMSD is calculated by adding a value pre-calculated by the structure including the ligand pose to the energy value, calculating the final energy value with the penalty, and reporting it as the correct value and inputting it into the loss function).

[0169] 7. A step of optimizing the parameters of the model by executing regression analysis to minimize the above loss function.

[0170] Steps 1 to 7 above are repeated as training data is added.

[0171] Equations 7-8 below represent the loss function of the Akscore2_DockC model for the native protein-ligand complex, Decoys_Cross, and Decoys_Random output samples, which is identical to the loss function for the energy prediction part of Akscore2_DockS.

[0172] For the output samples of Decoys_Dock, the loss using the mean squared error (MSE) loss function is calculated by Equation 9, where e pred is the predicted energy by the Akscore2_DockC model, and e true is the observed binding free energy of the original native protein-ligand complex before redocking, and r trueis the RMSD value calculated by comparing it with the native pose of the original native protein-ligand complex before performing redocking. Since the samples in Decoys_Dock do not have ground truth energy values, we intend to penalize them by the distance from the true binding pose, which is called e true + r true was imitated by

[0173] Final total loss L of the Akscore2_DockC model DC_total L listed in Equations 7-9 DC_E_native , L DC_E_cross_random , L DC_E_dock is the sum of the three losses (Equation 10). L DC_E_cross_random Note that we added weights by multiplying the factor by 5.

[0174]

[0175]

[0176] During the training phase, the native protein-ligand complexes, Decoys_Dock, Decoys_Cross, and Decoys_Random, were fed to the model in a batch ratio of 1:4:2:2. Unlike the Akscore2_DockS model, Akscore2_DockC does not need to ignore specific losses on specific datasets.

[0177]

[0178] (5) Training Details and Ensemble Strategy for Evaluation

[0179] All the above models were built using PyTorch (PyTorch, https: / / pytorch.org / ) and the PyTorch Geometric library, and trained for approximately 1000 epochs using the Adam optimizer (Adam: A Method for Stochastic Optimization, https: / / arxiv.org / abs / 1412.6980) with a learning rate of 0.0001 and dropout rates of 0.1 for all dropout layers. The batch size was set to 30 for Akscore2_NonDock, and 36 for Akscore2_DockS and Akscore2_DockC.

[0180] An ensemble strategy was used to obtain the numerical values. Akscore2_NonDock is a non-docking-based model that outputs the interaction probability between the protein binding pocket and the ligand, ranging from 0 to 1. The other two models are docking-based models, and their final output is a real-valued energy. The non-docking and docking-based models are ensembled by multiplying the outputs of these models.

[0181] Meanwhile, for the ensemble strategy, the physical scoring term Autodock Score was combined with the deep learning model. Autodock Score(d) is the output (energy value) of the AutoDock GPU model. The AutoDock GPU model is a tool that considers and calculates intermolecular interactions (van der Waals energy, hydrogen bonding, electrostatic attraction, and solvation effects) and outputs binding energy values. For details, see Accelerating AutoDock4 with GPUs and Gradient-Based Local Search, J. Chem. Theory Comput. 2021, 10.1021 / acs.jctc.0c01006. This model uses a physics-based scoring function.

[0182] A scoring function is a function used to evaluate conformations generated through a sampling algorithm during the docking process. It is applied to compare non-binding and binding compounds, or to identify relatively accurate poses among multiple poses. Since the purpose of a scoring function is not to calculate precise binding affinity, but rather to enable relative comparisons through estimated values, absolute scoring function values ​​are not particularly meaningful and can only be meaningfully used when comparing calculated values ​​between structures.

[0183] The ensemble strategy is described in detail with reference to Figures 7 and 8 below.

[0184] According to Fig. 7, the final ensemble score is Akscore2_NonDock x Akscore2_DockS + αx Autodock_Score. Specifically, it is as follows: Equation 11.

[0185] ---- Equation 11

[0186] p is the protein-ligand interaction probability, which is the output of Akscore2_NonDock, e is the energy, which is the output of the Akscore2_DockS model, r is the RMSD, which is the output of the Akscore2_DockS model, and d is the Autodock_Score, which is the output (energy value) of the AutoDock GPU model. Here, α is a weighting parameter and is empirically set to 0.65, but is not limited thereto.

[0187] According to FIG. 8, the ensemble score S' according to another embodiment of the present invention can be described as Akscore2_NonDock × Akscore2_DockC + α × Autodock_Score, and specifically, as shown in the following Equation 12.

[0188] (12)

[0189] Here, α is a weighting parameter, empirically set to 0.65 but is not limited thereto. p is the protein-ligand interaction probability output from Akscore2_NonDock, c is the output from the Akscore2_DockC model, and d is the Autodock_Score, output from the AutoDock GPU model.

[0190] AutoDock GPU is widely known; see Accelerating AutoDock4 with GPUs and Gradient-Based Local Search, J. Chem. Theory Comput. 2021, 10.1021 / acs.jctc.0c01006. This model utilizes a physics-based scoring function.

[0191] The prediction method and learning method according to the above-described embodiment of the present invention can be implemented by a computer program executable in an information processing device.

[0192] 3. Test

[0193] Benchmark data set

[0194] As benchmark data, we used the CASF2016 benchmark dataset consisting of the core set of PDBbind v2016 (see Minyi Su, Qifan Yang, Yu Du, Guoqin Feng, Zhihai Liu, Yan Li, Renxiao Wang, "Comparative Assessment of Scoring Functions: The CASF-2016 Update", J. Chem. Inf. Model, 2019, Vol. 59: pp 895-913) and the LIT-PCBA dataset (LIT-PCBA: An Unbiased Data Set for Machine Learning and Virtual Screening).

[0195]

[0196] CASF2016 Benchmark Test Results

[0197] The CASF2016 benchmark dataset consists of 285 docking decoys of protein-ligand complexes, including docking decoys (i.e., conformational decoys) to assess the accuracy of identifying near-native conformations, and decoys screening data (i.e., cross-docked complexes) to assess the accuracy of distinguishing true binders from non-binders.

[0198] Conformational decoy data have RMSDs in the range of 0 to 10 Å relative to poses obtained by redocking ligands from 285 complexes.

[0199] Screening decoys were generated by docking 285 ligands onto 57 proteins, each with 100 poses. For evaluation, each of the 57 proteins had 5 of the 285 active ligands. 285 samples were used for scoring and ranking, 21,465 samples for docking, and 1,589,263 samples for the screening benchmark. For each target, a graphical representation of the protein pocket was generated using the benchmark data, and samples not recognized by the RDKit toolkit were excluded. When evaluating docking data, an RMSD of less than 2 Å was considered active.

[0200] The benchmark test results using the CASF-2016 dataset are shown in Table 4.

[0201] 모델(Model)ScoringRankingDockingForward ScreeningReverse ScreeningPearson's rSpearman's ρSuccess RateTop 1Average EFTop 1%Success RateTop 1%Success RateTop 1%Autodock GPU0.5970.44775.8%6.5922.8%18.6%AutoDock Vina0.6040.52890.2%7.729.8%13.7%PIGNet (single)0.7490.66885.8%(without native pose)18.550.0%---PIGNet (ensemble)0.7610.68287.0%(without native pose)19.655.4%---RTMScore (average)0.4550.52997.3%28.066.7%37.6%Akscore2_NonDock0.0540.12345.3%17.0521.1%29.5%Akscore2_DockS0.6500.57963.9%6.4115.8%11.2%Akscore2_DockC0.7060.52367.7%12.3028.1%20.0%Akscore2_NonDock × Akscore2_DockS0.3540.40575.1%24.3447.4%34.0%Akscore2_NonDock × Akscore2_DockC0.3900.33770.9%29.2949.1%42.1%Akscore2_NonDock × Akscore2_DockS+ Autodock_GPU0.5510.45582.8%29.4857.9%41.8%Akscore2_NonDock × Akscore2_DockC + Autodock_GPU0.5390.46482.8%34.1973.7%43.5%

[0202] In Table 4, AutoDock GPU is widely known, and is referred to as Accelerating AutoDock4 with GPUs and Gradient-Based Local Search, J. Chem. Theory Comput. 2021, 10.1021 / acs.jctc.0c01006. This model utilizes a physics-based scoring function.

[0203] The results of AutoDock Vina in Table 4 are taken from CASF2016 by Minyi Su, Qifan Yang, Yu Du, Guoqin Feng, Zhihai Liu, Yan Li, Renxiao Wang, "Comparative Assessment of Scoring Functions: The CASF-2016 Update", J. Chem. Inf. Model, 2019, Vol. 59: pp 895-913. AutoDock Vina and others are the results tested using open source.

[0204] The results of PIGNet are from the paper (PIGNet: a physics-informed deep learning model toward generalized drug-target interaction predictions), and the results of RTMScore are from the paper (Shen C, Zhang X, Deng Y, Gao J, Wang D, Xu L, Pan P, Hou T, Kang Y. Boosting Protein-Ligand Binding Pose Prediction and Virtual Screening Based on Residue-Atom Distance Likelihood Potential and Graph Transformer. J Med Chem. 2022 Aug 11;65(15):10691-10706. doi: 10.1021 / acs.jmedchem.2c00991. Epub 2022 Aug 2. PMID: 35917397.).

[0205] The best results for each column representing the evaluation factors are shown in bold.

[0206] Particularly important factors in Table 4 are the average EF and success rate in forward screening. The average EF (average enrichment factor) is an indicator of how much the active ratio among molecules selected by the model increases relative to the interaction activation ratio of all randomly selected molecules for the target protein. An EF greater than 1 indicates a high selectivity, a value near 1 indicates a lack of selectivity, and conversely, a value near 0 indicates that inactive molecules were selected. An EF of 1% is the EF when only the top 1% of molecules based on score are selected. The Akscore2_NonDock model shows a good result of 17.05, and the ensemble combining the outputs of the models of the present invention all have an EF of 24 or higher.

[0207] The top 1% success rate is the probability that a molecule binds to the target protein among the top 1% of molecules in terms of scores. In the case of the ensemble method, all methods show excellent results of over 40%. In particular, Akscore2_NonDock x Akscore2_DockC + Autodock_GPU shows the best results with an average EF 1% of 34.19 and a top 1% success rate of 73.7%.

[0208]

[0209] LIT-PCBA Benchmark Test Results

[0210] For another benchmark test, we used the LIT-PCBA dataset, generated from high-confidence PubChem bioassay data. The LIT-PCBA dataset contains hit data for 15 proteins with 9,780 active and 407,839 unique inactive ligands.

[0211] First, the lowest energy 3D conformer (3D conformer, shape, stereostructure) of each ligand was generated using the RDKit tool.

[0212] Second, since LIT-PCBA contains multiple receptor PDBs per target, the PDB with the highest resolution was selected. By redocking the active and inactive ligands for each protein, 50 poses were generated for 10,030 active and 2,795,002 inactive ligands.

[0213] The screening results for the LIT-PCBA dataset are shown in Table 5. The screening results for the LIT-PCBA dataset are the average enrichment factor (EF) obtained by considering the top 0.5%, 1%, and 5% candidates for each target protein. Although PIGNet did not disclose the trained weights, we trained PIGNet for 2,300 epochs and confirmed that its performance was similar to the PIGNet paper, which achieved an EF of 18.86 for 1% forward screening on CASF-2016. For RTMScore, RTMScore1, the best-performing model on CASF-2016, was used in this test. The best results in each column are indicated in bold. The results of AutoDock GPU were used from the paper (Accelerating AutoDock4 with GPUs and Gradient-Based Local Search, J. Chem. Theory Comput. 2021, 10.1021 / acs.jctc.0c01006).

[0214]

[0215] ModelForward ScreeningAverage EF 0.5%Average EF 1%Average EF 5%Autodock GPU2.582.121.55PIGNet (single)1.681.431.60RTMScore (RTMScore1)2.873.261.70Akscore2_NonDock1.060.960.83Akscore2_DockS1.121.411.53Akscore2_DockC3.272.601.59Akscore2_NonDock × Akscore2_DockS1.040.941.37Akscore2_NonDock × Akscore2_DockC1.681.531.35Akscore2_NonDock × Akscore2_DockS +Autodock_GPU2.752.261.59Akscore2_NonDock × Akscore2_DockC +Autodock_GPU1.531.661.89

[0216] 4. Review

[0217] In an embodiment of the present invention, for model training, a model was created using not only a data set of known binding affinity values ​​for actual protein-ligand complexes, but also a negative data set, where binding does not occur. This solves the problem of false positives in other conventional models due to training data limitations.

[0218] Models 2 and 3 improve prediction accuracy by designing a loss function that predicts binding affinity while penalizing input structures that do not actually bind or are incorrectly bound. The loss function is designed to determine whether the input protein-ligand structure is similar to the actual crystal structure based on the root mean square deviation (RMSD). That is, the further away from the actual crystal structure the higher the RMSD value, resulting in a penalty.

[0219] When there is a structure in which a protein and a ligand are bound, the first model predicts whether binding occurs, and if binding occurs, the second and third models are proposed that consider not only the energy but also the accuracy of the binding free energy and conformation.

[0220] The three models can be used independently to predict interaction probabilities, binding energies, and RSMD, binding energies, respectively.

[0221] When the output of the classification model (model 1) and the output of the regression model (model 2 or model 3) are multiplied, particularly good results can be obtained in screening.

[0222] The best prediction results are obtained by multiplying the output of the binary classification model by the output of the regression model (the second or third model) and adding the output of the Autodock_Score, multiplied by the weights. This can be confirmed by the validation factor values ​​in Tables 4 and 5.

[0223] In addition, according to an embodiment of the present invention, compared to graph construction using conventional 10 to 20-dimensional features, by constructing a node graph of 73-dimensional vectors, atomic characteristics can be better reflected, thereby improving the accuracy of the prediction model.

[0224] According to the present invention, the overfitting problem was solved by using both a data set with known binding affinity and a data set of protein ligands that do not actually bind as training data.

[0225] The second and third models are models that consider both binding energy and conformational accuracy, and by designing a loss function that can predict binding energy while penalizing structures that do not actually bind or are incorrectly bound, the problems of models that only consider energy or binding affinity can be solved. That is, the loss function was designed to evaluate whether a given protein-ligand structure is close to the actual (crystal) structure by RMSD (the lower the value, the more similar it is to the actual structure) and apply a penalty accordingly. In other words, the more different the structure is from the actual one, the higher the RMSD value, resulting in a penalty.

[0226] Additionally, excellent prediction results can be obtained by multiplying the joint probability predicted by the binary classification model and the output of the second and third models (regression model) to obtain the final result.

[0227] The present invention can be used in the pharmaceutical industry, such as drug discovery.

Claims

1. A method for predicting protein-ligand interactions executed by an information processing system including a processor and memory, A step of predicting a protein-ligand interaction probability by a first model based on GNN, wherein a protein graph and a ligand graph, each of which is structured as a separate graph for the protein pocket and ligand of a target protein-ligand pair, are input into the first model to output a protein-ligand interaction probability; A step of predicting protein-ligand binding free energy by a second or third model based on GNN, wherein a protein-ligand graph, which is a structure of a target protein-ligand complex structured as a graph, is input into the second or third model to output a protein-ligand binding free energy prediction value; and A step of generating a protein-ligand interaction prediction value by multiplying the protein-ligand interaction probability by the protein-ligand binding free energy prediction value, A protein-ligand interaction prediction method characterized in that the first, second and third models are GNN-based neural network models.

2. In paragraph 1, To the predicted protein-ligand interaction, a physical scoring term multiplied by a parameter is added to generate an ensemble score. A protein-ligand interaction prediction method characterized in that the above physical scoring term is an output value of an AutoDock GPU model.

3. In paragraph 1, The input data of the first, second and third models can be described as an undirected graph G=(V, E), where V represents a group of nodes, and one node represents one atom of a protein or ligand, and E is a group of edges representing the characteristics of interactions between nodes. The protein graph and ligand graph, which are input values ​​of the first model, have a 73-dimensional vector as a node feature and a 12-dimensional vector as an edge feature. The protein-ligand graph, which is the input value for the second and third models, has a 73-dimensional vector as a node feature and a 24-dimensional vector as an edge feature. The first model uses the protein binding pocket and the ligand separately as input data to predict the interaction probability of the target protein-ligand pair without considering the binding pose of the target protein-ligand pair. The second and third models are protein-ligand interaction prediction methods characterized in that they use protein-ligand complexes having binding poses as input data to predict binding free energy.

4. In paragraph 1, The step of predicting the protein-ligand interaction probability by the first model is: A step of constructing a protein graph and a ligand graph by structuring the protein pockets and ligands of a target protein and ligand pair into individual graphs; A step of inputting the protein graph and ligand graph into an embedding layer and outputting an embedding vector; A step of independently inputting the embedding vector of the protein graph and the embedding vector of the ligand graph into the protein GNN and the ligand GNN, respectively; A step of pooling the output values ​​of the protein GNN and the output values ​​of the ligand GNN and then fusing them together; A step of inputting the output value of the above fusion step into a multi-layer perceptron and outputting the output value of the binary classification; and A method for predicting protein-ligand interaction, characterized by including a step of predicting protein-ligand interaction probability by applying a sigmoid function to the output value of the above binary classification.

5. In paragraph 1, The above first model is, It includes a protein GNN into which a protein graph is input, a ligand GNN into which a ligand graph is input, a pooling layer for pooling the output values ​​of the protein GNN and the ligand GNN, a concatenation layer for concatenating the output values ​​of the pooling layer, and a multilayer perceptron for inputting the output values ​​of the concatenation layer and outputting the output values ​​of the binary classification, The above protein GNN and ligand GNN each consist of one or more GAT layers, A protein-ligand interaction prediction method characterized by optimizing the parameters of a neural network model by inputting the output value of the above binary classification into a binary cross entropy (BCE) loss function.

6. In paragraph 5, The training data set of the above first model includes a native data set and an augmented data set, The native data set above contains protein pockets and ligand pairs of real native protein-ligand complexes with observed binding affinities, The augmented dataset includes a cross decoy dataset and a random decoy dataset, A method for predicting protein-ligand interactions, wherein the augmented data set comprises pairs of protein pockets and ligands generated by randomly selecting ligands other than the native protein-ligand complex's own ligands and redocking them with the protein pockets.

7. In paragraph 1, The step of predicting the protein-ligand binding free energy by the second model is: A step of constructing a protein-ligand complex graph by structuring the structure of a protein-ligand complex of a target protein-ligand pair into a graph; A step of inputting the protein-ligand complex graph into an embedding layer and outputting a protein-ligand embedding vector; A step of inputting the above protein-ligand embedding vector into a protein-ligand GNN; A step of inputting the output value of the protein-ligand GNN to a pooling layer and pooling; and A method for predicting protein-ligand interactions, characterized by including a step of inputting the output value of a global mean pooling layer into an energy multi-layer perceptron and an RMSD multi-layer perceptron, respectively, to output energy and RMSD, respectively.

8. In paragraph 1, The above second model is, It includes a protein-ligand GNN, a pooling layer for pooling the output values ​​of the protein-ligand GNN, an energy multilayer perceptron that inputs the output values ​​of the pooling layer and outputs energy, and an RMSD multilayer perceptron that inputs the output values ​​of the pooling layer and outputs RMSD. A protein-ligand interaction prediction method characterized in that the above protein-ligand GNN is composed of 5 layers of GATv2.

9. In paragraph 1, The step of predicting protein-ligand binding free energy using the third model is: A step of constructing a protein-ligand graph by structuring the structure of a protein-ligand complex of a target protein-ligand pair into a graph; Step of inputting the protein-ligand graph into the protein-ligand GNN; A step of inputting the output of the protein-ligand GNN to a pooling layer and pooling; and A method for predicting protein-ligand interactions, characterized by comprising a step of inputting the output of the above pooling layer into a multi-layer perceptron to output binding free energy.

10. In paragraph 1, The third model is, It includes a protein-ligand GNN, a pooling layer for pooling the output values ​​of the protein-ligand GNN, and a multilayer perceptron into which the output values ​​of the pooling layer are input. A method for predicting protein-ligand interactions, characterized in that the parameters of a neural network are optimized through regression analysis.

11. In paragraph 1, The training data set for training the second model is Contains native and augmented data sets of protein-ligand complexes, A method for predicting protein-ligand interactions, wherein the augmented data set comprises a structural decoy data set, a Decoys_Cross data set, and a Decoys_Random data set.

12. In paragraph 2, Ensemble score S is, And, Here, p is the protein-ligand interaction probability which is the output value of the first model, e is the energy which is the output value of the second model, r is the RMSD which is the output value of the second model, d is the Autodock_Score which is the output value of the AutoDock GPU model, and α is a weighting parameter. A protein-ligand interaction prediction method characterized in that:

13. In paragraph 2, Ensemble score S is, And, A protein-ligand interaction prediction method characterized in that p is the protein-ligand interaction probability, which is the output value of the first model, c is the output value of the third model, d is the Autodock_Score, which is the output value of the AutoDock GPU model, and α is a weight parameter.

14. A protein-ligand interaction prediction method executed by a protein-ligand interaction prediction system including a processor and memory, A step of storing the protein pocket structure, ligand structure, and protein pocket-ligand complex structure of the target protein-ligand pair in memory by configuring them into a protein pocket graph, a ligand graph, and a protein pocket-ligand complex graph, respectively; A step of predicting a protein-ligand interaction probability using a first model, wherein a target protein pocket graph and a ligand graph are input into a protein-GNN and a ligand-GNN, respectively, and a step of predicting a protein-ligand interaction probability is output; A step of predicting protein-ligand binding free energy using a second model, wherein a protein-ligand graph, which is a structure of a target protein-ligand complex structured as a graph, is input into a protein-ligand GNN to output a protein-ligand binding free energy prediction value; and A method for predicting protein-ligand interaction, comprising the step of generating a protein-ligand interaction prediction value by multiplying the protein-ligand interaction probability by the protein-ligand binding free energy prediction value.

Images