CN113129997A - MD/QM/CSM method for extracting beta-CD host-object system representative conformation - Google Patents
MD/QM/CSM method for extracting beta-CD host-object system representative conformation Download PDFInfo
- Publication number
- CN113129997A CN113129997A CN202110337800.2A CN202110337800A CN113129997A CN 113129997 A CN113129997 A CN 113129997A CN 202110337800 A CN202110337800 A CN 202110337800A CN 113129997 A CN113129997 A CN 113129997A
- Authority
- CN
- China
- Prior art keywords
- conformation
- host
- guest
- molecule
- dominant
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- WHGYBXFWUBPSRW-FOUAGVGXSA-N beta-cyclodextrin Chemical compound OC[C@H]([C@H]([C@@H]([C@H]1O)O)O[C@H]2O[C@@H]([C@@H](O[C@H]3O[C@H](CO)[C@H]([C@@H]([C@H]3O)O)O[C@H]3O[C@H](CO)[C@H]([C@@H]([C@H]3O)O)O[C@H]3O[C@H](CO)[C@H]([C@@H]([C@H]3O)O)O[C@H]3O[C@H](CO)[C@H]([C@@H]([C@H]3O)O)O3)[C@H](O)[C@H]2O)CO)O[C@@H]1O[C@H]1[C@H](O)[C@@H](O)[C@@H]3O[C@@H]1CO WHGYBXFWUBPSRW-FOUAGVGXSA-N 0.000 title claims abstract description 64
- 238000000034 method Methods 0.000 title claims abstract description 49
- 238000004364 calculation method Methods 0.000 claims abstract description 46
- 150000001875 compounds Chemical class 0.000 claims abstract description 25
- 239000002904 solvent Substances 0.000 claims abstract description 21
- 230000000694 effects Effects 0.000 claims abstract description 13
- 238000012216 screening Methods 0.000 claims abstract description 13
- 238000002474 experimental method Methods 0.000 claims abstract description 7
- 238000010219 correlation analysis Methods 0.000 claims abstract description 5
- 238000000329 molecular dynamics simulation Methods 0.000 claims description 44
- 238000005457 optimization Methods 0.000 claims description 39
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 25
- 230000003993 interaction Effects 0.000 claims description 17
- XEGGRYVFLWGFHI-UHFFFAOYSA-N bendiocarb Chemical compound CNC(=O)OC1=CC=CC2=C1OC(C)(C)O2 XEGGRYVFLWGFHI-UHFFFAOYSA-N 0.000 claims description 14
- DUEPRVBVGDRKAG-UHFFFAOYSA-N carbofuran Chemical group CNC(=O)OC1=CC=CC2=C1OC(C)(C)C2 DUEPRVBVGDRKAG-UHFFFAOYSA-N 0.000 claims description 14
- 230000000737 periodic effect Effects 0.000 claims description 14
- CVXBEEMKQHEXEN-UHFFFAOYSA-N carbaryl Chemical compound C1=CC=C2C(OC(=O)NC)=CC=CC2=C1 CVXBEEMKQHEXEN-UHFFFAOYSA-N 0.000 claims description 13
- 229960005286 carbaryl Drugs 0.000 claims description 12
- 238000007621 cluster analysis Methods 0.000 claims description 12
- 230000002159 abnormal effect Effects 0.000 claims description 8
- 230000007704 transition Effects 0.000 claims description 8
- 238000004458 analytical method Methods 0.000 claims description 7
- 239000007864 aqueous solution Substances 0.000 claims description 7
- 238000004422 calculation algorithm Methods 0.000 claims description 7
- 230000009881 electrostatic interaction Effects 0.000 claims description 7
- 239000000126 substance Substances 0.000 claims description 7
- 238000004057 DFT-B3LYP calculation Methods 0.000 claims description 5
- 238000003775 Density Functional Theory Methods 0.000 claims description 5
- 230000010354 integration Effects 0.000 claims description 5
- 238000010438 heat treatment Methods 0.000 claims description 4
- 230000000452 restraining effect Effects 0.000 claims description 4
- 239000000243 solution Substances 0.000 claims description 4
- 238000001514 detection method Methods 0.000 claims description 3
- 230000005274 electronic transitions Effects 0.000 claims description 2
- 230000005281 excited state Effects 0.000 claims description 2
- UFEJKYYYVXYMMS-UHFFFAOYSA-N methylcarbamic acid Chemical compound CNC(O)=O UFEJKYYYVXYMMS-UHFFFAOYSA-N 0.000 claims description 2
- 229920000858 Cyclodextrin Polymers 0.000 description 54
- 239000001116 FEMA 4028 Substances 0.000 description 53
- 235000011175 beta-cyclodextrine Nutrition 0.000 description 53
- 229960004853 betadex Drugs 0.000 description 53
- 238000002519 electronic circular dichroism spectroscopy Methods 0.000 description 15
- 230000008859 change Effects 0.000 description 4
- 239000000575 pesticide Substances 0.000 description 3
- 230000015572 biosynthetic process Effects 0.000 description 2
- 238000002983 circular dichroism Methods 0.000 description 2
- 239000002131 composite material Substances 0.000 description 2
- 230000003340 mental effect Effects 0.000 description 2
- 241000607479 Yersinia pestis Species 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000009471 action Effects 0.000 description 1
- 238000005284 basis set Methods 0.000 description 1
- 239000000073 carbamate insecticide Substances 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 230000000536 complexating effect Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000002212 electronic circular dichroism spectrum Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000036541 health Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 239000003446 ligand Substances 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 239000000447 pesticide residue Substances 0.000 description 1
- 102000004169 proteins and genes Human genes 0.000 description 1
- 108090000623 proteins and genes Proteins 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
- G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
- G16B15/00—ICT specially adapted for analysing two-dimensional or three-dimensional molecular structures, e.g. structural or functional relations or structure alignment
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/23—Clustering techniques
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Life Sciences & Earth Sciences (AREA)
- Theoretical Computer Science (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- General Engineering & Computer Science (AREA)
- Artificial Intelligence (AREA)
- General Physics & Mathematics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Chemical & Material Sciences (AREA)
- Biophysics (AREA)
- Crystallography & Structural Chemistry (AREA)
- Evolutionary Computation (AREA)
- Health & Medical Sciences (AREA)
- Biotechnology (AREA)
- General Health & Medical Sciences (AREA)
- Medical Informatics (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Organic Low-Molecular-Weight Compounds And Preparation Thereof (AREA)
Abstract
The invention discloses an MD/QM/CSM method for extracting a representative conformation of a beta-CD host-object system. The method comprises the following steps: (1) calculating the topological structure and RESP charge of the host molecule; (2) calculating the topological structure and RESP charge of the guest molecule; (3) constructing an initial conformation of the host-guest complex; (4) screening dominant conformations of the host-guest complex, the host molecule and the guest molecule; (5) performing ECD calculation on the dominant conformation of the host-guest compound, and comparing with an ICD experiment to determine the accuracy of a screening structure; (6) and calculating the inclusion free energy of the dominant conformations of the host-guest complex, the host molecule and the guest molecule by adopting a QM/CSM method, and performing linear correlation analysis on the calculated value and the experimental value. The invention can accurately extract the representative conformation of the host-guest complex system, provides an accurate structural model for the calculation of the inclusion free energy, and simultaneously calculates the inclusion free energy and the solvent effect in vacuum, so that the calculation of the inclusion free energy of the host-guest system is more perfect, and the invention provides guidance for experiments.
Description
Technical Field
The invention belongs to the field of theoretical calculation of a supramolecular host-object system, and relates to an MD/QM/CSM method for extracting a representative conformation of a beta-CD host-object system.
Background
N-methyl carbamate insecticides are widely used in agriculture to kill various pests that may damage crops and improve the quality of agricultural products. But poses a threat to human health during use, and therefore a simple, rapid and sensitive method is needed to quantify the pesticide residue. Because the pesticide molecules added into the organic medium have strong fluorescent signals, when the residues are measured, the content of the pesticide molecules can be rapidly measured by adding the pesticide molecules into a beta-cyclodextrin (beta-CD) medium.
The beta-CD and the derivatives thereof have the properties of internal hydrophobicity and external hydrophilicity, and can be included with different guest molecules, so that the solubility, the bioavailability and the like of the guest molecules are improved. The analysis of the structure and the complexing mechanism of the beta-CD host-guest complex is of great significance for deeply understanding the molecular recognition phenomenon in the system. With the development of computers and the continuous exploration of people on calculation methods, theoretical calculation plays an increasingly important role, reliable theoretical support is provided for experimental data, and a new visual angle is provided for designing novel supramolecular materials.
There are many theoretical calculation methods for predicting the inclusion capacity of a host and an object. In recent years, Gilson et al have organized a series of SAMPL (statistical Assessment of the Modeling of Proteins and ligands) and have conducted a comprehensive evaluation of the theoretical methods commonly used to calculate the inclusion free energy of host and guest. They have pointed out in the summary of SAMPL research results that most methods can calculate the large trend of host-guest inclusion free energy change, but they still have the problem of accuracy, the theoretical calculation value has low linear correlation with the experimental measurement value or has larger root mean square error, and meanwhile, the universality and portability of the calculation method still have larger problems (J Computaided Mol Des,2017,31(1): 1-19). Kessler et al calculated the inclusion free energies of 5 host molecules for the same guest molecule using DFT, PMF (positional of mean force) and MD-WHAM (weighted histogram analysis method with molecular dynamics relationships) methods. They found that the linear correlation between the calculation result of the DFT method and the experimental value is high, but the root mean square error is large; for the PM6-DH2X method, the root mean square error of the calculated and experimental values is less than that of the DFT method, but there is an error in the ranking of the relative stabilities of the different host-guest complexes; however, when the MD-WHAM is used for calculation, the root mean square error of the calculated value and the experimental value is the smallest, but the relative stability of different host-guest complexes cannot be correctly sequenced (Journal of Computational Chemistry,2012,33(29): 2310-2317).
Disclosure of Invention
The invention aims to provide an MD/QM/CSM method for extracting a representative conformation of a beta-CD host-object system.
The technical solution for realizing the purpose of the invention is as follows:
an MD/QM/CSM method for extracting a representative conformation of a beta-CD host-object system comprises the following steps:
(1) adopting a q4md-CD force field for the beta-CD of a main molecule to obtain a topological structure of the beta-CD and calculating RESP charge at the level of HF/6-31G;
(2) adopting a GAFF force field for N-methyl carbamate guest molecules to obtain a topological structure of the guest molecules and calculating at an HF/6-31G level to obtain RESP charges;
(3) adjusting the orientations of the host molecules subjected to charge and topology calculation in the step (1) and the guest molecules subjected to charge and topology calculation in the step (2) to construct an initial conformation of the host-guest complex;
(4) performing MD simulation on the initial conformation of the host-guest compound to obtain a compound conformation, adjusting a threshold value, performing cluster analysis and structure optimization, and screening out the dominant conformation of the compound; simultaneously, respectively carrying out MD simulation on the host molecules in the step (1) and the guest molecules in the step (2), adjusting a threshold value, carrying out cluster analysis and structure optimization, and respectively screening out dominant conformations of the host molecules and the guest molecules;
(5) performing Electronic Circular Dichroism (ECD) calculation on the dominant conformation of the compound obtained in the step (4), comparing with an Induced Circular Dichroism (ICD) experiment, determining the accuracy of the screened structure, performing the step (6) if the dominant conformation of the compound is consistent, and returning to the step 3 if the dominant conformation of the compound is inconsistent;
(6) calculating the dominant conformation of the compound subjected to ECD calculation and detection in the step (5) and the dominant conformation of the host molecule and the guest molecule screened in the step (4) by using a QM/CSM method to calculate inclusion free energy, and performing linear correlation analysis on the calculated inclusion free energy and an experimental value; if the calculated value has good linear correlation with the experimental value, the dominant conformation is the exact representative conformation of the host-guest complex, and if the calculated value has no good linear correlation with the experimental value, the step (3) is returned to.
In step (1) of the present invention, the structural formula of β -CD is:
in step (2) of the present invention, the N-methyl carbamate guest molecule is furacarb (PC), Bendiocarb (BC), Carbaryl (CY) or Carbofuran (CF), and the structural formula is as follows:
further, in the step (3), the specific method is as follows: and (3) opening the host molecule subjected to charge and topology calculation in the step (1) and the guest molecule subjected to charge and topology calculation in the step (2) on a gussview, and constructing the initial conformation of the host-guest complex by the ratio of beta-CD to the guest molecule of 1: 1.
Further, in step (4), the specific method for screening the dominant conformation of the host-guest complex is as follows:
a) dissolving the initial structure of the host-guest complex in TIP3P water, and setting the buffer distance from the host-guest complex to the boundary of the periodic aqueous solution box
Long range electrostatic interactions were treated using the PME method with the cut-off distance for non-covalent interactions set to
MD simulation is carried out under the periodic boundary condition, the SHAKE algorithm is used for restraining the chemical bond containing H, the restraint force is firstly exerted on the main body molecule, the first-step structural optimization is carried out on the system to eliminate the possible abnormal interaction between the solvent water molecule and the main body molecule and the compound thereof, then the first-step structural optimization is carried out under the unconstrained condition, and the temperature rise of the system is completed after the optimizationHeating the system to 300K through 400ps by adopting Langevin dynamics, performing MD simulation under 1atm and 300K by adopting NPT ensemble after the heating is finished, setting the integral step length as 2fs, storing energy and structure once every 2 ps, wherein the total MD time is 4ns, and 2000 conformations are obtained in the MD track;
b) after MD simulation, checking whether the object molecule enters a beta-CD cavity or not, if the object molecule does not enter the beta-CD cavity, readjusting the initial orientation of the beta-CD and the object molecule, and performing MD simulation again until the object molecule can enter the beta-CD cavity; if the object molecule can enter a beta-CD cavity, performing cluster analysis on conformations generated by MD simulation, and classifying the conformations into 6-8 classes by adjusting RMSD threshold;
c) extracting the structures of the first two largest clusters, carrying out unconstrained structure optimization on the MM level, after the optimization is completed, comprehensively considering the energy of the obtained conformation and the size of the cluster where the conformation is located, selecting an advantageous conformation, preferentially selecting a minimum value point from the largest clusters as the advantageous conformation if the difference between the conformation numbers in the two clusters is larger, and selecting a global minimum value point from the two conformations as the advantageous conformation if the conformation numbers in the two clusters are closer.
Further, in step (4), the specific method for screening the dominant conformation of the host molecule is as follows:
a) dissolving the initial structure of the host molecule in TIP3P water, and setting the buffer distance of the host molecule to the boundary of the periodic aqueous solution box
Long range electrostatic interactions were treated using the PME method with the cut-off distance for non-covalent interactions set to
MD simulation is carried out under the condition of periodic boundary, the SHAKE algorithm is used for restraining the chemical bond containing H, the restraint force is firstly exerted on the main molecules, the structure optimization of the system is carried out in the first step, so as to eliminate the possible abnormal interaction between the solvent water molecules and the main molecules, and then the abnormal interaction is carried out under the condition of no restraintPerforming structure optimization again, finishing the temperature rise of the system, adopting Langevin dynamics to heat the system to 300K through 400ps, adopting NPT ensemble after the temperature rise is finished, performing MD simulation under 1atm and 300K, setting the integration step length to 2fs, storing energy and structure once every 2 ps, setting the total MD time to 4ns, and obtaining 2000 conformations in the MD track;
b) performing clustering analysis on conformations generated by MD simulation, and dividing the conformations into 6-8 classes by adjusting RMSD threshold;
c) extracting the structures of the first two largest clusters, carrying out unconstrained structure optimization on the MM level, after the optimization is completed, comprehensively considering the energy of the obtained conformation and the size of the cluster where the conformation is located, selecting an advantageous conformation, preferentially selecting a minimum value point from the largest clusters as the advantageous conformation if the difference between the conformation numbers in the two clusters is larger, and selecting a global minimum value point from the two conformations as the advantageous conformation if the conformation numbers in the two clusters are closer.
Further, in the step (4), a specific method for screening out the dominant conformation of the guest molecule is as follows:
a) dissolving the initial structure of the host molecule in TIP3P water, and setting the buffer distance of the host molecule to the boundary of the periodic aqueous solution box
Long range electrostatic interactions were treated using the PME method with the cut-off distance for non-covalent interactions set to
MD simulation is carried out under the periodic boundary condition, a SHAKE algorithm is used for constraining chemical bonds containing H, constraint force is applied to main molecules, first-step structural optimization is carried out on a system to eliminate abnormal interaction possibly existing between solvent water molecules and the main molecules, then, structural optimization is carried out again under the condition of no constraint, the system is heated up after optimization is completed, Langevin dynamics is adopted to heat the system to 300K after 400ps, an NPT ensemble is adopted after the temperature rise is completed, MD simulation is carried out under 1atm and 300K, and integration step length is integratedSetting the energy and the structure as 2fs, storing the energy and the structure once every 2 ps, wherein the total MD time is 2ns, and 2000 conformations are obtained in the MD track;
b) performing clustering analysis on conformations generated by MD simulation, and dividing the conformations into 6-8 classes by adjusting RMSD threshold;
c) extracting the structures of the first two largest clusters, carrying out unconstrained structure optimization on the MM level, after the optimization is completed, comprehensively considering the energy of the obtained conformation and the size of the cluster where the conformation is located, selecting an advantageous conformation, preferentially selecting a minimum value point from the largest clusters as the advantageous conformation if the difference between the conformation numbers in the two clusters is larger, and selecting a global minimum value point from the two conformations as the advantageous conformation if the conformation numbers in the two clusters are closer.
Further, in the step (5), the ECD calculation method specifically includes:
removing water molecules in dominant conformations of the host-guest complex, performing TD-DFT calculation on a B3LYP/6-31G (d) group to obtain the wavelength, the vibrator intensity and the rotor intensity of an excited state, constructing an ECD spectrogram, comparing the ECD spectrogram with an ICD experimental spectrogram, simultaneously performing transition dipole moment calculation by using the B3LYP/6-31G (d) group to obtain the direction of electronic transition, comparing with the ICD experiment, and checking the accuracy of the dominant conformations of the host-guest complex.
Further, in step (5), the QM/CSM calculation method is specifically as follows:
a) removing water molecules in the dominant conformations of the host-guest complex, the host molecules and the guest molecules, and calculating the vacuum inclusion free energy by adopting a QM/CSM method;
b) optimizing the dominant conformation with water molecules removed in vacuum at the PM3 level and calculating the frequency, wherein after the frequency is calculated, if no virtual frequency appears, the screened dominant conformation is a stable representative conformation;
c) calculating the solvent effect by using a PCM method, and directly calculating the dominant conformation of the removed water molecules on a PM3 base group to obtain the solvent effect;
d) the vacuum inclusion free energy is added to the solvent effect to yield the total inclusion free energy in solution.
Compared with the prior art, the invention has the advantages that:
(1) according to the invention, ECD theoretical calculation and ICD experimental characterization in the circular dichroism technology are unified, so that a composite structure of a subject-object system can be accurately screened out from a complex MD conformational space, and accurate structural data is provided for subsequent energy calculation;
(2) the calculation model established by the improved method can calculate the inclusion free energy in vacuum, and also considers the solvent effect, so that the calculation of the inclusion free energy of the host-guest system is more complete.
Drawings
FIG. 1 is a flow chart of the MD/QM/CSM improvement method of the subject-object system of the present invention.
FIG. 2 is a graph showing the change of the centroid distance of the MD simulated subjects and objects of four subject and object systems of beta-CD/PC, beta-CD/BC, beta-CD/CY and beta-CD/CF along with time.
FIG. 3 is the ECD spectrum of dominant conformations of four host-guest complexes, beta-CD/PC, beta-CD/BC, beta-CD/CY and beta-CD/CF.
FIG. 4 is a diagram showing transition dipole moments of dominant conformations of four host-guest complexes, β -CD/PC, β -CD/BC, β -CD/CY and β -CD/CF.
FIG. 5 is a graph of analysis of the linear dependence of calculated inclusion free energy of four host-guest systems, namely, beta-CD/PC, beta-CD/BC, beta-CD/CY and beta-CD/CF, on experimental values.
Detailed Description
The invention is further illustrated with reference to the following specific embodiments and the accompanying drawings.
Example 1
(1) For beta-CD the topology was derived using a q4md-CD force field and RESP charge calculations were performed at HF/6-31G-level. The q4md-CD force field gives the topology of the cyclodextrin glucose units and the RESP charge calculated at HF/6-31G-level, with the GLYCAM04 force field and Amber99SB as force field parameters.
(2) The four N-methylcarbamate guest molecules PC, BC, CY, and CF were geometrically optimized using the B3LYP/6-31G (d) method and basal group, respectively, on Gussian 09. Then, using the programs from antechammber, tleap and parmchk2 of Amber software, topology was obtained in the GAFF force field and RESP charge calculations for guest molecules were performed at HF/6-31G levels.
(3) Opening the beta-CD and N-methyl carbamate guest molecules of the host molecules after charge and topology calculation on the guaassview, and constructing the initial structure of the beta-CD host-guest compound according to the proportion of 1:1 of the beta-CD and the guest molecules.
(4) Performing MD simulation on the initial conformation of the host-guest compound to obtain a compound conformation, adjusting a threshold value, performing cluster analysis and structure optimization, and screening out the dominant conformation of the compound; simultaneously, respectively carrying out MD simulation on the host molecules in the step (1) and the guest molecules in the step (2), adjusting a threshold value, carrying out cluster analysis and structure optimization, and respectively screening out dominant conformations of the host molecules and the guest molecules, specifically:
(4.1) screening of dominant conformations of the complex:
a) the initial structure of the beta-CD host-guest complex was MD simulated on Amber software. Dissolving the initial structure of the host-guest complex in TIP3P water, and setting the buffer distance from the host-guest complex to the boundary of the periodic aqueous solution box
Long range electrostatic interactions were treated using the PME method with the cut-off distance for non-covalent interactions set to
MD simulations were performed under periodic boundary conditions and using the shift algorithm to constrain H-containing chemical bonds. Firstly, applying constraint force to main molecules, carrying out first-step structural optimization on the system to eliminate abnormal interaction possibly existing between solvent water molecules and the main molecules and compounds thereof, and then carrying out structural optimization again under the unconstrained condition to finish heating the system. The system was warmed to 300K over 400ps using Langevin kinetics. Performing MD simulation under 1atm and 300K by adopting NPT ensemble after temperature rise is completed, setting the integration step length to 2fs, and storing energy once every 2 psAnd structure, total MD time 4ns, totaling 2000 conformations obtained in the MD trace.
b) After MD simulation, checking whether the object molecule enters a beta-CD cavity or not, if the object molecule does not enter the beta-CD cavity, readjusting the initial orientation of the beta-CD and the object molecule, and performing MD simulation again until the object molecule can enter the beta-CD cavity; and if the guest molecules can enter the beta-CD cavity, performing cluster analysis on the conformations generated by the MD simulation, and classifying the conformations into 6-8 classes by adjusting the RMSD threshold value.
c) And extracting the structures of the first two largest clusters, and performing unconstrained structure optimization on the MM level. After the optimization is completed, the energy of the conformation obtained by optimization and the size of cluster in which the conformation is positioned are comprehensively considered, and the dominant conformation is selected. If the conformational numbers in the two clusters are greatly different, preferentially selecting the minimum value point from the largest cluster as the dominant conformation; if the number of conformations in two clusters are relatively close, the global minimum point in the two classes of conformations is selected as the dominant conformation. Water molecules in the selected dominant conformation were removed for subsequent ECD and QM/CSM calculations.
(4.2) screening of dominant conformations of host and guest molecules:
in order to ensure the conformational quality of the host molecule and the guest molecule, the MD simulation is respectively carried out on the host molecule and the guest molecule, and the calculation parameters are basically consistent with the MD simulation of the compound. The guest molecule system is small, and the truncation distance of the non-covalent bond action is set as
The conformational space is also small, setting the time for MD simulation to 2 ns. Meanwhile, the conformations generated by the MD simulation of the host molecules and the guest molecules are subjected to cluster analysis, which is consistent with the cluster analysis method of the compound, and the dominant conformations of the host molecules and the guest molecules are screened out. Water molecules in the selected dominant conformation were removed for subsequent ECD and QM/CSM calculations.
Fig. 2 is a graph showing the variation of the centroid distance of the subject and object over time during the MD simulation. As can be seen from FIG. 2, the PC, CY and CF are placed at the large opening end of the beta-CD, and the BC is placed at the small opening end of the beta-CDWhen terminated, the resulting MD trace yields the best dominant conformation. Along the optimal MD trajectory, the four subject-object centroid distances vary with time as shown in fig. 2. Around t ≈ 200ps, the guest molecule enters the β -CD cavity and is then included in the cavity. After tending to balance, the physical and mental distances of the beta-CD/BC host and guest are all in
The principal and the guest have physical and mental distances of other three composite systems
Fluctuating up and down.
(5) The dominant conformations of the host-guest complexes extracted from the MD simulations were subjected to ECD and transition dipole moment calculations, as shown in fig. 3 and 4. The difference delta epsilon between the absorption coefficients of the left-handed circularly polarized light and the right-handed circularly polarized light of the beta-CD/BC is positive, and the delta epsilon of other three systems is negative, which is consistent with ICD experimental data. In the calculation results of the transition dipole moment, the direction of the transition dipole moment of the beta-CD/BC system is found to be parallel to the central axis of the beta-CD, and the directions of the transition dipole moments of the beta-CD/PC, beta-CD/CY and beta-CD/CF systems are all perpendicular to the central axis of the beta-CD, which is consistent with the calculation results and experimental data of the ECD, and the accuracy of the selected dominant conformation is proved.
And (3) if the ECD calculation and transition dipole moment calculation results of the extracted dominant conformations are inconsistent with ICD experimental data, returning to the step (3), readjusting the initial orientation between the subject and the object by using Guassview software, constructing a new MD initial conformation, extracting the dominant conformations through MD simulation and cluster analysis, and calculating the ECD.
(6) And (3) calculating the inclusion free energy of the beta-CD and N-methyl carbamate object molecules extracted from the MD simulation and the dominant conformation of the compound subjected to ECD calculation and detection in the step 5 by adopting a QM/CSM method. And (3) optimizing the extracted dominant conformation in vacuum at the PM3 level and calculating the frequency, wherein after the frequency is calculated, no virtual frequency appears, and the screened dominant conformation is a stable representative conformation. The solvent effect was calculated using the PCM method and the dominant conformation from which the MD simulation was extracted was calculated directly on the PM3 basis set to obtain the solvent effect. The vacuum inclusion free energy is added to the solvent effect to yield the total inclusion free energy in solution. And carrying out linear correlation analysis on the calculated inclusion free energy and the experimental value.
And (3) if the two systems have no good linear correlation or the relative stability sequence of each system is inconsistent with the experimental value, returning to the step (3), readjusting the initial orientation between the subject and the object by using the guassiviview software, constructing a new MD initial conformation, and repeating the steps until the inclusion free energy of the dominant conformation of the extracted complex has good linear correlation with the experimental value and the relative stability sequence of each system is consistent with the experimental value.
The calculated value of the dominant conformation inclusion free energy extracted by the steps has good linear correlation with the experimental value, the change trend of the calculated value is consistent with the experimental value, and the result is shown in table 1 (the last two columns are the calculated result and the experimental value of the inclusion free energy), so that the dominant conformation extracted by MD simulation is considered to be reasonable.
TABLE 1
FIG. 5 is a graph showing the linear correlation analysis between the calculated inclusion free energy and the experimental value of four subject-object systems of beta-CD/PC, beta-CD/BC, beta-CD/CY and beta-CD/CF, which are calculated by an improved MD/QM/CSM method, wherein a good linear correlation exists between the calculated inclusion free energy and the experimental value, and R is20.99. The vacuum inclusion free energy and solvent effect reflect the strength of the host and solvent-guest molecule interactions, respectively, and the data in Table 1 indicate that the former favors the formation of β -CD/PC and β -CD/CY complexes, while the latter favors the formation of β -CD/BC and β -CD/CF complexes. From the trend of change, the stability trend of other three systems except the beta-CD/CY system is mainly determined by the interaction between a host and an object, and the solvent effect also plays a certain role. Therefore, the results of calculating the four host-object systems by using the improved MD/QM/CSM method based on ECD calculation and ICD experimental data are reliable, and theoretical guidance can be performed on the experiments.
Claims (8)
1. The MD/QM/CSM method for extracting the representative conformation of the beta-CD host-object system is characterized by comprising the following steps:
(1) adopting a q4md-CD force field for the beta-CD of a main molecule to obtain a topological structure of the beta-CD and calculating RESP charge at the level of HF/6-31G;
(2) adopting a GAFF force field for N-methyl carbamate guest molecules to obtain a topological structure of the guest molecules and calculating at an HF/6-31G level to obtain RESP charges;
(3) adjusting the orientations of the host molecules subjected to charge and topology calculation in the step (1) and the guest molecules subjected to charge and topology calculation in the step (2) to construct an initial conformation of the host-guest complex;
(4) performing MD simulation on the initial conformation of the host-guest compound to obtain a compound conformation, adjusting a threshold value, performing cluster analysis and structure optimization, and screening out the dominant conformation of the compound; simultaneously, respectively carrying out MD simulation on the host molecules in the step (1) and the guest molecules in the step (2), adjusting a threshold value, carrying out cluster analysis and structure optimization, and respectively screening out dominant conformations of the host molecules and the guest molecules;
(5) performing ECD calculation on the dominant conformation of the compound obtained in the step (4), comparing with an ICD experiment, determining the accuracy of the screened structure, performing the step (6) if the dominant conformation of the compound is consistent, and returning to the step 3 if the dominant conformation of the compound is inconsistent;
(6) calculating the dominant conformation of the compound subjected to ECD calculation and detection in the step (5) and the dominant conformation of the host molecule and the guest molecule screened in the step (4) by using a QM/CSM method to calculate inclusion free energy, and performing linear correlation analysis on the calculated inclusion free energy and an experimental value; if the calculated value has good linear correlation with the experimental value, the dominant conformation is the exact representative conformation of the host-guest complex, and if the calculated value has no good linear correlation with the experimental value, the step (3) is returned to.
2. The MD/QM/CSM method of claim 1, wherein in step (2), the N-methylcarbamate guest molecule is furacarb, bendiocarb, carbaryl or carbofuran, having the following structural formula:
3. the MD/QM/CSM method according to claim 1, wherein in step (3), the specific method is as follows: and (3) opening the host molecule subjected to charge and topology calculation in the step (1) and the guest molecule subjected to charge and topology calculation in the step (2) on a gussview, and constructing the initial conformation of the host-guest complex by the ratio of beta-CD to the guest molecule of 1: 1.
4. The MD/QM/CSM method of claim 1, wherein the dominant conformation of the host-guest complex is selected in step (4) by the following steps:
a) dissolving the initial structure of the host-guest complex in TIP3P water, and setting the buffer distance from the host-guest complex to the boundary of the periodic aqueous solution box
Long range electrostatic interactions were treated using the PME method with the cut-off distance for non-covalent interactions set to
MD simulation is carried out under the periodic boundary condition, a SHAKE algorithm is used for constraining chemical bonds containing H, constraint force is applied to a main molecule, first-step structural optimization is carried out on a system to eliminate abnormal interaction possibly existing between solvent water molecules, the main molecule and compounds of the solvent water molecules, then, structural optimization is carried out again under the condition of no constraint, the system is heated up after optimization is completed, Langevin dynamics is adopted to heat the system to 300K through 400ps, NPT ensemble is adopted after heating up is completed, MD simulation is carried out under 1atm and 300K, the integration step length is set to 2fs, energy and structure are stored every 2 ps, and MD total time is setInterval is 4ns, and 2000 conformations are obtained in the MD track;
b) after MD simulation, checking whether the object molecule enters a beta-CD cavity or not, if the object molecule does not enter the beta-CD cavity, readjusting the initial orientation of the beta-CD and the object molecule, and performing MD simulation again until the object molecule can enter the beta-CD cavity; if the object molecule can enter a beta-CD cavity, performing cluster analysis on conformations generated by MD simulation, and classifying the conformations into 6-8 classes by adjusting RMSD threshold;
c) extracting the structures of the first two largest clusters, carrying out unconstrained structure optimization on the MM level, after the optimization is completed, comprehensively considering the energy of the obtained conformation and the size of the cluster where the conformation is located, selecting an advantageous conformation, preferentially selecting a minimum value point from the largest clusters as the advantageous conformation if the difference between the conformation numbers in the two clusters is larger, and selecting a global minimum value point from the two conformations as the advantageous conformation if the conformation numbers in the two clusters are closer.
5. The MD/QM/CSM method according to claim 1, wherein the dominant conformation of the host molecule is selected in step (4) by the following method:
a) dissolving the initial structure of the host molecule in TIP3P water, and setting the buffer distance of the host molecule to the boundary of the periodic aqueous solution box
Long range electrostatic interactions were treated using the PME method with the cut-off distance for non-covalent interactions set to
MD simulation is carried out under the periodic boundary condition, a SHAKE algorithm is used for restraining chemical bonds containing H, restraining force is firstly applied to main molecules, first-step structural optimization is carried out on a system to eliminate abnormal interaction possibly existing between solvent water molecules and the main molecules, then, structural optimization is carried out again under the condition of no restraint, the system is heated up after optimization is completed, and L is adoptedThe angevin dynamics heats the system to 300K through 400ps, an NPT ensemble is adopted after the temperature rise is finished, MD simulation is carried out under 1atm and 300K, the integral step length is set to be 2fs, energy and structure are stored every 2 ps, the total MD time is 4ns, and 2000 conformations are obtained in the MD track;
b) performing clustering analysis on conformations generated by MD simulation, and dividing the conformations into 6-8 classes by adjusting RMSD threshold;
c) extracting the structures of the first two largest clusters, carrying out unconstrained structure optimization on the MM level, after the optimization is completed, comprehensively considering the energy of the obtained conformation and the size of the cluster where the conformation is located, selecting an advantageous conformation, preferentially selecting a minimum value point from the largest clusters as the advantageous conformation if the difference between the conformation numbers in the two clusters is larger, and selecting a global minimum value point from the two conformations as the advantageous conformation if the conformation numbers in the two clusters are closer.
6. The MD/QM/CSM method according to claim 1, wherein the dominant conformation of the guest molecule is selected in step (4) by the following method:
a) dissolving the initial structure of the host molecule in TIP3P water, and setting the buffer distance of the host molecule to the boundary of the periodic aqueous solution box
Long range electrostatic interactions were treated using the PME method with the cut-off distance for non-covalent interactions set to
MD simulation is carried out under the periodic boundary condition, a SHAKE algorithm is used for constraining chemical bonds containing H, constraint force is firstly applied to main molecules, the first-step structural optimization is carried out on the system to eliminate abnormal interaction possibly existing between solvent water molecules and the main molecules, then, the structural optimization is carried out again under the condition of no constraint, the system is heated up after the optimization is completed, Langevin dynamics is adopted to heat the system to 300K after 400ps, NPT ensemble is adopted after the heating up is completed,performing MD simulation under 1atm and 300K, setting the integration step length as 2fs, storing energy and structure once every 2 ps, setting the total MD time as 2ns, and obtaining 2000 conformations in the MD track;
b) performing clustering analysis on conformations generated by MD simulation, and dividing the conformations into 6-8 classes by adjusting RMSD threshold;
c) extracting the structures of the first two largest clusters, carrying out unconstrained structure optimization on the MM level, after the optimization is completed, comprehensively considering the energy of the obtained conformation and the size of the cluster where the conformation is located, selecting an advantageous conformation, preferentially selecting a minimum value point from the largest clusters as the advantageous conformation if the difference between the conformation numbers in the two clusters is larger, and selecting a global minimum value point from the two conformations as the advantageous conformation if the conformation numbers in the two clusters are closer.
7. The MD/QM/CSM method according to claim 1, wherein in step (5), the ECD calculation method is as follows:
removing water molecules in dominant conformations of the host-guest complex, performing TD-DFT calculation on a B3LYP/6-31G (d) group to obtain the wavelength, the vibrator intensity and the rotor intensity of an excited state, constructing an ECD spectrogram, comparing the ECD spectrogram with an ICD experimental spectrogram, simultaneously performing transition dipole moment calculation by using the B3LYP/6-31G (d) group to obtain the direction of electronic transition, comparing with the ICD experiment, and checking the accuracy of the dominant conformations of the host-guest complex.
8. The MD/QM/CSM method according to claim 1, wherein in step (5), the QM/CSM calculation method is specifically as follows:
a) removing water molecules in the dominant conformations of the host-guest complex, the host molecules and the guest molecules, and calculating the vacuum inclusion free energy by adopting a QM/CSM method;
b) optimizing the dominant conformation with water molecules removed in vacuum at the PM3 level and calculating the frequency, wherein after the frequency is calculated, if no virtual frequency appears, the screened dominant conformation is a stable representative conformation;
c) calculating the solvent effect by using a PCM method, and directly calculating the dominant conformation of the removed water molecules on a PM3 base group to obtain the solvent effect;
d) the vacuum inclusion free energy is added to the solvent effect to yield the total inclusion free energy in solution.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110337800.2A CN113129997B (en) | 2021-03-30 | 2021-03-30 | MD/QM/CSM method for extracting beta-CD host-object system representative conformation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110337800.2A CN113129997B (en) | 2021-03-30 | 2021-03-30 | MD/QM/CSM method for extracting beta-CD host-object system representative conformation |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113129997A true CN113129997A (en) | 2021-07-16 |
| CN113129997B CN113129997B (en) | 2022-09-13 |
Family
ID=76774932
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110337800.2A Expired - Fee Related CN113129997B (en) | 2021-03-30 | 2021-03-30 | MD/QM/CSM method for extracting beta-CD host-object system representative conformation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113129997B (en) |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN106407740A (en) * | 2016-09-05 | 2017-02-15 | 南京大学 | Method for screening anti-androgen activity of flavonoid compounds based on molecular dynamics simulation |
| US20170285007A1 (en) * | 2016-04-01 | 2017-10-05 | Nanjing University | Screening methods for thyroid hormone disruptors based on co-regulator involved simulations |
-
2021
- 2021-03-30 CN CN202110337800.2A patent/CN113129997B/en not_active Expired - Fee Related
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20170285007A1 (en) * | 2016-04-01 | 2017-10-05 | Nanjing University | Screening methods for thyroid hormone disruptors based on co-regulator involved simulations |
| CN106407740A (en) * | 2016-09-05 | 2017-02-15 | 南京大学 | Method for screening anti-androgen activity of flavonoid compounds based on molecular dynamics simulation |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113129997B (en) | 2022-09-13 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Vassetti et al. | Assessment of GAFF2 and OPLS-AA general force fields in combination with the water models TIP3P, SPCE, and OPC3 for the solvation free energy of druglike organic molecules | |
| JP6942741B2 (en) | Identification for spectroscopic classification with reduced false positives | |
| Wang et al. | Predicting the impacts of mutations on protein-ligand binding affinity based on molecular dynamics simulations and machine learning methods | |
| CN110850020A (en) | Traditional Chinese medicine identification method based on artificial intelligence | |
| CN102279213A (en) | Method for rapid diagnosis of crop disease by volatile matter | |
| US12002551B2 (en) | Method to distinguish peroxisome proliferator-activated receptor gamma full agonist, partial agonist and antagonist with different activities and identification thereof | |
| US20210142038A1 (en) | Outlier detection for spectroscopic classification | |
| CN104792652A (en) | Multi-index rapid detection method for radix astragali | |
| CN103258244A (en) | Method for predicting inhibiting concentration of pyridazine HCV NS5B polymerase inhibitor based on particle swarm optimization support vector machine | |
| CN105372202B (en) | Transgene cotton variety ecotype method | |
| CN113624694A (en) | Inversion method and device for atmospheric methane concentration | |
| CN113129997B (en) | MD/QM/CSM method for extracting beta-CD host-object system representative conformation | |
| Del Galdo et al. | Cpl spectra of camphor derivatives in solution by an integrated QM/MD approach | |
| CN112837743A (en) | Medicine repositioning method based on machine learning | |
| CN107255646A (en) | A kind of method of fast quantification Predicting Stability of Drugs | |
| Zhang et al. | Machine learning-assisted quantitative prediction of thermal decomposition temperatures of energetic materials and their thermal stability analysis | |
| CN111781163B (en) | Method for eliminating influence of soil granularity on soil parameter detection of discrete near-infrared band | |
| Hu et al. | Interpretable prediction of protein-ligand interaction by convolutional neural network | |
| Liu et al. | ATR‐FTIR Spectroscopy Preprocessing Technique Selection for Identification of Geographical Origins of Gastrodia elata Blume | |
| CN107077532B (en) | Similarity evaluation method for determining structural effect of solvent reactivity, and system using the same | |
| Kitchen et al. | Computational techniques to support hit triage | |
| Jacek et al. | The log P parameter as a molecular descriptor in the computer-aided drug design–an overview | |
| CN113744815A (en) | MD/QM/CSM verification method for self-assembled supramolecular material | |
| Villar et al. | Substructural analysis in drug discovery | |
| NL2029908B1 (en) | Method for evaluating cold-hot nature of chinese herbal medicine based on fingerprint |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20220913 |
|
| CF01 | Termination of patent right due to non-payment of annual fee |