A method for generating stable tautomers and protonation states of drug molecules

Abstract

The present invention discloses a method for generating stable tautomers and protonation states of drug molecules. This method traverses all possible isomeric states of the molecule through an isomer conversion template, scores each isomeric state using a machine learning potential energy function, and determines the structure of the low-energy tautomer of the molecule. Then, based on the small molecule dissociation constant prediction method, the dissociation constants of each potential dissociation site in each low-energy tautomer are predicted to generate all possible protonation states of each isomer. The present invention integrates four aspects: tautomer conversion, scoring, ranking, and protonation of drug molecules, and can efficiently process small molecule compounds.

Description

Technical Field

[0001] The present invention relates to the field of computer-aided drug design, and is a method capable of rapidly generating stable tautomers and protonation states of drug molecules in aqueous solution. It can be applied in the fields of drug virtual screening and drug molecule design, can shorten the calculation time of drug virtual screening, and can guide chemists to design small molecule drugs and optimize the activity of small molecule drugs. Background Art

[0002] Tautomerism is an important issue in drug discovery. Many types of drug molecules can exist in multiple tautomeric states under physiological conditions. Since drug-like molecules generally contain a large number of heterocycles, hydrogen protons can transfer in the heterocycles, affecting the interaction between the drug molecule and the target protein.

[0003] The tautomeric state distribution of drug molecules in aqueous solution affects many properties of drugs, such as membrane permeability, reactivity, toxicity, partition coefficient, bioavailability, and solubility. Hydrogen bond donors and acceptors are important pharmacophore features and play an important role in protein-ligand binding. The interconversion of tautomeric states can change a hydrogen bond donor to an acceptor and vice versa, and this conversion will seriously affect the binding activity between the drug molecule and the target. Efficient and rapid processing of the tautomeric states of drug molecules is crucial for virtual screening and the construction of chemical information databases. Therefore, developing a rapid and effective method to obtain low-energy tautomers of drug-like molecules in aqueous solution is of great value for computer-aided drug discovery.

[0004] Identifying stable tautomers in aqueous solution can be accomplished by tautomer enumeration and tautomer ranking. Currently, two computational methods are widely used for tautomer ranking. One type is the empirical scoring method, which ranks based on empirical rules summarized from experimental results and computational data. The other is the energy-based scoring method, which ranks based on the internal energy and solvation energy of tautomers calculated by quantum chemistry.

[0005] Compared with the empirical-based scoring method, the energy-based scoring method has higher accuracy and robustness. However, the energy-based scoring method requires a large amount of computing power and calculation time and is difficult to process large compound libraries. In recent years, machine learning potential energy functions have developed rapidly and been widely applied. Therefore, using machine learning to develop a tautomer scoring method can not only speed up the calculation but also achieve high calculation accuracy. Summary of the Invention

[0006] The object of the present invention is to provide a method for generating stable tautomers and protonation states of drug molecules. In this method, first, the molecule is cut into fragments of different sizes using the designed cutting rules, and all possible tautomeric structures are generated for each fragment based on the tautomer conversion template. Then, a machine learning-based energy scoring method is used to score and rank each tautomer to obtain the stable tautomeric structures of each fragment. Next, the stable-state tautomeric structures of all fragments are combined to obtain all stable-state tautomers of the complete molecule. Finally, the dissociation constant prediction method MolGpKa is used to calculate the dissociation constants of potential dissociation centers in each tautomer to generate all possible protonation states. This method can quickly generate stable tautomeric structures and possible protonation states in aqueous solution for large compound libraries, improve the efficiency of virtual screening of drug molecules, and can also guide medicinal chemists to optimize the structures of drug molecules.

[0007] The object of the present invention is achieved as follows:

[0008] A method for generating stable tautomers and protonation states of drug molecules, the method comprising the following specific steps:

[0009] Step 1: Extract one million compound structures from the pubchem database, optimize the molecular conformations of each compound using quantum chemical methods, and calculate the solvation energy of each compound in aqueous solution;

[0010] Step 2: Based on the molecular conformation and solvation energy data obtained in Step 1, construct a deep learning model MolSolv. The model MolSolv is a solvation energy prediction model that can quickly predict the solvation energy of small molecule drugs;

[0011] Step 3: When generating stable tautomers and protonation states of small molecule compounds in solution, first, based on the cutting rules, cut the small molecule compounds into molecular fragments, and find all molecular fragments that can undergo tautomer conversion as variable regions based on the variable region recognition template;

[0012] Step 4: For each variable region found in Step 3, based on the collected SMIRKS conversion templates of tautomers, use RDKit to traverse the tautomeric states of each variable region;

[0013] Step 5: For each tautomeric structure of the same variable region, use RDKit to generate up to 50 molecular conformations;

[0014] Step 6: Use the ANI-2x model to optimize each molecular conformation generated by RDKit in Step 5, and the optimized optimal molecular conformation is used as the input for the calculation of the internal energy and hydration solvation energy of the next molecule;

[0015] Step 7: For the molecular conformations optimized in Step 6, use the deep learning-based solvation energy model MolSolv to predict the solvation energy of the molecular conformations in aqueous solution, use ANI-2x to predict the internal energy of the molecular conformations, and sum the internal energy and solvation energy of the molecular conformations as the final energy of the optimized molecular conformations. To speed up the calculation, the thermodynamic correction term is ignored; the specific calculation formula is as follows: ΔG = ΔG internal + ΔG solvation , where ΔG internal represents the internal energy of this molecular conformation, and ΔG solvation represents the solvation energy of this molecular conformation in aqueous solution;

[0016] Step 8: Sort the energies of all molecular conformations of each tautomer calculated in Step 7, and take the lowest energy as the final energy of the tautomer;

[0017] Step 9: Sort the tautomers based on the final energy of each tautomer obtained in Step 8, and take the tautomers with energy not exceeding 2.76 kcal / mol as the stable tautomers of the variable region in aqueous solution;

[0018] Step 10: Based on the connection information between fragments, use RDKit to combine the stable tautomers of all fragments to obtain the stable tautomer structure of the complete molecule;

[0019] Step 11: Use MolGpKa to predict the dissociation constants of all dissociation centers in the stable tautomer structure obtained in Step 10, and then use RDKit to generate all protonation states of each stable tautomer.

[0020] Among them, in Step 1, the quantum chemical method used to optimize the molecular conformations of each compound is B3LYP / 6-31G*; the method used to calculate the solvation energy of each compound in aqueous solution is M062X / 6-31G*SMD.

[0021] In Step 2, the deep learning model MolSolv is constructed based on the molecular conformation and solvation energy data calculated in Step 1. The model is a graph convolutional neural network, and the node features describe the basic properties and surrounding chemical environment of the atom represented by the node, including atom type, hybridization state, aromaticity, state of the ring where the atom is located, and atom environment fingerprint.

[0022] In Step 3, the cutting rules are as follows:

[0023] 1) The chemical bond to be cut can only be a single bond outside the ring;

[0024] 2) The atoms at both ends of the cleaved chemical bond cannot both be aromatic atoms;

[0025] 3) The atom at one end of the cleaved chemical bond cannot be an atom of the types NH, NH2, SH, or OH attached to a ring;

[0026] 4) The atom at one end of the cleaved chemical bond cannot be connected to a double bond or a triple bond. The variable region recognition template is composed of the chemical informatics language SMARTS. The RDKit package is used to parse SMARTS and identify the variable region.

[0027] In step 4, the tautomer conversion template is described in the SMIRKS chemical informatics language. The RDKit package is used to parse SMRIKKS and generate tautomers based on the template.

[0028] Compared with the traditional method, the beneficial effects of the present invention are as follows:

[0029] (1) The present invention uses more comprehensive and reliable tautomer conversion rules, which can generate more diverse tautomers for small molecule compounds. The tautomer conversion rules are represented by the molecular reaction language SMIRKS. Based on the molecular reaction module of the chemical informatics software RDKit and these tautomer conversion rules, all possible tautomer structures of each small molecule compound can be generated quickly.

[0030] (2) The present invention uses customized rules to cleave compounds into fragments, generates tautomers for each fragment, and uses a customized scoring method to calculate the energy and rank the tautomers. Here, cleaving the molecule into fragments can effectively reduce the calculation time and improve the accuracy of the machine learning energy function.

[0031] (3) The present invention trains a machine learning solvation energy model based on a solventization energy database of millions of levels. Based on an advanced graph convolutional neural network and introducing molecular conformation information, the solvation energy model has better robustness and accuracy.

[0032] (4) The present invention uses a machine learning energy function to score and rank each tautomer, which has higher calculation accuracy compared with the empirical scoring method; compared with using quantum chemical methods, it can save a large amount of computing power, and can process large compound libraries, improving the efficiency of virtual drug screening.

[0033] (5) The present invention can also analyze the influence of substituent effects on tautomer conversion, helping medicinal chemists design and modify drug molecules, and optimizing the physicochemical properties and binding activities of molecules.

[0034] (6) The present invention integrates four aspects: drug molecule tautomer conversion, scoring, ranking, and protonation, and can efficiently process small molecule compounds.

[0035] (7) The results show that by using the present invention, stable tautomers and protonation states of small molecule compounds can be generated more accurately and quickly. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] Figure 1 is a flowchart of the present invention;

[0037] Figure 2 is a performance evaluation graph of the machine learning scoring method on the test set;

[0038] Figure 3 is a performance evaluation graph of MolTaut for predicting the relative energy between tautomers. DETAILED DESCRIPTION OF THE INVENTION

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

[0040] Refer to Figure 1 , a method for generating stable tautomers and protonation states of drug molecules according to the present invention, which comprises the following specific steps:

[0041] Step 1: Extract 1 million neutral compound structures containing only elements C, H, O, N, S, F, and Cl from the pubchem database. For each compound, first generate a molecular conformation using RDKit, and then optimize the geometric conformation of each compound using the quantum chemistry software Gaussian16 in vacuum at the B3LYP / 6-31G* level. Then, based on the optimized geometric structure of the compound, calculate the solvation energy of each compound using the continuum medium model SMD in Gaussian16 at the M062X / 6-31G* level;

[0042] Step 2: Build a prediction model, MolSolv, that can quickly predict the solvation energy of small molecule compounds based on the molecular geometric structure and solvation energy dataset optimized in Step 1. Here, build a solvation energy prediction model based on a graph convolutional neural network, and the node features describe the basic properties and surrounding chemical environment of the atom represented by the node, including atom type, hybridization state, aromaticity, ring state, and atom environment fingerprint. Atom types include C, H, O, N, S, F, Cl, hybridization states include sp, sp2, sp3, sp3d, sp3d2, sp3d3, aromaticity indicates whether the atom is in an aromatic ring, and the atom environment fingerprint describes the chemical environment around the atom, introducing the geometric conformation information of the molecule. The node features are all calculated using RDKit.

[0043] Step 3: When small molecule compounds need to generate stable tautomers and protonation states, first find all variable regions in the small molecule compounds where tautomer conversions can occur, and extract the substructures of these regions;

[0044] Step 4: For each variable region, generate all possible tautomer structures of the variable region based on the tautomer conversion template. The conversion template is described by the chemical reaction language SMIRKS, and RDKit is used for tautomer conversion;

[0045] Step 5: For each tautomer structure of the same variable region, use RDKit to generate up to 50 molecular conformations;

[0046] Step 6: Use the ANI-2x model to optimize each molecular conformation as the input for the next energy calculation;

[0047] Step 7: Combine the solvation energy prediction model MolSolv based on deep learning and the molecular internal energy prediction model ANI-2x to calculate the energy of each molecular conformation;

[0048] Step 8: Take the lowest energy among all conformations calculated in Step 7 as the final energy of this tautomer;

[0049] Step 9: Sort all tautomers of this variable region based on the energy obtained in Step 8, and take the tautomers with energy less than 2.8 kcal / mol as the stable tautomers of this variable region;

[0050] Step 10: Combine the stable tautomer structures of each variable region to obtain the stable tautomer structure of the final complete molecule;

[0051] Step 11: Use the small molecule compound dissociation constant prediction software to predict the dissociation constants of all potential dissociation centers in the molecule, and generate the potential protonation states of each stable tautomer based on the dissociation constants.

[0052] In Step 2, a solvation energy prediction model is constructed and trained based on Pytorch, and the graph convolutional neural network is implemented by Pytorch Geometric. The model uses the gradient descent method to optimize the model parameters, uses the Adam optimizer, sets the learning rate to 0.001, the training batch size to 512, and trains for 1000 steps on an NVIDIA 2080Ti graphics card. The five-fold cross-validation method is used for training and validating the model. The size ratio of the training set to the validation set is 8:2. The model parameters are saved when the RMSE of the model on the validation set is the smallest, and finally 5 models are obtained. When the model is predicting, the average value of the prediction results of the 5 models is taken as the final result. The performance of the model on the test set is referred to Figure 2 : Figure 2A is the performance of the solvation energy prediction model in predicting the solvation energy of molecules, with RMSE being 0.31 kcal / mol and MAE being 0.15 kcal / mol; Figure 2 B is the performance of the solvation energy prediction model MolSolv in predicting the relative solvation energy between tautomer pairs, with RMSE being 0.90 kcal / mol and MAE being 0.59 kcal / mol.

[0053] The method for extracting the variable region in a molecule described in step 3 is specifically as follows: First, based on specific cleavage rules, use the substructure matching method in RDKit to find all the chemical bonds in the molecule that can be cleaved, then use RDKit to cleave all the cleavable chemical bonds and mark the connection information, such as Figure 1 shown. Perform tautomer conversion on all molecular fragments. The fragments that can undergo effective tautomer conversion are the variable regions, and those that cannot undergo effective tautomer conversion are the immutable regions. Only some special exocyclic single bonds in the molecule can be cleaved, which are defined by the cheminformatics language SMARTS. The specific SMARTS cleavage rules are as follows: [#6+0;!$(*=,#[!#6])]!@!=!#[!#0;!#1;!$([NH,OH,SH]-[*;a]);!$(*=,#[*;!R])]. The following is a detailed description of the rules:

[0054] 1) The chemical bond to be cleaved can only be an exocyclic single bond;

[0055] 2) The atoms at both ends of the chemical bond to be cleaved cannot be aromatic atoms simultaneously;

[0056] 3) The atom at one end of the chemical bond to be cleaved cannot be an atom of the NH, NH2, SH, OH type connected to the ring;

[0057] 4) The atom at one end of the chemical bond to be cleaved cannot be connected to a double bond and a triple bond. The variable region recognition template is composed of the cheminformatics language SMARTS. Use the RDKit package to parse the SMARTS and identify the variable region.

[0058] In step 4, based on the tautomer conversion templates (SMIRKS) collected from the literature, use RDKit to traverse all possible tautomer states of each variable region. The tautomer conversion templates only include the rules for proton transfer and do not consider the changes in the ring. There are a total of 54 reaction rules. When performing tautomer conversion on molecular fragments, first let RDKit read each SMIRKS conversion template, and predict all possible tautomer conversions of the fragment based on the RDKit reaction prediction module. The templates are as follows:

[0059]

[0060]

[0061]

[0062]

[0063]

[0064] In Step 5, since the energy of the tautomer is related to the conformation, each tautomer needs to generate multiple molecular conformations. Here, the number of conformations generated can be selected according to the size of the molecule. To balance the calculation accuracy and calculation time, the final number of conformations cannot be greater than 50. The molecular conformations are generated by the ETKDG method integrated in RDKit.

[0065] In Step 6, based on the ANI-2x deep learning potential energy function and the BFGS algorithm, the ASE software package is used to optimize the molecular conformations generated in Step 5. The maximum number of optimization steps is 25,000, and the optimization stops when the force on a single atom is less than 0.005.

[0066] In Step 7, a deep learning-based solvation energy model is used to predict the SMD solvation energy of the conformation, and ANI-2x is used to predict the molecular internal energy of the conformation. The sum of the molecular internal energy and the solvation energy of the conformation is used as the final energy of the conformation. To speed up the calculation, the thermochemical correction term is ignored. The formula is as follows:

[0067] ΔG t = 0.72 × ΔE ANI-2x + ΔG MolSolv

[0068] Based on the energy calculation formula, the calculation performance of MolTaut for the relative energy between tautomers in aqueous solution was evaluated on the experimental dataset Tautobase. For the calculation results of the DFT (wB97X / 6-31G* / / M062X / 6-31G* / SMD) method, the specific results are referred to Figure 3 , Figure 3 A is the performance of the MolTaut method on the experimental dataset, Figure 3 B is the performance of the DFT method on the experimental dataset. The RMSE of MolTaut = 3.15 kcal / mol, and the MAE = 2.39 kcal / mol, which is similar to the calculation structure of the DFT method.

[0069] In Step 8, the energies of all the conformations calculated in Step 7 are sorted, and the conformation with the lowest energy is taken as the final energy of this tautomer.

[0070] In step 9, the tautomers are sorted based on the energy of each tautomer obtained in step 8, and the tautomers with energy lower than 2.8 kcal / mol are taken as the stable tautomer structures.

[0071] In step 10, based on the connection information between the fragments, the stable tautomers of all fragments are combined using RDKit to obtain the stable tautomer structure of the complete molecule.

[0072] In step 11, MolGpKa is used to predict the dissociation constants of all potential dissociation centers in the molecule, and then RDKit is used to generate all possible protonation states of each stable tautomer.

Claims

1. A method for generating stable tautomers and protonation states of drug molecules, characterized in that, The method includes the following specific steps: Step 1: Extract 1 million compound structures from the PubChem database, optimize the molecular conformation of each compound using quantum chemical methods, and calculate the solvation energy of each compound in aqueous solution; Step 2: Based on the molecular conformation and solvation energy data obtained in Step 1, construct a deep learning model MolSolv. The model MolSolv is a solvation energy prediction model that can quickly predict the solvation energy of small molecule drugs; Step 3: When generating the stable tautomers and protonation states of small molecule compounds in solution, first, based on the cleavage rules, cleave the small molecule compounds into molecular fragments, and find all molecular fragments that can undergo tautomerization conversion as variable regions based on the variable region recognition template; Step 4: For each variable region found in Step 3, based on the collected SMIRKS conversion templates of tautomers, use RDKit to traverse the tautomer states of each variable region; Step 5: For each tautomer structure of the same variable region, use RDKit to generate up to 50 molecular conformations; Step 6: Use the ANI-2x model to optimize each molecular conformation generated by RDKit in Step 5, and use the optimized optimal molecular conformation as the input for the calculation of the next molecular internal energy and hydration solvation energy; Step 7: For the molecular conformation optimized in Step 6, use the deep learning-based solvation energy model MolSolv to predict the solvation energy of the molecular conformation in aqueous solution, and use ANI-2x to predict the internal energy of the molecular conformation. The sum of the internal energy and solvation energy of the molecular conformation is used as the final energy of the optimized molecular conformation. To accelerate the calculation speed, the thermodynamic correction term is ignored. The specific calculation formula is as follows: ΔG = ΔG internal + ΔG solvation , where ΔG internal represents the internal energy of this molecular conformation, and ΔG solvation represents the solvation energy of this molecular conformation in aqueous solution; Step 8: Sort the energies of all molecular conformations of each tautomer calculated in Step 7, and take the lowest energy as the final energy of the tautomer; Step 9: Sort the tautomers based on the final energy of each tautomer obtained in Step 8, and take the tautomers with an energy not exceeding 2.76 kcal / mol as the stable tautomers of the variable region in aqueous solution; Step 10: Based on the connection information between the fragments, use RDKit to combine the stable tautomers of all fragments to obtain the stable tautomer structure of the complete molecule; Step 11: Use MolGpKa to predict the dissociation constants of all dissociation centers in the stable tautomer structure obtained in Step 10, and then use RDKit to generate all protonation states of each stable tautomer.

2. The method according to claim 1, characterized in that In Step 1, the quantum chemical method used to optimize the molecular conformation of each compound is B3LYP / 6-31G*; the method used to calculate the solvation energy of each compound in aqueous solution is M062X / 6-31G*SMD.

3. The method according to claim 1, characterized in that, In Step 2, the deep learning model MolSolv constructed based on the molecular conformation and solvation energy data calculated in Step 1 is a graph convolutional neural network. The node features describe the basic properties and surrounding chemical environment of the atoms represented by the nodes, including atom type, hybridization state, aromaticity, the state of the ring where the atom is located, and the atom environment fingerprint.

4. The method according to claim 1, wherein In Step 3, the cleavage rules are as follows: 1) The chemical bond to be cleaved can only be a single bond outside the ring; 2) The atoms at both ends of the chemical bond to be cleaved cannot be aromatic atoms at the same time; 3) The atom at one end of the chemical bond to be cleaved cannot be an atom of the NH, NH2, SH, or OH type connected to the ring; 4) The atom at one end of the cleaved chemical bond cannot be connected to double bonds and triple bonds. The variable region recognition template is composed of the chemical informatics language SMARTS. The RDKit package is used to parse SMARTS and identify the variable region.

5. The method according to claim 1, wherein In step 4, the tautomer conversion template is described in the SMIRKS chemical informatics language. The RDKit package is used to parse SMRIKKS and generate tautomers based on the template.

Images