Efficient implementation of constant pH molecular dynamics on modern graphics processors.

The treatment of pH sensitive ionization states for titratable residues in proteins is often omitted from molecular dynamics (MD) simulations. While static charge models can answer many questions regarding protein conformational equilibrium and protein-ligand interactions, pH-sensitive phenomena such as acid-activated chaperones and amyloidogenic protein aggregation are inaccessible to such models. Constant pH molecular dynamics (CPHMD) coupled with the Generalized Born with a Simple sWitching function (GBSW) implicit solvent model provide an accurate framework for simulating pH sensitive processes in biological systems. Although this combination has demonstrated success in predicting pKa values of protein structures, and in exploring dynamics of ionizable side-chains, its speed has been an impediment to routine application. The recent availability of low-cost graphics processing unit (GPU) chipsets with thousands of processing cores, together with the implementation of the accurate GBSW implicit solvent model on those chipsets (Arthur and Brooks, J. Comput. Chem. 2016, 37, 927), provide an opportunity to improve the speed of CPHMD and ionization modeling greatly. Here, we present a first implementation of GPU-enabled CPHMD within the CHARMM-OpenMM simulation package interface. Depending on the system size and nonbonded force cutoff parameters, we find speed increases of between one and three orders of magnitude. Additionally, the algorithm scales better with system size than the CPU-based algorithm, thus allowing for larger systems to be modeled in a cost effective manner. We anticipate that the improved performance of this methodology will open the door for broad-spread application of CPHMD in its modeling pH-mediated biological processes. © 2016 Wiley Periodicals, Inc.

Parallelization and improvements of the generalized born model with a simple sWitching function for modern graphics processors.

Two fundamental challenges of simulating biologically relevant systems are the rapid calculation of the energy of solvation and the trajectory length of a given simulation. The Generalized Born model with a Simple sWitching function (GBSW) addresses these issues by using an efficient approximation of Poisson-Boltzmann (PB) theory to calculate each solute atom's free energy of solvation, the gradient of this potential, and the subsequent forces of solvation without the need for explicit solvent molecules. This study presents a parallel refactoring of the original GBSW algorithm and its implementation on newly available, low cost graphics chips with thousands of processing cores. Depending on the system size and nonbonded force cutoffs, the new GBSW algorithm offers speed increases of between one and two orders of magnitude over previous implementations while maintaining similar levels of accuracy. We find that much of the algorithm scales linearly with an increase of system size, which makes this water model cost effective for solvating large systems. Additionally, we utilize our GPU-accelerated GBSW model to fold the model system chignolin, and in doing so we demonstrate that these speed enhancements now make accessible folding studies of peptides and potentially small proteins.

Crowding induced collective hydration of biological macromolecules over extended distances.

Ultrafast two-dimensional infrared (2D-IR) spectroscopy reveals picosecond protein and hydration dynamics of crowded hen egg white lysozyme (HEWL) labeled with a metal-carbonyl vibrational probe covalently attached to a solvent accessible His residue. HEWL is systematically crowded alternatively with polyethylene glycol (PEG) or excess lysozyme in order to distinguish the chemically inert polymer from the complex electrostatic profile of the protein crowder. The results are threefold: (1) A sharp dynamical jamming-like transition is observed in the picosecond protein and hydration dynamics that is attributed to an independent-to-collective hydration transition induced by macromolecular crowding that slows the hydration dynamics up to an order of magnitude relative to bulk water. (2) The interprotein distance at which the transition occurs suggests collective hydration of proteins over distances of 30-40 Å. (3) Comparing the crowding effects of PEG400 to our previously reported experiments using glycerol exposes fundamental differences between small and macromolecular crowding agents.

Predicting extreme pKa shifts in staphylococcal nuclease mutants with constant pH molecular dynamics.

Accurate computational methods of determining protein and nucleic acid pK(a) values are vital to understanding pH-dependent processes in biological systems. In this article, we use the recently developed method constant pH molecular dynamics (CPHMD) to explore the calculation of highly perturbed pK(a) values in variants of staphylococcal nuclease (SNase). Simulations were performed using the replica exchange (REX) protocol for improved conformational sampling with eight temperature windows, and yielded converged proton populations in a total sampling time of 4 ns. Our REX-CPHMD simulations resulted in calculated pK(a) values with an average unsigned error (AUE) of 0.75 pK units for the acidic residues in Δ + PHS, a hyperstable variant of SNase. For highly pK(a)-perturbed SNase mutants with known crystal structures, our calculations yielded an AUE of 1.5 pK units and for those mutants based on modeled structures an AUE of 1.4 pK units was found. Although a systematic underestimate of pK shifts was observed in most of the cases for the highly perturbed pK mutants, correlations between conformational rearrangement and plasticity associated with the mutation and error in pK(a) prediction was not evident in the data. This study further extends the scope of electrostatic environments explored using the REX-CPHMD methodology and suggests that it is a reliable tool for rapidly characterizing ionizable amino acids within proteins even when modeled structures are employed.

Substituted imidazole of 5-fluoro-2-[4-[(2-phenyl-1H-imidazol-5-yl)methyl]-1-piperazinyl]pyrimidine Inactivates cytochrome P450 2D6 by protein adduction.

5-Fluoro-2-[4-[(2-phenyl-1H-imidazol-5-yl)methyl]-1-piperazinyl]pyrimidine (SCH 66712) is a potent mechanism-based inactivator of human cytochrome P450 2D6 that displays type I binding spectra with a K(s) of 0.39 ± 0.10 µM. The partition ratio is ~3, indicating potent inactivation that addition of exogenous nucleophiles does not prevent. Within 15 min of incubation with SCH 66712 and NADPH, ∼90% of CYP2D6 activity is lost with only ~20% loss in ability to bind CO and ~25% loss of native heme over the same time. The stoichiometry of binding to the protein was 1.2:1. SDS-polyacrylamide gel electrophoresis with Western blotting and autoradiography analyses of CYP2D6 after incubations with radiolabeled SCH 66712 further support the presence of a protein adduct. Metabolites of SCH 66712 detected by mass spectrometry indicate that the phenyl group on the imidazole ring of SCH 66712 is one site of oxidation by CYP2D6 and could lead to methylene quinone formation. Three other metabolites were also observed. For understanding the metabolic pathway that leads to CYP2D6 inactivation, metabolism studies with CYP2C9 and CYP2C19 were performed because neither of these enzymes is significantly inhibited by SCH 66712. The metabolites formed by CYP2C9 and CYP2C19 are the same as those seen with CYP2D6, although in different abundance. Modeling studies with CYP2D6 revealed potential roles of various active site residues in the oxidation of SCH 66712 and inactivation of CYP2D6 and showed that the phenyl group of SCH 66712 is positioned at 2.2 Å from the heme iron.

Predicting Ligand-Dissociation Energies of 3d Coordination Complexes with Auxiliary-Field Quantum Monte Carlo.

Transition-metal complexes are ubiquitous in biology and chemical catalysis, yet they remain difficult to accurately describe with ab initio methods because of the presence of a large degree of dynamic electron correlation, and, in some cases, strong static correlation which results from a manifold of low-lying states. Progress has been hindered by a scarcity of high-quality gas-phase experimental data, while exact ab initio predictions are usually computationally unaffordable because of the large size of the relevant complexes. In this work, we present a data set of 34 tetrahedral, square planar, and octahedral 3d metal-containing complexes with gas-phase ligand-dissociation energies that have reported uncertainties of ≤2 kcal/mol. We perform all-electron phaseless auxiliary-field quantum Monte Carlo (ph-AFQMC) calculations utilizing multideterminant trial wave functions selected by a black box procedure. We compare the results with those from the density functional theory (DFT) with the B3LYP, B97, M06, PBE0, ωB97X-V, and DSD-PBEP86/2013 functionals and a localized orbital variant of the coupled cluster theory with single, double, and perturbative triple excitations (DLPNO-CCSD(T)). We find mean averaged errors of 1.07 ± 0.27 kcal/mol for our most sophisticated ph-AFQMC approach versus 2.81 kcal/mol for DLPNO-CCSD(T) and 1.49-3.78 kcal/mol for DFT. We find maximum errors of 2.96 ± 1.71 kcal/mol for our best ph-AFQMC method versus 9.15 kcal/mol for DLPNO-CCSD(T) and 5.98-13.69 kcal/mol for DFT. The reasonable performance of a number of DFT functionals is in stark contrast to the much poorer accuracy previously demonstrated for diatomic species, suggesting a moderation in electron correlation because of ligand coordination in most cases. However, the unpredictably large errors for a small subset of cases with both DFT and DLPNO-CCSD(T) methods leave cause for concern, especially in light of the unreliability of common multireference indicators. In contrast, the robust and, in principle, systematically improvable results of ph-AFQMC for these realistic complexes establish the method as a useful tool for elucidating the electronic structure of transition-metal-containing complexes and predicting their gas-phase properties.

In silico prediction of annihilators for triplet-triplet annihilation upconversion via auxiliary-field quantum Monte Carlo.

The energy of the lowest-lying triplet state (T1) relative to the ground and first-excited singlet states (S0, S1) plays a critical role in optical multiexcitonic processes of organic chromophores. Focusing on triplet-triplet annihilation (TTA) upconversion, the S0 to T1 energy gap, known as the triplet energy, is difficult to measure experimentally for most molecules of interest. Ab initio predictions can provide a useful alternative, however low-scaling electronic structure methods such as the Kohn-Sham and time-dependent variants of Density Functional Theory (DFT) rely heavily on the fraction of exact exchange chosen for a given functional, and tend to be unreliable when strong electronic correlation is present. Here, we use auxiliary-field quantum Monte Carlo (AFQMC), a scalable electronic structure method capable of accurately describing even strongly correlated molecules, to predict the triplet energies for a series of candidate annihilators for TTA upconversion, including 9,10 substituted anthracenes and substituted benzothiadiazole (BTD) and benzoselenodiazole (BSeD) compounds. We compare our results to predictions from a number of commonly used DFT functionals, as well as DLPNO-CCSD(T0), a localized approximation to coupled cluster with singles, doubles, and perturbative triples. Together with S1 estimates from absorption/emission spectra, which are well-reproduced by TD-DFT calculations employing the range-corrected hybrid functional CAM-B3LYP, we provide predictions regarding the thermodynamic feasibility of upconversion by requiring (a) the measured T1 of the sensitizer exceeds that of the calculated T1 of the candidate annihilator, and (b) twice the T1 of the annihilator exceeds its S1 energetic value. We demonstrate a successful example of in silico discovery of a novel annihilator, phenyl-substituted BTD, and present experimental validation via low temperature phosphorescence and the presence of upconverted blue light emission when coupled to a platinum octaethylporphyrin (PtOEP) sensitizer. The BTD framework thus represents a new class of annihilators for TTA upconversion. Its chemical functionalization, guided by the computational tools utilized herein, provides a promising route towards high energy (violet to near-UV) emission.

Singlet-Triplet Energy Gaps of Organic Biradicals and Polyacenes with Auxiliary-Field Quantum Monte Carlo.

The energy gap between the lowest-lying singlet and triplet states is an important quantity in chemical photocatalysis, with relevant applications ranging from triplet fusion in optical upconversion to the design of organic light-emitting devices. The ab initio prediction of singlet-triplet (ST) gaps is challenging due to the potentially biradical nature of the involved states, combined with the potentially large size of relevant molecules. In this work, we show that phaseless auxiliary-field quantum Monte Carlo (ph-AFQMC) can accurately predict ST gaps for chemical systems with singlet states of highly biradical nature, including a set of 13 small molecules and the ortho-, meta-, and para- isomers of benzyne. With respect to gas-phase experiments, ph-AFQMC using CASSCF trial wave functions achieves a mean averaged error of ∼1 kcal/mol. Furthermore, we find that in the context of a spin-projection technique, ph-AFQMC using unrestricted single-determinant trial wave functions, which can be readily obtained for even very large systems, produces equivalently high accuracy. We proceed to show that this scalable methodology is capable of yielding accurate ST gaps for all linear polyacenes for which experimental measurements exist, that is, naphthalene, anthracene, tetracene, and pentacene. Our results suggest a protocol for selecting either unrestricted Hartree-Fock or Kohn-Sham orbitals for the single-determinant trial wave function, based on the extent of spin-contamination. These findings pave the way for future investigations of specific photochemical processes involving large molecules with substantial biradical character.

On Achieving High Accuracy in Quantum Chemical Calculations of 3 d Transition Metal-Containing Systems: A Comparison of Auxiliary-Field Quantum Monte Carlo with Coupled Cluster, Density Functional Theory, and Experiment for Diatomic Molecules.

The bond dissociation energies of a set of 44 3 d transition metal-containing diatomics are computed with phaseless auxiliary-field quantum Monte Carlo (ph-AFQMC) utilizing a correlated sampling technique. We investigate molecules with H, N, O, F, Cl, and S ligands, including those in the 3dMLBE20 database first compiled by Truhlar and co-workers with calculated and experimental values that have since been revised by various groups. In order to make a direct comparison of the accuracy of our ph-AFQMC calculations with previously published results from 10 DFT functionals, CCSD(T), and icMR-CCSD(T), we establish an objective selection protocol which utilizes the most recent experimental results except for a few cases with well-specified discrepancies. With the remaining set of 41 molecules, we find that ph-AFQMC gives robust agreement with experiment superior to that of all other methods, with a mean absolute error (MAE) of 1.4(4) kcal/mol and maximum error of 3(3) kcal/mol (parentheses account for reported experimental uncertainties and the statistical errors of our ph-AFQMC calculations). In comparison, CCSD(T) and B97, the best performing DFT functional considered here, have MAEs of 2.8 and 3.7 kcal/mol, respectively, and maximum errors in excess of 17 kcal/mol (for the CoS diatomic). While a larger and more diverse data set would be required to demonstrate that ph-AFQMC is truly a benchmark method for transition metal systems, our results indicate that the method has tremendous potential, exhibiting unprecedented consistency and accuracy compared to other approximate quantum chemical approaches.

Phaseless Auxiliary-Field Quantum Monte Carlo on Graphical Processing Units.

We present an implementation of phaseless Auxiliary-Field Quantum Monte Carlo (ph-AFQMC) utilizing graphical processing units (GPUs). The AFQMC method is recast in terms of matrix operations which are spread across thousands of processing cores and are executed in batches using custom Compute Unified Device Architecture kernels and the GPU-optimized cuBLAS matrix library. Algorithmic advances include a batched Sherman-Morrison-Woodbury algorithm to quickly update matrix determinants and inverses, density-fitting of the two-electron integrals, an energy algorithm involving a high-dimensional precomputed tensor, and the use of single-precision floating point arithmetic. These strategies accelerate ph-AFQMC calculations with both single- and multideterminant trial wave functions, though particularly dramatic wall-time reductions are achieved for the latter. For typical calculations we find speed-ups of roughly 2 orders of magnitude using just a single GPU card compared to a single modern CPU core. Furthermore, we achieve near-unity parallel efficiency using 8 GPU cards on a single node and can reach moderate system sizes via a local memory-slicing approach. We illustrate the robustness of our implementation on hydrogen chains of increasing length and through the calculation of all-electron ionization potentials of the first-row transition metal atoms. We compare long imaginary-time calculations utilizing a population control algorithm with our previously published correlated sampling approach and show that the latter improves not only the efficiency but also the accuracy of the computed ionization potentials. Taken together, the GPU implementation combined with correlated sampling provides a compelling computational method that will broaden the application of ph-AFQMC to the description of realistic correlated electronic systems.

Ultrafast 2D-IR and simulation investigations of preferential solvation and cosolvent exchange dynamics.

Using a derivative of the vitamin biotin labeled with a transition-metal carbonyl vibrational probe in a series of aqueous N,N-dimethylformamide (DMF) solutions, we observe a striking slowdown in spectral diffusion dynamics with decreased DMF concentration. Equilibrium solvation dynamics, measured with the rapidly acquired spectral diffusion (RASD) technique, a variant of heterodyne-detected photon-echo peak shift experiments, range from 1 ps in neat DMF to ∼3 ps in 0.07 mole fraction DMF/water solution. Molecular dynamics simulations of the biotin-metal carbonyl solute in explicit aqueous DMF solutions show marked preferential solvation by DMF, which becomes more pronounced at lower DMF concentrations. The simulations and the experimental data are consistent with an interpretation where the slowdown in spectral diffusion is due to solvent exchange involving distinct cosolvent species. A simple two-component model reproduces the observed spectral dynamics as well as the DMF concentration dependence, enabling the extraction of the solvent exchange time scale of 8 ps. This time scale corresponds to the diffusive motion of a few Å, consistent with a solvent-exchange mechanism. Unlike most previous studies of solvation dynamics in binary mixtures of polar solvents, our work highlights the ability of vibrational probes to sense solvent exchange as a new, slow component in the spectral diffusion dynamics.

Heterogeneous preferential solvation of water and trifluoroethanol in homologous lysozymes.

Cytoplasmic osmolytes can significantly alter the thermodynamic and kinetic properties of proteins relative to those under dilute solution conditions. Spectroscopic experiments of lysozymes in cosolvents indicate that such changes may arise from the heterogeneous, site-specific hydrophobic interactions between protein surface residues and individual solvent molecules. In pursuit of an accurate and predictive model for explaining biomolecular interactions, we study the averaged structural characteristics of mixed solvents with homologous lysozyme solutes using all-atom molecular dynamics. By observing the time-averaged densities of different aqueous solutions of trifluoroethanol, we deduce trends in the heterogeneous solvent interactions over each protein's surface, and investigate how the homology of protein structure does not necessarily translate to similarities in solvent structure and composition-even when observing identical side chains.

Site-specific hydration dynamics of globular proteins and the role of constrained water in solvent exchange with amphiphilic cosolvents.

The thermodynamic driving forces for protein folding, association, and function are often determined by protein-water interactions. With a novel covalently bound labeling approach, we have used sensitive vibrational probes, site-selectively conjugated to two lysozyme variants-in conjunction with ultrafast two-dimensional infrared (2D-IR) spectroscopy-to investigate directly the protein-water interface. By probing alternatively a topologically flat, rigid domain and a flexible domain, we find direct experimental evidence for spatially heterogeneous hydration dynamics. The hydration environment around globular proteins can vary from exhibiting bulk-like hydration dynamics to dynamically constrained water, which results from stifled hydrogen bond switching dynamics near extended hydrophobic surfaces. Furthermore, we leverage preferential solvation exchange to demonstrate that the liberation of dynamically constrained water is a sufficient driving force for protein-surface association reactions. These results provide an intuitive picture of the dynamic aspects of hydrophobic hydration of proteins, illustrating an essential function of water in biological processes.

Molecular analysis and modeling of inactivation of human CYP2D6 by four mechanism based inactivators.

Human cytochrome P450 2D6 (CYP2D6) is involved in metabolism of approximately 25% of pharmaceutical drugs. Inactivation of CYP2D6 can lead to adverse drug interactions. Four inactivators of CYP2D6 have previously been identified: 5-Fluoro-2-[4-[(2-phenyl-1H-imidazol-5-yl)methyl]-1-piperazinyl]pyrimidine(SCH66712), (1-[(2-ethyl- 4-methyl-1H-imidazol-5-yl)-methyl]-4-[4-(trifluoromethyl)-2-pyridinyl]piperazine(EMTPP), paroxetine, and 3,4- methylenedioxymethamphetamine (MDMA). All four contain planar, aromatic groups as well as basic nitrogens common to CYP2D6 substrates. SCH66712 and EMTPP also contain piperazine groups and substituted imidazole rings that are common in pharmaceutical agents, though neither of these compounds is clinically relevant. Paroxetine and MDMA contain methylenedioxyphenyls. SCH66712 and EMTPP are both known protein adductors while paroxetine and MDMA are probable heme modifiers. The current study shows that each inactivator displays Type I binding with Ks values that vary by 2-orders of magnitude with lower Ks values associated with greater inactivation. Comparison of KI, kinact, and partition ratio values shows SCH66712 is the most potent inactivator. Molecular modeling experiments using AutoDock identify Phe120 as a key interaction for all four inactivators with face-to-face and edge-to-face pi interactions apparent. Distance between the ligand and heme iron correlates with potency of inhibition. Ligand conformations were scored according to their binding energies as calculated by AutoDock and correlation was observed between molecular models and Ks values.