CN113129997B - 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
- CN113129997B CN113129997B CN202110337800.2A CN202110337800A CN113129997B CN 113129997 B CN113129997 B CN 113129997B CN 202110337800 A CN202110337800 A CN 202110337800A CN 113129997 B CN113129997 B CN 113129997B
- 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.)
- Expired - Fee Related
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 63
- 238000000034 method Methods 0.000 title claims abstract description 48
- 238000004364 calculation method Methods 0.000 claims abstract description 47
- 150000001875 compounds Chemical class 0.000 claims abstract description 27
- 239000002904 solvent Substances 0.000 claims abstract description 21
- 238000012216 screening Methods 0.000 claims abstract description 16
- 230000000694 effects Effects 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 37
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 26
- 230000003993 interaction Effects 0.000 claims description 15
- 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
- 238000007621 cluster analysis Methods 0.000 claims description 13
- 229960005286 carbaryl Drugs 0.000 claims description 12
- 238000010438 heat treatment Methods 0.000 claims description 11
- 230000007704 transition Effects 0.000 claims description 8
- 230000002159 abnormal effect Effects 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
- 238000004458 analytical method Methods 0.000 claims description 6
- 239000007864 aqueous solution Substances 0.000 claims description 6
- 238000004057 DFT-B3LYP calculation Methods 0.000 claims description 5
- 238000003775 Density Functional Theory Methods 0.000 claims description 5
- 239000000243 solution Substances 0.000 claims description 5
- 238000001514 detection method Methods 0.000 claims description 3
- 230000009471 action Effects 0.000 claims description 2
- 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 53
- 239000001116 FEMA 4028 Substances 0.000 description 52
- 235000011175 beta-cyclodextrine Nutrition 0.000 description 52
- 229960004853 betadex Drugs 0.000 description 52
- 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
- 238000010586 diagram Methods 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
- 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
- 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
- 230000010354 integration 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 fluorescence signals, when the residues are measured, the pesticide molecules are added into a beta-cyclodextrin (beta-CD) medium, so that the content of the pesticide molecules can be measured quickly.
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 Assembly of the Modeling of Proteins and ligands) studies, and have conducted a comprehensive evaluation of the theoretical methods commonly used to calculate the inclusion free energy of hosts and guests. 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-guest 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 inclusion free energy by using 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) together by adopting a QM/CSM method, 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
Performing MD simulation under a periodic boundary condition, constraining chemical bonds containing H by using a SHAKE algorithm, applying constraint force to a main molecule, performing first-step structural optimization on a system to eliminate abnormal interaction possibly existing between solvent water molecules, the main molecule and a compound of the solvent water molecules and the main molecule, performing primary structural optimization under an unconstrained condition, completing temperature rise of the system, heating the system to 300K through 400ps by adopting Langevin dynamics, performing MD simulation under 1atm and 300K by adopting NPT ensemble after temperature rise is completed, setting an integral step length to be 2fs, storing energy and structure every 2 ps, keeping the total MD time to be 4ns, and obtaining 2000 conformations in an 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
Performing MD simulation under a periodic boundary condition, constraining chemical bonds containing H by using an SHAKE algorithm, applying a constraint force to a main molecule, performing first-step structural optimization on a system to eliminate abnormal interaction possibly existing between solvent water molecules and the main molecule, performing primary structural optimization under an unconstrained condition, finishing optimization to heat the system, heating the system to 300K through 400ps by adopting Langevin dynamics, performing MD simulation under 1atm and 300K by adopting an NPT (non-positive-negative-pressure) ensemble after finishing heating, setting an integral step length to 2fs, storing energy and structure once every 2 ps, and obtaining 2000 conformations in a total MD time of 4MD and ns track;
b) performing cluster 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 water solution box as
Long range electrostatic interactions were treated using the PME method with a cutoff distance for non-covalent interactions set to
Performing MD simulation under a periodic boundary condition, constraining chemical bonds containing H by using an SHAKE algorithm, applying a constraint force to a main molecule, performing first-step structural optimization on a system to eliminate abnormal interaction possibly existing between solvent water molecules and the main molecule, performing primary structural optimization under an unconstrained condition, finishing optimization to heat the system, heating the system to 300K through 400ps by adopting Langevin dynamics, performing MD simulation under 1atm and 300K by adopting an NPT (non-positive-negative-pressure) ensemble after finishing heating, setting an integral step length to 2fs, storing energy and structure once every 2 ps, and obtaining 2000 conformations in a MD total time of 2 MD and ns 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 in which the conformation is located, selecting a dominant conformation, preferentially selecting a minimum point from the largest clusters as the dominant conformation if the conformational numbers in the two clusters are greatly different, and selecting a global minimum point from the two conformations as the dominant conformation if the conformational numbers in the two clusters are relatively close.
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) base group to obtain the wavelength, the oscillator intensity and the rotor intensity of an excited state, constructing an ECD spectrogram, comparing the ECD spectrogram with an ICD experimental spectrogram, performing transition dipole moment calculation by using the B3LYP/6-31G (d) base 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 the step (5), the QM/CSM calculating method specifically includes:
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, and if no virtual frequency appears after the frequency calculation, indicating that 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 conformation 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 diagram of the change of the centroid distance of a simulated subject and an object in four subject and object systems MD with time, wherein the distances of the centroid of the simulated subject and the object are beta-CD/PC, beta-CD/BC, beta-CD/CY and beta-CD/CF.
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 TThe buffer distance from the host-guest complex to the boundary of the periodic aqueous solution box in IP3P water is set as
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 constraint-free condition to finish heating the system. The system was warmed to 300K over 400ps using Langevin kinetics. And performing MD simulation under 1atm and 300K by adopting an NPT ensemble after the temperature rise is finished, setting the integration step length to be 2fs, storing energy and structure once every 2 ps, wherein the total MD time is 4ns, and obtaining 2000 conformations in an 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; and if the object molecule can enter a beta-CD cavity, performing cluster analysis on the conformation generated by MD simulation, and classifying the conformation into 6-8 classes by adjusting RMSD threshold.
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 MD traces obtained when placing PC, CY, CF at the large mouth end of the β -CD and BC at the small mouth end of the β -CD produce 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 (4) 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 guasview software, constructing a new MD initial conformation, extracting the dominant conformations through MD simulation and cluster analysis, and performing ECD calculation.
(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 is 2 0.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 (5)
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 to the main molecule beta-CD to obtain the topological structure of the beta-CD and calculating the 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;
the specific method for screening out the dominant conformation of the host-guest complex is as follows:
a) dissolving an initial structure of a host-guest compound in TIP3P water, setting a buffer distance from the host-guest compound to a box boundary of a periodic aqueous solution to be 8A, processing long-range electrostatic interaction by using a PME method, setting a truncation distance of non-covalent bond interaction to be 10A, performing MD simulation under a periodic boundary condition, constraining chemical bonds containing H by using an SHAKE algorithm, applying constraint force to host molecules, performing first-step structural optimization on a system to eliminate abnormal interaction possibly existing between solvent water molecules and the host molecules and the compound thereof, performing primary structural optimization under an unconstrained condition, finishing temperature rise of the system, heating the system to 300K through 400ps by using Langevin dynamics, performing NPT ensemble after temperature rise is finished, performing MD simulation under 1atm and 300K, and setting an integral step length to be 2fs, the energy and structure are saved once every 2 ps, 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 the conformation generated by MD simulation, and dividing the conformation 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;
specific methods for screening the dominant conformation of the host molecule are as follows:
a) dissolving an initial structure of a main body molecule in TIP3P water, setting a buffer distance from the main body molecule to a box boundary of a periodic aqueous solution to be 8A, processing long-range electrostatic interaction by using a PME method, setting a truncation distance of non-covalent bond action to be 10A, performing MD simulation under a periodic boundary condition, constraining chemical bonds containing H by using an SHAKE algorithm, applying constraint force to the main body molecule, performing first-step structural optimization on a system to eliminate abnormal interaction possibly existing between solvent water molecules and the main body molecule, performing structural optimization again under an unconstrained condition, finishing temperature rise of the system by optimizing, heating the system to 300K through 400ps by adopting Langevin dynamics, heating up and then adopting an NPT system, performing MD simulation under 1atm and 300K, setting an integral step length to be 2fs, storing energy and structure once every 2 ps, total time of MD is 4ns, and a total of 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;
the specific method for screening out the dominant conformation of the guest molecule is as follows:
a) dissolving an initial structure of a main body molecule in TIP3P water, setting a buffer distance from the main body molecule to a box boundary of a periodic aqueous solution to be 8A, processing long-range electrostatic interaction by using a PME method, setting a truncation distance of non-covalent bond interaction to be 6A, performing MD simulation under a periodic boundary condition, constraining chemical bonds containing H by using an SHAKE algorithm, applying constraint force to the main body molecule, performing first-step structural optimization on a system to eliminate abnormal interaction possibly existing between solvent water molecules and the main body molecule, performing structural optimization again under an unconstrained condition, finishing temperature rise of the system by optimizing, heating the system to 300K through 400ps by adopting Langevin dynamics, heating up and then adopting an NPT system, performing MD simulation under 1atm and 300K, setting an integral step length to be 2fs, storing energy and structure once every 2 ps, total time of MD is 2ns, and a total of 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 in which the conformation is located, selecting a dominant conformation, preferentially selecting a minimum point from the largest clusters as the dominant conformation if the conformational numbers in the two clusters are greatly different, and selecting a global minimum point from the two classes of conformations as the dominant conformation if the conformational numbers in the two clusters are relatively close;
(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 inclusion free energy by using 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) together by adopting a QM/CSM method, 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 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.
5. 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 a solvent effect by using a PCM method, and directly calculating the dominant conformation of removed water molecules on a PM3 basis 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 CN113129997A (en) | 2021-07-16 |
| CN113129997B true 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) |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105893759B (en) * | 2016-04-01 | 2018-08-24 | 南京大学 | A kind of thyroid hormone replacement therapy virtual screening and its active quantitative calculation method of interference being total to regulatory factor based on nuclear receptor |
| CN106407740A (en) * | 2016-09-05 | 2017-02-15 | 南京大学 | Method for screening anti-androgen activity of flavonoid compounds based on molecular dynamics simulation |
-
2021
- 2021-03-30 CN CN202110337800.2A patent/CN113129997B/en not_active Expired - Fee Related
Also Published As
| Publication number | Publication date |
|---|---|
| CN113129997A (en) | 2021-07-16 |
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 | |
| US20230273122A1 (en) | Reduced false positive identification for spectroscopic classification | |
| Roy et al. | On various metrics used for validation of predictive QSAR models with applications in virtual screening and focused library design | |
| CN102222178B (en) | Method for screening and/or designing medicines aiming at multiple targets | |
| Xiao et al. | GPCR-2L: predicting G protein-coupled receptors and their types by hybridizing two different modes of pseudo amino acid compositions | |
| Duffy et al. | Early phase drug discovery: cheminformatics and computational techniques in identifying lead series | |
| CN102279213B (en) | Method for rapid diagnosis of crop disease by volatile matter | |
| CN103258244A (en) | Method for predicting inhibiting concentration of pyridazine HCV NS5B polymerase inhibitor based on particle swarm optimization support vector machine | |
| CN104792652A (en) | Multi-index rapid detection method for radix astragali | |
| CN113129997B (en) | MD/QM/CSM method for extracting beta-CD host-object system representative conformation | |
| CN107526939B (en) | Rapid alignment method for small molecular structure | |
| Del Galdo et al. | Cpl spectra of camphor derivatives in solution by an integrated QM/MD approach | |
| CN103353445A (en) | Technical method for quickly identifying drought resistance of wheat by using near-infrared spectroscopy | |
| Li et al. | Prediction of kinase–substrate relations based on heterogeneous networks | |
| CA2477459C (en) | Comparative field analysis (comfa) utilizing topomeric alignment of molecular fragments | |
| Liu et al. | ATR‐FTIR Spectroscopy Preprocessing Technique Selection for Identification of Geographical Origins of Gastrodia elata Blume | |
| Kitchen et al. | Computational techniques to support hit triage | |
| JP6395276B2 (en) | Similarity evaluation method for structural effects determining solvent reactivity and system using the same | |
| CN113744815A (en) | MD/QM/CSM verification method for self-assembled supramolecular material | |
| Villar et al. | Substructural analysis in drug discovery | |
| CN103424160B (en) | A kind of method measuring rice field water holding layer water level depth | |
| CN118351928A (en) | Improved MD/QM/CSM method for screening representative conformations of beta-CD or host-guest system derived therefrom | |
| WO2020147668A1 (en) | Absolute binding free energy perturbation method for predicting drug target binding strength | |
| CN110929804B (en) | Method, device, equipment and medium for identifying production area of cultivated product | |
| Gangwal et al. | Overview and recent advances in Qsar studies |
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 |