Method for calculating absolute binding free energy without precise complex structure

Abstract

The application discloses a kind of absolute binding free energy calculation method without accurate complex structure, comprising: step one, based on acceptor protein and ligand structure, using software to obtain rough complex structure, and with Wang Landau simulation method, with the molecular structure system of ligand and protein, the potential energy function of intermolecular interaction as main simulation parameter as input to carry out molecular dynamics simulation, obtain the trajectory file of complex system molecular dynamics;Step two, the high-dimensional vector information corresponding to the trajectory file data establishes training set, validation set and test set, trains the flow model to find the lowest free energy conformation, i.e. the structure with the highest probability density in training result is selected as the optimal structure;Step three, the optimal structure generated is used as input file, and binding free energy is calculated.The application proposes a kind of absolute binding free energy calculation method without accurate complex structure, solves the problem that accurate complex structure is difficult to obtain.

Description

Technical Field

[0001] This invention relates to the field of computational biology, specifically to a method for calculating absolute binding free energy without requiring precise complex structure. Background Technology

[0002] The current new drug development market mainly relies on affinity experiments between drug molecules and receptor proteins to preliminarily determine the drug's efficacy. However, simple virtual screening combined with molecular docking is too approximate and has a low target hit rate. Existing methods for obtaining the structure of protein-ligand complexes commonly use molecular docking software, all based on structure searches using empirical scoring functions. Current methods for calculating binding free energy require precise protein-ligand complex structures, which are very difficult to obtain. For example, patent application CN201911135246.9 discloses a method for predicting the binding free energy of proteins and ligands based on a progressive neural network. This method transforms the structural information of proteins and ligands into one-dimensional tensors, establishes training, validation, and test sets, uses the training set to train the progressive neural network, optimizes and finds hyperparameters for prediction, and calculates the binding free energy for comparison with molecular docking results. This solves the technical problem of converting the three-dimensional structure of protein and ligand molecules into a computationally comprehensible tensor and inputting it into the progressive neural network for training and optimization, further accelerating the calculation speed and improving prediction accuracy. These methods use molecular docking software to obtain the complex structure, and the calculation of the binding free energy is highly dependent on the initial structure. An accurate protein-ligand complex structure is required; otherwise, the calculation results will be poor. Summary of the Invention

[0003] One object of the present invention is to solve at least the above-mentioned problems and to provide at least the advantages that will be described later.

[0004] Another objective of this invention is to provide a method for calculating the absolute binding free energy without requiring a precise complex structure. This method uses molecular docking software to obtain a rough complex structure of the receptor protein and ligand. The WangLandau simulation method is then used, employing the molecular structure of the ligand and protein, the potential energy function of intermolecular interactions, temperature, and pressure as the main simulation parameters for molecular dynamics simulation. This yields a trajectory file of the complex system's molecular dynamics. The high-dimensional vector information corresponding to the trajectory file data is used to divide the data into training, validation, and test sets. The flow model is then trained to find the conformation with the lowest free energy. The structure with the highest probability density in the training results is selected as the optimal structure. This optimal structure is then used as the input file for calculating the binding free energy.

[0005] To achieve these and other advantages according to the present invention, a method for calculating the absolute binding free energy without requiring precise composite structure is provided, comprising:

[0006] Step 1: Based on the individual structures of the receptor protein and ligand, a rough complex structure is obtained using molecular docking software. The WangLandau simulation method is then used, with the molecular structure system of the ligand and protein, the potential energy function of intermolecular interactions, and temperature and pressure as the main simulation parameters, to perform molecular dynamics simulation and obtain the trajectory file of the molecular dynamics of the complex system.

[0007] Step 2: Based on the high-dimensional vector information corresponding to the trajectory file data, establish a training set, a validation set, and a test set. Use a flow model to perform density estimation. Use the atomic coordinates of the complex structure in the trajectory file as the training set in the flow model training process. Train the flow model to find the lowest free energy conformation. That is, select the structure with the highest probability density in the training results of the training process as the optimal structure.

[0008] Step 3: Use the generated optimal structure as the input file to perform the binding free energy calculation.

[0009] Preferably, in the WangLandau simulation method in step one, a rough docking structure of the receptor protein and ligand is generated by docking software, and the docking structure is enhanced by sampling to obtain as many complex structures as possible.

[0010] Preferably, the enhanced sampling method weakens the protein-ligand interaction through an algorithm, the main formula of which is Equation 1;

[0011]

[0012] Among them, U SS U represents the potential energy of the interaction between a protein and its ligand. PP U represents the potential energy inside a protein. PW U represents the potential energy between proteins and small molecules. WW The potential energy of the small molecule is represented by β, which is an adjustment parameter that can be updated within the range of (0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1).

[0013] Preferably, after obtaining the molecular dynamics trajectory in step one, the traj.dcd file and parameter file are obtained. The mdanalysis module is used to read the trajectory file and align the Cartesian coordinates of all conformations to the reference conformation.

[0014] Preferably, the process of training the streaming model consists of two steps. The first step is to optimize the Autoencoder based on the reconstruction loss function, reducing the high-dimensional vector to 2-10 dimensions. The second step is to maximize the logarithmic sum of the probability density of the training set data based on the principle of maximum likelihood estimation, and determine the probability density estimate based on the stability of the negative logarithmic sum of the probability density during the training process, thereby optimizing the streaming model.

[0015] Preferably, in step two, an additional flow model module is added to the Autoencoder to estimate the density of the high-dimensional vector corresponding to the complex structure and find the optimal structure based on the density distribution.

[0016] Preferably, the flow model module employs multiple Affinecoupling layers, which use the probability distribution x before the change as the input layer x and the probability distribution z after the change as the output layer z, and the calculation process is repeated 8 times.

[0017] Preferably, the training set is used to train the model, obtain set variables, and calculate the change in free energy along the set variables.

[0018] The present invention has at least the following beneficial effects:

[0019] 1. The flow model described above can be used to transform the complex distribution of the complex structure into an invertible function with an analytical distribution. Therefore, the probability density of the complex distribution can be inversely derived from the probability density of the transformed analytical distribution. By using the probability density estimation, the structure with the highest probability of occurrence can be found, thereby obtaining the optimal structure.

[0020] 2. A method for calculating the absolute binding free energy without requiring a precise complex structure is proposed, which solves the problem of difficulty in obtaining a precise complex structure.

[0021] Other advantages, objectives and features of the present invention will become apparent in part from the following description, and in part from those skilled in the art through study and practice of the invention. Attached Figure Description

[0022] Figure 1 This is a flowchart of a method for calculating the absolute binding free energy without requiring a precise complex structure, as described in one of the technical solutions of this invention.

[0023] Figure 2 This is a rough structural diagram of the composite structure in one of the technical solutions of the present invention;

[0024] Figure 3 This is a structural diagram of a deep learning model in one of the technical solutions of the present invention;

[0025] Figure 4This is a graph showing the relationship between the changes in variables of an invertible function. Detailed Implementation

[0026] The present invention will now be described in further detail with reference to the accompanying drawings, so that those skilled in the art can implement it based on the description.

[0027] like Figure 1 As shown, this invention provides a method for calculating the absolute binding free energy without requiring precise complex structure, including:

[0028] Step 1: Complex structure sampling based on Wang-LandauSolute scaling

[0029] S1 uses the open-source software Autodocvina to dock with the receptor protein structure and the ligand structure. It performs a preliminary search for possible binding sites and binding structures, evaluates each complex structure, and generates a score for each complex structure. The structure with the highest score is selected as the input file for the next step. The docking results generally include multiple conformations and their corresponding energy scores. Users can select the most likely binding conformation based on the scores to further analyze the interaction between the protein and the ligand.

[0030] S2 employs the Wang-Landau Solute simulation method to perform a 2-microsecond Wang-Landau Solute scaling molecular dynamics simulation on the complex system to obtain a molecular dynamics trajectory file. This trajectory file contains a vast amount of complex structures, such as... Figure 2 As shown, in step two, the optimal structure will be searched among these structures to serve as the input for the combined free energy calculation.

[0031] Step Two: As Figure 3 As shown, the optimal structure is found by combining a deep learning model. This invention estimates the density of complex high-dimensional vectors based on a flow model. The principle is as follows: for two spaces before and after the transformation of an invertible function, the integral of the probability density element remains unchanged. Therefore, if we can construct a distribution that transforms a complex distribution into one with an analytical expression (assuming a Gaussian distribution), the probability density of the complex distribution can be inversely derived from the probability density of the transformed Gaussian distribution. This property is also known as the "change of variable theorem." Since the invertible function is monotonic and the integral of the probability density is equal to 1, the integral on the unit element is equal before and after the change in f, as shown in Equation 1 below:

[0032] p(x)dx = π(z)dz Equation 1

[0033] Where p(x) and π(z) represent the probability density distributions of x and z, respectively, dx is a infinitesimal element in the variable space before the change by the invertible function f, and dz is a infinitesimal element in the variable space after the change by f. Since the invertible function is monotonic and the integral of the probability density is equal to 1, that is, the integral on the unit infinitesimal element is equal before and after the change by f, Equation 1 can also be transformed into Equation 2:

[0034] p(x) = π(z)dz / dx = π(z)|J xz Equation 2

[0035] Where |J xz | That is, the Jacobian determinant of the function f, representing the partial derivative of the multivariate function. The subscript xz indicates that f is used to transform x to z. If z follows a standard normal distribution, then its probability density has an analytical expression π(z), and p(x) can be calculated using the above expression. The prerequisite for calculating p(x) is to construct an invertible function f. The flow model aims to use the neural network as f to meet the above requirements. It has the function of transforming the complex data distribution p(x) into a standard normal distribution and is convenient for solving the Jacobian determinant. The flow model uses special structures, such as orthogonal matrices or triangular matrices, to simplify the calculation of the Jacobian matrix and make the training of the entire model more efficient.

[0036] Step 3: Using the structure generated in the previous step as the input file, use the software BFEE2 to calculate the binding free energy. This software can automatically generate the preparation files required for the simulation and perform unified processing on the generated files to calculate the binding free energy.

[0037] The absolute binding free energy calculation method provided by this technical solution, which does not require a precise complex structure, obtains a rough complex structure of the receptor protein and ligand using molecular docking software. The molecular structure system of the ligand and protein, and the potential energy function of intermolecular interactions are used as the main simulation parameters. Molecular dynamics simulation is performed using the WangLandau simulation method to obtain a trajectory file of the complex system's molecular dynamics. The high-dimensional vector information corresponding to the trajectory file data is used to divide the data into training, validation, and test sets. The flow model is trained to find the conformation with the lowest free energy; that is, the structure with the highest probability density in the training results is selected as the optimal structure. The generated optimal structure is used as the input file for binding free energy calculation.

[0038] The present invention has at least the following beneficial effects:

[0039] 1. The flow model described above can be used to transform the complex distribution of the complex structure into an invertible function with an analytical distribution. Therefore, the probability density of the complex distribution can be inversely deduced from the probability density of the transformed analytical distribution. By using the probability density estimation, the structure with the highest probability of occurrence can be found, thereby obtaining the optimal structure.

[0040] 2. A method for calculating the absolute binding free energy without requiring a precise complex structure is proposed, which solves the problem of difficulty in obtaining a precise complex structure.

[0041] In another technical solution, such as Figure 2 As shown, in the absolute binding free energy calculation method that does not require precise complex structure, in the WangLandau simulation method in step one, a rough docking structure of the receptor protein and ligand is generated using CHARMM-GUI docking software. The docking structure is then enhanced by sampling to obtain as many complex structures as possible for modeling and input file preparation. A new algorithm is designed to weaken the protein-ligand interaction. This algorithm is named Wang-Landau Solute scaling. The main formula of the algorithm is shown in Equation 3 below, where U... SS U represents the potential energy of the interaction between a protein and its ligand. PP U represents the potential energy inside a protein. PW U represents the potential energy between proteins and small molecules. WW Potential energy representing small molecules;

[0042]

[0043] In another technical solution, the method for calculating the absolute binding free energy without requiring precise complex structure involves an enhanced sampling method that weakens the protein-ligand interaction through an algorithm. An adjustment parameter β is introduced, as shown in Equation 3. This parameter β can be updated within the range of (0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1). By adjusting parameter β, the relative motion of the protein and ligand can be locally heated, equivalent to a solute temperature between 300K and 900K, with a temperature change interval of 60K. The weight of β is updated every 10,000 steps, meaning that each β update results in a different potential energy for the interaction between the ligand and protein, leading to the collection of different structures.

[0044] In another technical solution, after obtaining the molecular dynamics trajectory in step one, the traj.dcd file and parameter file are obtained. The mdanalysis module is used to read the trajectory file. Since the structures of multiple complexes need to be compared, one of them is selected as the reference conformation, and the Cartesian coordinates of all conformations are aligned to the reference conformation to analyze their structural similarities or differences, so as to eliminate the influence of translation and rotation on data analysis.

[0045] In another technical solution, such as Figure 3 As shown, in the method for calculating the absolute binding free energy without requiring a precise complex structure, the training process of the streaming model consists of two steps. The first step is to optimize Autoencoder8 based on the reconstruction loss function, reducing the high-dimensional vector to 2-10 dimensions. The second step is to maximize the logarithmic sum of the probability densities of the training set data based on the principle of maximum likelihood estimation. The probability density estimate is determined based on the stability of the negative logarithmic sum of the probability densities during training. By reducing the loss, the streaming model is optimized to more accurately capture the probability distribution of the data, thereby improving the model's generation and inference performance. Furthermore, when training the streaming model, the backpropagation algorithm is typically used to calculate the gradient of the Jacobian matrix. This requires calculating the Jacobian matrix and its gradient for each invertible function, and then propagating this information back to the entire model using the chain rule to update parameters and optimize the model.

[0046] In another technical solution, such as Figure 3 As shown, in the method for calculating the absolute binding free energy without requiring a precise complex structure, step two adds a flow model module to Autoencoder8 to estimate the density of the high-dimensional vector corresponding to the complex structure. Based on the density distribution, the structure with the highest probability is found. Since Autoencoder itself does not model the probability distribution of the data, adding a flow model module can introduce more powerful probability modeling capabilities and flexibility in generating samples.

[0047] In another technical solution, such as Figure 3 As shown, in the method for calculating the absolute binding free energy without requiring a precise complex structure, the flow model module employs multiple Affinecoupling layers. The probability distribution x before the change is used as the input layer x. After copying and calculating the input layer x, the first layer output is obtained. After copying and calculating again, the probability distribution z is used as the output layer z. The calculation process is repeated 8 times, which ensures the nonlinearity of the neural network and facilitates the calculation of probability density.

[0048] In another technical solution, such as Figure 3As shown, in the method for calculating the absolute binding free energy without requiring a precise complex structure, the training set is used during the model training process, and the set variables are obtained during the model optimization process. The free energy change along the set variables is calculated, and the free energy change along the set variables is plotted so as to find the structure corresponding to the point of lowest free energy in the graph, which is the optimal structure.

[0049] The specific implementation method of this invention first uses Autodock-Vina for docking to obtain a rough structure, and then performs simulation based on this rough structure. First, this structure is minimized and optimized. Then, WangLandau simulation is performed on the bond lengths and bond angles to obtain a 2µs trajectory file containing all complex structures. The trajectory file type is a DCD file. The subsequent detailed operation steps are as follows:

[0050] 1. Extract Cartesian coordinates from the trajectory: First, extract the Cartesian coordinates of the protein backbone and small molecules in each frame from the DCD trajectory file, and save the coordinates in the form of a three-dimensional array with the shape of [simulated frame number, total number of atoms, 3];

[0051] 2. Alignment Trajectory: Selecting the first frame of the array as the reference conformation, and calculating the displacement vector and rotation matrix of the remaining frames based on minimizing the mean square deviation of the protein backbone. After calculation, the remaining frames are displacementd and rotated to align the proteins, resulting in an array with a consistent shape.

[0052] 3. Extract small molecule coordinates: Extract the coordinates of small molecules from the aligned array. The array shape is [simulation frame number, number of small molecule atoms, 3].

[0053] 4. Input the extracted small molecule coordinates into the neural network model: Taking the shape of the small molecule coordinate array as [100000, 33, 3] as an example, in the first layer of the Autoencoder, the Cartesian coordinates of each frame are completely flattened into [100000, 99]. 5. The flattened vectors are passed through the Dense layers in the Autoencoder. The dimension of the output is determined by the number of neurons in the Dense layer. For example, the first Dense layer has 100 neurons, so the shape of the output array is [100000, 100]. The other Dense layers are: 80, 60, 30, 60, 80, 100, 99. Finally, 99 is arranged into [33, 3] to calculate the mean squared error with the original input to optimize dimensionality reduction.

[0054] 6. The intermediate layer of the Autoencoder is also connected to the flow model for density estimation. Since density estimation requires the calculation of the Jacobian determinant, the dimensionality does not change in the intermediate process and remains constant at the dimension of the Autoencoder bottleneck layer, which is 30.

[0055] 7. After the model training is completed, the density can be estimated by the theory of the flow model to obtain the probability density array of all small molecules. The array shape is [100000,1], where 100000 represents each small molecule and 1 represents that the probability density is a constant. 8. Extract the structure corresponding to the highest probability density and use BFEE2 to calculate the binding free energy.

[0056] The number of devices and processing scale described herein are for the purpose of simplifying the description of the invention. Applications, modifications, and variations of the invention will be readily apparent to those skilled in the art.

[0057] Although embodiments of the present invention have been disclosed above, they are not limited to the applications listed in the specification and embodiments. They can be applied to various fields suitable for the present invention. For those skilled in the art, other modifications can be easily made. Therefore, without departing from the general concept defined by the claims and their equivalents, the present invention is not limited to the specific details and illustrations shown and described herein.

Claims

1. A method for calculating absolute binding free energy without the need for precise complex structures, characterized in that, The application relates to a method for predicting the structure of a protein-ligand complex. Step one: based on the structures of a receptor protein and a ligand, a coarse complex structure is obtained by using a molecular docking software; a molecular dynamics simulation is performed by using a WangLandau simulation method, taking the molecular structure system of the ligand and the protein, the potential energy function of the intermolecular interaction, and the temperature and pressure as input simulation parameters, so as to obtain a trajectory file of the molecular dynamics of the complex system; Step two: a training set, a verification set and a test set are established based on high-dimensional vector information corresponding to the trajectory file data; a flow model is used for density estimation; the atomic coordinates of the complex structure in the trajectory file are taken as the training set in the training process of the flow model, and the flow model is trained to find the lowest free energy conformation, that is, the structure with the highest probability density in the training result of the training process is selected as the optimal structure; Step three: the generated optimal structure is taken as an input file, and binding free energy calculation is performed.

2. The method of claim 1, wherein the absolute binding free energy of a complex is calculated without the need for an accurate complex structure. In the WangLandau simulation method process in step one, a coarse structure docking structure of a receptor protein and a ligand is generated by using a docking software, and the docking structure is subjected to enhanced sampling to obtain a complex structure.

3. The method of claim 2, wherein the absolute binding free energy of a complex is calculated without the need for an accurate complex structure. The enhanced sampling method weakens the interaction between the protein and the ligand by using an algorithm, and the formula of the algorithm comprises formula 1. Formula 1 wherein, represents the potential energy of the interaction of the protein with the ligand, represents the potential energy within the protein, represents the potential energy between the protein and the small molecule, represents the potential energy of the small molecule, is a tuning parameter, said tuning parameter being updated in the range (0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1).

4. The method of claim 3, wherein the absolute binding free energy of a complex is calculated without the need for an accurate complex structure. After obtaining the molecular dynamics trajectory in step one, a traj.dcd file and a parameter file are obtained; the mdanalysis module is used to read the trajectory file, and the Cartesian coordinates of all conformations are aligned to a reference conformation.

5. The method of claim 1, wherein the absolute binding free energy of a complex structure is calculated without the need for an accurate complex structure, and wherein the absolute binding free energy of a complex structure is calculated without the need for an accurate complex structure. The training process of the flow model comprises two steps: step one is to optimize an Autoencoder based on a reconstruction loss function, so as to reduce the high-dimensional vector to 2-10 dimensions; and step two is to maximize the sum of the logarithm of the probability density of the training set data based on the principle of maximum likelihood estimation, and the probability density estimation is determined according to the stability of the negative logarithm sum of the probability density in the training process, so as to optimize the flow model.

6. The method of claim 1, wherein the absolute binding free energy of a complex structure is calculated without the need for an accurate complex structure, and wherein the absolute binding free energy of a complex structure is calculated without the need for an accurate complex structure. In step two, a flow model module is additionally added on the basis of the Autoencoder, density estimation is performed on the high-dimensional vector corresponding to the complex structure, and the optimal structure is found according to the density distribution.

7. The method of claim 6, wherein the absolute binding free energy of a complex is calculated without the need for an accurate complex structure. The flow model module employs a plurality of Affine coupling layers which transform the probability distribution before the change x as input layer x , the probability distribution after the change z as output layer z The computation process is repeated 8 times.

8. The method of claim 1, wherein the absolute binding free energy of a complex structure is calculated without the need for an accurate complex structure, and wherein the absolute binding free energy of a complex structure is calculated without the need for an accurate complex structure. In the training model process, the training set is used for training, a set variable is obtained, and the free energy change along the set variable is calculated.

Images