Intelligent molecular design method based on mathematical programming method and skeleton-group coupling

Abstract

The invention discloses an intelligent molecular design method based on a mathematical programming method and skeleton-group coupling. According to the algorithm, a Beis-Murcko skeleton and a group fragment are taken as structural members, a potential skeleton subset is obtained through skeleton similarity retrieval, and a candidate structure space is constructed; on the basis, fragment selection and counting of candidate molecules are modeled as a mixed integer nonlinear programming (MINLP) problem, and structural feasibility constraints, property threshold constraints and scores / probabilities output by a deep learning model are uniformly introduced as target functions or constraint conditions. Meanwhile, an improved SMILES-based structure generation mechanism is adopted in the algorithm to establish consistency bridging between a fragment set and structure character string representation in the solving process, a fragment feasible solution is generated in combination with a decomposition solving strategy, and then structure representation generation and nonlinear property / model evaluation screening are carried out. And finally, outputting a molecular structure which meets the constraint and is optimal or approximately optimal in target.

Description

Technical Field

[0001] This invention belongs to the fields of computer-aided molecular design, molecular informatics and optimization computation, and specifically relates to a molecular intelligent design algorithm and its system implementation method that combines skeleton-group fragmentation representation, machine learning / deep learning property prediction and mathematical programming (mixed integer nonlinear programming MINLP). Background Technology

[0002] Molecular structure design is a typical problem at the intersection of chemistry and information technology. Its goal is to screen or generate molecular structures that meet specified objectives from a vast chemical structure space within acceptable time and computational resources. These objectives may be reflected in property indicators (such as stability, solubility, hydrophobicity, polarity correlation indicators, synthetic feasibility indicators, etc.), structural constraints (such as connection sites, valence rules, structural complexity restrictions, etc.), or scores / probabilities output by data-driven models (such as activity / performance probabilities or comprehensive scores under a specific task).

[0003] In multiple scenarios such as pharmaceuticals, materials, catalysis, solvents and functional chemicals, molecular structure design generally has two common characteristics: (1) the structure space is discrete and the combination is explosive: when the combination of "skeleton + group / substituent" is adopted, the number of candidates increases exponentially with the scale of the skeleton set and the group set; (2) the constraints and objectives are diverse and coupled: there is often more than one design objective, and there are multiple hard constraints (must be met) and soft constraints (try to meet them as much as possible), which are strongly coupled with the structure selection process.

[0004] In existing technologies, a common approach is to generate a large number of candidate molecules based on rules or enumeration, and then filter and rank them through property calculations or predictions. This approach is relatively straightforward, but when the structure space is large and there are many constraints, it generates a large number of candidates with invalid structures or that clearly do not meet the property thresholds, causing a large amount of computational resources to be consumed on invalid samples. Furthermore, structural and property constraints can often only be verified after generation, making it difficult to guarantee that "structure is compliant" during the generation process. Another approach is to use heuristic search or evolutionary algorithms to find better structures by iteratively improving the objective function. Although this can reduce the enumeration scale to some extent, it often requires complex penalty terms or repair strategies to maintain feasibility in scenarios with multiple constraints. The results are sensitive to hyperparameters and initial conditions, lack stability and repeatability, and are difficult to provide a "optimal or near-optimal" interpretation with clear solution semantics. In recent years, data-driven generative models (such as latent variable-based models, adversarial generative models, reinforcement learning, or diffusion models) have been used for molecular generation and optimization. They can learn data distributions and generate structures. However, in the case of multiple objectives and multiple hard constraints, there are still common problems such as difficulty in guaranteeing hard constraints, unstable out-of-distribution generalization, and difficulty in unifying the two-stage process of "generate first and then verify".

[0005] In engineering practice, to enable controlled searches within the structural space of combinatorial explosion, candidate molecules are often represented as "the selection and counting of skeleton and group fragments," that is, using discrete decision variables to characterize structural combinations. This fragmented representation facilitates the introduction of structural constraints such as connection sites, valence rules, skeleton uniqueness, and complexity thresholds, and is naturally compatible with discrete optimization frameworks such as integer programming. However, most property prediction models and task scoring models (especially deep learning models) typically require structural representations such as SMILES, molecular diagrams, or 3D conformations as input. Therefore, a natural disconnect exists between the "fragment decision variable space" and the "structural representation space": the optimization model can select combinations in the fragment space, but model evaluation often relies on structural strings or graph representations to calculate property or probability scores. If a bridging mechanism that can strictly constrain and automatically generate consistent structural representations during the optimization process is lacking, the optimization solution cannot stably call property prediction or deep models, and problems such as the fragment combinations obtained from the solution failing to generate valid structural representations, isomer ambiguity leading to unreproducible results, and structural feasibility relying solely on post-processing checks weakening overall consistency may arise.

[0006] Furthermore, when it is necessary to simultaneously satisfy the objective maximization or constraint satisfaction of structural rules, property thresholds, and model prediction output, a single sampling or post-filtering method is insufficient to express all constraints within a unified framework. Mathematical programming models such as mixed-integer nonlinear programming can uniformly express discrete structure selection, linear property constraints, and nonlinear property / model prediction constraints, but direct solutions often become difficult due to the large scale of discrete variables, strong coupling between nonlinearity and discreteness, and the nested relationship between structural representation generation and model evaluation, resulting in low solution efficiency, a high proportion of infeasible solutions, and unstable engineering implementation.

[0007] Therefore, there is an urgent need for a molecular intelligent design algorithm with an algorithm at its core: based on the fragmented representation of the skeleton-group, it integrates structural feasibility constraints and property constraints into a solvable mathematical programming framework, and uses a constrained structural string generation bridging mechanism to consistently associate fragment variables with structural representations. At the same time, it uses a decompositional solution strategy to first use linear and discrete constraints to quickly screen out large-scale infeasible combinations, and then performs nonlinear property prediction and deep model evaluation on a small number of candidates, thereby significantly improving search efficiency and result stability while ensuring feasibility. Summary of the Invention

[0008] The purpose of this invention is to provide a molecular intelligent design algorithm and implementation process with an algorithm as its core. This algorithm defines a controllable candidate structure space using a skeleton-group fragmentation representation. It integrates structural feasibility rules, property threshold requirements, and scores or probability assessments given by machine learning / deep learning models into a unified mathematical programming framework. By constructing a mixed integer nonlinear programming (MINLP) model and introducing consistency bridging constraints between fragment variables and structural string representations, it enables the stable generation of legal structural representations and the invocation of prediction models to complete target evaluation during the optimization process. Furthermore, by combining a decompositional solution strategy, it efficiently filters out infeasible solutions in a large-scale combinatorial space and outputs molecular structures that satisfy the constraints and are optimal or near-optimal to the target. This achieves generalized, interpretable, and implementable intelligent optimization design for multi-constraint, multi-objective molecular design tasks.

[0009] To achieve the above objectives, the present invention adopts the following technical solution, including the following specific steps: A molecular intelligent design method based on mathematical programming and skeleton-group coupling includes the following specific steps: Step 1: Establish a molecular database; Step 2: Establish a database containing the Bemis-Murcko skeleton: Extract the Bemis-Murcko skeleton from molecular structures in the molecular database using the Bemis–Murcko algorithm in RDKit; Step 3: Based on the skeletal structure of the target molecule, use a skeletal similarity algorithm to search for a subset of skeletal structures similar to the target molecule's skeletal structure from the skeletal structure database. At the same time, a set of commonly used groups is selected to form a collection. The skeleton subset G1 and the group set G2 together define the fragment space of the candidate molecule, and the candidate molecule is generated by the combination of "skeleton fragment + group fragment". Step 4: Designing promising candidate molecules: The problem of selecting fragment combinations for candidate molecules is formulated as a mixed-integer nonlinear programming (MINLP) model. The MINLP model uses fragment selection and fragment counting as core decision variables, structural and property constraints as feasibility criteria, and an objective function to drive a global optimization search for candidate molecule fragment combinations. By solving the MINLP model, the optimal candidate molecule structures with high binding affinity probabilities are identified from the feasible solutions. These feasible solutions are generated by optimizing the molecular skeleton and functional groups, and the combination process is constrained by the MINLP model. The objective function is to maximize the molecular properties to be optimized. 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 used to predict the molecular properties to be optimized;

[0010] Step 6: Molecular structure constraints: General equation (2) represents the octet rule Price key rules and chemical complexity The structural constraints allow for the generation of structurally rational molecules through the combination of the skeleton and groups.

[0011] Step 7: Molecular property constraints: General equations (3) and (4) represent the properties predicted by the groups and SMILES, respectively;

[0012] Step 8: Other constraints: General equation (5) represents the improved SMILES-based isomer generation algorithm, which is used to automatically convert the fragment set of candidate molecules into the corresponding molecular SMILES string. The fragments of the candidate molecules include the backbone and groups.

[0013] In the above general equation, F obj It is the objective function. Property These are molecular properties that need to be optimized. i The fragments involved in the candidate molecules, n i This represents the number of fragments involved in the candidate molecule. s It is represented by SMILES for molecules. m It is the type of structural constraint. and It is a structural constraint m upper and lower boundaries, p It is a molecular property. k It is a property type. and yes p k The upper and lower bounds; Step 9: Use a decomposition algorithm to solve the MINLP model; if there is no best candidate molecule that meets all constraints, return to step 4 to relax the constraint range.

[0014] Furthermore, step 1 specifically includes: Step 1.1: Collect small molecule structures with unique identifiers from public or private data sources; retrieve corresponding structure records in the structure database based on the unique identifiers; delete molecules for which no valid structure records can be obtained, and the remaining molecules can be further screened; Step 1.2: Screen molecules for properties based on a preset set of property thresholds, and obtain the property information of the screened molecules and the SMILES string of isomers; After applying the standards in steps 1.1 to 12 above, a molecular database containing small molecules and their identifiers, structural representations, SMILES strings of isomers, and property information is established.

[0015] Furthermore, in step 3, the skeleton-based similarity algorithm combines six similarity algorithms and four molecular representation methods to form a total of 24 combinations. Each combination is used to identify skeletons similar to the target molecule skeleton. The three most similar skeletons obtained from each combination are selected, and duplicates are removed to obtain the final subset of similar skeletons. The six similarity algorithms include Tanimoto, Dice, Cosine, Sokal, Russell, and Kulczynski, and the four molecular representation methods include topological fingerprinting, MACCS key, ECFP fingerprinting, and FCFP fingerprinting.

[0016] Furthermore, step 6 specifically includes: The specific molecular structure constraints of the MINLP model are given by equations (6) to (12); structural constraints This means that a candidate molecule selects only one skeleton from the skeleton set; Among them, the structural constraints of chemical complexity are:

[0017] In equations (6) to (12), It is the fragment involved in the candidate molecule. Quantity, It is a fragment valence bond number, It is a collection of fragments. It is a subset of the skeleton. It is a group set. yes a subset of ={-CH3, -CH2, -CH, CH2=CH-, -CH=CH-, CH2=C<, -CH=C<}; yes The number of segments, yes The number of segments, yes The number of valence bonds in a fragment.

[0018] Furthermore, step 7 specifically includes: in the MINLP model, the property set used to screen candidate molecules is entered into the decomposition solution process in two forms, linear and nonlinear; Furthermore, step 8 specifically includes: The SMILES string of the candidate molecule is a prediction Property The model input information has the same properties; in the structural constraints, a "bridge" needs to be established to associate the SMILES strings with the fragment sets, otherwise the MINLP model cannot be solved successfully; therefore, by using an improved SMILES-based isomer generation algorithm, the fragment sets of candidate molecules can be automatically converted into their corresponding SMILES strings, enabling the MINLP model to solve successfully. This algorithm can also identify isomers of candidate molecules containing the same fragment set. The improved SMILES-based isomer generation algorithm is executed through the RWMol module in the RDKit library: a molecule is represented by a set of fragments, including the backbone and functional groups, and this information is fed into the isomer generation algorithm; the specific algorithm execution flow is as follows: (1) Select the skeleton as the "seed", and then add the groups as "leaves" to the skeleton according to the pre-set group addition order; check whether there are any "seeds" in the seed-leaf structure. If it meets the requirements, add another group to the seed-leaf structure according to the pre-defined group addition order; if not, proceed to the next step. (2) Check whether all groups have been added to the seed-leaf structure; if not all have been added, it means that there are extra groups that have not been added to the seed-leaf structure. At this time, delete the seed-leaf structure and repeat step (1) according to the next preset group addition order; if all have been added, proceed to the next step. (3) Check whether the generated seed-leaf structure can generate SMILES structure through the RDKit library; if not, repeat steps (1)-(2) according to the next preset group addition order; if it can, save the SMILES result and repeat steps (1)-(2) according to the next preset group addition order; after trying all preset group addition orders, obtain all possible SMILES structures.

[0019] Furthermore, step 9 specifically includes: Step 9.1: Solve the MINLP optimization model using a decomposition algorithm: 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 structure by adhering to the octet rule, valence bond rule, and chemical complexity, as well as linear property constraints. Then, use the BARON solver in GAMS to generate... N One feasible solution, which is a candidate molecule represented by a fragment set; a parameter needs to be set in GAMS. N 1 max To define the maximum number of solutions for the MILP model; Step 9.3: Subproblem 2 – NLP: Using constraints from an improved SMILES-based heterogeneity generation algorithm, based on... N 1. Fragment set generation N The SMILES strings of the two candidate molecules, N 2≥ N 1; Step 9.4: Subproblem 3 – NLP: Considering nonlinear property constraints, use the nonlinear property prediction model to calculate N The corresponding properties of the SMILES strings of the two candidate molecules are analyzed, and those that do not meet the constraints are eliminated; the remaining molecules are selected. N The SMILES strings of the three candidate molecules were further evaluated; Step 9.5: Sub-problem 4 – NLP: Considering the constraints of the deep learning model, use a non-linear deep learning model to calculate the objective function. Property The SMILES strings of candidate molecules are sorted according to the objective function.

[0020] Compared with the prior art, the present invention has the following main advantages: (1) First, this invention extends the traditional CAMD (Computer-Aided Molecular Design) concept to a more complex structural space and establishes a solvable optimization framework based on mathematical programming (MINLP). Compared to traditional CAMD tasks targeting relatively simple structures such as small molecule solvents, candidate molecular structures with skeleton-group combinations often have higher structural complexity and a larger combinatorial space. This significantly increases the size of the MINLP model and reduces the solution efficiency, especially in cases involving multi-ring structures or multiple connection sites, requiring more structural fragments and counting variables to form MILP subproblems (subproblem 1 in the decomposition algorithm), making it more difficult for the model to obtain a high-quality feasible solution within engineering time. To address this challenge, this invention improves the solvability and search efficiency of the model through structural space reconstruction and constraint design at the algorithm level.

[0021] (2) Secondly, even if the MILP subproblem can generate structurally feasible combinations of fragments, traditional CAMD is still prone to generating structures that formally satisfy the rules but are unreasonable in terms of chemical intuition or engineering usability during the combination process. This is typically manifested as abnormal stacking of ring structures or substitution patterns, resulting in anomalous structures that satisfy local constraints. To solve this type of "abnormal but feasible structure" problem, this invention introduces a skeleton representation extracted based on the Bemis–Murcko algorithm and embeds the skeleton as a structural prior and constraint carrier into the CAMD optimization process, so that the candidate molecules maintain rationality and interpretability at the global structural level. At the same time, in order to avoid incorporating the entire skeleton into the model and causing a solution burden, this invention further adopts a skeleton-based similarity retrieval strategy to identify a subset of skeletons that are similar to the target structure before modeling, and perform combination optimization only within the potential skeleton range, thereby significantly compressing the candidate space and reducing the MILP scale, and improving the overall computational efficiency.

[0022] (3) Furthermore, this invention deeply couples the data-driven prediction model as an embeddable algorithm component with the optimization framework. By introducing deep learning model constraints / objectives, the scores or probabilities of candidate molecules can be uniformly expressed and optimized in the mathematical programming solution process. This invention can use advanced deep learning technology to construct a prediction model and use its output as part of the objective function or constraint conditions of the MINLP model, thereby forming a closed-loop molecular intelligent design algorithm process of "prediction-optimization-output".

[0023] In summary, this invention establishes a molecular intelligent design algorithm framework that is centered on mathematical programming, based on skeleton-group fragmentation representation, uses similarity retrieval to compress the search space, and is coupled with a deep learning prediction model. This framework can efficiently generate candidate molecules that meet structural and property requirements under multiple constraints and large-scale combinatorial spaces, reducing invalid searches and redundant calculations, improving the stability, interpretability, and engineering feasibility of optimization solutions, and facilitating the discovery of candidate molecules with novel structural features, thereby improving design efficiency and diversity. Attached Figure Description

[0024] Figure 1 This is a flowchart of the molecular intelligent design framework based on mathematical programming and skeleton-group in this invention; Figure 2 This is a schematic diagram of the molecular structure of the present invention, consisting of a backbone (ring structure and connecting structure) and side chains.

[0025] Figure 3 This is the chemical space of the candidate molecule design results plotted using ECFP fingerprinting and principal component analysis (PCA) methods according to the present invention. (Integers (0~14) in the right-hand legend represent 15 skeletons, and dot "7" represents a commercial molecule). Detailed Implementation

[0026] The process and effects of the present invention will be explained and described below with reference to the accompanying drawings and embodiments. The embodiments of the present invention are implemented under the premise of the technical solution of the present invention, and detailed implementation methods and specific operation processes are given. The embodiments use drug molecule design as an application example, and the database sources, target selection, molecular docking and molecular dynamics verification involved are all exemplary settings, but the protection scope of the present invention is not limited to the following embodiments.

[0027] See Figure 1 This invention specifically discloses a molecular intelligent design method based on mathematical programming and skeleton-group coupling, comprising the following steps: Step 1: A drug database was created using the DrugBank V5.0 database (https: / / go.drugbank.com / ) for subsequent drug skeleton generation.

[0028] Step 2: Generate the drug backbone for subsequent MINLP-based drug design models. The backbone (e.g., ...) is used in the CAMD method. Figure 2 By replacing traditional cyclic groups, more rationally structured drug candidates can be designed. Therefore, a database containing 2,898 Bemis-Murcko skeletons was established for subsequent drug candidate design. These skeletons were generated using the Bemis–Murcko algorithm in RDKit from 4,781 drugs in the drug database.

[0029] Step 3: To reduce unnecessary search space and save computational time costs in candidate drug design, a skeleton-based similarity algorithm is used to identify a subset of skeletons similar to those of axitinib (an FDA-approved anti-tumor drug) before building the MINLP model. A total of 72 skeletons were obtained, and after removing duplicates, 19 skeletons were obtained. Commonly used group sets were then selected to form a set. .

[0030] Step 4: High-throughput, intelligent design of promising drug candidates using optimized mathematical programming methods. The ab initio 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 candidate with high binding affinity among the feasible solutions can be identified. These feasible solutions are generated by optimizing the combination of drug skeletons and group sets under model constraints. Note that each MINLP model is applicable to one target and can be used to design promising drug candidates for it. The MINLP model formula is as follows: Objective function: Subject to steps 5-8: Step 5: Deep learning constraints: The general equation (Equation (1)) represents the deep learning model used to identify target-ligand complexes with high binding affinity.

[0031] Step 6: Drug structure constraints: The general equation (equation (2)) represents the octet rule ( Price key rules ) and chemical complexity ( The structural constraints of the skeleton and groups allow for the generation of structurally rational molecules.

[0032] Step 7: Drug property constraints: The general equation (Equation (3~4)) represents the "Ribinsky Five Rules" property (relative molecular mass) MW The number of hydrogen bond acceptors HBA The number of hydrogen bond donors HBD Octyl alcohol-water partition coefficient log P Rotatable angles ROT ( ROT frag )), composite feasibility score ( SA ) and compositional complexity score ( SC ).

[0033] Step 8: Other constraints: The general equation (Equation (5)) represents the improved SMILES-based isomer generation algorithm, which is used to automatically convert the fragment (backbone and group) set of candidate drugs into the corresponding drug SMILES string.

[0034] In the above general equation, F obj It is the objective function. Prob bind It is the probability of high binding affinity of the target-ligand complex. i This represents the fragments (skeleton and groups) involved in the candidate drug. n i This represents the number of fragments involved in the candidate drug. s It is represented by SMILES for molecules. m It is the type of structural constraint. and It is a structural constraint m upper and lower boundaries, p It is a drug property. kIt is a property type. and yes p k The upper and lower boundaries.

[0035] Step 9: Since the MINLP model involves a large number of nonlinear constraints, a decomposition algorithm is used to solve the MINLP model. If no optimal candidate drug meets all constraints, return to step 4 to relax the constraint range.

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

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

[0038] Step 1.2: Then, using a web crawler script and CAS numbers, search the PubChem database for the CID numbers (Unique Identifiers of Chemicals in the PubChem database, https: / / pubchem.ncbi.nlm.nih.gov / search / ) of 7,746 drugs. After removing drugs without CID numbers, 7,474 drugs remain for further screening.

[0039] Step 1.3: Use the "Ribinsky Five Rules" property (relative molecular mass ( MW ≤500, the number of hydrogen bond donors ( HBD ≤5, the number of hydrogen bond acceptors ( HBA Octyl alcohol-water partition coefficient ≤10 (log P )≤5, number of rotatable angles ( ROT Drugs with good pharmacokinetic properties were screened using a criterion of ≤10. The Ribinski Five Rules properties and isomer SMILES strings for all 7,474 drugs were obtained via the official PubChem database interface based on their CID numbers.

[0040] Finally, after applying the above standards, a drug database containing 4,781 small molecule drugs and their CAS numbers, CID numbers, isomer SMILES strings, and "Ribinski's Five Rules" properties was established.

[0041] Furthermore, step 3 specifically includes: Step 3.1: If all 2,898 skeletons are used in a single MINLP-based drug design model, the excessively large model size will hinder the efficient solution of the MINLP model. Therefore, we propose a skeleton-based similarity algorithm to identify skeletons similar to the target drug skeleton from the 2,898 skeletons. The subset of similar skeletons (determined by case studies) is obtained by combining six similarity algorithms and four molecular representation methods to form 24 possible combinations. Each combination identifies skeletons similar to the target drug skeleton, and the three most similar skeletons from each combination are selected, with duplicates removed to obtain the final subset of similar skeletons. The six similarity algorithms include Tanimoto, Dice, Cosine, Sokal, Russell, and Kulczynski, and the four molecular representation methods include topological fingerprinting, MACCS key, ECFP fingerprinting, and FCFP fingerprinting.

[0042] Step 3.2: After deleting duplicate skeletons, a certain number of similar skeletons can be obtained depending on the specific problem. In addition, 29 commonly used group sets were selected for our MINLP model ( (See Table 1).

[0043] Table 1 Selected group groups

[0044] Furthermore, step 6 specifically includes: The specific drug structure constraints of the MINLP model are given by equations (6-12). This indicates that a candidate drug selects only one skeleton from the skeleton group.

[0045] In equations (6~12), These are the fragments (skeleton and groups) involved in the candidate drug. Quantity, It is a fragment valence bond number, It is a collection of fragments. It is a skeleton set. It is a group set. yes subset of ={-CH3, -CH2, -CH, CH2=CH-, -CH=CH-, CH2=C<, -CH=C<}); yes The number of segments, yes The number of segments, yes The number of valence bonds in a fragment.

[0046] Furthermore, step 7 specifically includes: Step 7.1: In the MINLP model, the properties of the "Ribinski Five Rules" ( MW , HBA , HBD , log P , ROT The result was calculated using the RDKit-based quantitative drug similarity estimation method (QED); the upper limit of the Ribinski Five Rules property is... MW ≤500, HBD ≤5, HBA ≤10, log P ≤5 and ROT ( ROT frag )≤10; in equation (3) ROT frag It is the linear sum of the number of rotatable angles of the fragments involved in the candidate drug, while in equation (4) ROT The rotatable angles of the entire molecule are calculated based on the SMILES string of the candidate drug, which is a non-linear property; considering ROT frag ≤ ROT , ROT frag Drugs with a value >10 will not meet the requirements. ROT The constraint is ≤10; therefore, when using a decomposition algorithm to solve the MINLP model, the following is introduced: ROT frag Capable of calculating ROT Previously, some infeasible solutions were removed, which helps improve the solution efficiency of the MINLP model.

[0047] Step 7.2: SA (1~10) and SC (1~5) Properties are used to ensure that the designed drug candidate is easy to synthesize, among which low SA and SC The value indicates that the molecule is easy to synthesize. Based on our empirical knowledge, their constraint range is set to... SA ≤6 and SC ≤3.4.

[0048] Furthermore, step 8 specifically includes: The SMILES string of the candidate drug is a prediction Prob bind , ROT , log P , SA and SC The model input information is of a certain nature; in the structural constraints, a "bridge" needs to be established to associate the SMILES string with the fragment set, otherwise the MINLP model cannot be solved successfully; therefore, by using an improved SMILES-based isomer generation algorithm, the fragment set of candidate drugs can be automatically converted into its corresponding SMILES string, so that the MINLP model can be solved successfully. At the same time, the algorithm can identify isomers of candidate drugs that contain the same fragment set.

[0049] The improved SMILES-based isomer generation algorithm is executed through the RWMol module in the RDKit library: a molecule is represented by a set of fragments, including the backbone and functional groups, and this information is fed into the isomer generation algorithm; the specific algorithm execution flow is as follows: (1) Select the skeleton as the "seed", and then add the groups as "leaves" to the skeleton according to the pre-set group addition order; check whether there are any "seeds" in the seed-leaf structure. If it meets the requirements, add another group to the seed-leaf structure according to the pre-defined group addition order; if not, proceed to the next step. (2) Check whether all groups have been added to the seed-leaf structure; if not all have been added, it means that there are extra groups that have not been added to the seed-leaf structure. At this time, delete the seed-leaf structure and repeat step (1) according to the next preset group addition order; if all have been added, proceed to the next step. (3) Check whether the generated seed-leaf structure can generate SMILES structure through the RDKit library; if not, repeat steps (1)-(2) according to the next preset group addition order; if it can, save the SMILES result and repeat steps (1)-(2) according to the next preset group addition order; after trying all preset group addition orders, obtain all possible SMILES structures.

[0050] Furthermore, step 9 specifically includes: Step 9.1: Generally speaking, if the nonlinear equations in the MINLP model are not particularly complex, the MINLP model can be solved directly using the BARON solver in the GAMS software (https: / / www.gams.com / ).

[0051] Step 9.2: However, our MINLP model, based on SMILES-based property prediction models (e.g., deep learning models), contains extremely complex nonlinear equations, making direct solution of the MINLP model extremely difficult. Therefore, a decomposition-based solution algorithm is employed 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.

[0052] Step 9.3: Subproblem 1 (MILP): First, constrain the structure by adhering to the octet rule, valence bond rule, and chemical complexity. MW , HBA , HBD and ROT frag The linear property constraints are used in GAMS to generate a certain number of ( ) N 1) Feasible solutions (candidate drugs represented by a set of fragments). Note that a parameter needs to be set in GAMS. N 1 max This defines the maximum number of solutions for the MILP model. Using mathematical programming methods, under constraints of structural and linear properties, we obtain... N 1 = 33,759 feasible solutions (represented by a set of fragments), taking 210 seconds on our desktop computer (Intel(R) Core(TM) i7-10700F CPU @2.90GHz 24.0 GB RAM).

[0053] Step 9.4: Subproblem 2 (NLP): Using the constraints of the improved SMILES-based heterogeneity generation algorithm, respectively based on... N Generation of a set of fragments N 2 = 262,741 candidate drug SMILES strings (3,113 seconds).

[0054] Step 9.5: Sub-problem 3 (NLP): Consider ROT , log P , SA and SC The nonlinear property constraints are used to calculate the nonlinear property prediction model. N Determine the corresponding properties of the SMILES strings for the two candidate drugs and remove those that do not meet the constraints. Remaining... N The SMILES strings of 3 = 105,164 candidate drugs were used for further analysis.

[0055] 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, indicating that our MINLP-based drug design model can not only find existing drug candidates, but also design new drug candidates (95.36%).

[0056] Secondly, chemical spaces were created using ECFP fingerprinting and principal component analysis (PCA) methods (see...). Figure 3 To characterize the structural diversity of 105,164 designed drug candidates, among which x shaft and y The axes are PC1 and PC2, two principal components. In the legend on the right, integers (0~14) represent 15 skeletons, and point "7" represents axitinib. Figure 3 This demonstrates the wide distribution of the designed drug candidates in chemical space, showcasing the powerful capability of our MINLP-based drug design model in designing a variety of structurally diverse drug candidates similar to axitinib.

[0057] Step 9.6: Sub-problem 4 (NLP): Considering the constraints of the deep learning model, use a non-linear deep learning model to compute... N 3 target-ligand complexes Prob bind (Objective function) was used, and the SMILES strings of candidate drugs were sorted according to the objective function. The ranking results showed that 433 designed candidate drugs outperformed axitinib (97.96%) in terms of high binding affinity. Some of the top-ranked candidate drugs were further validated by other physics-based methods such as molecular docking and molecular dynamics simulations.

[0058] The foregoing description of specific exemplary embodiments of the invention is for illustrative and explanatory purposes. These descriptions are not intended to limit the invention to the precise forms disclosed, and it will be apparent that many changes and variations can be made in accordance with the foregoing teachings. The exemplary embodiments were chosen and described in order to explain the specific principles of the invention and its practical application, thereby enabling those skilled in the art to implement and utilize various different exemplary embodiments of the invention, as well as various different choices and variations. The scope of the invention is intended to be defined by the claims and their equivalents.

Claims

1. A molecular intelligent design method based on mathematical programming and skeleton-group coupling, characterized in that, The specific steps include the following: Step 1: Establish a molecular database; Step 2: Establish a database containing the Bemis-Murcko skeleton: Extract the Bemis-Murcko skeleton from molecular structures in the molecular database using the Bemis–Murcko algorithm in RDKit; Step 3: Based on the skeletal structure of the target molecule, use a skeletal similarity algorithm to search for a subset of skeletal structures similar to the target molecule's skeletal structure from the skeletal structure database. At the same time, a set of commonly used groups is selected to form a collection. The skeleton subset G1 and the group set G2 together define the fragment space of the candidate molecule, and the candidate molecule is generated by the combination of "skeleton fragment + group fragment". Step 4: Designing promising candidate molecules: The problem of selecting fragment combinations for candidate molecules is formulated as a mixed-integer nonlinear programming (MINLP) model. The MINLP model uses fragment selection and fragment counting as core decision variables, structural and property constraints as feasibility criteria, and an objective function to drive a global optimization search for candidate molecule fragment combinations. By solving the MINLP model, the optimal candidate molecule structures with high binding affinity probabilities are identified from the feasible solutions. These feasible solutions are generated by optimizing the molecular skeleton and functional groups, and the combination process is constrained by the MINLP model. The objective function is to maximize the molecular properties to be optimized. 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 used to predict the molecular properties to be optimized; Step 6: Molecular structure constraints: General equation (2) represents the octet rule Price key rules and chemical complexity The structural constraints allow for the generation of structurally rational molecules through the combination of the skeleton and groups. Step 7: Molecular property constraints: General equations (3) and (4) represent the properties predicted by the groups and SMILES, respectively; Step 8: Other constraints: General equation (5) represents the improved SMILES-based isomer generation algorithm, which is used to automatically convert the fragment set of candidate molecules into the corresponding molecular SMILES string. The fragments of the candidate molecules include the backbone and groups. In the above general equation, F obj It is the objective function. Property These are molecular properties that need to be optimized. i The fragments involved in the candidate molecules, n i This represents the number of fragments involved in the candidate molecule. s It is represented by SMILES for molecules. m It is the type of structural constraint. and It is a structural constraint m upper and lower boundaries, p It is a molecular property. k It is a property type. and yes p k The upper and lower bounds; Step 9: Use a decomposition algorithm to solve the MINLP model; if there is no best candidate molecule that meets all constraints, return to step 4 to relax the constraint range.

2. The molecular intelligent design method based on mathematical programming and skeleton-group coupling according to claim 1, characterized in that, Step 1 specifically includes: Step 1.1: Collect small molecule structures with unique identifiers from public or private data sources; retrieve corresponding structure records in the structure database based on the unique identifiers; delete molecules for which no valid structure records can be obtained, and the remaining molecules can be further screened; Step 1.2: Screen molecules for properties based on a preset set of property thresholds, and obtain the property information of the screened molecules and the SMILES string of isomers; After applying the standards in steps 1.1 to 12 above, a molecular database containing small molecules and their identifiers, structural representations, SMILES strings of isomers, and property information is established.

3. The molecular intelligent design method based on mathematical programming and skeleton-group coupling 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. Each combination is used to identify skeletons similar to the target molecule skeleton. The three most similar skeletons obtained from each combination are selected, and duplicates are removed to obtain the final subset of similar skeletons. The six similarity algorithms include Tanimoto, Dice, Cosine, Sokal, Russell, and Kulczynski, and the four molecular representation methods include topological fingerprinting, MACCS key, ECFP fingerprinting, and FCFP fingerprinting.

4. The molecular intelligent design method based on mathematical programming and skeleton-group coupling according to claim 1, characterized in that, Step 6 specifically includes: The specific molecular structure constraints of the MINLP model are given by equations (6) to (12); structural constraints This means that a candidate molecule selects only one skeleton from the skeleton set; Among them, the structural constraints of chemical complexity are: In equations (6) to (12), It is the fragment involved in the candidate molecule. Quantity, It is a fragment The number of valence bonds, It is a collection of fragments. It is a subset of the skeleton. It is a group set. yes a subset of ={-CH3, -CH2, -CH, CH2=CH-, -CH=CH-, CH2=C<, -CH=C<}; yes The number of segments, yes The number of segments, yes The number of valence bonds in a fragment.

5. The molecular intelligent design method based on mathematical programming and skeleton-group coupling according to claim 1, characterized in that, Step 7 specifically includes: In the MINLP model, the property set used to screen candidate molecules is entered into the decomposition solution process in two forms, linear and nonlinear.

6. The molecular intelligent design method based on mathematical programming and skeleton-group coupling according to claim 1, characterized in that, Step 8 specifically includes: The SMILES string of the candidate molecule is a prediction Property The model input information has the same properties; in the structural constraints, a "bridge" needs to be established to associate the SMILES strings with the fragment sets, otherwise the MINLP model cannot be solved successfully; therefore, by using an improved SMILES-based isomer generation algorithm, the fragment sets of candidate molecules can be automatically converted into their corresponding SMILES strings, enabling the MINLP model to solve successfully. This algorithm can also identify isomers of candidate molecules containing the same fragment set. The improved SMILES-based isomer generation algorithm is executed through the RWMol module in the RDKit library: a molecule is represented by a set of fragments, including the backbone and functional groups, and this information is fed into the isomer generation algorithm; the specific algorithm execution flow is as follows: (1) Select the skeleton as the "seed", and then add the groups as "leaves" to the skeleton according to the pre-set group addition order; check whether there are any "seeds" in the seed-leaf structure. If it meets the requirements, add another group to the seed-leaf structure according to the pre-defined group addition order; if not, proceed to the next step. (2) Check whether all groups have been added to the seed-leaf structure; if not all have been added, it means that there are extra groups that have not been added to the seed-leaf structure. At this time, delete the seed-leaf structure and repeat step (1) according to the next preset group addition order; if all have been added, proceed to the next step. (3) Check whether the generated seed-leaf structure can generate SMILES structure through the RDKit library; if not, repeat steps (1)-(2) according to the next preset group addition order; if it can, save the SMILES result and repeat steps (1)-(2) according to the next preset group addition order; after trying all preset group addition orders, obtain all possible SMILES structures.

7. The molecular intelligent design method based on mathematical programming and skeleton-group coupling according to claim 1, characterized in that, Step 9 specifically includes: Step 9.1: Solve the MINLP optimization model using a decomposition algorithm: 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 structure by adhering to the octet rule, valence bond rule, and chemical complexity, as well as linear property constraints. Then, use the BARON solver in GAMS to generate... N One feasible solution, which is a candidate molecule represented by a fragment set; a parameter needs to be set in GAMS. N 1 max To define the maximum number of solutions for the MILP model; Step 9.3: Subproblem 2 – NLP: Using constraints from an improved SMILES-based heterogeneity generation algorithm, based on... N 1. Fragment set generation N The SMILES strings of the two candidate molecules, N 2≥ N 1; Step 9.4: Subproblem 3 – NLP: Considering nonlinear property constraints, use the nonlinear property prediction model to calculate N The corresponding properties of the SMILES strings of the two candidate molecules are analyzed, and those that do not meet the constraints are eliminated; the remaining molecules are selected. N The SMILES strings of the three candidate molecules were further evaluated; Step 9.5: Sub-problem 4 – NLP: Considering the constraints of the deep learning model, use a non-linear deep learning model to calculate the objective function. Property The SMILES strings of candidate molecules are sorted according to the objective function.

Images