Target scoring function optimization method and system

Abstract

The present invention relates to a method for optimizing a target scoring function, comprising: a. selecting the corresponding organic small molecule type and quantity and performing mixed solvent molecular dynamics simulations; b. constructing a database of target-ligand complex structures and binding free energies for the corresponding targets based on the mixed solvent molecular dynamics simulation results; and c. using the constructed database to train the weights of the scoring function, assigning different weights to atoms near the target binding pocket. The present invention also relates to a target scoring function optimization system. The present invention fully considers the individual characteristics of the target and improves the screening performance of molecular docking.

Description

Technical Field

[0001] The present invention relates to the field of molecular dynamics, and in particular to a target scoring function optimization method and system. Background Art

[0002] Predicting the affinity of drugs (including small molecules and macromolecules such as peptides, antibodies, and proteins) for their targets plays a crucial role in drug screening. Molecular docking is a key virtual drug screening technique that consists of two parts: conformational search and conformational evaluation. First, a conformational generation algorithm is used to generate a series of ligand conformations. A scoring function is then used to calculate the binding free energy between the compound's conformation and the receptor. Finally, potential active molecules are selected based on the binding free energy ranking, or the binding conformation between the active molecule and the target is predicted.

[0003] In addition to the machine learning-based scoring functions that have emerged in recent years, traditional scoring functions can be divided into three categories: physics-based, empirical, and statistical. Physics-based scoring functions are based on force fields, so the potential energy function and parameters of the force field determine the magnitude of the predicted value. The prediction mechanism is to accumulate van der Waals and electrostatic interactions between receptor-ligand atomic pairs, solvation effects, and the torsional energy of the ligand. Physics-based scoring functions include DOCK and AutoDock4. Empirical scoring functions calculate binding free energy by accumulating important energy factors, such as hydrogen bonds, hydrophobic effects, and steric clashes. When constructing empirical scoring functions, a training library with known activity is used to optimize the parameters of each energy term through linear regression. Representative methods include X-Score, AutoDock Vina, and Vinardo. Statistical scoring functions assume that the frequency of each atom pair in a large number of receptor-ligand complexes is related to the interaction between the two atoms. This frequency is then converted into a distance-dependent mean force potential. The potential energy of each receptor-ligand atom pair is then accumulated to produce a predicted binding free energy. While the principles of these three traditional scoring functions differ, their specific composition is the sum of the product of each weight and the corresponding term, representing the universal law of drug-target interactions understood by the scoring function.

[0004] Molecular dynamics simulation (MDS) is a computational framework that uses classical molecular mechanics, which treats atoms as solid spheres obeying Newton's laws of motion. As time changes, the forces acting on the atoms and their spatial coordinates change accordingly, allowing the entire system to be simulated. Molecular mechanics involves a set of parameters, known as force fields, that approximate the molecular potential energy surface based on the molecular topology. The force field consists of terms that describe both intramolecular and intermolecular interactions. Intramolecular interaction terms include bonding interactions between neighboring atoms (such as chemical bonds, bond angles, and dihedral angles) and non-bonding interactions within the molecule (such as Coulomb and van der Waals forces). Intermolecular interaction terms consist of non-bonded interactions between molecules. The total potential energy of the system is the sum of these energy terms. The force field divides atoms into different types based on their surrounding environments (e.g., carbon atoms in different chemical environments, such as alkyl, ether, or ester). Each atom type has corresponding force field parameters, allowing the vast majority of molecules in chemical space to be described with a small number of parameters. Compared to the classic protein-water solution molecular dynamics simulation, the solvent environment of the mixed solvent molecular dynamics simulation is an aqueous solution containing multiple organic solvents (probes). After the simulation, the grid points are used to count the number of times the probe appears on the protein surface. The binding site of the protein can be detected by the grid point with the highest value. The binding free energy of the target-probe can also be calculated using the grid point value and the inverse Boltzmann distribution formula (Formula 1). The physical meaning of this formula is the measure of the free energy change of moving an atom from the solvent to the grid point i. Where: N i is the actual frequency of the probe at grid point i; N0 is the expected frequency of the probe at grid point i if there is no interaction with the protein; R is the gas constant; T represents the absolute temperature.

[0005] ΔG i =-RT ln(N i / N0) Formula 1

[0006] The training library of the existing scoring function, namely the target-ligand complex structure and affinity data set, all comes from wet experiments, and the parameters of the scoring function are constant when applied to different targets.

[0007] When using a scoring function to predict affinity for a target, the prediction results will be more accurate if the weights can be adjusted based on the properties of the target itself. According to the general process, parameter adjustment requires a large amount of activity data for the corresponding target and compound. This creates a paradox: because the scoring function is applied to screen for active compounds, in this case there is little or no activity data for the corresponding target and compound, and parameter adjustment is impossible; if there is more activity data for the corresponding target and compound, parameter adjustment can be made, but this also indicates that the drug development level of this target is very sufficient, and there is no need to use a scoring function for initial screening. As a result, there is currently no method or system that can adjust the weights of the scoring function according to the target structure. Summary of the Invention

[0008] In view of this, it is necessary to provide a target scoring function optimization method and system.

[0009] The present invention provides a target scoring function optimization method, which includes the following steps: a. selecting corresponding organic small molecule types and quantities and performing mixed solvent molecular dynamics simulation; b. constructing a target-ligand complex structure and binding free energy database of the corresponding target based on the mixed solvent molecular dynamics simulation results; c. using the constructed database to train the weights of the scoring function and assign different weights to atoms near the target binding pocket.

[0010] Preferably, the step a comprises:

[0011] First, the target structure is processed before simulation to fill in missing atoms. Then, the complete structure of the target is used to perform pre-simulation preparation for mixed solvent molecular dynamics simulation, which includes adding water, probes, and ions. Finally, the topology file, parameter file, and configuration file for molecular dynamics simulation are obtained.

[0012] Preferably, the mixed solvent comprises an organic solvent; the organic solvent comprises: propanol, acetamide, isopropylamine, dimethyl sulfoxide, trifluoroethane, tert-butane, acetic acid, guanidine salt, imidazole and benzene.

[0013] Preferably, the ratio of the number of probes to the number of water molecules is 1:20.

[0014] Preferably, the step b specifically includes:

[0015] First, the trajectory is shifted and superimposed based on the protein's center of mass to generate a simulation trajectory containing only the protein and the probe. Then, the probe frequency is counted based on the grid point. Based on the receptor-probe interaction frequency and Equation 1, the binding free energy of each probe with the receptor at different positions is calculated.

[0016] ΔG i=-RT ln(N i / N0) Formula 1

[0017] Where: N i is the actual frequency of the probe at grid point i; N0 is the expected frequency of the probe at grid point i if there is no interaction with the protein; R is the gas constant; T represents the absolute temperature.

[0018] Preferably, the step b further comprises:

[0019] Screening of probe-receptor complexes:

[0020] The screening process involves two steps: first, calculating the binding free energy of the corresponding probe-receptor based on the grid point values ​​while estimating the binding free energy of the complex using the average values ​​of the surrounding grid points. When both values ​​are less than -1.2 kcal / mol, the complex is extracted from the trajectory. Second, VMD clusters the conformations of the probes based on RMSD, with a cutoff of 0.8. After the clustering results are obtained, conformations that are not clustered or have less than 4 clusters are removed. For conformations with more than 6 clusters, up to 6 conformations are randomly retained, and the remaining conformations are saved in PDB format.

[0021] Preferably, the step c specifically includes:

[0022] First, the topology files corresponding to the receptor and probe in each complex were converted from PDB format to PDBQT format. Then, the scores of each atom of the target on the probe molecule in hydrogen bonding, van der Waals, solvent effect and Coulomb interaction were calculated to obtain the input features. The MSMD binding free energy prediction value was used as the label to obtain the training set. Then, JAX was used to build the training model and train the parameters according to the training set.

[0023] Preferably, the initial values ​​of the parameters are not randomly generated, but are weights of the scoring function.

[0024] Preferably, the step c further comprises:

[0025] The optimized parameters were checked to see the degree of change relative to the initial values ​​and the number of changes in each degree. If the change was large or negative, the parameters were re-optimized. If the degree of change and the number of changes were within a reasonable range, the optimized parameters were obtained and used for subsequent molecular docking.

[0026] The present invention provides a target scoring function optimization system, which includes a simulation module, a construction module, and an allocation module, wherein: the simulation module is used to select corresponding organic small molecule types and quantities and perform mixed solvent molecular dynamics simulation; the construction module is used to construct a target-ligand complex structure and binding free energy database of the corresponding target based on the mixed solvent molecular dynamics simulation results; and the allocation module is used to use the constructed database to train the weights of the scoring function and assign different weights to atoms near the target binding pocket.

[0027] Instead of relying on real-world activity data, the present invention adjusts the scoring function based on simulated data derived from the corresponding target structure. Using a database of target-organic small molecule affinities derived from molecular dynamics simulations, weights are adjusted at the atomic level to optimize the scoring function. This fully considers the individual characteristics of the target, improving the screening performance of molecular docking. BRIEF DESCRIPTION OF THE DRAWINGS

[0028] Figure 1 Flowchart of the target scoring function optimization method of the present invention;

[0029] Figure 2 This is a hardware architecture diagram of the target scoring function optimization system of the present invention. DETAILED DESCRIPTION

[0030] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0031] See Figure 1 , which is a flowchart of a preferred embodiment of the target scoring function optimization method of the present invention.

[0032] Step S1: Select the corresponding organic small molecule type and quantity and perform mixed solvent molecular dynamics simulation. Specifically:

[0033] First, the target structure is processed before simulation to fill in missing atoms. Then, the complete structure of the target is used to perform pre-simulation preparation for mixed solvent molecular dynamics simulation, which includes adding water, probes, and ions. Finally, the topology file, parameter file, and configuration file for molecular dynamics simulation are obtained.

[0034] Adjust the selection and ratio of probe types in mixed solvent molecular dynamics simulation within an appropriate range. In this embodiment, 10 organic solvents were selected; the organic solvents included propanol, acetamide, isopropylamine, dimethyl sulfoxide, trifluoroethane, tert-butane, acetic acid, guanidine salt, imidazole and benzene, and the ratio of the number of probes to water molecules was 1:20. Since hydrophobic probes such as tert-butane, trifluoroethane, and benzene will aggregate during the simulation, the hydrophobic probes will collectively aggregate on one side of the solution and cannot freely explore the protein surface and its pockets. By affecting the frequency of occurrence N i This affects the calculation results of the inverse Boltzmann distribution formula. To improve the exploration capability of the hydrophobic probes, this embodiment reduces the concentration of the hydrophobic probes in the solution while applying an additional force to the aggregated hydrophobic probes to prevent the aggregation of the hydrophobic probes.

[0035] Step S2: Based on the results of the mixed solvent molecular dynamics simulation, a database of target-ligand complex structures and binding free energies corresponding to the target is constructed. Specifically:

[0036] The greatest highlight and uniqueness of the present invention is that it uses mixed solvent molecular dynamics simulation to provide training samples for parameter optimization of the scoring function.

[0037] After the mixed-solvent molecular dynamics simulation, the trajectory is first shifted and superimposed based on the protein's center of mass to generate a simulation trajectory containing only the protein and probe. This is because excluding unnecessary water molecules can speed up subsequent calculations. A grid-based probe frequency statistics are then performed. The center of mass of the ligand in the protein crystal structure is used as the central coordinate of the statistical grid. The statistical length of each dimension is 24 angstroms, and the length of each grid point is set to 0.25 to facilitate accurate statistics. Based on the receptor-probe interaction frequency and Equation 1, the binding free energy of each probe at different positions on the receptor is calculated.

[0038] The number of probe-protein complexes extracted from the simulation trajectory is huge. In order to keep the training set at a reasonable size and representative, the probe-receptor complex is screened in this embodiment. The screening includes two operations: the first part is to calculate the binding free energy of the corresponding probe-receptor based on the grid value, and estimate the binding free energy of the complex using the numerical average of the surrounding grid points. When the above two values ​​are both less than -1.2kcal / mol, the complex is extracted from the trajectory; the second part is to use VMD to cluster the conformations of the probe according to RMSD, with a cutoff value of 0.8. After obtaining the clustering results, the conformations that are not clustered or have less than 4 in each class are removed. For the conformations with more than 6 in each class, 6 conformations are randomly retained, and the remaining conformations are saved in PDB format.

[0039] Step S3: Use the constructed database to train the weights of the scoring function and assign different weights to atoms near the target binding pocket. Specifically:

[0040] This embodiment selects the scoring function of AutoDock4 as the optimization object. The optimization method of MSMD for scoring functions is universal. Although this embodiment only takes AutoDock4 as an example, it should be noted that the optimization processes of other scoring functions are similar.

[0041] First, the topology files corresponding to the receptor and probe in each complex were converted from PDB format to PDBQT format. Then, the scores of each atom of the target on the probe molecule in hydrogen bonding, van der Waals, solvent effect and Coulomb interaction were calculated to obtain the input features. The MSMD binding free energy prediction value was used as the label to obtain the training set. Then, JAX was used to build the training model and train the parameters according to the training set.

[0042] It is worth noting that the initial values ​​of the parameters are not randomly generated, but the weights of the AutoDock4 scoring function. The loss function consists of two parts: the first part is the sum of the absolute values ​​of the differences between the predicted values ​​and the label values; the second part is the sum of the change ratios of each parameter compared to the initial values. The role of the second part of the loss function is to prevent the parameters from undergoing huge changes that do not conform to physical meaning. Such a combination can enable the scoring function to maintain the physical meaning of each scoring item while improving the prediction performance. Compared to AutoDock4, the parameters of the van der Waals, hydrogen bond, dissociation and Coulomb force terms for all atoms are the same. The parameter optimization operation of this embodiment assigns independent parameters to each atom in the protein PDBQT file, so that the atoms near the binding pocket are given different weights, which better reflects the binding tendency of the receptor-ligand.

[0043] The optimized parameters were checked to see the degree of change relative to the initial values ​​and the number of changes in each degree. If the change was large or negative, the parameters were re-optimized. If the degree of change and the number of changes were within a reasonable range, the optimized parameters were obtained and used for subsequent molecular docking.

[0044] See Figure 2 , which is a hardware architecture diagram of the target scoring function optimization system 10 of the present invention.

[0045] The system includes: a simulation module 101, a construction module 102, and an allocation module 103.

[0046] The simulation module 101 is used to select the corresponding organic small molecule type and quantity to perform mixed solvent molecular dynamics simulation. Specifically:

[0047] First, the target structure is processed before simulation to fill in missing atoms. Then, the complete structure of the target is used to perform pre-simulation preparation for mixed solvent molecular dynamics simulation, which includes adding water, probes, and ions. Finally, the topology file, parameter file, and configuration file for molecular dynamics simulation are obtained.

[0048] Adjust the selection and ratio of probe types in mixed solvent molecular dynamics simulation within an appropriate range. In this embodiment, 10 organic solvents were selected; the organic solvents included propanol, acetamide, isopropylamine, dimethyl sulfoxide, trifluoroethane, tert-butane, acetic acid, guanidine salt, imidazole and benzene, and the ratio of the number of probes to water molecules was 1:20. Since hydrophobic probes such as tert-butane, trifluoroethane, and benzene will aggregate during the simulation, the hydrophobic probes will collectively aggregate on one side of the solution and cannot freely explore the protein surface and its pockets. By affecting the frequency of occurrence N i This affects the calculation results of the inverse Boltzmann distribution formula. To improve the exploration capability of the hydrophobic probes, this embodiment reduces the concentration of the hydrophobic probes in the solution while applying an additional force to the aggregated hydrophobic probes to prevent the aggregation of the hydrophobic probes.

[0049] The construction module 102 is used to construct a target-ligand complex structure and binding free energy database of the corresponding target based on the mixed solvent molecular dynamics simulation results. Specifically:

[0050] The greatest highlight and uniqueness of the present invention is that it uses mixed solvent molecular dynamics simulation to provide training samples for parameter optimization of the scoring function.

[0051] After the mixed-solvent molecular dynamics simulation, the trajectory is first shifted and superimposed based on the protein's center of mass to generate a simulation trajectory containing only the protein and probe. This is because excluding unnecessary water molecules can speed up subsequent calculations. A grid-based probe frequency statistics are then performed. The center of mass of the ligand in the protein crystal structure is used as the central coordinate of the statistical grid. The statistical length of each dimension is 24 angstroms, and the length of each grid point is set to 0.25 to facilitate accurate statistics. Based on the receptor-probe interaction frequency and Equation 1, the binding free energy of each probe at different positions on the receptor is calculated.

[0052] The number of probe-protein complexes extracted from the simulation trajectory is huge. In order to keep the training set at a reasonable size and representative, the probe-receptor complex is screened in this embodiment. The screening includes two operations: the first part is to calculate the binding free energy of the corresponding probe-receptor based on the grid value, and estimate the binding free energy of the complex using the numerical average of the surrounding grid points. When the above two values ​​are both less than -1.2kcal / mol, the complex is extracted from the trajectory; the second part is to use VMD to cluster the conformations of the probe according to RMSD, with a cutoff value of 0.8. After obtaining the clustering results, the conformations that are not clustered or have less than 4 in each class are removed. For the conformations with more than 6 in each class, 6 conformations are randomly retained, and the remaining conformations are saved in PDB format.

[0053] The allocation module 103 is used to train the weights of the scoring function using the constructed database, and to allocate different weights to atoms near the target binding pocket.

[0054] Specifically:

[0055] This embodiment selects the scoring function of AutoDock4 as the optimization object. The optimization method of MSMD for scoring functions is universal. Although this embodiment only takes AutoDock4 as an example, it should be noted that the optimization processes of other scoring functions are similar.

[0056] First, the topology files corresponding to the receptor and probe in each complex were converted from PDB format to PDBQT format. Then, the scores of each atom of the target on the probe molecule in hydrogen bonding, van der Waals, solvent effect and Coulomb interaction were calculated to obtain the input features. The MSMD binding free energy prediction value was used as the label to obtain the training set. Then, JAX was used to build the training model and train the parameters according to the training set.

[0057] It is worth noting that the initial values ​​of the parameters are not randomly generated, but the weights of the AutoDock4 scoring function. The loss function consists of two parts: the first part is the sum of the absolute values ​​of the differences between the predicted values ​​and the label values; the second part is the sum of the change ratios of each parameter compared to the initial values. The role of the second part of the loss function is to prevent the parameters from undergoing huge changes that do not conform to physical meaning. Such a combination can enable the scoring function to maintain the physical meaning of each scoring item while improving the prediction performance. Compared to AutoDock4, the parameters of the van der Waals, hydrogen bond, dissociation and Coulomb force terms for all atoms are the same. The parameter optimization operation of this embodiment assigns independent parameters to each atom in the protein PDBQT file, so that the atoms near the binding pocket are given different weights, which better reflects the binding tendency of the receptor-ligand.

[0058] The optimized parameters were checked to see the degree of change relative to the initial values ​​and the number of changes in each degree. If the change was large or negative, the parameters were re-optimized. If the degree of change and the number of changes were within a reasonable range, the optimized parameters were obtained and used for subsequent molecular docking.

[0059] This application proposes a scoring function optimization method that uses molecular dynamics simulation data to construct a target-organic small molecule affinity database, thereby adjusting the weights of the scoring function and improving the molecular docking performance. This application is applicable to the optimization of various classic scoring functions.

[0060] Although the present invention has been described with reference to the current preferred embodiments, those skilled in the art should understand that the above-mentioned preferred embodiments are only used to illustrate the present invention and are not used to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the scope of protection of the present invention.

Claims

1. A target scoring function optimization method, characterized in that: The method comprises the following steps: a. Select the corresponding organic small molecule type and quantity and perform mixed solvent molecular dynamics simulation; b. Construct a database of target-ligand complex structures and binding free energies for the corresponding targets based on the results of mixed solvent molecular dynamics simulations; c. Use the constructed database to train the weights of the scoring function and assign different weights to atoms near the target binding pocket.

2. The method according to claim 1, wherein Described step a comprises: First, the target structure is processed before simulation to fill in missing atoms; then, the complete structure of the target is used to perform pre-simulation preparation for mixed solvent molecular dynamics simulation, which includes adding water, probes and ions; finally, the topology file, parameter file and configuration file for molecular dynamics simulation are obtained.

3. The method according to claim 2, wherein The mixed solvent includes an organic solvent; the organic solvent includes: propanol, acetamide, isopropylamine, dimethyl sulfoxide, trifluoroethane, tert-butane, acetic acid, guanidine salt, imidazole and benzene.

4. The method according to claim 3, wherein The ratio of the number of probes to water molecules is 1:

20.

5. The method according to claim 4, wherein The step b specifically includes: First, the trajectory is shifted and superimposed based on the protein's center of mass to generate a simulation trajectory containing only the protein and the probe. Then, the probe frequency is counted based on the grid point. Based on the receptor-probe interaction frequency and Equation 1, the binding free energy of each probe with the receptor at different positions is calculated. ΔG i =-RT ln(N i / N0) Formula 1 Where: N i is the actual frequency of the probe appearing at grid point i; N0 is the expected frequency of the probe appearing at grid point i if there is no interaction with the protein; R is the gas constant; and T represents the absolute temperature.

6. The method according to claim 5, wherein Described step b also comprises: Screening of probe-receptor complexes: The screening process involves two steps: first, calculating the binding free energy of the corresponding probe-receptor based on the grid point values ​​while estimating the binding free energy of the complex using the average values ​​of the surrounding grid points. When both values ​​are less than -1.2 kcal / mol, the complex is extracted from the trajectory. Second, VMD clusters the conformations of the probes based on RMSD, with a cutoff of 0.

8. After the clustering results are obtained, conformations that are not clustered or have less than 4 clusters are removed. For conformations with more than 6 clusters, up to 6 conformations are randomly retained, and the remaining conformations are saved in PDB format.

7. The method according to claim 6, wherein The step c specifically includes: First, the topology files corresponding to the receptor and probe in each complex were converted from PDB format to PDBQT format. Then, the scores of each atom of the target on the probe molecule in hydrogen bonding, van der Waals, solvent effect and Coulomb interaction were calculated to obtain the input features. The MSMD binding free energy prediction value was used as the label to obtain the training set. Then, JAX was used to build the training model and train the parameters according to the training set.

8. The method according to claim 7, wherein The initial values ​​of the parameters are not randomly generated, but are the weights of the scoring function.

9. The method according to claim 8, wherein Described step c also comprises: The optimized parameters were checked to see the degree of change relative to the initial values ​​and the number of changes in each degree. If the change was large or negative, the parameters were re-optimized. If the degree of change and the number of changes were within a reasonable range, the optimized parameters were obtained and used for subsequent molecular docking.

10. A target scoring function optimization system, characterized in that: The system includes a simulation module, a construction module, and an allocation module, wherein: The simulation module is used to select the corresponding organic small molecule type and quantity to perform mixed solvent molecular dynamics simulation; The building module is used to construct a target-ligand complex structure and binding free energy database of the corresponding target based on the mixed solvent molecular dynamics simulation results; The assignment module is used to train the weights of the scoring function using the constructed database, and assign different weights to atoms near the target binding pocket.

Images