Fragment density functional theory calculation of NMR chemical shifts for proteins with implicit solvation

PW-SMD: A Plane-Wave Implicit Solvation Model Based on Electron Density for Surface Chemistry and Crystalline Systems in Aqueous Solution

Electron density-based implicit solvation models are a class of techniques for quantifying solvation effects and calculating free energies of solvation without an explicit representation of solvent molecules. Integral to the accuracy of solvation modeling is the proper definition of the solvation shell separating the solute molecule from the solvent environment, allowing for a physical partitioning of the free energies of solvation. Unlike state-of-the-art implicit solvation models for molecular quantum chemistry calculations, e.g., the solvation model based on solute electron density (SMD), solvation models for systems under periodic boundary conditions with plane-wave (PW) basis sets have been limited in their accuracy. Furthermore, a unified implicit solvation model with both homogeneous solution-phase and heterogeneous interfacial structures treated on equal footing is needed. In order to address this challenge, we developed a high-accuracy solvation model for periodic PW calculations that is applicable to molecular, ionic, interfacial, and bulk-phase chemistry. Our model, PW-SMD, is an extension of the SMD molecular solvation model to periodic systems in water. The free energy of solvation is partitioned into the electrostatic and cavity-dispersion-solvent structure (CDS) contributions. The electrostatic contributions of the solvation shell surrounding solute structures are parametrized based on their geometric and physical properties. In addition, the nonelectrostatic contribution to the solvation energy is accounted for by extending the CDS formalism of SMD to incorporate periodic boundary conditions. We validate the accuracy and robustness of our solvation model by comparing predicted solvation free energies against experimental data for molecular and ionic systems, carved-cluster composite energetic models of solvated reaction energies and barriers on surface systems, and deep-learning-accelerated ab initio molecular dynamics (AIMD). Our developed periodic implicit solvation model shows significantly improved accuracy compared to previous work (namely, solvation models in aqueous solution) and can be applied to simulate solvent effects in a wide range of surface and crystalline materials.

A Structure Independent Molecular Fragment Interfuse Model for Mesoscale Dissipative Particle Dynamics Simulation of Peptides

There is a need to develop robust computational models for mesoscale simulation of the structure of peptides over large length scales toward the discovery of novel peptides for medical applications to address the issues of peptide aggregation, enzymatic degradation, and short half-life. The primary objective was to predict the structure and conformation of peptides whose native structures are not known. This work presents a new model for computation of interaction parameters between the beads in coarse-grained dissipative particle dynamics (DPD) simulation that is properly calibrated for amino acids, supports compressibility requirement of water molecules, and accounts for subtle differences in the structure of amino acids and the charge in the side chain of charged amino acids. This new model is referred to as Structure Independent Molecular Fragment Interfuse Model, abbreviated as SIMFIM, because it accounts for specific interactions between different beads, which represent molecular fragments of the amino acids, in calculating nonbonded interaction parameters in the absence of knowing the actual peptide structure. The electrostatic interactions are incorporated in this model by using a normal distribution of charges around the center of the beads to prevent the collapse of oppositely charged soft beads. The uniquely parameterized DPD force field in the SIMFIM model is optimized for a given peptide with respect to the degree of coarse-grained graining for simulating the peptide over long times and length scales. The SIMFIM model was tested in this work using four peptides, namely, TrpZip2, Rubrivinodin, Lihuanodin, and IC3-CB1/Gai peptides, whose structures were sourced from the Protein Data Bank. The SIMFIM model predicted radius of gyration (Rg) values for the peptides closer to the actual structures as compared to the conventional model, and there was less deviation between the predicted and actual structures of the peptides.

Predicted and Experimental NMR Chemical Shifts at Variable Temperatures: The Effect of Protein Conformational Dynamics

NMR chemical shifts provide a sensitive probe of protein structure and dynamics but remain challenging to predict and interpret. We examine the effect of protein conformational distributions on 15N chemical shifts for dihydrofolate reductase (DHFR), comparing QM/MM predicted shifts with experimental shifts in solution as well as frozen distributions. Representative snapshots from MD trajectories exhibit variation in predicted 15N chemical shifts of up to 25 ppm. The average over the fluctuations is in significantly better agreement with room temperature solution experimental values than the prediction for any single optimal conformations. Meanwhile, solid-state NMR (SSNMR) measurements of frozen solutions at 105 K exhibit broad lines whose widths agree well with the widths of distributions of predicted shifts for samples from the trajectory. The backbone torsion angle ψi-1 varies over 60° on the picosecond time scale, compensated by φi. These fluctuations can explain much of the shift variation.

Automated Fragmentation Quantum Mechanical Calculation of 15N and 13C Chemical Shifts in a Membrane Protein

In this work, we developed an accurate and cost-effective automated fragmentation quantum mechanics/molecular mechanics (AF-QM/MM) method to calculate the chemical shifts of 15N and 13C of membrane proteins. The convergence of the AF-QM/MM method was tested using Krokinobacter eikastus rhodopsin 2 as a test case. When the distance threshold of the QM region is equal to or larger than 4.0 Å, the results of the AF-QM/MM calculations are close to convergence. In addition, the effects of selected density functionals, basis sets, and local chemical environment of target atoms on the chemical shift calculations were systematically investigated. Our results demonstrate that the predicted chemical shifts are more accurate when important environmental factors including cross-protomer interactions, lipid molecules, and solvent molecules are taken into consideration, especially for the 15N chemical shift prediction. Furthermore, with the presence of sodium ions in the environment, the chemical shift of residues, retinal, and retinal Schiff base are affected, which is consistent with the results of the solid-state nuclear magnetic resonance (NMR) experiment. Upon comparing the performance of various density functionals (namely, B3LYP, B3PW91, M06-2X, M06-L, mPW1PW91, OB95, and OPBE), the results show that mPW1PW91 is a suitable functional for the 15N and 13C chemical shift prediction of the membrane proteins. Meanwhile, we find that the improved accuracy of the 13Cβ chemical shift calculations can be achieved by the employment of the triple-ζ basis set. However, the employment of the triple-ζ basis set does not improve the accuracy of the 15N and 13Cα chemical shift calculations nor does the addition of a diffuse function improve the overall prediction accuracy of the chemical shifts. Our study also underscores that the AF-QM/MM method has significant advantages in predicting the chemical shifts of key ligands and nonstandard residues in membrane proteins than most widely used empirical models; therefore, it could be an accurate computational tool for chemical shift calculations on various types of biological systems.

Predicted and Experimental NMR Chemical Shifts at Variable Temperatures: The Effect of Protein Conformational Dynamics

NMR chemical shifts provide a sensitive probe of protein structure and dynamics. Prediction of shifts, and therefore interpretation of shifts, particularly for the frequently measured amidic 15 N sites, remains a tall challenge. We demonstrate that protein 15 N chemical shift prediction from QM/MM predictions can be improved if conformational variation is included via MD sampling, focusing on the antibiotic target, E. coli Dihydrofolate reductase (DHFR). Variations of up to 25 ppm in predicted 15 N chemical shifts are observed over the trajectory. For solution shifts the average of fluctuations on the low picosecond timescale results in a superior prediction to a single optimal conformation. For low temperature solid state measurements, the histogram of predicted shifts for locally minimized snapshots with specific solvent arrangements sampled from the trajectory explains the heterogeneous linewidths; in other words, the conformations and associated solvent are 'frozen out' at low temperatures and result in inhomogeneously broadened NMR peaks. We identified conformational degrees of freedom that contribute to chemical shift variation. Backbone torsion angles show high amplitude fluctuations during the trajectory on the low picosecond timescale. For a number of residues, including I60, ψ varies by up to 60º within a conformational basin during the MD simulations, despite the fact that I60 (and other sites studied) are in a secondary structure element and remain well folded during the trajectory. Fluctuations in ψ appear to be compensated by other degrees of freedom in the protein, including φ of the succeeding residue, resulting in "rocking" of the amide plane with changes in hydrogen bonding interactions. Good agreement for both room temperature and low temperature NMR spectra provides strong support for the specific approach to conformational averaging of computed chemical shifts.

Improving IDP theoretical chemical shift accuracy and efficiency through a combined MD/ADMA/DFT and machine learning approach

This work extends the multi-scale computational scheme for the quantum mechanics (QM) calculations of Nuclear Magnetic Resonance (NMR) chemical shifts (CSs) in proteins that lack a well-defined 3D structure. The scheme couples the sampling of an intrinsically disordered protein (IDP) by classical molecular dynamics (MD) with protein fragmentation using the adjustable density matrix assembler (ADMA) and density functional theory (DFT) calculations. In contrast to our early investigation on IDPs (Pavlíková Přecechtělová et al., J. Chem. Theory Comput., 2019, 15, 5642-5658) and the state-of-the art NMR calculations for structured proteins, a partial re-optimization was implemented on the raw MD geometries in vibrational normal mode coordinates to enhance the accuracy of the MD/ADMA/DFT computational scheme. In addition, machine-learning based cluster analysis was performed on the scheme to explore its potential in producing protein structure ensembles (CLUSTER ensembles) that yield accurate CSs at a reduced computational cost. The performance of the cluster-based calculations is validated against results obtained with conventional structural ensembles consisting of MD snapshots extracted from the MD trajectory at regular time intervals (REGULAR ensembles). CS calculations performed with the refined MD/ADMA/DFT framework employed the 6-311++G(d,p) basis set that outperformed IGLO-III calculations with the same density functional approximation (B3LYP) and both explicit and implicit solvation. The partial geometry optimization did not universally improve the agreement of computed CSs with the experiment but substantially decreased errors associated with the ensemble averaging. A CLUSTER ensemble with 50 structures yielded ensemble averages close to those obtained with a REGULAR ensemble consisting of 500 MD frames. The cluster based calculations thus required only a fraction of the computational time.

Dynamic Coupling of Tyrosine 185 with the Bacteriorhodopsin Photocycle, as Revealed by Chemical Shifts, Assisted AF-QM/MM Calculations and Molecular Dynamic Simulations

Aromatic residues are highly conserved in microbial photoreceptors and play crucial roles in the dynamic regulation of receptor functions. However, little is known about the dynamic mechanism of the functional role of those highly conserved aromatic residues during the receptor photocycle. Tyrosine 185 (Y185) is one of the highly conserved aromatic residues within the retinal binding pocket of bacteriorhodopsin (bR). In this study, we explored the molecular mechanism of its dynamic coupling with the bR photocycle by automated fragmentation quantum mechanics/molecular mechanics (AF-QM/MM) calculations and molecular dynamic (MD) simulations based on chemical shifts obtained by 2D solid-state NMR correlation experiments. We observed that Y185 plays a significant role in regulating the retinal cis–trans thermal equilibrium, stabilizing the pentagonal H-bond network, participating in the orientation switch of Schiff Base (SB) nitrogen, and opening the F42 gate by interacting with the retinal and several key residues along the proton translocation channel. Our findings provide a detailed molecular mechanism of the dynamic couplings of Y185 and the bR photocycle from a structural perspective. The method used in this paper may be applied to the study of other microbial photoreceptors.

Fragment-Based Ab Initio Molecular Dynamics Simulation for Combustion

We develop a fragment-based ab initio molecular dynamics (FB-AIMD) method for efficient dynamics simulation of the combustion process. In this method, the intermolecular interactions are treated by a fragment-based many-body expansion in which three- or higher body interactions are neglected, while two-body interactions are computed if the distance between the two fragments is smaller than a cutoff value. The accuracy of the method was verified by comparing FB-AIMD calculated energies and atomic forces of several different systems with those obtained by standard full system quantum calculations. The computational cost of the FB-AIMD method scales linearly with the size of the system, and the calculation is easily parallelizable. The method is applied to methane combustion as a benchmark. Detailed reaction network of methane reaction is analyzed, and important reaction species are tracked in real time. The current result of methane simulation is in excellent agreement with known experimental findings and with prior theoretical studies.

Automated fragmentation quantum mechanical calculation of 13C and 1H chemical shifts in molecular crystals

In this work, the automated fragmentation quantum mechanics/molecular mechanics (AF-QM/MM) approach was applied to calculate the ¹³C and ¹H nuclear magnetic resonance (NMR) chemical shifts in molecular crystals. Two benchmark sets of molecular crystals were selected to calculate the NMR chemical shifts. Systematic investigation was conducted to examine the convergence of AF-QM/MM calculations and the impact of various density functionals with different basis sets on the NMR chemical shift prediction. The result demonstrates that the calculated NMR chemical shifts are close to convergence when the distance threshold for the QM region is larger than 3.5 Å. For ¹³C chemical shift calculations, the mPW1PW91 functional is the best density functional among the functionals chosen in this study (namely, B3LYP, B3PW91, M06-2X, M06-L, mPW1PW91, OB98, and OPBE), while the OB98 functional is more suitable for the ¹H NMR chemical shift prediction of molecular crystals. Moreover, with the B3LYP functional, at least a triple-ζ basis set should be utilized to accurately reproduce the experimental ¹³C and ¹H chemical shifts. The employment of diffuse basis functions will further improve the accuracy for ¹³C chemical shift calculations, but not for the ¹H chemical shift prediction. We further proposed a fragmentation scheme of dividing the central molecule into smaller fragments. By comparing with the results of the fragmentation scheme using the entire central molecule as the core region, the AF-QM/MM calculations with the fragmented central molecule can not only achieve accurate results but also reduce the computational cost. Therefore, the AF-QM/MM approach is capable of predicting the ¹³C and ¹H NMR chemical shifts for molecular crystals accurately and effectively, and could be utilized for dealing with more complex periodic systems such as macromolecular polymers and biomacromolecules. The AF-QM/MM program for molecular crystals is available at https://github.com/shiman1995/NMR.

Using quantum chemistry to estimate chemical shifts in biomolecules

An automated fragmentation quantum mechanics/molecular mechanics approach (AFNMR) has shown promising results in chemical shift calculations for biomolecules. Sample results for ubiquitin, and an RNA hairpin and helix are presented, and used to recent directions in quantum calculations. Trends in chemical shift are stable with regards to change in density functional or basis sets, and the use of the small "pcSseg-0" basis, which was optimized for chemical shift prediction [1], opens the way to more extensive conformational averaging, which can often be necessary, even for fairly well-defined structures.

Polarizable continuum models provide an effective electrostatic embedding model for fragment-based chemical shift prediction in challenging systems

Ab initio nuclear magnetic resonance chemical shift prediction provides an important tool for interpreting and assigning experimental spectra, but it becomes computationally prohibitive in large systems. The computational costs can be reduced considerably by fragmentation of the large system into a series of contributions from many smaller subsystems. However, the presence of charged functional groups and the need to partition the system across covalent bonds create complications in biomolecules that typically require the use of large fragments and careful descriptions of the electrostatic environment. The present work shows how a model that combines chemical shielding contributions from non-overlapping monomer and dimer fragments embedded in a polarizable continuum model provides a simple, easy-to-implement, and computationally inexpensive approach for predicting chemical shifts in complex systems. The model's performance proves rather insensitive to the continuum dielectric constant, making the selection of the optimal embedding dielectric less critical. The PCM-embedded fragment model is demonstrated to perform well across systems ranging from molecular crystals to proteins.

Accurate Backbone 13 C and 15 N Chemical Shift Tensors in Galectin-3 Determined by MAS NMR and QM/MM: Details of Structure and Environment Matter

Chemical shift tensors obtained from solid-state NMR spectroscopy are very sensitive reporters of structure and dynamics in proteins. While accurate 13 C and 15 N chemical shift tensors are accessible by magic angle spinning (MAS) NMR, their quantum mechanical calculations remain challenging, particularly for 15 N atoms. Here we compare experimentally determined backbone 13 Cα and 15 NH chemical shift tensors by MAS NMR with hybrid quantum mechanics/molecular mechanics/molecular dynamics (MD-QM/MM) calculations for the carbohydrate-binding domain of galectin-3. Excellent agreement between experimental and computed 15 NH chemical shift anisotropy values was obtained using the Amber ff15ipq force field when solvent dynamics was taken into account in the calculation. Our results establish important benchmark conditions for improving the accuracy of chemical shift calculations in proteins and may aid in the validation of protein structure models derived by MAS NMR.

Revealing binding selectivity of inhibitors toward bromodomain-containing proteins 2 and 4 using multiple short molecular dynamics simulations and free energy analyses

Emerging evidences indicate bromodomain-containing proteins 2 and 4 (BRD2 and BRD4) play critical roles in cancers, inflammations, cardiovascular diseases and other pathologies. Multiple short molecular dynamics (MSMD) simulations combined with molecular mechanics generalized Born surface area (MM-GBSA) method were applied to investigate the binding selectivity of three inhibitors 87D, 88M and 89G towards BRD2 over BRD4. The root-mean-square fluctuation (RMSF) analysis indicates that the structural flexibility of BRD4 is stronger than that of BRD2. Moreover the calculated distances between the C_α atoms in the centres of the ZA_loop and BC_loop of BRD4 are also bigger than that of BRD2. The rank of binding free energies calculated using MM-GBSA method agrees well with that determined by experimental data. The results show that 87D can bind more favourably to BRD2 than BRD4, while 88M has better selectivity on BRD4 over BRD2. Residue-based free-energy decomposition method was utilized to estimate the inhibitor-residue interaction spectrum and the results not only identify the hot interaction spots of inhibitors with BRD2 and BRD4, but also demonstrate that several common residues, including (W370, W374), (P371, P375), (V376, V380) and (L381, L385) belonging to (BRD2, BRD4), generate significant binding difference of inhibitors to BRD2 and BRD4.

Molecular mechanism of inhibitor bindings to bromodomain-containing protein 9 explored based on molecular dynamics simulations and calculations of binding free energies

Recently, bromodomain-containing protein 9 (BRD9) has been a prospective therapeutic target for anticancer drug design. Molecular dynamics (MD) simulations combined with molecular mechanics generalized Born surface area (MM-GBSA) method were adopted to explore binding modes of three inhibitors (5SW, 5U2, and 5U6) to BRD9 and identify the hot spot of the inhibitor-BRD9 binding. The results indicate that the inhibitor 5SW has the strongest binding ability to BRD9 among the current three inhibitors. Furthermore, the rank of the binding free energies predicted by MM-GBSA approach agrees with that determined by the experimental values. In addition, inhibitor-residue interactions were computed by using residue-based free-energy decomposition method and the results suggest that residue His42 produces the CH-H interactions, residues Asn100, Ile53 and Val49 produce the CH-[Formula: see text] interactions with three inhibitors and Tyr106, Phe45 and Phe44 generate the π-π interactions with inhibitors. Notably, the residue Asn140 forms hydrogen bonding interactions with three inhibitors. This research is expected to provide useful molecular basis and dynamics information at atomic levels for the design of potent inhibitors inhibiting the activity of BRD9.

Fantasy versus reality in fragment-based quantum chemistry

Since the introduction of the fragment molecular orbital method 20 years ago, fragment-based approaches have occupied a small but growing niche in quantum chemistry. These methods decompose a large molecular system into subsystems small enough to be amenable to electronic structure calculations, following which the subsystem information is reassembled in order to approximate an otherwise intractable supersystem calculation. Fragmentation sidesteps the steep rise (with respect to system size) in the cost of ab initio calculations, replacing it with a distributed cost across numerous computer processors. Such methods are attractive, in part, because they are easily parallelizable and therefore readily amenable to exascale computing. As such, there has been hope that distributed computing might offer the proverbial "free lunch" in quantum chemistry, with the entrée being high-level calculations on very large systems. While fragment-based quantum chemistry can count many success stories, there also exists a seedy underbelly of rarely acknowledged problems. As these methods begin to mature, it is time to have a serious conversation about what they can and cannot be expected to accomplish in the near future. Both successes and challenges are highlighted in this Perspective.

A Fragment Quantum Mechanical Method for Metalloproteins

An accurate energy calculation of metalloprotein is of crucial importance and also a theoretical challenge. In this work, a metal molecular fractionation with conjugate caps (metal-MFCC) approach is developed for efficient linear-scaling quantum calculation of potential energy and atomic forces of metalloprotein. In this approach, the potential energy of a given protein is calculated by a linear combination of potential energies of the neighboring residues, two-body interaction energy between non-neighboring residues that are spatially in close contact and the potential energy of the metal binding group. The calculation of each fragment is embedded in a field of point charges representing the remaining protein environment. Numerical studies were carried out to check the performance of this method, and the calculated potential energies and atomic forces all show excellent agreement with the full system calculations at the M06-2X/6-31G(d) level. By combining the energy calculation with molecular dynamic simulation, we performed an ab initio structural optimization for a zinc finger protein with high efficiency. The present metal-MFCC approach is linear-scaling with a low prefactor and trivially parallelizable. The individual fragment typically contains about 50 atoms, and it is thus possible to be calculated at higher levels of the quantum chemistry method. This fragment method can be routinely applied to perform structural optimization and ab initio molecular dynamic simulation for metalloproteins of any size.

Functional roles of tyrosine 185 during the bacteriorhodopsin photocycle as revealed by in situ spectroscopic studies

Tyrosine 185 (Y185), one of the aromatic residues within the retinal (Ret) chromophore binding pocket in helix F of bacteriorhodopsin (bR), is highly conserved among the microbial rhodopsin family proteins. Many studies have investigated the functions of Y185, but its underlying mechanism during the bR photocycle remains unclear. To address this research gap, in situ two-dimensional (2D) magic-angle spinning (MAS) solid-state NMR (ssNMR) of specifically labelled bR, combined with light-induced transient absorption change measurements, dynamic light scattering (DLS) measurements, titration analysis and site-directed mutagenesis, was used to elucidate the functional roles of Y185 during the bR photocycle in the native membrane environment. Different interaction modes were identified between Y185 and the Ret chromophore in the dark-adapted (inactive) state and M (active) state, indicating that Y185 may serve as a rotamer switch maintaining the protein dynamics, and plays an important role in the efficient proton-pumping mechanism in the bR purple membrane.

Automated Fragmentation QM/MM Calculation of NMR Chemical Shifts for Protein-Ligand Complexes

In this study, the automated fragmentation quantum mechanics/molecular mechanics (AF-QM/MM) method was applied for NMR chemical shift calculations of protein-ligand complexes. In the AF-QM/MM approach, the protein binding pocket is automatically divided into capped fragments (within ~200 atoms) for density functional theory (DFT) calculations of NMR chemical shifts. Meanwhile, the solvent effect was also included using the Poission-Boltzmann (PB) model, which properly accounts for the electrostatic polarization effect from the solvent for protein-ligand complexes. The NMR chemical shifts of neocarzinostatin (NCS)-chromophore binding complex calculated by AF-QM/MM accurately reproduce the large-sized system results. The 1H chemical shift perturbations (CSP) between apo-NCS and holo-NCS predicted by AF-QM/MM are also in excellent agreement with experimental results. Furthermore, the DFT calculated chemical shifts of the chromophore and residues in the NCS binding pocket can be utilized as molecular probes to identify the correct ligand binding conformation. By combining the CSP of the atoms in the binding pocket with the Glide scoring function, the new scoring function can accurately distinguish the native ligand pose from decoy structures. Therefore, the AF-QM/MM approach provides an accurate and efficient platform for protein-ligand binding structure prediction based on NMR derived information.

Full QM Calculation of RNA Energy Using Electrostatically Embedded Generalized Molecular Fractionation with Conjugate Caps Method

In this study, the electrostatically embedded generalized molecular fractionation with conjugate caps (concaps) method (EE-GMFCC) was employed for efficient linear-scaling quantum mechanical (QM) calculation of total energies of RNAs. In the EE-GMFCC approach, the total energy of RNA is calculated by taking a proper combination of the QM energy of each nucleotide-centric fragment with large caps or small caps (termed EE-GMFCC-LC and EE-GMFCC-SC, respectively) deducted by the energies of concaps. The two-body QM interaction energy between non-neighboring ribonucleotides which are spatially in close contact are also taken into account for the energy calculation. Numerical studies were carried out to calculate the total energies of a number of RNAs using the EE-GMFCC-LC and EE-GMFCC-SC methods at levels of the Hartree-Fock (HF) method, density functional theory (DFT), and second-order many-body perturbation theory (MP2), respectively. The results show that the efficiency of the EE-GMFCC-SC method is about 3 times faster than the EE-GMFCC-LC method with minimal accuracy sacrifice. The EE-GMFCC-SC method is also applied for relative energy calculations of 20 different conformers of two RNA systems using HF and DFT, respectively. Both single-point and relative energy calculations demonstrate that the EE-GMFCC method has deviations from the full system results of only a few kcal/mol.

Automated Fragmentation Polarizable Embedding Density Functional Theory (PE-DFT) Calculations of Nuclear Magnetic Resonance (NMR) Shielding Constants of Proteins with Application to Chemical Shift Predictions

Full-protein nuclear magnetic resonance (NMR) shielding constants based on ab initio calculations are desirable, because they can assist in elucidating protein structures from NMR experiments. In this work, we present NMR shielding constants computed using a new automated fragmentation (J. Phys. Chem. B 2009, 113, 10380-10388) approach in the framework of polarizable embedding density functional theory. We extend our previous work to give both basis set recommendations and comment on how large the quantum mechanical region should be to successfully compute ¹³C NMR shielding constants that are comparable with experiment. The introduction of a probabilistic linear regression model allows us to substantially reduce the number of snapshots that are needed to make comparisons with experiment. This approach is further improved by augmenting snapshot selection with chemical shift predictions by which we can obtain a representative subset of snapshots that gives the smallest predicted error, compared to experiment. Finally, we use this subset of snapshots to calculate the NMR shielding constants at the PE-KT3/pcSseg-2 level of theory for all atoms in the protein GB3.

Protein structure refinement using a quantum mechanics-based chemical shielding predictor

We show that a QM-based predictor of a protein backbone and CB chemical shifts is of comparable accuracy to empirical chemical shift predictors after chemical shift-based structural refinement that removes small structural errors (errors in chemical shifts shown in red).

The accurate prediction of protein chemical shifts using a quantum mechanics (QM)-based method has been the subject of intense research for more than 20 years but so far empirical methods for chemical shift prediction have proven more accurate. In this paper we show that a QM-based predictor of a protein backbone and CB chemical shifts (ProCS15, PeerJ, 2016, 3, e1344) is of comparable accuracy to empirical chemical shift predictors after chemical shift-based structural refinement that removes small structural errors. We present a method by which quantum chemistry based predictions of isotropic chemical shielding values (ProCS15) can be used to refine protein structures using Markov Chain Monte Carlo (MCMC) simulations, relating the chemical shielding values to the experimental chemical shifts probabilistically. Two kinds of MCMC structural refinement simulations were performed using force field geometry optimized X-ray structures as starting points: simulated annealing of the starting structure and constant temperature MCMC simulation followed by simulated annealing of a representative ensemble structure. Annealing of the CHARMM structure changes the CA-RMSD by an average of 0.4 Å but lowers the chemical shift RMSD by 1.0 and 0.7 ppm for CA and N. Conformational averaging has a relatively small effect (0.1–0.2 ppm) on the overall agreement with carbon chemical shifts but lowers the error for nitrogen chemical shifts by 0.4 ppm. If an amino acid specific offset is included the ProCS15 predicted chemical shifts have RMSD values relative to experiments that are comparable to popular empirical chemical shift predictors. The annealed representative ensemble structures differ in CA-RMSD relative to the initial structures by an average of 2.0 Å, with >2.0 Å difference for six proteins. In four of the cases, the largest structural differences arise in structurally flexible regions of the protein as determined by NMR, and in the remaining two cases, the large structural change may be due to force field deficiencies. The overall accuracy of the empirical methods are slightly improved by annealing the CHARMM structure with ProCS15, which may suggest that the minor structural changes introduced by ProCS15-based annealing improves the accuracy of the protein structures. Having established that QM-based chemical shift prediction can deliver the same accuracy as empirical shift predictors we hope this can help increase the accuracy of related approaches such as QM/MM or linear scaling approaches or interpreting protein structural dynamics from QM-derived chemical shift.

Enhanced NMR Discrimination of Pharmaceutically Relevant Molecular Crystal Forms through Fragment-Based Ab Initio Chemical Shift Predictions

Chemical shift prediction plays an important role in the determination or validation of crystal structures with solid-state nuclear magnetic resonance (NMR) spectroscopy. One of the fundamental theoretical challenges lies in discriminating variations in chemical shifts resulting from different crystallographic environments. Fragment-based electronic structure methods provide an alternative to the widely used plane wave gauge-including projector augmented wave (GIPAW) density functional technique for chemical shift prediction. Fragment methods allow hybrid density functionals to be employed routinely in chemical shift prediction, and we have recently demonstrated appreciable improvements in the accuracy of the predicted shifts when using the hybrid PBE0 functional instead of generalized gradient approximation (GGA) functionals like PBE. Here, we investigate the solid-state ¹³C and ¹⁵N NMR spectra for multiple crystal forms of acetaminophen, phenobarbital, and testosterone. We demonstrate that the use of the hybrid density functional instead of a GGA provides both higher accuracy in the chemical shifts and increased discrimination among the different crystallographic environments. Finally, these results also provide compelling evidence for the transferability of the linear regression parameters mapping predicted chemical shieldings to chemical shifts that were derived in an earlier study.

Subsystem Analysis for the Fragment Molecular Orbital Method and Its Application to Protein-Ligand Binding in Solution

A subsystem analysis is derived incorporating interfragment interactions into the fragment properties, such as energies or charges. The relative stabilities of three alanine isomers, the α-helix, the β-turn, and the extended form are studied and the differences in fragment properties are elucidated. The analysis is further elaborated for studies of binding energies. The binding of the Trp-cage protein (PDB: 1L2Y ) to two ligands is studied in detail. Binding energies defined for each fragment can be used as a convenient descriptor for analyzing contributions to binding in solution.

Modeling Polymorphic Molecular Crystals with Electronic Structure Theory

Interest in molecular crystals has grown thanks to their relevance to pharmaceuticals, organic semiconductor materials, foods, and many other applications. Electronic structure methods have become an increasingly important tool for modeling molecular crystals and polymorphism. This article reviews electronic structure techniques used to model molecular crystals, including periodic density functional theory, periodic second-order Møller-Plesset perturbation theory, fragment-based electronic structure methods, and diffusion Monte Carlo. It also discusses the use of these models for predicting a variety of crystal properties that are relevant to the study of polymorphism, including lattice energies, structures, crystal structure prediction, polymorphism, phase diagrams, vibrational spectroscopies, and nuclear magnetic resonance spectroscopy. Finally, tools for analyzing crystal structures and intermolecular interactions are briefly discussed.

The influence of Mg2+ coordination on 13C and 15N chemical shifts in CKI1RD protein domain from experiment and molecular dynamics/density functional theory calculations

Sequence dependence of (13) C and (15) N chemical shifts in the receiver domain of CKI1 protein from Arabidopsis thaliana, CKI1RD , and its complexed form, CKI1RD •Mg(2+), was studied by means of MD/DFT calculations. MD simulations of a 20-ns production run length were performed. Nine explicitly hydrated structures of increasing complexity were explored, up to a 40-amino-acid structure. The size of the model necessary depended on the type of nucleus, the type of amino acid and its sequence neighbors, other spatially close amino acids, and the orientation of amino acid NH groups and their surface/interior position. Using models covering a 10 and a 15 Å environment of Mg(2+), a semi-quantitative agreement has been obtained between experiment and theory for the V67-I73 sequence. The influence of Mg(2+) binding was described better by the 15 Å as compared to the 10 Å model. Thirteen chemical shifts were analyzed in terms of the effect of Mg(2+) insertion and geometry preparation. The effect of geometry was significant and opposite in sign to the effect of Mg(2+) binding. The strongest individual effects were found for (15) N of D70, S74, and V68, where the electrostatics dominated; for (13) Cβ of D69 and (15) N of K76, where the influences were equal, and for (13) Cα of F72 and (13) Cβ of K76, where the geometry adjustment dominated. A partial correlation between dominant geometry influence and torsion angle shifts upon the coordination has been observed.

Quantum chemical studies on nucleophilic sites in calcium ion bound zwitterionic calmodulin loops

Quantum chemical calculation on Ca²⁺ bound zwitterionic calmodulin-loops shows terminal capping contributions localized near HOMO and LUMO, which decay exponentially and presence of nucleophilic site at the phenyl-ring O_h of Y99.

A systematic study on RNA NMR chemical shift calculation based on the automated fragmentation QM/MM approach

¹H, ¹³C and ¹⁵N NMR chemical shift calculations on RNAs were performed using the automated fragmentation quantum mechanics/molecular mechanics (AF-QM/MM) approach.

Systematic Study of Locally Dense Basis Sets for NMR Shielding Constants

This paper presents a systematic study of partitioning schemes for locally dense basis sets in the context of NMR shielding calculations. The partitionings explored were based exclusively on connectivity and utilized the basis sets from the pcS-n series. Deviations from pcS-4 shieldings were calculated for a set of 28 organic molecules at the HF, B3LYP, and KT3 levels of theory, with the primary goal being the determination of an efficient scheme that achieves maximal deviations of 0.1 ppm for (1)H and 1 ppm for (13)C. Both atom based and group based divisions of basis sets were examined, with the latter providing the most promising results. It is demonstrated that for the systems studied, at least pcS-1 is required for all parts of the molecule. This, coupled with pcS-3 on the group of interest and pcS-2 on the adjacent groups, is sufficient to achieve the desired level of accuracy at a minimal computational expense. In addition, the suitability of the pcS-n basis sets for post-SCF methods was confirmed through a comparison with other standard basis sets at the MP2 level.

Fragment-Based Electronic Structure Approach for Computing Nuclear Magnetic Resonance Chemical Shifts in Molecular Crystals

First-principles chemical shielding tensor predictions play a critical role in studying molecular crystal structures using nuclear magnetic resonance. Fragment-based electronic structure methods have dramatically improved the ability to model molecular crystal structures and energetics using high-level electronic structure methods. Here, a many-body expansion fragment approach is applied to the calculation of chemical shielding tensors in molecular crystals. First, the impact of truncating the many-body expansion at different orders and the role of electrostatic embedding are examined on a series of molecular clusters extracted from molecular crystals. Second, the ability of these techniques to assign three polymorphic forms of the drug sulfanilamide to the corresponding experimental (13)C spectra is assessed. This challenging example requires discriminating among spectra whose (13)C chemical shifts differ by only a few parts per million (ppm) across the different polymorphs. Fragment-based PBE0/6-311+G(2d,p) level chemical shielding predictions correctly assign these three polymorphs and reproduce the sulfanilamide experimental (13)C chemical shifts with 1 ppm accuracy. The results demonstrate that fragment approaches are competitive with the widely used gauge-invariant projector augmented wave (GIPAW) periodic density functional theory calculations.

Fragment-based (13)C nuclear magnetic resonance chemical shift predictions in molecular crystals: An alternative to planewave methods

We assess the quality of fragment-based ab initio isotropic (13)C chemical shift predictions for a collection of 25 molecular crystals with eight different density functionals. We explore the relative performance of cluster, two-body fragment, combined cluster/fragment, and the planewave gauge-including projector augmented wave (GIPAW) models relative to experiment. When electrostatic embedding is employed to capture many-body polarization effects, the simple and computationally inexpensive two-body fragment model predicts both isotropic (13)C chemical shifts and the chemical shielding tensors as well as both cluster models and the GIPAW approach. Unlike the GIPAW approach, hybrid density functionals can be used readily in a fragment model, and all four hybrid functionals tested here (PBE0, B3LYP, B3PW91, and B97-2) predict chemical shifts in noticeably better agreement with experiment than the four generalized gradient approximation (GGA) functionals considered (PBE, OPBE, BLYP, and BP86). A set of recommended linear regression parameters for mapping between calculated chemical shieldings and observed chemical shifts are provided based on these benchmark calculations. Statistical cross-validation procedures are used to demonstrate the robustness of these fits.

Probing Origin of Binding Difference of inhibitors to MDM2 and MDMX by Polarizable Molecular Dynamics Simulation and QM/MM-GBSA Calculation

Binding abilities of current inhibitors to MDMX are weaker than to MDM2. Polarizable molecular dynamics simulations (MD) followed by Quantum mechanics/molecular mechanics generalized Born surface area (QM//MM-GBSA) calculations were performed to investigate the binding difference of inhibitors to MDM2 and MDMX. The predicted binding free energies not only agree well with the experimental results, but also show that the decrease in van der Walls interactions of inhibitors with MDMX relative to MDM2 is a main factor of weaker bindings of inhibitors to MDMX. The analyses of dihedral angles based on MD trajectories suggest that the closed conformation formed by the residues M53 and Y99 in MDMX leads to a potential steric clash with inhibitors and prevents inhibitors from arriving in the deep of MDMX binding cleft, which reduces the van der Waals contacts of inhibitors with M53, V92, P95 and L98. The calculated results using the residue-based free energy decomposition method further prove that the interaction strength of inhibitors with M53, V92, P95 and L98 from MDMX are obviously reduced compared to MDM2. We expect that this study can provide significant theoretical guidance for designs of potent dual inhibitors to block the p53-MDM2/MDMX interactions.

Probing Difference in Binding Modes of Inhibitors to MDMX by Molecular Dynamics Simulations and Different Free Energy Methods

The p53-MDMX interaction has attracted extensive attention of anti-cancer drug development in recent years. This current work adopted molecular dynamics (MD) simulations and cross-correlation analysis to investigate conformation changes of MDMX caused by inhibitor bindings. The obtained information indicates that the binding cleft of MDMX undergoes a large conformational change and the dynamic behavior of residues obviously change by the presence of different structural inhibitors. Two different methods of binding free energy predictions were employed to carry out a comparable insight into binding mechanisms of four inhibitors PMI, pDI, WK23 and WW8 to MDMX. The data show that the main factor controlling the inhibitor bindings to MDMX arises from van der Waals interactions. The binding free energies were further divided into contribution of each residue and the derived information gives a conclusion that the hydrophobic interactions, such as CH-CH, CH-π and π-π interactions, are responsible for the inhibitor associations with MDMX.

ProCS15: a DFT-based chemical shift predictor for backbone and Cβ atoms in proteins

We present ProCS15: a program that computes the isotropic chemical shielding values of backbone and Cβ atoms given a protein structure in less than a second. ProCS15 is based on around 2.35 million OPBE/6-31G(d,p)//PM6 calculations on tripeptides and small structural models of hydrogen-bonding. The ProCS15-predicted chemical shielding values are compared to experimentally measured chemical shifts for Ubiquitin and the third IgG-binding domain of Protein G through linear regression and yield RMSD values of up to 2.2, 0.7, and 4.8 ppm for carbon, hydrogen, and nitrogen atoms. These RMSD values are very similar to corresponding RMSD values computed using OPBE/6-31G(d,p) for the entire structure for each proteins. These maximum RMSD values can be reduced by using NMR-derived structural ensembles of Ubiquitin. For example, for the largest ensemble the largest RMSD values are 1.7, 0.5, and 3.5 ppm for carbon, hydrogen, and nitrogen. The corresponding RMSD values predicted by several empirical chemical shift predictors range between 0.7–1.1, 0.2–0.4, and 1.8–2.8 ppm for carbon, hydrogen, and nitrogen atoms, respectively.

AFNMR: automated fragmentation quantum mechanical calculation of NMR chemical shifts for biomolecules

We evaluate the performance of the automated fragmentation quantum mechanics/molecular mechanics approach (AF-QM/MM) on the calculation of protein and nucleic acid NMR chemical shifts. The AF-QM/MM approach models solvent effects implicitly through a set of surface charges computed using the Poisson-Boltzmann equation, and it can also be combined with an explicit solvent model through the placement of water molecules in the first solvation shell around the solute; the latter substantially improves the accuracy of chemical shift prediction of protons involved in hydrogen bonding with solvent. We also compare the performance of AF-QM/MM on proteins and nucleic acids with two leading empirical chemical shift prediction programs SHIFTS and SHIFTX2. Although the empirical programs outperform AF-QM/MM in predicting chemical shifts, the differences are in some cases small, and the latter can be applied to chemical shifts on biomolecules which are outside the training set employed by the empirical programs, such as structures containing ligands, metal centers, and non-standard residues. The AF-QM/MM described here is implemented in version 5 of the SHIFTS software, and is fully automated, so that only a structure in PDB format is required as input.