Methods for structure prediction and affinity estimation of small molecule drug protein-ligand complexes

Abstract

The present invention discloses a method for predicting the structure and estimating the affinity of small molecule drug protein-ligand complexes. The method comprises the following steps: constructing a protein-ligand map, extracting and updating features from the map, and obtaining a feature update module; combining the feature update module with an affinity estimation module to perform affinity estimation and obtain protein-ligand affinity; combining the feature update module with a structure prediction module to perform complex structure prediction and obtain the protein-ligand three-dimensional structure; and constructing a protein-ligand network model, training the protein-ligand network model based on the protein-ligand affinity and the protein-ligand three-dimensional structure, adjusting network parameters, and completing structure prediction. The present invention improves prediction precision and accuracy.

Description

Technical Field

[0001] The present invention belongs to the technical field of small molecule complex structure prediction, and in particular relates to a method for predicting the structure and estimating the affinity of a small molecule drug protein-ligand complex. Background Art

[0002] According to statistics, small molecules accounted for nearly 70% of all approved drugs between 2015 and 2020. Generally speaking, small molecule ligands achieve their biological activity by binding to biological macromolecules. Proteins are the most common carriers of various functions in the body. Most active small molecules bind to protein pockets to form protein-ligand complexes, thereby inhibiting or activating the corresponding protein function. Therefore, understanding the structural details of these complexes is crucial for gaining a deeper understanding of the interaction between proteins and small molecule ligands, both for improving existing drugs and developing new ones. Current techniques for structural determination include cryo-electron microscopy, X-ray crystallography, and nuclear magnetic resonance (NMR) spectroscopy. These techniques have solved the structures of over 122,000 proteins, with X-ray crystallography accounting for over 85% of these, and also for elucidating a large number of protein-ligand complexes. However, both X-ray crystallography and cryo-electron microscopy are time- and resource-intensive.

[0003] To rapidly determine the molecular structure of protein-small molecule complexes, numerous computer-assisted methods have been developed, with molecular docking being the most commonly used computer prediction method. Generally speaking, molecular docking can be divided into two interrelated steps: sampling, which first samples the conformations of the ligand in the protein's active site; and scoring, which ranks these conformations using a scoring function. Ideally, the sampling algorithm should be able to reproduce the experimental binding mode, while the scoring function should rank the conformations closest to the native state highest among all generated conformations.

[0004] Currently, the mainstream method for combining machine learning with hybrid prediction methods is primarily based on molecular docking. Recent research and improvements have focused on using machine learning as a scoring function. This includes various mainstream machine learning frameworks such as CNNs, RNNs, and GNNs. In addition, there are also reports using machine learning as a guide for minimization algorithms. Summary of the Invention

[0005] The purpose of the present invention is to propose a method for predicting the structure and affinity of small molecule drug protein-ligand complexes to improve prediction accuracy and efficiency.

[0006] To achieve the above objectives, the present invention provides a method for predicting the structure and estimating the affinity of a small molecule drug protein-ligand complex, comprising the following steps:

[0007] constructing a protein-ligand graph, and extracting and updating features of the protein-ligand graph using a feature updating module;

[0008] Combining the feature update module with the structure prediction module to perform complex structure prediction and obtain the protein-ligand three-dimensional structure; combining the feature update module with the affinity estimation module to perform affinity estimation and obtain the protein-ligand affinity;

[0009] A protein-ligand network model is constructed, and the protein-ligand network model is trained based on the protein-ligand three-dimensional structure and the protein-ligand affinity, and the network parameters are adjusted to complete the prediction.

[0010] Optionally, the protein-ligand graph includes a complete graph and a heterogeneous graph;

[0011] The complete graph includes: constructing a bidirectional complete graph based on a given protein-ligand pair, wherein each node represents an atom, and each atom is interconnected with other atoms to represent the complete graph; the complete graph consists of a set of vectors storing the characteristics of each node, and a set of matrices storing the characteristics of the edges.

[0012] The isomeric graph includes: according to a given protein-ligand pair, representing the full connection of the protein and the ligand small molecule as the isomeric graph.

[0013] Optionally, the feature update module includes: information aggregation and information feedforward;

[0014] The information aggregation includes: updating the features of the central node, obtaining the aggregated features of the central node; aggregating the feedforward information received by each central node in the feature update module to obtain new features and coordinate information

[0015] The information feedforward includes: transmitting the features and coordinate information of the adjacent nodes around each central node in the feature updating module to the central node.

[0016] Optionally, the affinity estimation process includes:

[0017] The node feature update output is passed through a two-layer multilayer perceptron to obtain a single-dimensional affinity prediction for each protein-ligand pair, calculated as follows:

[0018]

[0019] where y aff is the affinity prediction for the protein-ligand pair, N fIndicates the number of graph nodes.

[0020] Optionally, the complex structure prediction process includes:

[0021] The attention mechanism calculated by the feature update module is converted into a single-dimensional distance gradient. The distance gradient represents the contribution of each adjacent node to the distance update of the central node. The updated coordinates of the central node are obtained by summing the distance gradients of each adjacent node, which is calculated as follows:

[0022]

[0023]

[0024] where x i represents the 3D coordinate of the center node i, x j Represents the 3D coordinates of its adjacent node j.

[0025] Optionally, the protein-ligand network model is trained based on the protein-ligand three-dimensional structure and the protein-ligand affinity. The training process includes: using a Monte Carlo strategy to train the model, only training one round at a time, not performing data calculation for the cycles after this round, and not performing gradient updates for the cycles before this round, to complete the training.

[0026] Technical effect of the present invention: The present invention discloses a method for predicting the structure and estimating the affinity of a small molecule drug protein-ligand complex, which can increase the prediction accuracy and efficiency. On the one hand, considering that the sampling space of molecular docking is limited within a limited time, the accuracy has an upper limit; on the other hand, since a large amount of sampling is not required, the network directly generates a structure close to the natural state, which is more efficient. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] The accompanying drawings, which constitute part of this application, are intended to provide a further understanding of this application. The exemplary embodiments and descriptions of this application are intended to explain this application and do not constitute an improper limitation on this application. In the accompanying drawings:

[0028] Figure 1 Schematic diagram of the process of predicting the structure and estimating the affinity of a small molecule drug protein-ligand complex according to an embodiment of the present invention;

[0029] Figure 2 Schematic diagram of the structure of a method for predicting the structure and estimating the affinity of a small molecule drug protein-ligand complex according to an embodiment of the present invention. DETAILED DESCRIPTION

[0030] It should be noted that, in the absence of conflict, the embodiments and features of the embodiments in this application can be combined with each other. The present application will be described in detail below with reference to the accompanying drawings and in combination with the embodiments.

[0031] It should be noted that the steps shown in the flowcharts of the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and that, although a logical order is shown in the flowcharts, in some cases, the steps shown or described can be executed in an order different from that shown here.

[0032] like Figure 1-2 As shown, this embodiment provides a method for predicting the structure and estimating the affinity of a small molecule drug protein-ligand complex, comprising the following steps:

[0033] constructing a protein-ligand graph, and extracting and updating features of the protein-ligand graph using a feature updating module;

[0034] Combining the feature update module with the structure prediction module to perform complex structure prediction and obtain the protein-ligand three-dimensional structure; combining the feature update module with the affinity estimation module to perform affinity estimation and obtain the protein-ligand affinity;

[0035] A protein-ligand network model is constructed, and the protein-ligand network model is trained based on the protein-ligand three-dimensional structure and the protein-ligand affinity, and the network parameters are adjusted to complete the prediction.

[0036] This paper presents a deep learning-based end-to-end feature learning framework for protein-ligand molecular complexes, and a method flow for completing two major drug development tasks, drug structure prediction and affinity estimation, by combining three different task-related modules based on this framework. The protein-ligand network framework is constructed, and the constructed complete graph of protein pockets and ligands is used for drug structure prediction and affinity estimation. The process is as follows: Figure 1 As shown. The framework is as Figure 2 As shown in FIG. 2 , the protein-ligand graph construction of the present invention can be performed in two graph representations, namely, a complete graph and a heterogeneous graph.

[0037] Complete Graph: For a given protein-ligand pair, a bidirectional complete graph is constructed, where each node represents an atom, and each node is connected to every other atom (including itself). The graph represents the protein-ligand complex, with each node encoding a single atom. Each node contains the chemical features of that atom. The feature dimensions for protein nodes are 79, and for ligands, 45. These chemical features are obtained from RDKit. In addition, each node contains its 3D coordinates, with the protein coordinates retained at their original positions (i.e., those in the input file). To enable the use of symmetry-sensitive loss functions and data augmentation, the coordinates of all ligand nodes are randomly initialized to a 4nm cube around the center of the protein pocket, defined as the geometric center of the native ligand structure. Nodes contain features and coordinate positions in space, and the nodes in the graph are connected by edges, which also encode their chemical features. The number of features is 7. It is worth noting that distance information is also included in the edge features, but the distances involved in complex prediction are masked by -1 (i.e., initialized to -1), including the edges between protein and ligand and some edges within the ligand, as described later. First, the ligand is divided into sub-fragments according to its rotatable bonds. (1) The distances of the edges within each fragment are initialized to the original value, (2) The distances of the edges connected by covalent bonds are initialized to the original value, and (3) The edges between different fragments are initialized to -1. The above three considerations avoid the leakage of information about the interaction between protein and ligand and the specific conformation of the ligand.

[0038] Heterogeneous graph: For a given protein-ligand pair, the protein and ligand small molecule are represented as a fully connected bidirectional complete graph as shown above. However, this is a heterogeneous graph, that is, considering two or more types of nodes and edges. Different nodes and edges will be learned by different model weights, that is, considering the different effects of each node or edge type on the information propagation process. The initialization and isomorphism of the node position and characteristics and edge characteristics are Figure 1 However, considering that the edge features between protein and ligand are not rich in information, the edge features between protein and ligand are not given.

[0039] Sampling feature update module

[0040] The present invention designs a sampling and looping method to make the model fully utilize a large number of protein atoms. A sampling strategy is used to generate a subgraph consisting of all "core atoms" and some randomly selected "background atoms", as well as their respective edges, to input into the network. The core atoms involve all ligand atoms, as well as the C α , Cβ atoms. Other protein atoms are considered background atoms that describe the details of the amino acid structure. The maximum number of nodes in a single sampling is 500, and all edges are assigned based on the actual relationships between atoms.

[0041] Since the nodes are not fully used in a single sampling, a loop strategy is used to improve the representability and robustness of the method. In each loop, the features and coordinates of the sampled subgraph are updated using the stacked feature and coordinate update modules N l times. In addition, except for the first round of sampling, after each sampling, the features finally output by the core atom node will be added to a newly sampled subgraph through an element-wise gated residual (denoted as Ψ, the same as defined in the feature update block), while the features finally output by the background atom are ignored for "partial update". The updating principle of the edge features is the same. The coordinates of the ligand node in the newly sampled subgraph are directly inherited from the updated values ​​predicted by the previous round subgraph. Then, the subgraph will be fed into the model again, and the cycle will be repeated N times. c times (including the first sampling, a total of N c times sampling, in this work N c =3).

[0042] Each cycle uses a feature update module based on a Graph Transformer architecture to extract and update the features of nodes and edges in the subgraph. Before being formally input into the feature update module, all features are transformed through a Linear layer (linear transformation, including bias) and embedded in a latent space with a fixed number of dimensions. Since the node feature dimensions of proteins and ligands are initialized with different numbers of dimensions, a Linear layer is used to change them respectively, and the output dimension size is d f Similarly, the edge feature also uses a Linear layer to change, and the output dimension size is d e Since the initialization dimension of the edge features of the protein and ligand are the same, the same Linear layer is used for embedding. Then the LeakyReLU function (with leaky linear rectifier function) is used for a nonlinear transformation. The node features after output are recorded as The edge feature is recorded as

[0043] The feature update module mainly includes two operation steps, information aggregation and information feedforward. In information aggregation, for the lth feature update module, where l∈{1,...,N l Each central node updates its own features by aggregating the features and coordinate information of itself and its neighboring nodes. The specific way of information aggregation is achieved through multi-head attention. The query of attention is denoted as q, the key is denoted as k, and the value is denoted as v. The detailed description is as follows:

[0044]

[0045]

[0046]

[0047]

[0048] in and They represent the coordinates of the i-th node (i usually represents the central node number in the information aggregation of this section) and its j-th neighbor node (j usually represents the neighbor node number of the i-th node in the information aggregation of this section) in the l-1th module. In this paper, W is represented as a learnable linear layer with bias. For simplicity, its bias is omitted in the formula. Concat(·) represents a tensor concatenation operation. μ is a learnable weight. h represents the number of hidden layer dimensions of the h-th attention head, where h∈{1,...,N h}(The number of multi-head attention is N h ). In this work, d h =d f / N h ⊙ represents element-wise multiplication.

[0049] After the attention mechanism's q, k, and v are calculated, these variables are used to aggregate information, as follows:

[0050]

[0051]

[0052]

[0053]

[0054] in Neighbor(i) represents the set of neighboring nodes of the i-th node (i.e., the central node).

[0055] Before the information feedforward of the lth feature update module, a new element-wise gated residual connection (denoted as Ψ) is defined to enhance the robustness and learning ability of the model. For example, the output of the residual connection is expressed as:

[0056]

[0057]

[0058] in Norm(·) represents layer normalization.

[0059] Then After information feedforward, the updated node features of the lth feature update module are finally as follows:

[0060] in Similarly, the edge information feedforward steps are as follows:

[0061]

[0062] in Finally, the updated edge features of the lth feature update module are obtained

[0063] 2) Complex structure prediction

[0064] This paper proposes an accurate complex prediction method based on coordinate updating. This method combines a structure prediction module with a feature update module to update the 3D coordinates of nodes. Specifically, the attention mechanism calculated by the feature update module is converted into a single-dimensional distance gradient. These distance gradients represent the contribution of each adjacent node to the distance update of the central node. The updated coordinates of the central node are ultimately obtained by summing the distance gradients of each adjacent node, as expressed as follows:

[0065]

[0066]

[0067] where x i represents the 3D coordinate of the center node i, x j Represents the 3D coordinates of its adjacent node j. It is calculated by the formula in the feature update module. pass Perform a linear transformation (where ), and then multiply it by the direction modulus vector in the 3D space between the center node i and the neighboring node j to get the distance gradient contribution of the neighboring node j to the center node i. By adding the distance gradient contributions of all neighboring nodes (Neighbor(i)), the coordinate update change of the center node i is obtained Where i is not equal to j. Finally, for all N h Multi-head attention, the center node coordinate update changes are determined by the learnable weights The weighted sum is obtained.

[0068] 3) Affinity estimation

[0069] The present invention achieves affinity prediction by adding an additional affinity prediction module and assists the model in predicting protein-ligand structure. In this case, the final node feature is updated and output After a two-layer multilayer perceptron (MLP), a single-dimensional affinity prediction is obtained for each protein-ligand pair, which is expressed as follows:

[0070]

[0071] where y aff It is the affinity prediction for the protein-ligand pair. N f Indicates the number of graph nodes.

[0072] 4) Model training and inference:

[0073] In order to take into account the depth of the model after multiple cycles, the Monte Carlo strategy is used to train the model in each round of training. Specifically, for a total of N c Round loop, only one round is trained each time, no data calculation is performed for the loop after this round, and no gradient update is performed for the loop before this round. For example, in N c =3, each training round is selected from [1,3] for training. When the second round is selected for training, the first two rounds of forward data calculation will be performed normally, and the third round of data calculation will not be performed. When the network parameter weight gradient is updated, the parameter weight of the network is not calculated and the weight is not updated in the first round. The parameter weight is updated normally in the second round.

[0074] A total of 150 epochs are trained. The learning rate for the last 5 epochs is set to 1×10 -6 The Adam optimizer was used to update the model's parameter weights. The learning rate decayed exponentially with a constant of 0.99 per epoch. To ensure stability and efficiency in model training, protein-ligand structures and affinities were processed internally within the model using 1 / 10 of their original values ​​(e.g., coordinates in PDB files are in angstroms, while the units used in the model calculations are in nm). Affinities were calculated using the negative logarithm.

[0075] The loss function for predicting the protein-ligand complex structure is given a set of N ligIn each training, a sequence index is used to match the predicted atomic positions to their natural real positions. This is to ensure that the loss of the predicted atomic positions to the real positions is calculated on atoms with the same chemical identity. Afterwards, the loss of the structure prediction is obtained by calculating the average of the differences between all ligand atoms and their real natural positions using the above index sequence (this loss is recorded as L coor ).

[0076]

[0077] Where Nodes(lig) represents all nodes (atoms) of the ligand in the protein-ligand pair. However, many ligands have symmetrical structures. In such cases, for a given symmetric small molecule, many equivalent index sequences can be obtained.

[0078] Considering that different index sequences will result in different losses, in order to solve this problem, a unique symmetric-aware loss function is designed, and the loss calculated is recorded as L sym , used to correctly calculate the loss and predict the true structure, as follows:

[0079] It is assumed that the ligand contains N s An equivalent index sequence, Calculate L for each index sequence coor Loss, and the one with the lowest loss value is selected as the index sequence used to calculate the structure coordinate loss, where the equivalent index sequence is obtained by the RDKit tool.

[0080] For affinity prediction, the mean square error (MSE) is used as the loss function, which is expressed as follows:

[0081]

[0082] Among them L aff represents the loss of affinity prediction.

[0083] We use a Monte Carlo method for training, with only one cycle trained at a time. The following representations are all cycles with gradient updates. Finally, for each protein-ligand pair, the loss is expressed as follows:

[0084]

[0085] It is obtained by weighted summation of three parts: (1) L calculated from the output of each coordinate update module symLoss averaged, (2) L calculated from the predicted structure output by the final coordinate update module sym loss, (3) affinity loss, where γ1 and γ2 are constant parameters used to balance the loss of structure and the loss of affinity.

[0086] The ensemble method is used to predict the structure and affinity of protein-ligand, where the number of ensembles is recorded as N ens (Optional N ens =10).

[0087] The above description is merely a preferred embodiment of the present application, but the scope of protection of the present application is not limited thereto. Any changes or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in the present application should be included in the scope of protection of the present application. Therefore, the scope of protection of the present application should be based on the scope of protection of the claims.

Claims

1. A method for predicting the structure and estimating the affinity of a small molecule drug protein-ligand complex, characterized in that: The following steps are involved: constructing a protein-ligand graph, and extracting and updating features of the protein-ligand graph using a feature updating module; Combining the feature update module with the structure prediction module to perform complex structure prediction and obtain the protein-ligand three-dimensional structure; combining the feature update module with the affinity estimation module to perform affinity estimation and obtain the protein-ligand affinity; Constructing a protein-ligand network model, training the protein-ligand network model based on the protein-ligand three-dimensional structure and the protein-ligand affinity, adjusting network parameters, and completing prediction; The protein-ligand graph includes a complete graph and a heterogeneous graph; The complete graph comprises: constructing a bidirectional complete graph according to a given protein-ligand pair, wherein each node represents an atom, and each atom is connected to other atoms to form the complete graph; the complete graph comprises a set of vectors storing features of each node, and a set of matrices storing features of edges; The isomeric graph comprises: according to a given protein-ligand pair, representing the full connection of the protein and the ligand small molecule as the isomeric graph; The feature update module includes: information aggregation and information feedforward; The information aggregation includes: updating the features of the central node to obtain the aggregated features of the central node; aggregating the feedforward information received by each central node in the feature update module to obtain new features and coordinate information; The information feedforward includes: transmitting the features and coordinate information of the adjacent nodes around each central node in the feature update module to the central node; The affinity estimation process includes: The node feature update output is passed through a two-layer multilayer perceptron to obtain a single-dimensional affinity prediction for each protein-ligand pair, calculated as follows: where y aff is the affinity prediction for the protein-ligand pair, N f Indicates the number of graph nodes.

2. The method for predicting the structure and estimating the affinity of a small molecule drug protein-ligand complex according to claim 1, wherein: The complex structure prediction process includes: The attention mechanism calculated by the feature update module is converted into a single-dimensional distance gradient. The distance gradient represents the contribution of each adjacent node to the distance update of the central node. The updated coordinates of the central node are obtained by summing the distance gradients of each adjacent node, which is calculated as follows: where x i represents the 3D coordinate of the center node i, x j Represents the 3D coordinates of its neighboring node j.

3. The method for predicting the structure and estimating the affinity of a small molecule drug protein-ligand complex according to claim 1, wherein: The training process includes: using a Monte Carlo strategy to train the model, only training one round at a time, not performing data calculations on the cycles after this round, and not performing gradient updates on the cycles before this round, to complete the training.

Images