Prediction method for forming crystal structures of two drugs

Abstract

The present invention provides a computer-aided drug design method for predicting the structure of a drug co-crystal or a co-amorphous substance. On the basis of obtaining drug structure information, whether the combination of the two drugs can form a hydrogen bond is calculated based on the electrostatic potential surface energy (MEPS) between the two drug molecules, and then a molecular dynamics simulation (MD) method and a molecular mechanics-Boltzmann Poisson surface area (MM-PBSA) free energy calculation are further adopted to establish a fast, efficient and accurate prediction system for the structure of the co-crystal or amorphous substance of the two drugs.

Description

Technical Field

[0001] The present invention relates to the field of computer-aided molecular design, and in particular to a method for predicting whether two drugs can form a co-crystal or a co-amorphous structure, and is particularly suitable for predicting drug crystal structures of psychotropic drugs. Background Art

[0002] Drug development is a long and costly task. For example, the development of drugs for the treatment of mental illness is extremely difficult to develop a safe and effective new drug in the absence of a relatively poor understanding of the complex pathology of mental illness. In order to solve the dilemma of new molecular entity drug development, scientists have come up with different solutions. Still taking drugs for the treatment of mental illness as an example, one of the solutions is to use a large number of existing drugs targeting mental illness in a multimodal combination in clinical practice to achieve the purpose of treatment. In fact, through drug repositioning (i.e., the process of reverse screening, combining or discovering new uses of existing drugs), and conventional experimental methods to screen whether two drugs can achieve enhanced effects when used simultaneously, research is widely underway. However, the problem that the dosage of the selected drug combination is different from that of a single drug is very serious in clinical practice, especially it is unrealistic to change the physicochemical properties of the drug and thus enhance the bioavailability of the drug. To solve these problems, many researchers have begun to study new drug crystal forms, drug co-crystals and drug co-amorphous forms. Because experimental methods for screening drug crystal structures are time-consuming and expensive, drug combination methods based on computer-assisted rational screening that can form two drug crystal structures are booming.

[0003] For drugs, there are usually many solid structures, among which cocrystals are formed by non-covalent bonding. Drug cocrystals are homogeneous crystalline solid substances composed of two or more small molecule compounds in a certain stoichiometric ratio. Drug cocrystals can be composed of active drug ingredients and comorphs, or can be composed of two active drug ingredients. Co-amorphous drugs are single-phase amorphous binary systems with a single glass transition temperature formed by the combination of active drug ingredients with other drugs or comorphs. Drug cocrystals or co-amorphous have the advantages of improving the physical and chemical properties of drugs such as solubility, dissolution rate, stability and bioavailability, and have become a new approach to drug research and development.

[0004] In recent years, great progress has been made in the screening of drug combinations that can form drug crystals using computational molecular electrostatic potential surface methods. The core of this method is to analyze the charge distribution between the active drug ingredients and the co-former, and to evaluate which drug combination is most likely to form crystals by ranking the strength of hydrogen bond donors and acceptors of different functional groups of the co-crystal components, thereby reducing the cost and workload of experimental screening. However, the molecular electrostatic potential surface method is limited in accurately predicting drug crystals composed of two active drug ingredients. This is mainly because the method only considers the strength of hydrogen bond forces that can be formed, and ignores other intermolecular forces, such as halogen bonds, π-π stacking, van der Waals forces, and electrostatic forces, resulting in current methods being limited to combinatorial screening in large databases. There is still an urgent need to establish a systematic, effective, and rapid computational method for screening and predicting the crystal structure of drug combinations. Summary of the invention

[0005] The purpose of the present invention is to use computer-aided drug design to establish a prediction method for the crystal structure composed of two active drug ingredients, so as to achieve the purpose of effectively and quickly predicting whether the combination of two drugs can form a co-crystal or a co-amorphous state by only using the structural information of the drugs.

[0006] The present invention is achieved through the following technical solutions:

[0007] When only drug information is known, the integration method based on the electrostatic potential surface (MEPS) between the two drug molecules is used to calculate whether the combination of the two drugs can form a hydrogen bond. The molecular dynamics simulation (MD) method and the molecular mechanics-Boltzmann surface area (MM-PBSA) free energy calculation are further used to establish a prediction system for the crystal structures of the two drugs.

[0008] Optionally, before performing the above prediction calculation or after obtaining the prediction calculation results, the information on drug combination applications can be combined to provide further screening or evaluation indicators for obtaining drug combinations and effective crystal structures that can produce synergistic effects.

[0009] A method for predicting whether two drug molecules can form a co-crystal or a co-amorphous structure comprises the following steps:

[0010] Obtain drug information for drug crystal structure prediction;

[0011] Based on the drug information, MEPS calculation of the combination of two drug molecules is performed to evaluate the strength of hydrogen bonds between drug molecules;

[0012] According to the drug information, the binding simulation calculation of the two drug molecules is performed to obtain the trajectory conformation set of the simulated binding between the two drug molecules;

[0013] Based on the acquired trajectory conformation information of the simulated binding between the two drug molecules, the intermolecular free energy calculation is performed to screen the binding conformation with the minimum free energy between the two drug molecules;

[0014] Based on the strength of the hydrogen bond between the two drug molecules and the strength of the non-covalent bond between the two drug molecules in the binding conformation with the minimum free energy, a structural prediction model for the formation of a co-crystal or a co-amorphous substance between the two drug molecules is screened.

[0015] According to the present invention, a step can be optionally added to the method, which is to obtain information on the combined application of drugs, combine the combined application information, and further screen out a prediction model for the crystal structure of the two drug molecules based on the calculation results (the strength of the hydrogen bond between the two drug molecules, and the strength of the non-covalent bond between the two drug molecules in the binding conformation with the minimum free energy).

[0016] According to the present invention, the information on the combined use of drugs may be the information on the combined use of two drugs reported in the prior art, or may be the molecular prediction calculation result from the reverse docking of the drug with the target protein.

[0017] In one embodiment of the present invention, a step of obtaining the interaction information between the drug and the target protein and performing molecular prediction calculations for reverse docking of the drug and the target protein is optionally added to the method to predict the possibility of a combined effect between the two drugs.

[0018] According to the present invention, the drug information used for drug crystal structure prediction includes, but is not limited to: structural information of drug molecules (e.g., in mol or sdf format), drug target protein information, etc. The drug information can be obtained from drug databases, such as the Kyoto Encyclopedia of Genes and Genomes (KEGG) and the DrugBank database, or by experimental determination.

[0019] According to the present invention, when calculating the MEPS of the combination of two drug molecules and evaluating the strength of the hydrogen bond between the two drug molecules, the drug MEPS calculation formula is:

[0020]

[0021] ΔE=E cc +nE1-mE2 (2)

[0022] Among them, α i is the hydrogen bond donor parameter, β j is the hydrogen bond acceptor parameter, E represents the surface potential energy of the small molecule, and E cc is the crystal energy, E1 is the energy of drug 1, E2 is the energy of drug 2, and n and m are the number of drug 1 and drug 2 in the crystal, respectively. When calculating the MEPS of a combination of two drug molecules, it can be set according to the number of two drugs that usually appear in the co-crystal or co-amorphous compound, and feedback adjustment can be performed through the calculation results. Usually n and m are independently 1 or 2. In an embodiment of the present invention, the values ​​of n and m are 1, respectively.

[0023] Preferably, the drug structure is optimized based on Gaussian B3LYP parameter density functional and 6-31+G* basis set, and the local maximum and minimum values ​​are obtained according to the molecular electrostatic potential of the electron density isosurface. Obtain local maxima and minima in the molecular electrostatic potential of the electron density isosurface.

[0024] According to the calculated local maximum and minimum values, they are converted into α i and β j The hydrogen bonding parameters of the two drug molecules were calculated, and the E1 and E2 energy values ​​of the single drug and the ΔE energy value of the crystal were calculated. The |ΔE| calculated between the two drug molecules was used to score the hydrogen bonding strength of the combination of the two drug molecules (S MEPS ). Those skilled in the art will appreciate that the scoring method can be adjusted and varied according to the different structures of the drug. In one embodiment of the present invention, the scoring method is: when |ΔE|≥3 kJ / mol, S MEPS =5 points; 3kJ / mol>|ΔE|≥1kJ / mol, S MEPS =3 points; when |ΔE|<1kJ / mol, S MEPS =1 point.

[0025] According to the present invention, the combination simulation of two drug molecules is performed based on drug information, which is to simulate the trajectory of the intermolecular interaction between the two drug molecules. The calculation process includes:

[0026] (1) The initial structure of the system is established through experimental data or theoretical models, and the Amberff99SB force field is used for simulation at room temperature and pressure. The simulated complex system is placed in a cubic and periodic solvent box filled with SPC water molecule models, and the shortest distance between the complex composed of two drug small molecules and the solvent box is set to Counterions were added to the solvent to make the whole system electrically neutral. The system was minimized for 1000 steps using the steepest descent method to avoid possible atomic collisions or unreasonable spatial structures of the system. Then, the system was subjected to 100ps of equilibrium treatment under the NVT and NPT systems, respectively, so that the system temperature rose to 300K and the pressure was stabilized at one atmosphere. The temperature and pressure were fitted to the expected platform using the V-rescale and Parrinello-Rahman algorithms, respectively. At the specified temperature and pressure, molecular dynamics simulations were performed for 10 to 50 ns for each small molecule complex system, with a time step of 2 fs. During the simulation, the bond lengths of all atoms were constrained using the LINCS algorithm, and the long-range electrostatic forces were calculated using the Particle-Mesh Ewald (PME) method with a grid width of 0.12 nm.

[0027] (2) assigning each atom an initial velocity, which can be chosen from a Maxwell-Boltzmann distribution at a certain temperature;

[0028] (3) After the system reaches equilibrium, the particle motion trajectory that changes with time is obtained, and the obtained motion trajectory is sampled and selected to obtain the optimal combination of the two drug small molecules.

[0029] The trajectory regards microscopic particles as Newtonian classical particles, and uses the differential equation of Newton's second law to obtain the change of the position and velocity of atoms in the system over time. The calculation method is as follows:

[0030]

[0031] Formula (3) is the mass m i Atom i is under force F i Under the action of i The equation of motion is , and t is the calculation time.

[0032] The algorithm established in the simulation process is the leap-frog algorithm. This algorithm first uses the acceleration at time t to calculate the velocity after half a time step:

[0033]

[0034] Then we can calculate the position after the next step:

[0035]

[0036] The leapfrog method assumes that the system is The speed of time and the position at time t is known, and then by Calculate the force and acceleration received by the particle The prediction system The speed of time and the position at time t+δt Then calculate the speed at time t:

[0037]

[0038] According to the present invention, the intermolecular free energy calculation is performed based on the acquired trajectory conformation information of the simulated binding between the two drug molecules to screen the binding conformation with the minimum free energy between the two drug molecules. The calculation method is as follows:

[0039] G=E MM -TS MM +G PB +G surface (7)

[0040] Where E MM is the energy of the gas phase, T is the temperature, S MM is the entropy of the weather, G PB is the molecular polar energy, G surface is the nonpolar energy of the molecule.

[0041] First, E MM Partial calculation:

[0042] E MM =E bonded +E nonbonded＝ E bonded +E vdw +E elec (8)

[0043] Where E bonded Including chemical bonds, bond angles, and dihedral interactions; E nonbonded Including electrostatic interaction and van der Waals interaction, the computational potential functions commonly used in the art can be used for modeling. In one embodiment of the present invention, Coulomb and Lennard-Jones (LJ) potential functions are used for modeling respectively.

[0044] Then to G PB Partial calculation:

[0045] By solving the Poisson-Boltzmann equation for G PB Perform the calculation:

[0046]

[0047] where φ(r) is the electrostatic potential, ε(r) is the dielectric constant and ρ f(r) is a fixed charge density. The constant κ 2 is the reciprocal of the Debye length and depends on the ionic strength of the solution.

[0048] Calculate G again surface part:

[0049] G surface is the non-polar solvation energy, using the Solvent Accessible Surface Area (SASA) method:

[0050] G nonpolar =G surface =γA+b (10)

[0051] γ is the surface tension coefficient of the solvent, A is SASA, and b is the fitting parameter. When calculating, you need to specify the spacing of the cubic grid points (take The dielectric constant of pure water is ∈=78.45, and the dielectric constant of the solute is ∈=2. The empirical constant for calculating the nonpolar solvation energy is b=0.

[0052] In one embodiment of the present invention, the surface area A is calculated by the Shrake-Rupley numerical method.

[0053] When calculating the entropy contribution, we first need to perform energy minimization, where no cutoff is used for non-bonded interactions, then perform a normal analysis, calculate the mass-weighted Hessian matrix, diagonalize it, and use the obtained frequencies to calculate the entropy contribution:

[0054]

[0055] Where h is Planck's constant, νi is the ith frequency, and β = 1 / k B T.

[0056] According to the present invention, the conformation of the binding simulation trajectory of two drug molecules and the intermolecular MM-PBSA free energy calculation are integrated, and the binding conformation set of the two drug molecules can be formed by screening. Generally, through drug molecule binding simulation, each pair of drug combinations will have thousands of conformations of the simulation trajectory. The G value of each pair of drug combinations is the binding free energy of the two drug molecules. According to the G value, a pair of conformations with the smallest G value is selected from the simulation trajectory conformation set of each pair of drug combinations.

[0057] According to the present invention, the calculated E nonbonded The strength of the drug intermolecular force is scored by the selected conformation with the smallest G value (S nonbonded ). In one embodiment of the present invention, the scoring method is: |E nonbonded|≥10kcal / mol, S nonbonded =5 points; 10kcal / mol>|E nonbonded |≥7kcal / mol, S nonbonded =3 points; |E nonbonded |<7kcal / mol, S nonbonded =1 point.

[0058] According to the present invention, the strength of hydrogen bonding between drug molecules and the strength of non-covalent bonding between two drug molecules in the binding conformation with the minimum free energy can be further integrated to establish a structural prediction model for the formation of co-crystals or co-amorphous substances by two drug molecules. The scoring method of the prediction model is:

[0059] S total =S MEPS +S nonbonded That is, the strength of the non-covalent bonds and hydrogen bonds between two drug molecules is scored, and based on the score, whether a co-crystal or co-amorphous structure can be formed between the two drugs is predicted, and the possible co-crystal or co-amorphous structure is screened out.

[0060] According to the present invention, the information on the combined use of drugs refers to whether two drugs may be used in combination. The information can be obtained by using drug combinations that have been publicly reported or used in the art, such as a combination drug group that has been widely used in clinical practice. Alternatively, the information can be obtained by calculating the interaction between the drug and the target protein using a calculation method known in the art.

[0061] In one embodiment of the present invention, drug reverse molecular docking calculation is used to obtain information on drug combination application. The drug reverse molecular docking calculation process includes: (1) pre-treatment of the drug and the target protein by dehydration, hydrogenation, etc.; (2) using software such as AutoDock Vina software, docking the drug molecule optimized by pretreatment in step (1) to the active site of its target protein, thereby obtaining the interaction information between the drug and the target protein; (3) the binding site of the drug is described by a cubic box, the center of the box is located near the active center of the target protein, and the box covers the residues in the entire active area; (4) using the Monte Carlo method to search for conformations, and at the same time obtain energy scores for different molecular docking conformations; (5) evaluating the binding of the drug to the target protein based on the molecular docking energy score.

[0062] According to the present invention, only S total The value can be used to determine whether the two drugs can form a co-crystal or a co-amorphous structure. totalBased on the predicted value, the crystal structures of two drugs that are predicted to form co-crystals or co-amorphous structures and may have a combined effect are screened out in combination with whether the two drugs have the possibility of being used in combination. After obtaining the predicted structure, pharmaceutical and / or pharmacodynamic experiments can be carried out to obtain drug pairs that can eventually form drug co-crystals or co-amorphous structures.

[0063] The present invention provides information on whether a drug combination forms a drug crystal structure by analyzing the interaction characteristics between two drugs, and integrates molecular mechanics and molecular dynamics simulation methods to predict drug crystal structures faster and more effectively. At present, the most commonly used cocrystal screening calculation method in this field is MEPS. Although this method has a short calculation time, it can quickly screen coformers that can form cocrystals with the target API from a large database. However, MEPS calculation mainly looks at the most conventional hydrogen bonding between API and coformer (for example: amine-carboxylate, acid-carboxylate, etc.), which is not accurate for predicting whether cocrystals are formed between two APIs. In general, in addition to hydrogen bonding, other intermolecular forces between two APIs, such as halogen bonds, π-π stacking, van der Waals forces and electrostatic forces, are also very important. The present invention uses molecular dynamics simulation and MM-PBSA to calculate these forces, which not only has a fast calculation speed, but also improves the prediction accuracy. It only takes 3 to 5 hours / 40 core pairs to calculate a group of cocrystals using the method of the present invention, which greatly improves work efficiency.

[0064] Studies have shown that in clinical practice, some drugs (such as psychotropic drugs) can achieve better results when used in combination than when used alone, and new drug crystal forms that do not change the physical and chemical properties of the drug have advantages in terms of improved bioavailability and improved stability. Therefore, building a drug combination that can produce synergistic effects through drug crystal structure screening can improve the value of clinical applications. The present invention uses a reverse molecular docking method to reveal the interaction between drugs and target proteins, analyzes the drug action mechanism and drug combinations that can produce synergistic effects, and combines the predicted results of combined drug use with the predicted results of co-crystal structure, which can further improve the screening efficiency of drug co-crystals or co-amorphous compounds.

[0065] With the increasing application of drug repositioning and drug crystal screening in the drug development process, and the discovery of more and more drug co-crystals and co-amorphous compounds, a more practical foundation has been laid for the present invention.

[0066] The prediction method of the present invention can provide potential drug combinations that can form drug crystal structures and have synergistic effects, and is of great value in accelerating the development of new crystal form research in drug combination therapy programs. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1 is a flowchart of a specific implementation process of the prediction method in the embodiment;

[0068] Figure 2 It is the X-ray powder diffraction pattern of the physical mixture of clozapine and quetiapine, quetiapine, and clozapine;

[0069] Figure 3 is the X-ray powder diffraction pattern of amorphous compound of clozapine and quetiapine in a molar ratio of 1:1;

[0070] Figure 4 is the mDSC graph of the co-amorphous product of clozapine and quetiapine in a molar ratio of 1:1;

[0071] Figure 5 It is a comparison chart of infrared spectra of amorphous substance of clozapine and quetiapine in a molar ratio of 1:1, quetiapine and clozapine. DETAILED DESCRIPTION

[0072] In order to make the content of the present invention more clearly understood, the present invention is further described in detail below based on specific examples and drawings, taking psychotropic drugs as an example. The embodiments are not intended to limit the scope of protection of the present invention.

[0073] According to the method of the present invention, this example adopts the following steps:

[0074] 1. Obtain basic information about the drug

[0075] 18 known psychotropic drugs were obtained from the drug database (DrugBank), as shown in Table 1. The 18 drugs were sorted from 1 to 18 for subsequent crystal structure screening. Drug No. 15 is the ester structure of drug No. 14, so drug No. 14, which has a pharmacological effect, was selected for investigation and subsequent reverse molecular docking analysis.

[0076] Table 1 Structures and names of 18 psychotropic drugs

[0077]

[0078]

[0079] 2. Obtain information on drug-target protein interactions

[0080] The target protein information of the above 17 drugs was obtained from the Kyoto Encyclopedia of Genes and Genomes (KEGG) and DrugBank, as shown in Table 2. From the statistical results in Table 2, except for drug No. 8 clobazam which only acts on γ-aminobutyric acid (γ-GABA) receptors, the other drugs are multi-target drugs, mainly serotonin-2A subunit (5HT2A) and dopamine 2 subunit (DRD2), and secondly, dopamine 3 subunit (DRD3), serotonin-1A subunit (5HT1A), dopamine 1 subunit (DRD1), dopamine 4 subunit (DRD4), serotonin- 2C subunit (5HT2C), serotonin-3 subunit (5HT3), serotonin-6 subunit (5HT6) and serotonin-7 subunit (5HT7) receptors, of which 14 drugs act on 5HT2A, accounting for 82.35%, 16 drugs act on DRD2, accounting for 94.18%, 7 drugs act on DRD3 and 5 drugs act on 5HT1A, accounting for 50% and 35.71%, respectively.

[0081] Table 2 Information on 17 psychotropic drugs and their target proteins

[0082] (“+” in the table indicates that the drug interacts with the target protein)

[0083]

[0084] Similarly, according to the survey of KEGG database, it is found that most psychotropic drugs are agonists or antagonists of dopamine and 5-hydroxytryptamine. The first-generation psychotropic drugs chlorpromazine and haloperidol are typical dopamine antagonists. Chlorpromazine binds but does not activate dopamine receptors, thereby blocking dopamine or exogenous agonists. Haloperidol competitively blocks DRD2 receptors in the limbic system, thereby eliminating dopamine neurotransmission.

[0085] The second-generation psychotropics risperidone, clozapine, olanzapine, and quetiapine act as atypical dopamine and 5HT2A receptor antagonists. Studies have found that risperidone can improve positive and negative symptoms in the treatment of schizophrenia. Paliperidone is the active metabolite of risperidone. Risperidone selectively antagonizes serotonin through cortical 5HT2 receptors and competes with limbic DRD2 receptors to a lesser extent. Clozapine can weakly block DRD2 receptors and relieve schizophrenia symptoms. Quetiapine is an antagonist of histamine receptors and can be used to treat major depressive disorder in addition to schizophrenia. The new second-generation psychotropics iloperidone, asenapine, and lurasidone are popular in clinical practice, but their pharmacological characteristics are not clear. Ziprasidone is an atypical psychotropic that acts as an antagonist at DRD2 and 5HT2A receptors.

[0086] The third-generation psychotropic drugs aripiprazole, blonanserin, and cariprazine have a mechanism of action that is mainly related to dopamine. Aripiprazole, as a DRD2 receptor agonist, is an atypical schizophrenia drug with better tolerability than typical psychotropic drugs. Blonanserin is a dopamine agonist that affects the life cycle of serotonin and the survival of neurons. Cariprazine is a partial agonist of DRD2 and DRD3 receptors, and also acts as an antagonist of 5HT2A receptors.

[0087] Other drugs: Amisulpride is a dopamine antagonist. Low doses of amisulpride preferentially block presynaptic DRD2 and DRD3 receptors, enhancing dopamine transmission, while higher doses of amisulpride block postsynaptic receptors, inhibiting dopamine overactivity. Clobazam is a partial agonist of the γ-GABA receptor. Unlike endogenous γ-GABA ligands, clobazam binds to the γ-GABA receptor in an allosteric manner, increasing the frequency of chloride channel opening and the membrane permeability to chloride ions.

[0088] 3. Calculate the electrostatic potential surface energy between two drug molecules

[0089] Any two of the above 17 drugs were combined to obtain 153 drug pairs. The MEPS of the 153 drug combinations were calculated, and the strength of the hydrogen bond between the two drugs in each drug pair was evaluated.

[0090] The MEPS calculation formula for a drug is:

[0091]

[0092] ΔE=E cc +nE1-mE2 (2)

[0093] Among them, α i is the hydrogen bond donor parameter, β j is the hydrogen bond acceptor parameter, E represents the surface potential energy of the small molecule, and E cc is the crystal energy, E1 is the energy of drug 1, E2 is the energy of drug 2, n and m are the number of drug 1 and drug 2 in the crystal, respectively. To calculate the MEPS of the combination of two drugs, the values ​​of n and m are 1 respectively.

[0094] The drug structure was optimized based on Gaussian B3LYP parameter density functional and 6-31+G* basis set, and the Obtain local maxima and minima in the molecular electrostatic potential of the electron density isosurface.

[0095] According to the calculated local maximum and minimum values, they are converted into α i and β jThe hydrogen bonding parameters of the two drugs were calculated and the energy values ​​E of the individual drugs and the crystal were calculated. The calculated |ΔE| between the two drugs was used to score the strength of the hydrogen bonding effect of the combination of the two drugs (S MEPS ). The scoring method is: when |ΔE|≥3kJ / mol, S MEPS =5 points; 3kJ / mol>|ΔE|≥1kJ / mol, S MEPS =3 points; when |ΔE|<1kJ / mol, S MEPS =1 point.

[0096] 4. Perform a simulation calculation of the binding of two drug molecules

[0097] A binding model of two drug molecules is established to simulate the trajectory of intermolecular interaction between the two drug molecules. The calculation process includes:

[0098] For each of the 153 drug pairs, GROMACS was introduced to establish a two-drug molecule binding simulation model to simulate the process of intermolecular interaction between the two drug molecules. The Amberff99SB force field parameters were used to describe the drug binding system, and the drug binding was immersed in a cubic box filled with SPC water model, and the shortest distance between the complex and the solvent box was set to Counterions were added to the solvent to make the whole system electrically neutral. Before the MD simulation, the system was minimized for 1000 steps using the steepest descent method to avoid possible atomic collisions or unreasonable spatial structures. Then, the temperature of the simulation system was maintained at 300 K and the pressure of the system was maintained at 1.01×10 5 The coupling time was set to 100 ps. The simulated isothermal temperature of the controlled water was set to 4.5×10 5 Pa -1 To speed up the simulation, the bond lengths of all atoms were constrained using the LINCS method, and the PME method was used to calculate the long-range electrostatic interactions, with a grid width of 0.12 nm. Subsequently, 30 ns of MD simulation sampling was performed, and the time step in the simulation was set to 2 fs. The molecular structure trajectories of each of the 153 drug combinations were obtained.

[0099] 5. Perform intermolecular free energy calculation and select the binding conformation between two drug molecules:

[0100] According to the binding trajectory of the two drug molecules, the MM-PBSA free energy was calculated to obtain the structure of the two drug crystals that can be formed in 153 pairs. The calculation method is: G = E MM -TS MM +G PB+G surface

[0101] Where E MM is the energy of the gas phase, T is the temperature, S MM is the entropy of the weather, G PB is the molecular polar energy, G surface is the nonpolar energy of the molecule.

[0102] Calculate E MM The method is: E MM =E bonded +E nonbonded＝ E bonded +E vdw +E elec

[0103] Where E bonded Including chemical bonds, bond angles, dihedral interactions, E nonbonded Including electrostatic interactions and van der Waals interactions.

[0104] The conformations of the drug molecule binding simulation trajectory and the calculated MM-PBSA free energy are integrated to screen the conformation set that can form a co-crystal or co-amorphous combination of two drug molecules. Each pair of drug molecules has 1000 to 2000 simulation trajectory conformations. The G value of each pair of drug molecules is the binding free energy of the two drug molecules. According to the G value, a pair of conformations with the smallest G value is selected from the 1000 to 2000 simulation trajectory conformations of each pair of combinations. Then use the calculated E nonbonded The strength of the interaction between the two drug molecules is scored by using the selected conformation with the smallest G value (S nonbonded ). The scoring method is: |E nonbonded |≥10kcal / mol, S nonbonded =5 points; 10kcal / mol>|E nonbonded |≥7kcal / mol, S nonbonded =3 points; |E nonbonded |<7kcal / mol, S nonbonded =1 point.

[0105] The information obtained from MEPS and MM-PBSA was integrated to establish a crystal structure prediction model for the two drug molecules. The prediction model was scored as follows: total =S MEPS +S nonbonded . S total The evaluation criteria for the magnitude of the intermolecular forces of the combination of two drugs are shown in Table 3. Based on the prediction model, the calculated S total The combination of two drug molecules with scores ≥4 is relatively total The probability of forming a crystal structure by combining two drug molecules with a score <4 is higher, so Stotal Drug combinations with scores <4 were not included in the following reverse molecular docking analysis.

[0106] Table 3 S total Scoring table for crystal structures of two drug components with ≥4 points

[0107] Two-drug combination <![CDATA[S MEPS ]]> <![CDATA[S nonbonded ]]> <![CDATA[S total ]]> Clobazam-Clozapine 5 3 8 Clobazam-Olanzapine 5 3 8 Clobazam-quetiapine 5 3 8 Iloperidone-ziprasidone 3 5 8 Olanzapine-ziprasidone 3 5 8 Paliperidone-quetiapine 5 3 8 Asenapine-Olanzapine 1 5 6 Asenapine-paliperidone 3 3 6 Asenapine-paliperidone palmitate 3 3 6 Chlorpromazine-Clozapine 1 5 6 Chlorpromazine-haloperidol 3 3 6 Chlorpromazine-Iloperidone 1 5 6 Clozapine-quetiapine 1 5 6 Haloperidol-iloperidone 3 3 6 Haloperidol-lurasidone 3 3 6 Olanzapine-paliperidone palmitate 3 3 6 Olanzapine-paliperidone 3 3 6 Olanzapine-quetiapine 3 3 6 Olanzapine-risperidone 3 3 6 Clobazam-risperidone 3 3 6 Paliperidone-ziprasidone 1 3 4 Aripiprazole-paliperidone 1 3 4 Aripiprazole-risperidone 1 3 4

[0108] 6. Perform drug reverse molecular docking

[0109] In order to study the binding of drugs to the dopamine receptor family, 5-hydroxytryptamine receptor family and γ-aminobutyric acid receptor family, 23 target proteins were reversely docked with 17 drugs (drug No. 15 is the ester structure of drug No. 14, so only drug No. 14 was selected for reverse molecular docking analysis) to provide optional reference target protein information for the analysis of potential targets that drugs may bind to, or their mechanism of action and possible side effects.

[0110] According to the drug reverse molecular docking calculation process in the invention content part, the molecular docking results of 17 drugs and 23 target proteins are shown in Table 4. The binding affinity of the 17 drugs to 5HT2A is stronger than that to DRD2, reflecting the pharmacological action combination of high antagonism to 5HT2A and weak antagonism to DRD2 of psychotropic drugs; in terms of the binding to the 5-hydroxytryptamine receptor family, the 17 drugs have good binding affinity with 5HT1B, 5HT2B, and 5HT2C in addition to the known target 5HT2A, but 5HT1B and 5HT2B receptors are not expressed in the brain, so they are not considered as pharmacological action targets of psychotropic drugs, but this may be the cause of side effects, while 5HT2C receptors are expressed in the brain and have similar binding affinity to 5HT2A. Studies have shown that 5HT2C receptor antagonists can inhibit dopamine activity through the mesencephalon limbic and nigrostriatal pathways, which can alleviate the symptoms of mental illness and reduce its extrapyramidal side effects.

[0111] Among the 17 drugs that bind to dopamine receptors D1-5, DRD4 shows the best affinity, while the other four have similar affinity. In the pathways regulated by dopamine, DRD3 and DRD4 regulate the same pathway, which may lead to side effects such as endocrine disorders, weight gain and sexual dysfunction.

[0112] Table 4 Molecular docking scores of 17 drugs and 23 target proteins

[0113]

[0114]

[0115] Note: Scores < -5.0 are replaced with spaces. Molecular docking only considers the binding affinity between the target protein and the drug.

[0116] Through reverse molecular docking, the targets of many drugs are not limited to the currently known 5HT2A and DRD2, but can be expanded to the 5-hydroxytryptamine receptor family and the dopamine receptor family, providing more options for their pharmacological effects. By integrating the binding structure information of two drugs screened based on the prediction model and the drug combination information obtained based on drug reverse molecular docking, 10 pairs of two-drug molecular combinations were finally selected from the aforementioned 23 pairs of drugs, as shown in Table 5.

[0117] Table 5 10 crystal structure prediction diagrams of two drug molecules (energy unit: kcal / mol)

[0118]

[0119]

[0120]

[0121] 7. Conduct efficacy experiments based on the model screening results to evaluate the accuracy of the prediction model of the present invention

[0122] Based on the predicted 10 pairs of two-drug combinations, 8 pairs were selected for efficacy testing to evaluate the effectiveness of the above prediction model.

[0123] The purpose of this study was to give mice an intraperitoneal injection of 0.3 mg / kg MK-801 ((+)-MK-801 hydrogen maleate) to induce a high activity model in mice, screen the efficacy of different formulations of schizophrenia drugs, and compare the inhibition rate of the model with the inhibition rate of the two drugs used alone to observe whether the effect strength (inhibition rate) of the combined drug application is greater than the effect strength of the single drug.

[0124] Solvent preparation method: Take a certain amount of CMC-Na, add pure water and stir to dissolve, and prepare a pure water mixed solution with a concentration of 0.5% CMC-Na.

[0125] Experimental Design

[0126] animal

[0127] All animals were quarantined and adapted for 4 days, during which they were observed daily. The animal room had a 12h / 12h light-dark cycle, ≥8 times / h ventilation, temperature of 16-26°C, and relative humidity of 40-70%. The animals were randomly divided into a blank group, a model group (only treated with the solvent), and each test group according to their weight, with 8 mice in each group.

[0128] Test methods

[0129] Refer to the standard operating procedure of the pharmacology laboratory SOP-YL-JF-001-01. That is, about 1 hour after the test sample or solvent is given, 0.3 mg / kg of MK-801 is injected intraperitoneally, and the mice are placed in a self-activity box (a black polyethylene box with a size of 29 cm × 29 cm × 30 cm) for video recording. The video recording time is 60 minutes. After the video recording, video analysis is performed to evaluate the activity of the mice.

[0130] Drug delivery design:

[0131] Test 1

[0132] Only solvent was added to the blank group and the model group; 0.5 (mg / kg, ig) aripiprazole was added to the drug 1 group; 0.05 (mg / kg, ig) paliperidone was added to the drug 2 group; 0.5 (mg / kg, ig) aripiprazole and 0.05 (mg / kg, ig) paliperidone were added to the drug 1+2 group (combination application group).

[0133] Test 2

[0134] Only solvent was added to the blank group and the model group; 0.5 (mg / kg, ig) aripiprazole was added to the drug 1 group; 0.03 (mg / kg, ig) risperidone was added to the drug 2 group; 0.5 (mg / kg, ig) aripiprazole and 0.03 (mg / kg, ig) risperidone were added to the drug 1+2 group (combination application group).

[0135] Test 3

[0136] Only solvent was added to the blank group and model group; 3 (mg / kg, ig) clozapine was added to drug group 1; 20 (mg / kg, ig) quetiapine was added to drug group 2; 3 (mg / kg, ig) clozapine and 20 (mg / kg, ig) quetiapine were added to drug 1+2 group (combination application group).

[0137] Test 4

[0138] Only solvent was added to the blank group and the model group; 0.5 (mg / kg, ig) olanzapine was added to the drug 1 group; 0.05 (mg / kg, ig) paliperidone was added to the drug 2 group; 0.5 (mg / kg, ig) olanzapine and 0.05 (mg / kg, ig) paliperidone were added to the drug 1+2 group (combination group).

[0139] Test 5

[0140] Only solvent was added to the blank group and the model group; 0.075 (mg / kg, ig) asenapine was added to the drug 1 group; 0.5 (mg / kg, ig) olanzapine was added to the drug 2 group; 0.075 (mg / kg, ig) asenapine and 0.5 (mg / kg, ig) olanzapine were added to the drug 1+2 group (combination application group).

[0141] Test 6

[0142] Only solvent was added to the blank group and model group; 4 (mg / kg, ig) chlorpromazine was added to drug group 1; 3 (mg / kg, ig) clozapine was added to drug group 2; 4 (mg / kg, ig) chlorpromazine and 3 (mg / kg, ig) clozapine were added to drug 1+2 group (combination application group).

[0143] Test 7

[0144] Only the solvent was added to the blank group and the model group; 0.05 (mg / kg, ig) paliperidone was added to the drug 1 group; 2 (mg / kg, ig) ziprasidone was added to the drug 2 group; 0.05 (mg / kg, ig) paliperidone and 2 (mg / kg, ig) ziprasidone were added to the drug 1+2 group (combination application group).

[0145] Test 8

[0146] Only solvent was added to the blank group and model group; 0.1 (mg / kg, ig) iloperidone was added to drug group 1; 2 (mg / kg, ig) ziprasidone was added to drug group 2; 0.1 (mg / kg, ig) iloperidone and 2 (mg / kg, ig) ziprasidone were added to drug 1+2 group (combination application group).

[0147] Evaluation of the inhibition of high activity in mice after drug administration

[0148] The high activity inhibition rate of mice in the model group was calculated by analyzing the activity distance of mice after drug administration: drug group inhibition rate = 100-(activity distance of mice in the drug group / activity distance of mice in the model group)×100.

[0149] result

[0150] The test was conducted for two days. The test result of the model group on the first morning was 260885.92mm, and the test result on the second afternoon was 250057.16mm. The error range between the two was about 4% (when calculating the inhibition rate, we combined the results of the two days). The error of the results of paliperidone 0.05mg / kg on the first and second days was less than 1%. Therefore, this test system is reliable. For some other repeated groups, we only conducted one test. The detailed results are as follows:

[0151] Experiment 1: The inhibition rate of drug group 1 (aripiprazole) was 19.07%; the inhibition rate of drug group 2 (paliperidone) was 21.03%; the inhibition rate of drug group 1+2 (combination group) was 42.07%.

[0152] Experiment 2: The inhibition rate of drug group 1 (aripiprazole) was 19.07%; the inhibition rate of drug group 2 (risperidone) was 17.10%; the inhibition rate of drug group 1+2 (combination group) was 42.22%.

[0153] Experiment 3: The inhibition rate of drug group 1 (clozapine) was 28.77%; the inhibition rate of drug group 2 (quetiapine) was 43.77%; the inhibition rate of drug group 1+2 (combination group) was 84.39%.

[0154] Experiment 4: The inhibition rate of drug group 1 (olanzapine) was 5.97%; the inhibition rate of drug group 2 (paliperidone) was 21.03%; the inhibition rate of drug group 1+2 (combination group) was 41.01%.

[0155] Experiment 5: The inhibition rate of drug group 1 (asenapine) was 28.28%; the inhibition rate of drug group 2 (olanzapine) was 5.97%; the inhibition rate of drug group 1+2 (combination group) was 33.15%.

[0156] Experiment 6: The inhibition rate of drug 1 (chlorpromazine) group was 54.89%; the inhibition rate of drug 2 (clozapine) group was 28.77%; the inhibition rate of drug 1+2 group (combination group) was 76.26%.

[0157] Experiment 7: The inhibition rate of drug group 1 (paliperidone) was 21.03%; the inhibition rate of drug group 2 (ziprasidone) was 14.47%; the inhibition rate of drug group 1+2 (combination group) was 30.66%.

[0158] Experiment 8: The inhibition rate of drug group 1 (iloperidone) was 11.75%; the inhibition rate of drug group 2 (ziprasidone) was 14.47%; the inhibition rate of drug group 1+2 (combination group) was 41.96%.

[0159] In summary, it can be concluded that among the eight drug combinations predicted through screening, the inhibition rates of five groups (aripiprazole-paliperidone; aripiprazole-risperidone; clozapine-quetiapine; olanzapine-paliperidone; iloperidone-ziprasidone) after combined use were greater than the inhibition rates of the drugs alone.

[0160] 8. Evaluate the effectiveness of the prediction model of the present invention based on screening experiments

[0161] The effectiveness evaluation of the present invention conducted the above-mentioned 5 pairs of co-crystal screening experiments for combined use, among which clozapine-quetiapine co-amorphous and olanzapine-paliperidone co-crystals were screened out.

[0162] The following uses clozapine-quetiapine co-amorphous compound as an example to illustrate the preparation and confirmation process of cocrystals or co-amorphous compounds.

[0163] The two drugs are from Enhua Research Institute, clozapine batch number SC05-20170802, quetiapine (fumaric acid) batch number 0903-20170904.

[0164] According to the chemical molar ratio of clozapine and quetiapine of 1:1, grind for 2 minutes and mix on a shaker for 20 minutes to make them uniform. Tested by X-ray diffractometer, the results are shown in Figure 2 .

[0165] According to the chemical molar ratio of clozapine and quetiapine of 1:1, grind and mix them evenly, heat them on a hot plate at 185℃ for 0.5-1min to melt them, then take out the sample and place it on an aluminum plate for quenching to obtain a light yellow block product, which is ground into powder. Tested by X-ray diffractometer, the results are shown in Figure 3 .

[0166] The X-ray powder diffraction pattern of the screening experiment of the present invention was collected on a D8Advance X-ray powder diffractometer, and the test temperature was a normal ambient temperature. The parameters of the X-ray powder diffraction method of the present invention are as follows:

[0167] Target: Cu-Kα radiation

[0168] wavelength:

[0169] Tube voltage: 40kV

[0170] Tube current: 40mA

[0171] Step size: 0.02°

[0172] Scanning speed: 0.3° / s

[0173] The results showed that the structural spectrum of the co-amorphous product of clozapine and quetiapine was broad peak-shaped and had no sharp diffraction peaks.

[0174] The differential scanning calorimetry (DSC) diagram of the present invention was collected on TA Q2000. Figure 4 , the method parameters of differential scanning calorimetry (DSC) described in the present invention are as follows:

[0175] Range: -40-400℃

[0176] Heating rate: 10℃ / min

[0177] The glass transition temperature of the amorphous form of clozapine and quetiapine is 50.44°C. The Tg value is between the Tg values ​​of amorphous clozapine and quetiapine, and is close to the Tg value (49.67°C) calculated by the Gordon–Taylor empirical equation.

[0178] The infrared spectrum (IR) of the present invention was collected on a Bruker Tensor 27 Fourier transform infrared spectrometer. Figure 5 As shown, the method parameters of infrared spectroscopy (IR) described in the present invention are as follows:

[0179] Tableting conditions: Potassium bromide tableting

[0180] Infrared spectra of clozapine and quetiapine co-amorphous structures (cm -1 ): 3278.3±2cm -1 、3052.5±2cm -1 、2916.9±2cm -1 、2849.5±2cm -1 、1714.8±2cm -1 、1600.8±2cm -1 、1573.9±2cm -1 、1458.0±2cm -1 、1376.9±2cm -1 、1303.8±2cm -1 、1243.9±2cm -1 、1120.5±2cm -1 、1082.6±2cm -1 、1062.7±2cm -1 、1003.3±2cm -1 、981.9±2cm -1 、953.2±2cm -1 、894.7±2cm -1 、875.5±2cm -1 、814.6±2cm -1 、762.2±2cm -1 、742.1±2cm -1 、694.9±2cm -1 、668.4±2cm -1 、597.0±2cm -1 、568.9±2cm -1 、546.3±2cm -1 、522.8±2cm -1 and 462.4±2cm -1.

[0181] Clozapine and quetiapine have similar structures and basically the same groups, and their IR spectra are also basically similar. The IR spectra of the amorphous form of clozapine and quetiapine are between 3300 and 3500 cm -1 The peak shape is relatively blunt, which may be due to the formation of chemical bonds between the hydroxyl or secondary amine groups in the amorphous co-polymer, resulting in a change in the chemical environment; at 1715cm -1 The carbonyl peak appears at , indicating that the group connected to the carbonyl group is disconnected or the conjugation is reduced compared with a single API molecule.

[0182] According to screening experiments, the inventors of the present invention found that clozapine and quetiapine, two active ingredients of psychotropic drugs, can form a co-amorphous substance with a molar ratio of 1:1.

[0183] Taking psychotropic drugs as an example, the effectiveness of the method of the present invention is demonstrated, especially the crystal structure prediction method based on the information of drugs with combined effects. It provides a simple, effective and low-research-cost reliable model for rationally predicting whether two drugs with combined effects can form a crystal structure, which is of great value for future screening of other effective combination drug crystal structures.

Claims

1. A method for predicting whether two drug molecules can form a co-crystal or a co-amorphous structure, comprising the following steps: Obtain drug information for drug crystal structure prediction; Based on the drug information, the MEPS calculation of the combination of two drug molecules is performed to evaluate the strength of hydrogen bonds between drug molecules. The drug MEPS calculation formula is: E＝-∑ ij α i β j (1) ΔE = E cc + nE1 - mE2 (2) in, α i is the hydrogen bond donor parameter, β j is the hydrogen bond acceptor parameter, E represents the surface potential energy of the small molecule, and E cc is the crystal energy, E1 is the energy of drug 1, E2 is the energy of drug 2, n and m are the number of drug 1 and drug 2 in the crystal respectively; the local maximum and minimum values ​​calculated are converted into α i and β j The hydrogen bonding parameters of the two drug molecules are calculated, and the E1 and E2 energy values ​​of the single drug and the ΔE energy value of the crystal are calculated; the strength of the hydrogen bonding effect of the combination of the two drug molecules is scored based on |ΔE| MEPS ; According to the drug information, the binding simulation calculation of the two drug molecules is performed to obtain the trajectory conformation set of the simulated binding between the two drug molecules; According to the acquired trajectory conformation information of the simulated binding between the two drug molecules, the intermolecular free energy calculation is performed to screen the binding conformation with the minimum free energy between the two drug molecules. The intermolecular free energy calculation method is as follows: G=E MM -TS MM +G PB +G surface (7) Where E MM is the energy of the gas phase, T is the temperature, S MM is the entropy of the weather, G PB is the molecular polar energy, G surface is the nonpolar energy of the molecule; E MM The calculation of the part is: AND MM =And bonded +E nonbonded＝ AND bonded +E vdw +E elec (8) Where E bonded Including chemical bonds, bond angles, and dihedral interactions; E nonbonded Includes electrostatic interactions and van der Waals interactions; G PB The calculation of the part is: where φ(r) is the electrostatic potential, ε(r) is the dielectric constant and ρ f (r) is a fixed charge density; the constant κ 2 is the reciprocal of the Debye length and depends on the ionic strength of the solution; The solvent accessible surface area (SASA) method was used to calculate G surface : G nonpolar =G surface =γA+b (10) γ is the surface tension coefficient of the solvent, A is the SASA, and b is the fitting parameter; The entropy contribution is calculated by first performing energy minimization without using a cutoff for non-bonded interactions, then performing a normal analysis, calculating the mass-weighted Hessian matrix, diagonalizing it, and using the resulting frequencies to calculate the entropy contribution: Where h is Planck's constant, νi is the ith frequency, and β = 1 / k B T; Based on the strength of the hydrogen bond between the two drug molecules and the strength of the non-covalent bond between the two drug molecules in the binding conformation with the minimum free energy, a structural prediction model for the formation of a co-crystal or a co-amorphous substance between the two drug molecules is screened.

2. The method according to claim 1, characterized in that The MEPS of the combination of two drug molecules is calculated to evaluate the strength of the hydrogen bond between the two drug molecules. In the drug MEPS calculation formula, the drug structure is optimized based on Gaussian B3LYP parameter density functional and 6-31+G* basis set, and the local maximum and minimum values ​​are obtained according to the molecular electrostatic potential of the electron density isosurface.

3. The method according to claim 2, characterized in that From Obtain local maxima and minima in the molecular electrostatic potential of the electron density isosurface.

4. The method according to any one of claims 1 to 3, characterized in that: The calculation process of performing the binding simulation calculation of the two drug molecules to obtain the trajectory conformation set of the simulated binding between the two drug molecules includes: (1) Establish the initial structure of the system through experimental data or theoretical model, and use Amberff99SB force field to simulate at room temperature and pressure. Place the simulated complex system in a cubic and periodic solvent box filled with SPC water molecule model, and add counter ions to the solvent to make the whole system electrically neutral. Use the steepest descent method to minimize the system's energy for 1000 steps to avoid possible atomic collisions or unreasonable spatial structures of the system. Then, perform 100ps equilibrium treatment on the system under NVT and NPT systems respectively, so that the system temperature rises to 300K and the pressure is stabilized at one atmosphere. Use V-rescale and Parrinello-Rahman algorithms to fit the temperature and pressure to the expected platform respectively. Perform molecular dynamics simulation on each small molecule complex system at the specified temperature and pressure. (2) Selecting the initial velocity of each atom from the Maxwell-Boltzmann distribution at a certain temperature; (3) After the system reaches equilibrium, the particle motion trajectory that changes with time is obtained, and the obtained motion trajectory is sampled and selected to obtain the optimal combination of the two drug molecules.

5. The method according to claim 4, characterized in that In step (1) of the calculation process of performing a binding simulation calculation of two drug molecules to obtain a set of trajectory conformations of the simulated binding between the two drug molecules, the simulated complex system is placed in a cubic solvent box filled with SPC water molecule models and having a periodicity, and the shortest distance between the complex composed of the two drug small molecules and the solvent box is set to 6. The method according to claim 4, characterized in that In step (1) of the calculation process of performing a binding simulation calculation of two drug molecules to obtain a set of trajectory conformations of the simulated binding between the two drug molecules, a molecular dynamics simulation of 10 to 50 ns is performed on each small molecule complex system, and the time step is set to 2 fs.

7. The method according to claim 4, characterized in that In step (1) of the calculation process of performing a binding simulation calculation of two drug molecules to obtain a set of trajectory conformations for the simulated binding between the two drug molecules, during the molecular dynamics simulation, the bond lengths of all atoms were constrained using the LINCS algorithm, and the long-range electrostatic forces were calculated using the Particle-Mesh Ewald method with a grid width of 0.12 nm.

8. The method according to claim 4, characterized in that In step (3) of the calculation process of performing a binding simulation calculation of two drug molecules to obtain a set of trajectory conformations of the simulated binding between the two drug molecules, the position and velocity of atoms in the system are obtained over time through the differential equation of Newton's second law: Formula (3) is the mass m i Atom i is under force F i Under the action of i The motion equation at time t is the calculation time; Using the leap-frog algorithm, the acceleration at time t is used to calculate the velocity after half a time step: Then calculate the position after the next step: The calculation process assumes that the system is The speed of time and the position at time t is known, and then by Calculate the force and acceleration received by the particle The prediction system The speed of time and the position at time t+δt Then calculate the speed at time t:

9. The method according to any one of claims 1 to 3, characterized in that: For E nonbonded The electrostatic interaction and van der Waals interaction are modeled using Coulomb and Lennard-Jones (LJ) potential functions, respectively.

10. The method according to any one of claims 1 to 3, characterized in that: In the calculation of formula (10), the spacing of the cubic grid points is specified as The dielectric constant of pure water is ∈=78.45; the dielectric constant of the solute is ∈=2; b=0.

11. The method according to claim 10, characterized in that In the calculation of formula (10), A is calculated by the Shrake-Rupley numerical method.

12. The method according to any one of claims 1 to 3, characterized in that: The conformations of the binding simulation trajectory of the two drug molecules and the intermolecular MM-PBSA free energy are integrated to screen the binding conformation set that may form a co-crystal or co-amorphous product of the two drug molecules: according to the binding free energy G value of each pair of drug combinations, a pair of conformations with the smallest G value is selected from the conformation set of each pair of drug combinations, and then the calculated E nonbonded The value S is used to score the strength of the drug intermolecular force for the selected conformation with the smallest G value. nonbonded .

13. The method according to claim 12, characterized in that Based on the strength of the non-covalent bonds and hydrogen bonds between the two drug molecules, a score is given on whether the two drugs can form a co-crystal or co-amorphous structure: S total =S MEPS +S nonbonded .

14. The method according to any one of claims 1 to 3, characterized in that: A step is added to the method, which is to obtain information on the combined application of drugs, combine the combined application information, the strength of the hydrogen bond between the two drug molecules, and the strength of the non-covalent bond between the two drug molecules in the binding conformation with the minimum free energy, and further screen out a prediction model for the crystal structure of the two drug molecules.

15. The method according to claim 14, characterized in that Drug reverse molecular docking calculation is used to obtain information on drug combination application. The drug reverse molecular docking calculation process includes: (1) pre-processing the drug and target protein respectively; (2) using software to dock the drug molecule optimized by pre-processing in step (1) to the active site of its target protein, thereby obtaining the interaction information between the drug and the target protein; (3) the drug binding site is described by a cubic box, the center of the box is located near the active center of the target protein, and the box covers the residues in the entire active area; (4) using the Monte Carlo method to perform conformational search and obtain energy scores for different molecular docking conformations; (5) evaluating the binding of the drug to the target protein based on the molecular docking energy score.

Images