An optimized de novo drug design approach

Abstract

The present invention relates to an optimized de novo drug design method, and belongs to the fields of computer-aided drug design technology and biological and pharmaceutical informatics. The computer-aided molecular design method based on the mixed integer nonlinear programming model of the present invention is used to design candidate drugs, and is coupled with a deep learning model for predicting target-ligand binding affinity. As one of the most important properties in the field of drug design, the binding affinity property is regarded as an objective function in the subsequent mixed integer nonlinear programming model. By using the Bemis-Murcko algorithm to convert complex drug structures into drug skeletons, the skeleton structure is considered in the computer-aided molecular design method to ensure the rationality of the designed candidate drugs. In addition, a skeleton-based similarity algorithm is proposed to significantly reduce the scale of drug design problems, so that the mixed integer nonlinear programming model can effectively design candidate drugs.

Description

Technical Field

[0001] The present invention relates to computer-aided drug design technology and the field of biological and pharmaceutical informatics, and in particular to a de novo drug design method. Background Art

[0002] As a high-value-added chemical product, small molecule drugs are often used to prevent diseases and protect public health. However, many challenges hinder the development of innovative drugs, including high listing costs, low success rates in clinical trials, and long development cycles. The average cost of developing a new drug is estimated to be more than $2.5 billion, and the development cycle is approximately 10 to 15 years. Therefore, in order to reduce the time and cost associated with experimental analysis and improve the hit rate of subsequent clinical studies, it is worth considering using computer-aided technologies and models to accelerate drug discovery in the early stages of drug development.

[0003] Virtual screening methods usually use molecular docking methods or molecular similarity algorithms to identify new potential drugs from large chemical libraries containing thousands to hundreds of millions of candidate compounds. Despite the high efficiency of virtual screening methods, the cost of traversing and screening promising drug candidates through enumeration methods is still high, considering that the number of molecules that can be considered as potential drug-like molecules is estimated to be 10 30 to 10 60 In the meantime, it is computationally infeasible to conduct a complete search in this chemical space. To avoid this problem, optimization-based de novo drug design methods have attracted widespread attention from scholars.

[0004] In contrast, de novo drug design methods can design molecular structures from scratch according to given requirements, and have the advantage of exploring a wider chemical space beyond the existing compound library. Generally speaking, optimization-based de novo drug design methods consist of three parts, including a generative model for creating a chemical space for candidate drugs, a property prediction model for evaluating the properties of candidate drugs, and an optimization algorithm for finding the best candidate drugs with ideal properties. For example, the drug design method based on genetic algorithms creates a chemical space by randomly combining groups through heuristic rules, and then optimizes the population (candidate drugs) according to the fitness function (objective function) through a gradient-free optimization algorithm. This drug design method can quickly identify candidate drugs that meet the property requirements in a vast chemical space outside the existing compound library. However, the drug design method based on genetic algorithms is prone to fall into local optimal solutions, which limits the exploration of more promising candidate drugs.

[0005] Recently, deep learning-based generative models have been developed for de novo drug design. Compared with genetic algorithm-based methods, deep generative models do not require heuristic rules. They are trained and predicted in a completely data-driven and high-throughput manner, with minimal requirement for expert knowledge. Widely used deep learning architectures include variational autoencoders, generative adversarial networks, etc. Nevertheless, deep generative models are still in their infancy because structurally unfeasible molecules are often observed in the generated molecules. In addition, the extrapolation ability of deep generative models is still unsatisfactory. For example, it has been shown that generative adversarial networks are unstable and prone to produce small chemical spaces (i.e., the model only generates a single molecular structure or a small number of similar molecular structures). In addition, deep generative models are often presented as "black boxes" with low interpretability. If the computational cost of "white box" models is as low as that of deep generative models, it is worth looking forward to developing "white box" models for drug design by using mathematical formulas with clear interpretability.

[0006] Computer-aided molecular design (CAMD) technology has long been used to design optimal molecular structures that meet target properties. It is often used to design small molecule solvents such as absorbents and extractants. The CAMD problem is usually formulated as a mixed integer nonlinear programming (MINLP) model consisting of an objective function, molecular structure constraints, and molecular property constraints. The MINLP model can optimize the combination of groups by using a highly interpretable mathematical formula to create a chemical space containing reasonable molecules. By solving the MINLP model with a gradient-based optimization algorithm, all feasible solutions (i.e., molecules represented by a set of groups) that satisfy the structural and property constraints can be obtained, and the optimal solution with the maximum or minimum target property can be determined. Therefore, mathematical programming models can generate molecules with the desired properties more efficiently than deep generative models, because deep generative models usually use fingerprint-based similarity indicators to generate molecules, and then traverse and evaluate the properties of the generated molecules by enumeration. However, when the nonlinear equations in the MINLP model are too complex, it is actually not feasible to solve the MINLP model directly. In order to solve this problem, some scholars have proposed a decomposition solution algorithm to solve the MINLP model with extremely strong non-convexity. Although the CAMD method has achieved great success in designing small molecule solvents, there are still two challenges in directly applying the traditional CAMD method to candidate drug design: First, the structure of the drug (especially the ring structure) is more complex and larger than that of the small molecule solvent. Therefore, a larger number of ring groups are required to establish the MINLP model for candidate drug design, which will increase the problem scale and difficulty of solving the MINLP model. Even if the MINLP model is successfully solved, the second challenge is that the traditional CAMD method is prone to produce some structurally feasible but abnormal ring structures. For example, if the cyclic group "aC-C#CH" is selected for the MINLP model, the traditional CAMD method will design a molecule similar to "CC(C)c1c(C#C)c(C#C)c(C#C)c(C#C)c1C#C" (SMILES representation). Summary of the invention

[0007] In view of the above problems existing in the prior art and in order to address the above two challenges, the purpose of the present invention is to establish an optimization-based de novo drug design method, in which a CAMD method based on a MINLP model is proposed for candidate drug design and coupled with a deep learning model for predicting target-ligand binding affinity. As one of the most important properties in the field of drug design, the binding affinity property is considered as the objective function in the subsequent MINLP model. The Bemis-Murcko algorithm is used to convert complex drug structures into drug skeletons, which are considered in our customized CAMD method to ensure the rationality of the designed candidate drugs. In addition, a skeleton-based similarity algorithm is proposed to greatly reduce the scale of the drug design problem and enable the MINLP model to effectively design candidate drugs. Finally, a case study involving the design of anti-tumor drug candidates is introduced to highlight the versatility and effectiveness of the established optimization-based de novo drug design framework.

[0008] To achieve the above object, the present invention adopts the following technical solution, including the following specific steps:

[0009] An optimized de novo drug design method comprising the following specific steps:

[0010] Step 1: Establish a drug database;

[0011] Step 2: Building a database containing Bemis-Murcko skeletons: The Bemis-Murcko skeletons were extracted from the drug structures in the drug database by using the Bemis–Murcko algorithm in RDKit.

[0012] Step 3: For the skeleton structure of the target drug, a skeleton-based similarity algorithm is used to search for a skeleton subset G1 similar to the skeleton of the target drug from the skeleton database, and a group of commonly used group groups is selected to form a set G2.

[0013] Step 4: Design promising drug candidates: The de novo drug design problem is formulated as a mixed integer nonlinear programming (MINLP) model consisting of an objective function, drug structure constraints, and drug property constraints (including deep learning models); by solving the MINLP model, the optimal candidate drug structure with a high probability of binding affinity is identified from the feasible solutions, which are generated by optimizing the drug skeleton and groups, and the combination process is restricted by the MINLP model constraints; the objective function is to maximize the probability of high binding affinity of the target-ligand complex. The MINLP model formula is as follows:

[0014] Objective function: Subject to steps 5-8:

[0015] Step 5: Deep learning constraints: The general equation (1) represents a deep learning model for identifying target-ligand complexes with high binding affinity.

[0016] Prob bjnd =f deep (s) (1)

[0017] Step 6: Drug structure constraints: The general equation (2) represents the structural constraints of the octet rule m1, the valence bond rule m2 and the chemical complexity m3. The combination of the skeleton and the groups can generate molecules with reasonable structures.

[0018]

[0019] Step 7: Constraints on drug properties: General equations (3) and (4) represent the “Lipinski’s Five Rules” properties: relative molecular mass MW, number of hydrogen bond acceptors HBA, number of hydrogen bond donors HBD, octanol-water partition coefficient logP, rotatable angle ROT (ROT frag ), synthesis feasibility score SA and synthesis complexity score SC.

[0020]

[0021]

[0022] Step 8: Other constraints: General equation (5) represents an improved SMILES-based isomer generation algorithm, which is used to automatically convert a set of fragments of a candidate drug into a corresponding drug SMILES string, where the fragments of the candidate drug include a skeleton and a group.

[0023] f SMILES (n i ,s)=0 (5)

[0024] In the general equation above, F obj is the objective function, Prob bind is the probability of high binding affinity of the target-ligand complex, i represents the fragment involved in the candidate drug, and n i represents the number of fragments involved in the drug candidate, s is the SMILES representation of the molecule, m is the type of structural constraint, and are the upper and lower bounds of the structural constraint m, p is the drug property, k is the property type, and Yes k The upper and lower bounds of .

[0025] Step 9: Use a decomposition-based solution algorithm to solve the MINLP model. If there is no optimal candidate drug that meets all constraints, return to step 4 to relax the constraints.

[0026] Step 10: The optimal solution of the MINLP model is further verified by physics-based molecular docking and molecular dynamics simulation methods.

[0027] Furthermore, step 1 specifically includes:

[0028] Step 1.1: Collect small molecule drugs with CAS numbers from the DrugBank database.

[0029] Step 1.2: Search the CID number of the drug in PubChem by using the crawler script and the drug CAS number; after removing the drugs without CID number, the remaining drugs are available for further screening.

[0030] Step 1.3: Use the “Lipinski’s Rule of Five” properties: relative molecular mass MW ≤ 500, the number of hydrogen bond donors HBD ≤ 5, the number of hydrogen bond acceptors HBA ≤ 10, the octanol-water partition coefficient logP ≤ 5, and the rotatable angle ROT ≤ 10 to screen drugs with good pharmacokinetic properties; the “Lipinski’s Rule of Five” properties and isomer SMILES strings of all drugs were obtained through the official website interface of the PubChem database according to their CID numbers.

[0031] After applying the above-mentioned criteria from step 1.1 to step 1.3, a drug database containing small molecule drugs and their CAS numbers, CID numbers, isomer SMILES strings, and “Lipinski’s Rule of Five” properties is established.

[0032] Furthermore, in step 3, the skeleton-based similarity algorithm is a combination of six similarity algorithms and four molecular representation methods, forming a total of 24 combinations, using each combination to identify a skeleton similar to the target drug skeleton, taking the three most similar skeletons obtained from each combination, removing duplicates, and obtaining the final similar skeleton subset G1. The six similarity algorithms include Tanimoto, Dice, Cosine, Sokal, Russel, and Kulczynski, and the four molecular representation methods include topological fingerprints, MACCS keys, ECFP fingerprints, and FCFP fingerprints.

[0033] Furthermore, step 6 specifically includes:

[0034] The specific drug structure constraints of the MINLP model are given by equations (6) to (12); This means that a drug candidate selects only one skeleton from the skeleton set.

[0035] Among them, the structural constraints of chemical complexity are:

[0036]

[0037]

[0038] n i ≤3,i∈G2 (8)

[0039]

[0040]

[0041] The structural constraints of the octet rule are:

[0042]

[0043] The structural constraints of the valence bond rule are:

[0044]

[0045] In equations (6) to (12), n i is the number of fragments i involved in the candidate drug, v i is the number of valence bonds of fragment i, G is the fragment set, G1 is the skeleton subset, G2 is the group set, G3 is the subset of G2, G3={-CH3,-CH2, -CH,CH2=CH-,-CH=CH-,CH2=C<,-CH=C<}. is the number of i1 fragments, is the number of i2 segments, is the bond number of the i1 fragment.

[0046] Furthermore, step 7 specifically includes:

[0047] Step 7.1: In the MINLP model, the “Lipinski Rule of Five” properties (MW, HBA, HBD, logP, ROT) are calculated by the RDKit-based quantitative estimation method (QED); the upper limits of the “Lipinski Rule of Five” properties are MW ≤ 500, HBD ≤ 5, HBA ≤ 10, logP ≤ 5 and ROT (ROT frag )≤10; ROT in equation (3) frag is the linear sum of the rotatable angles of the fragments involved in the candidate drug, while ROT in equation (4) is the overall rotatable angle of the molecule calculated based on the SMILES string of the candidate drug, which is a nonlinear property; considering ROT frag ≤ROT,ROT frag >10 drugs will not satisfy the constraint of ROT≤10; therefore, when using the decomposition algorithm to solve the MINLP model, ROT is introduced fragSome infeasible solutions can be eliminated before calculating ROT, which helps to improve the solution efficiency of the MINLP model.

[0048] Step 7.2: SA and SC properties are used to ensure that the designed candidate drug molecule is easy to synthesize, where low SA and SC values ​​indicate that the molecule is easy to synthesize.

[0049] Furthermore, step 8 specifically includes:

[0050] The SMILES string of the candidate drug is the predicted Prob bind , ROT, logP, SA and SC properties of the model input information; a "bridge" needs to be established in the structural constraints to associate the SMILES string with the fragment set, otherwise the MINLP model cannot be successfully solved; therefore, by using the improved SMILES-based isomer generation algorithm, the fragment set of the candidate drug can be automatically converted into its corresponding SMILES string, so that the MINLP model can be successfully solved, and the algorithm can identify candidate drug isomers containing the same fragment set.

[0051] The improved SMILES-based isomer generation algorithm is implemented through the RWMol module in the RDKit library: a molecule is represented by a set of fragments, which include skeletons and groups, and this information will be fed into the isomer generation algorithm. The specific algorithm execution process is as follows:

[0052] (1) Select the skeleton as the "seed", and then combine the groups as "leaves" with the skeleton according to the pre-set group addition order. Check if there are "*" matches in the seed-leaf structure. If there are, add another group to the seed-leaf structure according to the pre-set group addition order. If not, proceed to the next step.

[0053] (2) Check whether all the groups have been added to the seed-leaf structure. If not, it means that there are extra groups that have not been added to the seed-leaf structure. In this case, delete the seed-leaf structure and repeat step (1) according to the next pre-set group addition sequence. If all the groups have been added, proceed to the next step.

[0054] (3) Check whether the generated seed-leaf structure can generate a SMILES structure through the RDKit library; if not, repeat steps (1)-(2) according to the next pre-set group addition sequence; if yes, save the SMILES result and repeat steps (1)-(2) according to the next pre-set group addition sequence; after trying all pre-set group addition sequences, all possible SMILES structures are obtained.

[0055] Furthermore, step 9 specifically includes:

[0056] Step 9.1 uses a decomposition solution algorithm to solve the MINLP optimization model: decompose the MINLP model into a mixed integer linear programming (MILP) subproblem and three nonlinear programming (NLP) subproblems.

[0057] Step 9.2: Subproblem 1 (MILP): First constrain the structure constraints of the octet rule, valence bond rule, and chemical complexity, as well as MW, HBA, HBD, and ROT frag The linear property constraints are used in GAMS to generate N1 feasible solutions using the BARON solver. The feasible solutions are candidate drugs represented by the fragment set. A parameter N1 needs to be set in GAMS. max To specify the maximum number of solutions for the MILP model.

[0058] Step 9.3: Subproblem 2 (NLP): Using the constraints of the improved SMILES-based isomer generation algorithm, generate SMILES strings of N2 candidate drugs based on the N1 fragment set, N2 ≥ N1.

[0059] Step 9.4: Subproblem 3 (NLP): Considering the nonlinear property constraints of ROT, logP, SA, and SC, use the nonlinear property prediction model to calculate the corresponding properties of the SMILES strings of the N2 candidate drugs, and eliminate those that do not meet the constraints; select the remaining N3 candidate drug SMILES strings for further evaluation.

[0060] Step 9.5: Subproblem 4 (NLP): Considering the constraints of the deep learning model, use the nonlinear deep learning model to calculate the objective function Prob of N3 target-ligand complexes bind , and sort the SMILES strings of candidate drugs according to the objective function.

[0061] Compared with the prior art, the present invention has the following beneficial effects:

[0062] (1) First, the present invention applies the CAMD method to the design of candidate drugs. In fact, it is not easy to directly apply the traditional CAMD method to candidate drugs. The first challenge is that the structure of the drug (especially the ring structure) is more complex and larger than the structure of the small molecule solvent, which will make the MINLP model difficult to solve. For example, the number of ring structures in the solvent is generally 0 to 2, while the number of ring structures in the drug is usually 2 to 5. This means that we need a larger number of ring groups to establish the MILP model (sub-problem 1 after using the decomposition algorithm). If the CAMD method is applied to candidate drugs, this will increase the size of the MILP model and the difficulty of solving it.

[0063] (2) Secondly, even if the MILP model is successfully solved, the second challenge is that the traditional CAMD method is prone to produce some structurally feasible but anomalous ring structures. For example, if the cyclic group "aC-C#CH" is selected for the MILP model, the traditional CAMD method will design a molecule similar to "CC(C)c1c(C#C)c(C#C)c(C#C)c(C#C)c1C#C" (SMILES representation).

[0064] To address the above two challenges, we used the Bemis–Murcko algorithm to obtain drug skeletons and introduced the skeletons into the traditional CAMD method to ensure the rationality of the designed drug candidates. These skeletons were generated from approved drugs or drug candidates in the DrugBank database. However, if we use all 2,898 generated drug skeletons to design drug candidates, the computational efficiency of the MILP model is still unacceptable. Considering that drug candidates with similar skeletons may have similar properties, before establishing the MILP model, we use a skeleton-based similarity algorithm to identify a subset of skeletons similar to the target marketed drug; in this way, we only need to focus on potential skeletons to design drug candidates, which can greatly reduce the size of the MILP model.

[0065] (3) Again, we used two state-of-the-art deep learning techniques, gate-enhanced attention mechanism and convolutional neural network, to establish a reasonable deep learning model for predicting the binding affinity of target-ligand complexes and integrated it into the CAMD method.

[0066] In summary, the present invention has established an optimization-based framework that combines a deep learning model with a customized MINLP-based CAMD method for candidate drug design. The present invention can significantly reduce the time and cost associated with experimental analysis, improve the hit rate of subsequent clinical studies, and design new chemical structures, thereby improving the efficiency and diversity of new drug design. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 A flowchart of the optimized de novo drug design framework of the present invention;

[0068] Figure 2 The schematic diagram of the drug structure of the present invention is composed of a skeleton (ring structure and connecting structure) and a side chain.

[0069] Figure 3 The chemical space of candidate drug design results drawn using the ECFP fingerprint and principal component analysis (PCA) method of the present invention. (The integers (0 to 14) in the legend on the right represent 15 skeletons, and the point "7" represents axitinib) DETAILED DESCRIPTION

[0070] The following is an explanation of the process and effect of the present invention in conjunction with the accompanying drawings and embodiments. The embodiments of the present invention are implemented on the premise of the technical solution of the present invention, and detailed implementation methods and specific operation processes are given, but the protection scope of the present invention is not limited to the following embodiments.

[0071] See also Figure 1 The embodiment of the present invention specifically discloses a de novo drug design framework based on optimization, comprising the following steps:

[0072] Step 1: A drug database was established with reference to the DrugBank V5.0 database (https: / / go.drugbank.com / ) for subsequent drug skeleton generation.

[0073] Step 2: Generate drug skeletons for subsequent MINLP-based drug design models. Figure 2 ) to replace the traditional cyclic groups, more structurally reasonable candidate drugs can be designed. Therefore, a database of 2,898 Bemis-Murcko skeletons was established for subsequent candidate drug design. These skeletons were generated by using the Bemis–Murcko algorithm in RDKit for 4,781 drugs in the drug database.

[0074] Step 3: In order to reduce unnecessary search space and save computational time cost for candidate drug design problems, a skeleton-based similarity algorithm was used to identify a subset of skeletons G1 similar to the skeleton of axitinib (an FDA-approved anti-tumor drug) before establishing the MINLP model. A total of 72 skeletons were obtained, and 19 skeletons were obtained after deleting duplicate values. At the same time, commonly used groups were selected to form the set G2.

[0075] Step 4: Design promising drug candidates in a high-throughput and intelligent manner based on optimized mathematical programming methods. The de novo drug design problem is formulated as a MINLP model consisting of an objective function (maximizing the probability of high binding affinity of the target-ligand complex), drug structure constraints, and drug property constraints (including deep learning models). By solving the MINLP model, the optimal drug candidates with high binding affinity among feasible solutions can be identified. These feasible solutions are generated by optimizing the drug skeleton and group groups under the constraints of the model. Note that each MINLP model is applicable to a target for which promising drug candidates can be designed. The MINLP model formula is as follows:

[0076] Objective function: Subject to steps 5-8:

[0077] Step 5: Deep learning constraints: The general equation (Eq. (1)) represents a deep learning model for identifying target-ligand complexes with high binding affinity.

[0078] Prob bind =f deep (s) (1)

[0079] Step 6: Drug structure constraints: The general equation (Eq. (2)) represents the structural constraints of the octet rule (m1), the valence bond rule (m2) and the chemical complexity (m3), and a structurally reasonable molecule can be generated through the combination of the skeleton and the groups.

[0080]

[0081] Step 7: Constraints on drug properties: The general equation (Equations (3-4)) represents the “Lipinski Rule of Five” properties (relative molecular weight MW, number of hydrogen bond acceptors HBA, number of hydrogen bond donors HBD, octanol-water partition coefficient logP, rotatable angle ROT (ROT frag )), synthesis feasibility score (SA) and synthesis complexity score (SC).

[0082]

[0083]

[0084] Step 8: Other constraints: The general equation (Eq. (5)) represents an improved SMILES-based isomer generation algorithm, which is used to automatically convert a set of fragments (skeleton and groups) of a candidate drug into a corresponding drug SMILES string.

[0085] f SMILES (n i ,s)=0 (5)

[0086] In the general equation above, F obj is the objective function, Prob bind is the probability of high binding affinity of the target-ligand complex, i represents the fragments (skeleton and groups) involved in the candidate drug, and n i represents the number of fragments involved in the drug candidate, s is the SMILES representation of the molecule, m is the type of structural constraint, and are the upper and lower bounds of the structural constraint m, p is the drug property, k is the property type, and Yes k The upper and lower bounds of .

[0087] Step 9: Since there are a large number of nonlinear constraints involved in the MINLP model, a decomposition-based solution algorithm is used to solve the MINLP model. If there is no optimal candidate drug that meets all the constraints, return to step 4 to relax the constraint range.

[0088] Step 10: The optimal solution of the MINLP model is further validated by physics-based molecular docking and molecular dynamics simulation methods. Finally, a case study involving the design of anti-tumor drug candidates is proposed to highlight the effectiveness of the de novo drug design framework based on the MINLP model.

[0089] Furthermore, step 1 specifically includes:

[0090] Step 1.1: First, 7,746 small molecule drugs with CAS numbers were collected from the DrugBank database V5.0 (https: / / go.drugbank.com / ).

[0091] Step 1.2: Then, the CID numbers (unique identifiers of chemicals in the PubChem database, https: / / pubchem.ncbi.nlm.nih.gov / search / ) of 7,746 drugs were searched in the PubChem database by using the crawler script and the CAS number. After removing drugs without CID numbers, 7,474 drugs were available for further screening.

[0092] Step 1.3: Screen drugs with good pharmacokinetic properties using the “Lipinski rule of five” properties (relative molecular weight (MW) ≤ 500, number of hydrogen bond donors (HBD) ≤ 5, number of hydrogen bond acceptors (HBA) ≤ 10, octanol-water partition coefficient (logP) ≤ 5, rotatable angle (ROT) ≤ 10). The “Lipinski rule of five” properties and isomer SMILES strings of all 7,474 drugs were obtained through the official website interface of the PubChem database according to their CID numbers.

[0093] Finally, after applying the above criteria, a drug database containing 4,781 small molecule drugs and their CAS numbers, CID numbers, isomer SMILES strings, and “Lipinski’s Rule of Five” properties was established.

[0094] Furthermore, step 3 specifically includes:

[0095] Step 3.1: If all 2,898 skeletons are used in a single MINLP-based drug design model, the excessive model size will hinder the effective solution of the MINLP model. Therefore, we propose a skeleton-based similarity algorithm to identify a subset of skeletons (G1) similar to the skeleton of the target drug from the 2,898 skeletons (determined by case studies). In this way, the skeleton-based similarity algorithm combines six similarity algorithms and four molecular representation methods to form a total of 24 combinations, using each combination to identify skeletons similar to the target drug skeleton, taking the three most similar skeletons obtained from each combination, removing duplicates, and obtaining the final similar skeleton subset G1; the six similarity algorithms include Tanimoto, Dice, Cosine, Sokal, Russel, Kulczynski, and the four molecular representation methods include topological fingerprint, MACCS key, ECFP fingerprint, and FCFP fingerprint.

[0096] Step 3.2: After removing duplicate skeletons, a certain number of similar skeletons can be obtained according to the specific problem. In addition, 29 commonly used group groups (G2) were selected for our MINLP model, as shown in Table 1.

[0097] Table 1 Selected groups

[0098]

[0099]

[0100] Furthermore, step 6 specifically includes:

[0101] The specific drug structure constraints of the MINLP model are given by equations (Equations (6-12)). Indicates that a candidate drug selects only one skeleton from the skeleton group.

[0102]

[0103]

[0104] n i ≤3,i∈G2 (8)

[0105]

[0106]

[0107]

[0108]

[0109] In equations (6 to 12), ni is the number of fragments (skeleton and groups) i involved in the candidate drug, v i is the number of valence bonds of fragment i, G is the fragment set, G1 is the skeleton set, G2 is the group set, and G3 is a subset of G2 (G3 = {-CH3, -CH2, -CH, CH2 = CH-, -CH = CH-, CH2 = C<, -CH = C<}); is the number of i1 fragments, is the number of i2 segments, is the bond number of the i1 fragment.

[0110] Furthermore, step 7 specifically includes:

[0111] Step 7.1: In the MINLP model, the “Lipinski Rule of Five” properties (MW, HBA, HBD, logP, ROT) are calculated by the RDKit-based quantitative estimation method (QED); the upper limits of the “Lipinski Rule of Five” properties are MW ≤ 500, HBD ≤ 5, HBA ≤ 10, logP ≤ 5 and ROT (ROT frag )≤10; ROT in equation (3) frag is the linear sum of the rotatable angles of the fragments involved in the candidate drug, while ROT in equation (4) is the overall rotatable angle of the molecule calculated based on the SMILES string of the candidate drug, which is a nonlinear property; considering ROT frag ≤ROT,ROT frag >10 drugs will not satisfy the constraint of ROT≤10; therefore, when using the decomposition algorithm to solve the MINLP model, ROT is introduced frag Some infeasible solutions can be eliminated before calculating ROT, which helps to improve the solution efficiency of the MINLP model.

[0112] Step 7.2: SA (1-10) and SC (1-5) properties are used to ensure that the designed drug candidate is easy to synthesize, where low SA and SC values ​​indicate that the molecule is easy to synthesize. Based on our empirical knowledge, their constraints are set to SA ≤ 6 and SC ≤ 3.4.

[0113] Furthermore, step 8 specifically includes:

[0114] The SMILES string of the candidate drug is the predicted Prob bind, ROT, logP, SA and SC properties of the model input information; a "bridge" needs to be established in the structural constraints to associate the SMILES string with the fragment set, otherwise the MINLP model cannot be successfully solved; therefore, by using the improved SMILES-based isomer generation algorithm, the fragment set of the candidate drug can be automatically converted into its corresponding SMILES string, so that the MINLP model can be successfully solved, and the algorithm can identify candidate drug isomers containing the same fragment set.

[0115] The improved SMILES-based isomer generation algorithm is implemented through the RWMol module in the RDKit library: a molecule is represented by a set of fragments, which include skeletons and groups. This information will be fed into the isomer generation algorithm; the specific algorithm execution process is as follows:

[0116] (1) Select the skeleton as the "seed", and then combine the groups as "leaves" with the skeleton according to the pre-set group addition order; check whether there is a "*" match in the seed-leaf structure; if there is, add another group to the seed-leaf structure according to the pre-set group addition order; if not, proceed to the next step;

[0117] (2) Check whether all the groups have been added to the seed-leaf structure; if not, it means that there are extra groups that have not been added to the seed-leaf structure. In this case, delete the seed-leaf structure and repeat step (1) according to the next pre-set group addition sequence; if all the groups have been added, proceed to the next step;

[0118] (3) Check whether the generated seed-leaf structure can generate a SMILES structure through the RDKit library; if not, repeat steps (1)-(2) according to the next pre-set group addition sequence; if yes, save the SMILES result and repeat steps (1)-(2) according to the next pre-set group addition sequence; after trying all pre-set group addition sequences, all possible SMILES structures are obtained.

[0119] Furthermore, step 9 specifically includes:

[0120] Step 9.1: Generally speaking, if the nonlinear equations in the MINLP model are not particularly complex, the BARON solver in the GAMS software (https: / / www.gams.com / ) can be used to solve the MINLP model directly.

[0121] Step 9.2: However, the SMILES-based property prediction model (e.g., deep learning model) in our MINLP model contains extremely complex nonlinear equations, which makes it extremely difficult to directly solve the MINLP model. To this end, a decomposition solution algorithm is used to solve our complex MINLP optimization model. It decomposes the MINLP model into a mixed integer linear programming (MILP) subproblem and three nonlinear programming (NLP) subproblems.

[0122] Step 9.3: Subproblem 1 (MILP): First restrict the structural constraints of the octet rule, valence bond rule, and chemical complexity, as well as MW, HBA, HBD, and ROT frag The linear property constraints are used in GAMS to generate a certain number (N1) of feasible solutions (candidate drugs represented by fragment sets). Note that a parameter N1 needs to be set in GAMS max To specify the maximum number of solutions for the MILP model, we used a mathematical programming method and obtained N1 = 33,759 feasible solutions (represented by a set of fragments) under the constraints of structural and linear properties, which took 210 seconds on our desktop computer (Intel(R) Core(TM) i7-10700F CPU@2.90GHz 24.0GB RAM).

[0123] Step 9.4: Sub-problem 2 (NLP): Using the constraints of the improved SMILES-based isomer generation algorithm, generate SMILES strings of N2=262,741 candidate drugs based on the N1 set of fragments (3,113 seconds).

[0124] Step 9.5: Subproblem 3 (NLP): Considering the nonlinear property constraints of ROT, logP, SA and SC, the nonlinear property prediction model is used to calculate the corresponding properties of the SMILES strings of N2 candidate drugs, and those that do not meet the constraints are eliminated. The remaining N3 = 105,164 SMILES strings of candidate drugs are used for further analysis.

[0125] First, we searched the SMILES strings of the 105,164 compounds we designed in the PubChem database. We found that 4,902 (4.64%) of the designed molecular structures already existed in PubChem, which shows that our MINLP-based drug design model can not only find existing drug candidates, but also design new drug candidates (95.36%).

[0126] Second, the chemical space was created using ECFP fingerprinting and principal component analysis (PCA) (see Figure 3) is used to characterize the structural diversity of 105,164 designed drug candidates, where the x-axis and y-axis represent the two principal components PC1 and PC2, respectively. The integers (0-14) in the legend on the right represent 15 skeletons, and the point "7" represents axitinib. Figure 3 The results show that the designed drug candidates are widely distributed in the chemical space, indicating the powerful ability of our MINLP-based drug design model in designing a variety of structurally diverse drug candidates similar to axitinib.

[0127] Step 9.6: Subproblem 4 (NLP): Consider the deep learning model constraints and use the nonlinear deep learning model to calculate the Prob of N3 target-ligand complexes bind (objective function), and the SMILES strings of the candidate drugs were ranked according to the objective function. The ranking results showed that 433 designed candidate drugs were superior to axitinib (97.96%) in terms of the probability of high binding affinity. Some of the top-ranked candidate drugs were further verified by other physics-based methods such as molecular docking and molecular dynamics simulation.

[0128] The foregoing description of specific exemplary embodiments of the present invention is for the purpose of illustration and demonstration. These descriptions are not intended to limit the present invention to the precise form disclosed, and it is clear that many changes and variations can be made based on the above teachings. The purpose of selecting and describing the exemplary embodiments is to explain the specific principles of the present invention and its practical application, so that those skilled in the art can realize and utilize various different exemplary embodiments of the present invention and various different selections and changes. The scope of the present invention is intended to be limited by the claims and their equivalents.

Claims

1. An optimized de novo drug design method, characterized in that: The specific steps include: Step 1: Establish a drug database; Step 2: Establish a database containing Bemis-Murcko skeletons: Extract the Bemis-Murcko skeletons from the drug structures in the drug database by using the Bemis–Murcko algorithm in RDKit; Step 3: For the skeleton structure of the target drug, a skeleton-based similarity algorithm is used to search for a skeleton subset G1 similar to the skeleton of the target drug from the skeleton database, and a set of commonly used group groups is selected to form a set G2; Step 4: Design promising drug candidates: The de novo drug design problem is formulated as a mixed integer nonlinear programming MINLP model consisting of an objective function, drug structure constraints, and drug property constraints including deep learning models; by solving the MINLP model, the optimal candidate drug structure with a high probability of binding affinity is identified from feasible solutions, which are generated by optimizing the drug skeleton and groups, and the combination process is restricted by the MINLP model constraints; the objective function is to maximize the probability of high binding affinity of the target-ligand complex; the MINLP model formula is as follows: Objective function: Subject to steps 5-8: Step 5: Deep learning constraints: The general equation (1) represents the deep learning model for identifying target-ligand complexes with high binding affinity; Probpind = fdeep(s)(1) Step 6: Drug structure constraints: General equation (2) represents the structural constraints of the octet rule m1, the valence bond rule m2 and the chemical complexity m3. The combination of the skeleton and the groups can generate molecules with reasonable structures. Step 7: Constraints on drug properties: General equations (3) and (4) represent the "Lipinski's Five Rules" properties: relative molecular mass MW, number of hydrogen bond acceptors HBA, number of hydrogen bond donors HBD, octanol-water partition coefficient logP, rotatable angle ROT or ROTfrag, synthetic feasibility score SA and synthetic complexity score SC; Step 8: Other constraints: General equation (5) represents an improved SMILES-based isomer generation algorithm, which is used to automatically convert the fragment set of the candidate drug into the corresponding drug SMILES string, where the fragment of the candidate drug includes the skeleton and the group; f SMILES (n i ,s)=0 (5)In the above general equation, F obj is the objective function, Prob bind is the probability of high binding affinity of the target-ligand complex, i represents the fragment involved in the candidate drug, and n i represents the number of fragments involved in the drug candidate, s is the SMILES representation of the molecule, m is the type of structural constraint, and are the upper and lower bounds of the structural constraint m, p is the drug property, k is the property type, and Yes k The upper and lower bounds of Step 9: Use a decomposition-based solution algorithm to solve the MINLP model; if there is no optimal candidate drug that meets all constraints, return to step 4 to relax the constraint range; Step 10: The optimal solution of the MINLP model is further verified by physics-based molecular docking and molecular dynamics simulation methods.

2. An optimized de novo drug design method according to claim 1, characterized in that: Step 1 specifically includes: Step 1.1: Collect small molecule drugs with CAS numbers from the DrugBank database; Step 1.2: Search the CID number of the drug in PubChem by using the crawler script and the drug CAS number; after removing the drugs without CID number, the remaining drugs are available for further screening; Step 1.3: Use the "Lipinski's Five Rules" properties: relative molecular mass MW ≤ 500, the number of hydrogen bond donors HBD ≤ 5, the number of hydrogen bond acceptors HBA ≤ 10, the octanol-water partition coefficient logP ≤ 5, and the number of rotatable angles ROT ≤ 10 to screen drugs with good pharmacokinetic properties; the "Lipinski's Five Rules" properties and isomer SMILES strings of all drugs are obtained through the official website interface of the PubChem database according to their CID numbers; After applying the above step 1.1 to step 1.3 criteria, a drug database containing small molecule drugs and their CAS numbers, CID numbers, isomer SMILES strings, and "Lipinski's Rule of Five" properties is established.

3. An optimized de novo drug design method according to claim 1, characterized in that: In step 3, the skeleton-based similarity algorithm combines six similarity algorithms and four molecular representation methods to form a total of 24 combinations, and uses each combination to identify a skeleton similar to the target drug skeleton. The three most similar skeletons obtained from each combination are taken, and the duplicates are eliminated to obtain the final similar skeleton subset G1; the six similarity algorithms include Tanimoto, Dice, Cosine, Sokal, Russel, and Kulczynski, and the four molecular representation methods include topological fingerprints, MACCS keys, ECFP fingerprints, and FCFP fingerprints.

4. The optimized de novo drug design method according to claim 1, characterized in that: Step 6 specifically includes: The specific drug structure constraints of the MINLP model are given by equations (6) to (12); This means that a drug candidate selects only one skeleton from the skeleton collection; Among them, the structural constraints of chemical complexity are: n≤3,i∈G2(8) The structural constraints of the octet rule are: The structural constraints of the valence bond rule are: In equations (6) to (12), n i is the number of fragments i involved in the candidate drug, v i is the number of valence bonds of fragment i, G is the fragment set, G1 is the backbone subset, G2 is the group set, G3 is the subset of G2, G3 = {-CH3, -CH2, -CH, CH2 = CH-, -CH = CH-, CH2 = C<, -CH = C<}; is the number of i1 fragments, is the number of i2 segments, is the bond number of the i1 fragment.

5. The optimized de novo drug design method according to claim 1, characterized in that: Step 7 specifically includes: Step 7.1: In the MINLP model, the "Lipinski Rule of Five" properties are calculated using the RDKit-based quantitative estimation method for drug similarity, QED; the upper limits of the "Lipinski Rule of Five" properties are MW ≤ 500, HBD ≤ 5, HBA ≤ 10, logP ≤ 5, ROT or ROT frag ≤10; ROT in equation (3) frag is the linear sum of the rotatable angles of the fragments involved in the candidate drug, while ROT in equation (4) is the overall rotatable angle of the molecule calculated based on the SMILES string of the candidate drug, which is a nonlinear property; considering ROT frag ≤ROT,ROT frag >10 drugs will not satisfy the constraint of ROT≤10; therefore, when using the decomposition algorithm to solve the MINLP model, ROT is introduced frag It can eliminate some infeasible solutions before calculating ROT, which helps to improve the solution efficiency of MINLP model; Step 7.2: SA and SC properties are used to ensure that the designed candidate drug molecule is easy to synthesize, where low SA and SC values ​​indicate that the molecule is easy to synthesize.

6. The optimized de novo drug design method according to claim 1, characterized in that: Step 8 specifically includes: The SMILES string of the candidate drug is the predicted Prob bind , ROT, logP, SA and SC properties of the model input information; a "bridge" needs to be established in the structural constraints to associate the SMILES string with the fragment set, otherwise the MINLP model cannot be successfully solved; therefore, by using the improved SMILES-based isomer generation algorithm, the fragment set of the candidate drug can be automatically converted into its corresponding SMILES string, so that the MINLP model can be successfully solved, and the algorithm can identify candidate drug isomers containing the same fragment set; The improved SMILES-based isomer generation algorithm is implemented through the RWMol module in the RDKit library: a molecule is represented by a set of fragments, which include skeletons and groups. This information will be fed into the isomer generation algorithm; the specific algorithm execution process is as follows: (1) Select the skeleton as the "seed", and then combine the groups as "leaves" with the skeleton according to the pre-set group addition order; check whether there is a "*" match in the seed-leaf structure; if there is, add another group to the seed-leaf structure according to the pre-set group addition order; if not, proceed to the next step; (2) Check whether all the groups have been added to the seed-leaf structure; if not, it means that there are extra groups that have not been added to the seed-leaf structure. In this case, delete the seed-leaf structure and repeat step (1) according to the next pre-set group addition sequence; if all the groups have been added, proceed to the next step; (3) Check whether the generated seed-leaf structure can generate a SMILES structure through the RDKit library; if not, repeat steps (1)-(2) according to the next pre-set group addition sequence; if yes, save the SMILES result and repeat steps (1)-(2) according to the next pre-set group addition sequence; after trying all pre-set group addition sequences, all possible SMILES structures are obtained.

7. The optimized de novo drug design method according to claim 1, characterized in that: Step 9 specifically includes: Step 9.1 uses a decomposition solution algorithm to solve the MINLP optimization model: decompose the MINLP model into a mixed integer linear programming MILP subproblem and three nonlinear programming NLP subproblems; Step 9.2: Subproblem 1 - MILP: First restrict the structural constraints of the octet rule, valence bond rule and chemical complexity, as well as MW, HBA, HBD and ROT frag The linear property constraints are used in GAMS to generate N1 feasible solutions using the BARON solver. The feasible solutions are candidate drugs represented by the fragment set. A parameter N1 needs to be set in GAMS. max To specify the maximum number of solutions of the MILP model; Step 9.3: Sub-problem 2 - NLP: Using the constraints of the improved SMILES-based isomer generation algorithm, generate SMILES strings of N2 candidate drugs based on the N1 fragment set, N2 ≥ N1; Step 9.4: Subproblem 3 - NLP: Consider the nonlinear property constraints of ROT, logP, SA and SC, use the nonlinear property prediction model to calculate the corresponding properties of the SMILES strings of N2 candidate drugs, and eliminate those that do not meet the constraints; select the remaining N3 candidate drug SMILES strings for further evaluation; Step 9.5: Subproblem 4 - NLP: Considering the constraints of the deep learning model, use the nonlinear deep learning model to calculate the objective function Prob of N3 target-ligand complexes bind , and sort the SMILES strings of candidate drugs according to the objective function.

Images