Structures in solutions from joint experimental-computational analysis: applications to cyclic molecules and studies of noncovalent interactions

The potential of an approach combining nuclear magnetic resonance (NMR) spectroscopy, molecular dynamics (MD) simulations, and quantum mechanical (QM) calculations for full structural characterizations in solution is assessed using cyclic organic compounds, namely, benzazocinone derivatives 1-3 with fused five- and eight-membered aliphatic rings, camphoric anhydride 4, and bullvalene 5. Various MD simulations were considered, using force field and semiempirical QM treatments, implicit and explicit solvation, and high-temperature MD calculations for selecting plausible molecular geometries for subsequent QM geometry optimizations using mainly B3LYP, M062X, and MP2 methods. The QM-predicted values of NMR parameters were compared to their experimental values for verification of the final structures derived from the MD/QM analysis. From these comparisons, initial estimates of quality thresholds (calculated as rms deviations) were 0.7-0.9 Hz for (3)J(HH) couplings, 0.07-0.11 ? for interproton distances, 0.05-0.08 ppm for (1)H chemical shifts, and 1.0-2.1 ppm for (13)C chemical shifts. The obtained results suggest that the accuracy of the MD analysis in predicting geometries and relative conformational energies is not critical and that the final geometry refinements of the structures selected from the MD simulations using QM methods are sufficient for correcting for the expected inaccuracy of the MD analysis. A unique example of C(sp(3))-H···N(sp(3)) intramolecular noncovalent interaction is also identified using the NMR/MD/QM and the natural bond orbital analyses. As the NMR/MD/QM approach relies on the final QM geometry optimization, comparisons of geometric characteristics predicted by different QM methods and those from X-ray and neutron diffraction measurements were undertaken using rigid and flexible cyclic systems. The joint analysis shows that intermolecular noncovalent interactions present in the solid state alter molecular geometries significantly compared to the geometries of isolated molecules from QM calculations.     

1H and 13C NMR chemical shift assignments and conformational analysis for the two diastereomers of the vitamin K epoxide reductase inhibitor brodifacoum

Proton and ¹³C NMR chemical shifts and ¹H? ¹H scalar couplings for the two diastereomers of the potent vitamin K epoxide reductase (VKOR) inhibitor brodifacoum have been determined at 293 K from acetone solutions containing both diastereomers. To facilitate difficult assignments, homo‐ and heteronuclear correlation spectra were acquired at 750 and 900 MHz over 268–303 K temperature range. Conformations of both diastereomers inferred from the scalar couplings and 1‐D NOE measurements reveal that one diastereomer (SS/RR) adopts a strained geometry in the cyclohexene ring system of the tetralin group. The NMR spectra also show evidence of line broadening due to conformational exchange at room temperature for the SR/RS diastereomer. These assignments and conformational analyses may be useful in studies of biomolecular interactions of brodifacoum with target proteins such as VKOR and in source determination of brodifacoum. Copyright © 2009 John Wiley & Sons, Ltd.     

NMR chemical shifts of the rhodopsin chromophore in the dark state and in bathorhodopsin: a hybrid QM/MM molecular dynamics study

We investigate nuclear magnetic resonance (NMR) parameters of the rhodopsin chromophore in the dark state of the protein and in the early photointermediate bathorhodopsin via first-principles molecular dynamics simulations and NMR chemical shift calculations in a hybrid quantum/classical (QM/MM) framework. NMR parameters are particularly sensitive to structural properties and to the chemical environment, which allows us to address different questions about the retinal chromophore in situ. Our calculations show that both the 13C and the 1H NMR chemical shifts are rather insensitive to the protonation state of Glu181, an ionizable amino acid side chain located in the vicinity of the isomerizing 11-cis bond. Thus, other techniques should be better suited to establish its protonation state. The calculated chemical shifts for bathorhodopsin further support our previously published theoretical structure, which is in very good agreement with more recent X-ray data.     

Computational protocols for prediction of solute NMR relative chemical shifts. a case study of L-tryptophan in aqueous solution

In this study, we have applied two different spanning protocols for obtaining the molecular conformations of L-tryptophan in aqueous solution, namely a molecular dynamics simulation and a molecular mechanics conformational search with subsequent geometry re-optimization of the stable conformers using a quantum mechanically based method. These spanning protocols represent standard ways of obtaining a set of conformations on which NMR calculations may be performed. The results stemming from the solute-solvent configurations extracted from the MD simulation at 300 K are found to be inferior to the results stemming from the conformations extracted from the MM conformational search in terms of replicating an experimental reference as well as in achieving the correct sequence of the NMR relative chemical shifts of L-tryptophan in aqueous solution. We find this to be due to missing conformations visited during the molecular dynamics run as well as inaccuracies in geometrical parameters generated from the classical molecular dynamics simulations.     

Effects of Quantum Nuclear Delocalisation on NMR Parameters from Path Integral Molecular Dynamics

The influence of nuclear delocalisation on NMR chemical shifts in molecular organic solids is explored using path integral molecular dynamics (PIMD) and density functional theory calculations of shielding tensors. Nuclear quantum effects are shown to explain previously observed systematic deviations in correlations between calculated and experimental chemical shifts, with particularly large PIMD‐induced changes (up to 23 ppm) observed for carbon atoms in methyl groups. The PIMD approach also enables isotope substitution effects on chemical shifts and J couplings to be predicted in excellent agreement with experiment for both isolated molecules and molecular crystals. An approach based on convoluting calculated shielding or coupling surfaces with probability distributions of selected bond distances and valence angles obtained from PIMD simulations is used to calculate isotope effects.     

Investigation of the differences in activity between hydroxycycloalkyl N1 substituted pyrazole derivatives as inhibitors of B-Raf kinase by using docking, molecular dynamics, QM/MM, and fragment-based de novo design: study of binding mode of diastereomer compounds

N1 substituted pyrazole derivatives show diverse B-Raf kinase inhibitory activities when different hydroxy-substituted cycloalkyl groups are placed at this position. Docking, molecular dynamics (MD) simulations, and hybrid calculation methods (Quantum Mechanics/Molecular Mechanics (QM/MM)) were performed on the complexes, in order to explain these differences. Docking of the inhibitors showed the same orientation that X-ray crystal structure of the analogous (1E)-5-[1-(4-piperidinyl)-3-(4-pyridinyl)-1H-pyrazol-4-yl]-2,3-dihydro-1H-inden-1-one oxime. MD simulations of the most active diastereomer compounds containing cis- and trans-3-hydroxycyclohexyl substituents showed stable interactions with residue Ile463 at the entrance of the B-Raf active site. On the other hand, the less active diastereomer compounds containing cis- and trans-2-hydroxycyclopentyl substituents showed interactions with inner residues Asn580 and Ser465. We found that the differences in activity can be explained by considering the dynamic interactions between the inhibitors and their surrounding residues within the B-Raf binding site. We also explained the activity trend by using a testing scoring function derived from more reliable QM/MM calculations. In addition, we search for new inhibitors from a virtual screening carried out by fragment-based de novo design. We generated a set of approximately 200 virtual compounds, which interact with Ile463 and fulfill druglikeness properties according to Lipinski, Veber, and Ghose rules.     

CHARMM36 all‐atom additive protein force field: Validation based on comparison to NMR data

Protein structure and dynamics can be characterized on the atomistic level with both nuclear magnetic resonance (NMR) experiments and molecular dynamics (MD) simulations. Here, we quantify the ability of the recently presented CHARMM36 (C36) force field (FF) to reproduce various NMR observables using MD simulations. The studied NMR properties include backbone scalar couplings across hydrogen bonds, residual dipolar couplings (RDCs) and relaxation order parameter, as well as scalar couplings, RDCs, and order parameters for side‐chain amino‐ and methyl‐containing groups. It is shown that the C36 FF leads to better correlation with experimental data compared to the CHARMM22/CMAP FF and suggest using C36 in protein simulations. Although both CHARMM FFs contains the same nonbond parameters, our results show how the changes in the internal parameters associated with the peptide backbone via CMAP and the χ₁ and χ₂ dihedral parameters leads to improved treatment of the analyzed nonbond interactions. This highlights the importance of proper treatment of the internal covalent components in modeling nonbond interactions with molecular mechanics FFs. © 2013 Wiley Periodicals, Inc.     

Pi stack structure and hole transfer couplings in DNA hairpins and DNA. A combined QM/MD study

Many key characteristics of hole transfer (HT) in DNA have been derived from spectroscopic studies of DNA hairpins. Because the capping groups in the hairpins can remarkably influence the structure and flexibility of the pi stack, and therefore, the charge transfer rate, the question arises of whether the HT parameters obtained for hairpins may be transferred to DNA oligomers. On the basis of large-time scale QM/MD simulations, we compare structural and electronic parameters of AT stacks in hairpins and DNA oligomers. We find that even in short hairpins, Sa-AA-Sd and Sa-AAA-Sd, the effects of the capping chromophores on the structure of the pi stack and the HT couplings properly averaged over MD trajectories are relatively small, and therefore, the hairpins are good models to study hole transfer through DNA. By contrast, the calculations of the electronic couplings based on the average structures of the systems lead to essential errors in the HT rate and the misleading statement that the charge transfer properties of DNA domains within hairpins are quite different from those of normal sequences.     

Computational characterization of reaction intermediates in the photocycle of the sensory domain of the AppA blue light photoreceptor

The AppA protein with the BLUF (blue light using flavin adenine dinucleotide) domain is a blue light photoreceptor that cycle between dark-adapted and light-induced functional states. We characterized possible reaction intermediates in the photocycle of AppA BLUF. Molecular dynamics (MD), quantum chemical and quantum mechanical-molecular mechanical (QM/MM) calculations were carried out to describe several stable structures of a molecular system modeling the protein. The coordinates of heavy atoms from the crystal structure (PDB code 2IYG) of the protein in the dark state served as starting point for 10 ns MD simulations. Representative MD frames were used in QM(B3LYP/cc-pVDZ)/MM(AMBER) calculations to locate minimum energy configurations of the model system. Vertical electronic excitation energies were estimated for the molecular clusters comprising the quantum subsystems of the QM/MM optimized structures using the SOS-CIS(D) quantum chemistry method. Computational results support the occurrence of photoreaction intermediates that are characterized by spectral absorption bands between those of the dark and light states. They agree with crystal structures of reaction intermediates (PDB code 2IYI) observed in the AppA BLUF domain. Transformations of the Gln63 side chain stimulated by photo-excitation and performed with the assistance of the chromophore and the Met106 side chain are responsible for these intermediates.     

Molecular dynamics simulations of the mononuclear zinc-beta-lactamase from Bacillus cereus.

Herein, we report molecular dynamics simulations of the mononuclear form of the Bacillus cereuszinc-beta-lactamase. We studied two different configurations which differ in the presence of a zinc-bound hydroxide or a zinc-bound water and in the protonation state of the essential His210 residue. Contacts of the catalytically important residues (Asp90, His210, Cys168, etc.) with the zinc center are characterized by the MD analyses. The nature of the Zn-OH(2) --> His210 proton transfer pathway connecting the two configurations was studied by means of QM calculations on cluster models while the relative stability of the two configurations was estimated from QM/MM calculations in the enzyme. From these results, a theoretical model for the kinetically active form of the B. cereus metalloenzyme is proposed. Some mechanistic implications and the influence of mutating the Cys168 residue are also discussed.     

Fluorescence spectra of organic dyes in solution: a time dependent multilevel approach

Classical all-atom molecular dynamics (MD) simulations and quantum mechanical (QM) time-dependent density functional theory (TD-DFT) calculations are employed to study the conformational and photophysical properties of the first emitter excited state of tetramethyl-rhodamine iso-thiocyanate fluorophore in aqueous solution. For this purpose, a specific and accurate force field has been parameterised from QM data to model the fluorophore's first bright excited state. During the MD simulations, the consequences of the π→π* electronic transition on the structure and microsolvation sphere of the dye has been analysed in some detail and compared to the ground state behaviour. Thereafter, fluorescence has been calculated at the TD-DFT level on configurations sampled from the simulated MD trajectories, allowing us to include time dependent solvent effects in the computed emission spectrum. The latter, when compared with the absorption spectrum, reproduces well the experimental Stokes shift, further validating the proposed multilevel computational procedure.     

High-resolution NMR spectra in inhomogeneous fields via IDEAL (intermolecular dipolar-interaction enhanced all lines) method

Intermolecular double-quantum technique is used to yield high-resolution NMR spectra in inhomogeneous magnetic fields. The method exploits the distant dipolar interactions between the solvent and solute nuclear spins. Chemical shifts, J couplings, multiplicity patterns, and relative areas are retained with the method. Except for a 1.5-fold change in the scale factor of J couplings, other parameters are consistent with those extracted from one-dimensional spectra obtained in a homogeneous field.     

Development of massive multilevel molecular dynamics simulation program,platypus (PLATform for dYnamic protein unified simulation), for the elucidation of protein functions

A massively parallel program for quantum mechanical‐molecular mechanical (QM/MM) molecular dynamics simulation, called Platypus (PLATform for dYnamic Protein Unified Simulation), was developed to elucidate protein functions. The speedup and the parallelization ratio of Platypus in the QM and QM/MM calculations were assessed for a bacteriochlorophyll dimer in the photosynthetic reaction center (DIMER) on the K computer, a massively parallel computer achieving 10 PetaFLOPs with 705,024 cores. Platypus exhibited the increase in speedup up to 20,000 core processors at the HF/cc‐pVDZ and B3LYP/cc‐pVDZ, and up to 10,000 core processors by the CASCI(16,16)/6‐31G** calculations. We also performed excited QM/MM‐MD simulations on the chromophore of Sirius (SIRIUS) in water. Sirius is a pH‐insensitive and photo‐stable ultramarine fluorescent protein. Platypus accelerated on‐the‐fly excited‐state QM/MM‐MD simulations for SIRIUS in water, using over 4000 core processors. In addition, it also succeeded in 50‐ps (200,000‐step) on‐the‐fly excited‐state QM/MM‐MD simulations for the SIRIUS in water. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

NMR studies of molecular conformations in alpha-cyclodextrin

A new approach for analysis of NMR parameters is proposed. The experimental data set includes scalar couplings, NOEs, and residual dipolar couplings. The method, which aims at construction of the conformational distribution function, is applied to alpha-cyclodextrin in isotropic solution and dissolved in a dilute liquid crystal. An attempt to analyze the experimental data using an average molecular conformation resulted in unacceptable errors. Our approach rests on the maximum entropy method (ME), which gives the flattest possible distribution, consistent with the experimental data. Very good agreement between experimental and calculated NMR parameters was observed. In fact, two conformational states were required in order to obtain a satisfactory agreement between calculated and experimental data. In addition, good agreement with Langevin dynamics computer simulations was obtained.     

QM/MM calculation of solvent effects on absorption spectra of guanine

Electronic spectra of guanine in the gas phase and in water were studied by quantum mechanical/molecular mechanical (QM/MM) methods. Geometries for the excited‐state calculations were extracted from ground‐state molecular dynamics (MD) simulations using the self‐consistent‐charge density functional tight binding (SCC‐DFTB) method for the QM region and the TIP3P force field for the water environment. Theoretical absorption spectra were generated from excitation energies and oscillator strengths calculated for 50 to 500 MD snapshots of guanine in the gas phase (QM) and in solution (QM/MM). The excited‐state calculations used time‐dependent density functional theory (TDDFT) and the DFT‐based multireference configuration interaction (DFT/MRCI) method of Grimme and Waletzke, in combination with two basis sets. Our investigation covered keto‐N7H and keto‐N9H guanine, with particular focus on solvent effects in the low‐energy spectrum of the keto‐N9H tautomer. When compared with the vertical excitation energies of gas‐phase guanine at the optimized DFT (B3LYP/TZVP) geometry, the maxima in the computed solution spectra are shifted by several tenths of an eV. Three effects contribute: the use of SCC‐DFTB‐based rather than B3LYP‐based geometries in the MD snapshots (red shift of ca. 0.1 eV), explicit inclusion of nuclear motion through the MD snapshots (red shift of ca. 0.1 eV), and intrinsic solvent effects (differences in the absorption maxima in the computed gas‐phase and solution spectra, typically ca. 0.1–0.3 eV). A detailed analysis of the results indicates that the intrinsic solvent effects arise both from solvent‐induced structural changes and from electrostatic solute–solvent interactions, the latter being dominant. © 2009 Wiley Periodicals, Inc. J Comput Chem 2010     

Force field parameters for rotation around chi torsion axis in nucleic acids

To raise the accuracy of the force field for nucleic acids, several parameters were elaborated, focusing on the rotation around chi torsion axis. The reliability of molecular dynamics (MD) simulation was significantly increased by improving the torsion parameters at C8--N9--C1'--X (X = H1', C2', O4') in A, G and those at C6--N1--C1'--X in C, T, and U. In this work, we constructed small models representing the chemical structure of A, G, C, T, and U, and estimated energy profile for chi-axis rotation by executing numerous quantum mechanical (QM) calculations. The parameters were derived by discrete Fourier transformation of the calculated QM data. A comparison in energy profile between molecular mechanical (MM) calculation and QM one shows that our presently derived parameters well reproduce the energy surface of QM calculation for all the above torsion terms. Furthermore, our parameters show a good performance in MD simulations of some nucleic acids. Hence, the present refinement of parameters will enable us to perform more accurate simulations for various types of nucleic acids.     

Towards Accurate Predictions of Proton NMR Spectroscopic Parameters in Molecular Solids

The factors contributing to the accuracy of quantum-chemical calculations for the prediction of proton NMR chemical shifts in molecular solids are systematically investigated. Proton chemical shifts of six solid amino acids with hydrogen atoms in various bonding environments (CH, CH₂, CH₃, OH, SH and NH₃) were determined experimentally using ultra-fast magic-angle spinning and proton-detected 2D NMR experiments. The standard DFT method commonly used for the calculations of NMR parameters of solids is shown to provide chemical shifts that deviate from experiment by up to 1.5 ppm. The effects of the computational level (hybrid DFT functional, coupled-cluster calculation, inclusion of relativistic spin-orbit coupling) are thoroughly discussed. The effect of molecular dynamics and nuclear quantum effects are investigated using path-integral molecular dynamics (PIMD) simulations. It is demonstrated that the accuracy of the calculated proton chemical shifts is significantly better when these effects are included in the calculations.     

Fragmentation-based QM/MM simulations: length dependence of chain dynamics and hydrogen bonding of polyethylene oxide and polyethylene in aqueous solutions

Molecular fragmentation quantum mechanics (QM) calculations have been combined with molecular mechanics (MM) to construct the fragmentation QM/MM method for simulations of dilute solutions of macromolecules. We adopt the electrostatics embedding QM/MM model, where the low-cost generalized energy-based fragmentation calculations are employed for the QM part. Conformation energy calculations, geometry optimizations, and Born-Oppenheimer molecular dynamics simulations of poly(ethylene oxide), PEO(n) (n = 6-20), and polyethylene, PE(n) ( n = 9-30), in aqueous solution have been performed within the framework of both fragmentation and conventional QM/MM methods. The intermolecular hydrogen bonding and chain configurations obtained from the fragmentation QM/MM simulations are consistent with the conventional QM/MM method. The length dependence of chain conformations and dynamics of PEO and PE oligomers in aqueous solutions is also investigated through the fragmentation QM/MM molecular dynamics simulations.     

Analytic energy gradient in combined second-order Møller-Plesset perturbation theory and polarizable force field calculation

Second-order M?ller-Plesset perturbation theory (MP2) is used to describe electronic correlation on the basis of Hartree-Fock (HF) variational calculations that incorporate induced dipole polarizable force fields (i.e., QM/MMpol style HF and MP2). The Z-vector equations for regular closed shell and open shell MP2 methods (RMP2, ZAPT2, and UMP2) are extended to include induced dipole contributions to determine the MP2 response density so that nuclear gradient and other properties can be efficiently evaluated. A better estimation of the induced dipole polarization energy can be obtained using the MP2 relaxed density. QM/MMpol style MP2 molecular dynamics simulations are performed for the ground state and first triplet state of acetone solvated by 1024 polarizable water molecules. A switching function is used to ensure energy conservation in QM/MM simulation under periodic boundary condition.