Free Energy Perturbation Calculations with Combined QM/MM Potentials Complications, Simplifications, and Applications to Redox Potential Calculations

Machine Learning Approach to Vertical Energy Gap in Redox Processes

A straightforward approach to calculating the free energy change (ΔG) and reorganization energy of a redox process is linear response approximation (LRA). However, accurate prediction of redox properties is still challenging due to difficulties in conformational sampling and vertical energy-gap sampling. Expensive hybrid quantum mechanical/molecular mechanical (QM/MM) calculations are typically employed in sampling energy gaps using conformations from simulations. To alleviate the computational cost associated with the expensive QM method in the QM/MM calculation, we propose machine learning (ML) methods to predict the vertical energy gaps (VEGs). We tested several ML models to predict the VEGs and observed that simple models like linear regression show excellent performance (mean absolute error ∼0.1 eV) in predicting VEGs in all test systems, even when using features extracted from cheaper semiempirical methods. Our best ML model (extra trees regressor) shows a mean absolute error of around 0.1 eV while using features from the cheapest QM method. We anticipate our approach can be generalized to larger macromolecular systems with more complex redox centers.

Flexibility of Binding Site is Essential to the Ca2+ Selectivity in EF-Hand Calcium-Binding Proteins

High binding affinity and selectivity of metal ions are essential to the function of metalloproteins. Thus, understanding the factors that determine these binding characteristics is of major interest for both fundamental mechanistic investigations and guiding of the design of novel metalloproteins. In this work, we perform QM cluster model calculations and quantum mechanics/molecular mechanics (QM/MM) free energy simulations to understand the binding selectivity of Ca2+ and Mg2+ in the wild-type carp parvalbumin and its mutant. While a nonpolarizable MM model (CHARMM36) does not lead to the correct experimental trend, treatment of the metal binding site with the DFTB3 model in a QM/MM framework leads to relative binding free energies (ΔΔGbind) comparable with experimental data. For the wild-type (WT) protein, the calculated ΔΔGbind is ∼6.6 kcal/mol in comparison with the experimental value of 5.6 kcal/mol. The good agreement highlights the value of a QM description of the metal binding site and supports the role of electronic polarization and charge transfer to metal binding selectivity. For the D51A/E101D/F102W mutant, different binding site models lead to considerable variations in computed binding affinities. With a coordination number of seven for Ca2+, which is shown by QM/MM metadynamics simulations to be the dominant coordination number for the mutant, the calculated relative binding affinity is ∼4.8 kcal/mol, in fair agreement with the experimental value of 1.6 kcal/mol. The WT protein is observed to feature a flexible binding site that accommodates a range of coordination numbers for Ca2+, which is essential to the high binding selectivity for Ca2+ over Mg2+. In the mutant, the E101D mutation reduces the flexibility of the binding site and limits the dominant coordination number of Ca2+ to be seven, thereby leading to reduced binding selectivity against Mg2+. Our results highlight that the binding selectivity of metal ions depends on both the structural and dynamical properties of the protein binding site.

Protein-protein interaction of RdRp with its co-factor NSP8 and NSP7 to decipher the interface hotspot residues for drug targeting: A comparison between SARS-CoV-2 and SARS-CoV

In this study we explored the molecular mechanism of RdRp (Non-Structural Protein, NSP12) interaction with its co-factors NSP7 and NSP8 which is the main toolbox for RNA replication and transcription of SARS-CoV-2 and SARS-CoV. The replication complex is a heterotetramer consists of one NSP12, one NSP7 and two NSP8. Extensive molecular dynamics (MD) simulations were applied on both the heterotetramer complexes to generate the conformations and were used to estimate the MMPBSA binding free energy (BFE) and per-residue energy decomposition of NSP12-NSP8 and NSP12-NSP7 and NSP7-NSP8 complexes. The BFE of SARS-CoV-2 heterotetramer complex with its corresponding partner protein was significantly higher as compared to SARS-CoV. Interface hotspot residues were predicted using different methods implemented in KFC (Knowledge-based FADA and Contracts), HotRegion and Robetta web servers. Per-residue energy decomposition analysis showed that the predicted interface hotspot residues contribute more energy towards the formation of complexes and most of the predicted hotspot residues are clustered together. However, there is a slight difference in the residue-wise energy contribution in the interface NSPs on heterotetramer viral replication complex of both coronaviruses. While the overall replication complex of SARS-CoV-2 was found to be slightly flexible as compared to SARS-CoV. This difference in terms of structural flexibility/stability and energetic characteristics of interface residues including hotspots at PPI interface in the viral replication complexes may be the reason of higher rate of RNA replication of SARS-CoV-2 as compared to SARS-CoV. Overall, the interaction profile at PPI interface such as, interface area, hotspot residues, nature of bonds and energies between NSPs, may provide valuable insights in designing of small molecules or peptide/peptidomimetic ligands which can fit into the PPI interface to disrupt the interaction.

DFT-based QM/MM with Particle-Mesh Ewald for Direct, Long-Range Electrostatic Embedding

We present a DFT-based, QM/MM implementation with long-range electrostatic embedding achieved by direct real-space integration of the particle mesh Ewald (PME) computed electrostatic potential. The key transformation is the interpolation of the electrostatic potential from the PME grid to the DFT quadrature grid, from which integrals are easily evaluated utilizing standard DFT machinery. We provide benchmarks of the numerical accuracy with choice of grid size and real-space corrections, and demonstrate that good convergence is achieved while introducing nominal computational overhead. Furthermore, the approach requires only small modification to existing software packages, as is demonstrated with our implementation in the OpenMM and Psi4 software. After presenting convergence benchmarks, we evaluate the importance of long-range electrostatic embedding in three solute/solvent systems modeled with QM/MM. Water and BMIM/BF4 ionic liquid were considered as ``simple' and ``complex' solvents respectively, with water and p-phenylenediamine (PPD) solute molecules treated at QM level of theory. While electrostatic embedding with standard real-space truncation may introduce negligible error for simple systems such as water solute in water solvent, errors become more significant when QM/MM is applied to complex solvents such as ionic liquids. An extreme example is the electrostatic embedding energy for oxidized PPD in BMIM/BF4 for which real-space truncation produces severe error even at 2-3 nm cutoff distances. This latter example illustrates that utilization of QM/MM to compute redox potentials within concentrated electrolytes/ionic media requires carefully chosen long-range electrostatic embedding algorithms, with our presented algorithm providing a general and robust approach.

Benchmark Study of Redox Potential Calculations for Iron-Sulfur Clusters in Proteins

Redox potentials have been calculated for 12 different iron–sulfur sites of 6 different types with 1–4 iron ions. Structures were optimized with combined quantum mechanical and molecular mechanical (QM/MM) methods, and the redox potentials were calculated using the QM/MM energies, single-point QM methods in a continuum solvent or by QM/MM thermodynamic cycle perturbations. We show that the best results are obtained with a large QM system (∼300 atoms, but a smaller QM system, ∼150 atoms, can be used for the QM/MM geometry optimization) and a large value of the dielectric constant (80). For absolute redox potentials, the B3LYP density functional method gives better results than TPSS, and the results are improved with a larger basis set. However, for relative redox potentials, the opposite is true. The results are insensitive to the force field (charges of the surroundings) used for the QM/MM calculations or whether the protein and solvent outside the QM system are relaxed or kept fixed at the crystal structure. With the best approach for relative potentials, mean absolute and maximum deviations of 0.17 and 0.44 V, respectively, are obtained after removing a systematic error of −0.55 V. Such an approach can be used to identify the correct oxidation states involved in a certain redox reaction.

We have studied redox potentials of 12 iron−sulfur sites of 6 types with 1−4 iron ions. Structures were optimized with combined quantum mechanical and molecular mechanical (QM/MM) methods, and the redox potentials were calculated with QM/MM, QM calculations in a continuum solvent or by QM/MM thermodynamic cycle perturbations. The best results are obtained with the second approach using ∼300 atoms in the QM model and a large dielectric constant.

Theoretical Modeling of Redox Potentials of Biomolecules

The estimation of the redox potentials of biologically relevant systems by means of theoretical-computational approaches still represents a challenge. In fact, the size of these systems typically does not allow a full quantum-mechanical treatment needed to describe electron loss/gain in such a complex environment, where the redox process takes place. Therefore, a number of different theoretical strategies have been developed so far to make the calculation of the redox free energy feasible with current computational resources. In this review, we provide a survey of such theoretical-computational approaches used in this context, highlighting their physical principles and discussing their advantages and limitations. Several examples of these approaches applied to the estimation of the redox potentials of both proteins and nucleic acids are described and critically discussed. Finally, general considerations on the most promising strategies are reported.

BENCHMARKS OF THE DENSITY FUNCTIONAL TIGHT-BINDING METHOD FOR REDOX, PROTONATION AND ELECTRONIC PROPERTIES OF QUINONES

Organic materials with controllable molecular design and sustainable resources are promising electrode materials. Crystalline quinones have been investigated in a variety of rechargeable battery chemistries due to their ubiquitous nature,...

The redox potential of a heme cofactor in Nitrosomonas europaea cytochrome c peroxidase: a polarizable QM/MM study

Redox reactions are crucial to biological processes that protect organisms against oxidative stress. Metalloenzymes, such as peroxidases which reduce excess reactive oxygen species into water, play a key role in detoxification mechanisms. Here we present the results of a polarizable QM/MM study of the reduction potential of the electron transfer heme in the cytochrome c peroxidase of Nitrosomonas europaea. We have found that environment polarization does not substantially affect the computed value of the redox potential. Particular attention has been given to analyzing the role of electrostatic interactions within the protein environment and the solvent on tuning the redox potential of the heme co-factor. We have found that the electrostatic interactions predominantly explain the fluctuations of the vertical ionization/attachment energies of the heme for the sampled configurations, and that the long range electrostatic interactions (up to 40 Å) contribute substantially to the absolute values of the vertical energy gaps.

Quantum Chemistry Calculations for Metabolomics

A primary goal of metabolomics studies is to fully characterize the small-molecule composition of complex biological and environmental samples. However, despite advances in analytical technologies over the past two decades, the majority of small molecules in complex samples are not readily identifiable due to the immense structural and chemical diversity present within the metabolome. Current gold-standard identification methods rely on reference libraries built using authentic chemical materials ("standards"), which are not available for most molecules. Computational quantum chemistry methods, which can be used to calculate chemical properties that are then measured by analytical platforms, offer an alternative route for building reference libraries, i.e., in silico libraries for "standards-free" identification. In this review, we cover the major roadblocks currently facing metabolomics and discuss applications where quantum chemistry calculations offer a solution. Several successful examples for nuclear magnetic resonance spectroscopy, ion mobility spectrometry, infrared spectroscopy, and mass spectrometry methods are reviewed. Finally, we consider current best practices, sources of error, and provide an outlook for quantum chemistry calculations in metabolomics studies. We expect this review will inspire researchers in the field of small-molecule identification to accelerate adoption of in silico methods for generation of reference libraries and to add quantum chemistry calculations as another tool at their disposal to characterize complex samples.

Exploring Coupled Redox and pH Processes with a Force-Field-Based Approach: Applications to Five Different Systems

Coupled redox and pH-driven processes are at the core of many important biological mechanisms. As the distribution of protonation and redox states in a system is associated with the pH and redox potential of the solution, having efficient computational tools that can simulate under these conditions becomes very important. Such tools have the potential to provide information that complement and drive experiments. In previous publications we have presented the implementation of the constant pH and redox potential molecular dynamics (C(pH,E)MD) method in AMBER and we have shown how multidimensional replica exchange can be used to significantly enhance the convergence efficiency of our simulations. In the current work, after an improvement in our C(pH,E)MD approach that allows a given residue to be simultaneously pH- and redox-active, we have employed our methodologies to study five different systems of interest in the literature. We present results for capped tyrosine dipeptide, two maquette systems containing one pH- and redox-active tyrosine (α₃Y and peptide A), and two proteins that contain multiple heme groups (diheme cytochrome c from Rhodobacter sphaeroides and Desulfovibrio vulgaris Hildenborough cytochrome c₃). We show that our results can provide new insights into previous theoretical and experimental findings by using a fully force-field-based and GPU-accelerated approach, which allows the simulations to be executed with high computational performance.

Development of a Robust Indirect Approach for MM → QM Free Energy Calculations That Combines Force-Matched Reference Potential and Bennett's Acceptance Ratio Methods

We use the PBE0/6-31G* density functional method to perform ab initio quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations under periodic boundary conditions with rigorous electrostatics using the ambient potential composite Ewald method in order to test the convergence of MM → QM/MM free energy corrections for the prediction of 17 small-molecule solvation free energies and eight ligand binding free energies to T4 lysozyme. The "indirect" thermodynamic cycle for calculating free energies is used to explore whether a series of reference potentials improve the statistical quality of the predictions. Specifically, we construct a series of reference potentials that optimize a molecular mechanical (MM) force field's parameters to reproduce the ab initio QM/MM forces from a QM/MM simulation. The optimizations form a systematic progression of successively expanded parameters that include bond, angle, dihedral, and charge parameters. For each reference potential, we calculate benchmark quality reference values for the MM → QM/MM correction by performing the mixed MM and QM/MM Hamiltonians at 11 intermediate states, each for 200 ps. We then compare forward and reverse application of Zwanzig's relation, thermodynamic integration (TI), and Bennett's acceptance ratio (BAR) methods as a function of reference potential, simulation time, and the number of simulated intermediate states. We find that Zwanzig's equation is inadequate unless a large number of intermediate states are explicitly simulated. The TI and BAR mean signed errors are very small even when only the end-state simulations are considered, and the standard deviations of the TI and BAR errors are decreased by choosing a reference potential that optimizes the bond and angle parameters. We find a robust approach for the data sets of fairly rigid molecules considered here is to use bond + angle reference potential together with the end-state-only BAR analysis. This requires QM/MM simulations to be performed in order to generate reference data to parametrize the bond + angle reference potential, and then this same simulation serves a dual purpose as the full QM/MM end state. The convergence of the results with respect to time suggests that computational resources may be used more efficiently by running multiple simulations for no more than 50 ps, rather than running one long simulation.

Parameterization of Molybdenum Disulfide Interacting with Water Using the Free Energy Perturbation Method

Water contact angles (WCA) are often used to parametrize force field parameters of novel 2D nanomaterials, such as molybdenum disulfide (MoS₂), which has emerged as a promising nanomaterial in many biomedical applications due to its unique and impressive properties. However, there is a wide range of water-MoS₂ contact angles in the literature depending on the aging process on the surface of a MoS₂ nanosheet and/or substrate material. In this study, we revisit and optimize existing parameters for the basal plane of MoS₂ with two popular water models, TIP3P and SPC/E, using the wide range of WCAs from various experiments. We develop and deploy the free energy perturbation method for parametrizing MoS₂ with experimentally determined WCAs for both fresh and aged surfaces. Energy decomposition analysis on the simulation trajectories reveals that MoS₂-water interaction is dominated by van der Waals interaction, which mainly comes from the top layer of MoS₂. We conclude that to describe both fresh and aged MoS₂ surfaces it is convenient to only adjust the Lennard-Jones parameter ε_S (the depth of the potential well of a sulfur atom), which displays a surprisingly linear correlation with WCAs.

Role of the Photosystem II as an Environment in the Oxidation Free Energy of the Mn Cluster from S1 to S2

The manganese cluster (CaMn₄O₅) in the photosystem II (PSII) is the reaction center of the light-driven oxidation reaction, which generates the molecular oxygen. In this paper, we address the issue of the effect of the environment on the free energy associated with the oxidation of the Mn cluster in S₁ state by conducting the large-scale quantum mechanical/molecular mechanical simulations, which involve the whole of the PSII monomer. It was found by the simulations at the level of the B3LYP functional that the environment surrounding the Mn cluster reduces the vertical oxidation free energy Δμ_vrt by 64.8 kcal/mol. A decomposition analysis of the free energy Δμ_vrt revealed that the system composed of peptide chains, ligands, lipids, and potassium ions contributes to lowering of Δμ_vrt by -98.0 kcal/mol, whereas the solvent water makes an opposite contribution of 38.9 kcal/mol. Reduction of the vertical oxidation free energy directly leads to the lowering of the activation free energy ΔG_ac for the electron transfer reaction from the Mn cluster in S₁ state to the neighboring Tyr_z⁺. Consequently, the electron transfer rate was found to be enhanced by a factor of 10¹² by virtue of the influence of the environment.

Accurate Quantum Chemical Free Energies at Affordable Cost

Free energy sampling methods allow studying the full dynamics of activated processes. Unfortunately, the affordable accuracy of the potential describing the energy and forces of the system is usually rather low. Here we introduce a new method that by combining metadynamics and free energy perturbation allows calculating accurate quantum chemical free energies for chemical reactions. To prove the effectiveness of this new approach we study the S_N2 reaction of CH₃F + Cl^- → CH₃Cl + F^- in vacuo and solvated by water. Comparisons are made with harmonic transition-state theory to show how this method could provide accurate equilibrium and rate constants for complex systems.

Polarizable embedding for simulating redox potentials of biomolecules

Redox reactions play a key role in various biological processes, including photosynthesis and respiration. Quantitative and predictive computational characterization of redox events is therefore highly desirable for enriching our knowledge on mechanistic features of biological redox-active macromolecules. Here, we present a computational protocol exploiting polarizable embedding hybrid quantum-classical approach and resulting in accurate estimates of redox potentials of biological macromolecules. A special attention is paid to fundamental aspects of the theoretical description such as the effects of environment polarization and of the long-range electrostatic interactions on the computed energetic parameters. Environment (protein and the solvent) polarization is shown to be crucial for accurate estimates of the redox potential: hybrid quantum-classical results with and without account for environment polarization differ by 1.4 V. Long-range electrostatic interactions are shown to contribute significantly to the computed redox potential value even at the distances far beyond the protein outer surface. The approach is tested on simulating reduction potential of cryptochrome 1 protein from Arabidopsis thaliana. The theoretical estimate (0.07 V) of the midpoint reduction potential is in good agreement with available experimental data (-0.15 V).

Neutron scattering in the biological sciences: progress and prospects

The scattering of neutrons can be used to provide information on the structure and dynamics of biological systems on multiple length and time scales. Pursuant to a National Science Foundation-funded workshop in February 2018, recent developments in this field are reviewed here, as well as future prospects that can be expected given recent advances in sources, instrumentation and computational power and methods. Crystallography, solution scattering, dynamics, membranes, labeling and imaging are examined. For the extraction of maximum information, the incorporation of judicious specific deuterium labeling, the integration of several types of experiment, and interpretation using high-performance computer simulation models are often found to be particularly powerful.

Accurate Quantum Mechanical/Molecular Mechanical Calculations of Reduction Potentials in Azurin Variants

Understanding the regulation mechanism and molecular determinants of the reduction potential of metalloprotein is a major challenge. An ab initio quantum mechanical/molecular mechanical (QM/MM) method combining the minimum free energy path (MFEP) and fractional number of electron (FNE) approaches has been developed in our group to simulate the redox processes of large systems. The FNE scheme provides an efficient unique description for the redox process, while the MFEP method provides improved conformational sampling on complex environments such as protein in the QM/MM calculations. The reduction potentials of wild-type and seven mutants of azurin, a type 1 copper metalloprotein, were simulated with the QM/MM-MFEP+FNE approach in this paper. A range of 350 mV for the variations of the reduction potentials of these azurin proteins was reproduced faithfully with relative errors around 20 mV. The correlation between structural interactions and reduction potentials observed in simulations provides in-depth insight into the regulation of reduction potentials, which potentially can also be very useful to the engineering of metalloprotein-based electrocatalysts in artificial photosynthesis. The excellent accuracy and efficiency of the QM/MM-MFEP+FNE approach demonstrate the potential for simulations of many electron transfer processes in condensed phases and biochemical systems.

A GPU-Accelerated Parameter Interpolation Thermodynamic Integration Free Energy Method

There has been a resurgence of interest in free energy methods motivated by the performance enhancements offered by molecular dynamics (MD) software written for specialized hardware, such as graphics processing units (GPUs). In this work, we exploit the properties of a parameter-interpolated thermodynamic integration (PI-TI) method to connect states by their molecular mechanical (MM) parameter values. This pathway is shown to be better behaved for Mg2+ → Ca2+ transformations than traditional linear alchemical pathways (with and without soft-core potentials). The PI-TI method has the practical advantage that no modification of the MD code is required to propagate the dynamics, and unlike with linear alchemical mixing, only one electrostatic evaluation is needed (e.g., single call to particle-mesh Ewald) leading to better performance. In the case of AMBER, this enables all the performance benefits of GPU-acceleration to be realized, in addition to unlocking the full spectrum of features available within the MD software, such as Hamiltonian replica exchange (HREM). The TI derivative evaluation can be accomplished efficiently in a post-processing step by reanalyzing the statistically independent trajectory frames in parallel for high throughput. We also show how one can evaluate the particle mesh Ewald contribution to the TI derivative evaluation without needing to perform two reciprocal space calculations. We apply the PI-TI method with HREM on GPUs in AMBER to predict p Ka values in double stranded RNA molecules and make comparison with experiments. Convergence to under 0.25 units for these systems required 100 ns or more of sampling per window and coupling of windows with HREM. We find that MM charges derived from ab initio QM/MM fragment calculations improve the agreement between calculation and experimental results.

Electrostatic Interactions in Protein Structure, Folding, Binding, and Condensation

Charged and polar groups, through forming ion pairs, hydrogen bonds, and other less specific electrostatic interactions, impart important properties to proteins. Modulation of the charges on the amino acids, e.g., by pH and by phosphorylation and dephosphorylation, have significant effects such as protein denaturation and switch-like response of signal transduction networks. This review aims to present a unifying theme among the various effects of protein charges and polar groups. Simple models will be used to illustrate basic ideas about electrostatic interactions in proteins, and these ideas in turn will be used to elucidate the roles of electrostatic interactions in protein structure, folding, binding, condensation, and related biological functions. In particular, we will examine how charged side chains are spatially distributed in various types of proteins and how electrostatic interactions affect thermodynamic and kinetic properties of proteins. Our hope is to capture both important historical developments and recent experimental and theoretical advances in quantifying electrostatic contributions of proteins.

Theoretical estimation of redox potential of biological quinone cofactors

Redox potentials are essential to understand biological cofactor reactivity and to predict their behavior in biological media. Experimental determination of redox potential in biological system is often difficult due to complexity of biological media but computational approaches can be used to estimate them. Nevertheless, the quality of the computational methodology remains a key issue to validate the results. Instead of looking to the best absolute results, we present here the calibration of theoretical redox potential for quinone derivatives in water coupling QM + MM or QM/MM scheme. Our approach allows using low computational cost theoretical level, ideal for long simulations in biological systems, and determination of the uncertainties linked to the calculations. © 2017 Wiley Periodicals, Inc.

Computing converged free energy differences between levels of theory via nonequilibrium work methods: Challenges and opportunities

We demonstrate that Jarzynski's equation can be used to reliably compute free energy differences between low and high level representations of systems. The need for such a calculation arises when employing the so-called "indirect" approach to free energy simulations with mixed quantum mechanical/molecular mechanical (QM/MM) Hamiltonians; a popular technique for circumventing extensive simulations involving quantum chemical computations. We have applied this methodology to several small and medium sized organic molecules, both in the gas phase and explicit solvent. Test cases include several systems for which the standard approach; that is, free energy perturbation between low and high level description, fails to converge. Finally, we identify three major areas in which the difference between low and high level representations make the calculation of ΔAlow→high difficult: bond stretching and angle bending, different preferred conformations, and the response of the MM region to the charge distribution of the QM region. © 2016 Wiley Periodicals, Inc.

Toward polarizable AMOEBA thermodynamics at fixed charge efficiency using a dual force field approach: application to organic crystals

First principles prediction of the structure, thermodynamics and solubility of organic molecular crystals, which play a central role in chemical, material, pharmaceutical and engineering sciences, challenges both potential energy functions and sampling methodologies. Here we calculate absolute crystal deposition thermodynamics using a novel dual force field approach whose goal is to maintain the accuracy of advanced multipole force fields (e.g. the polarizable AMOEBA model) while performing more than 95% of the sampling in an inexpensive fixed charge (FC) force field (e.g. OPLS-AA). Absolute crystal sublimation/deposition phase transition free energies were determined using an alchemical path that grows the crystalline state from a vapor reference state based on sampling with the OPLS-AA force field, followed by dual force field thermodynamic corrections to change between FC and AMOEBA resolutions at both end states (we denote the three step path as AMOEBA/FC). Importantly, whereas the phase transition requires on the order of 200 ns of sampling per compound, only 5 ns of sampling was needed for the dual force field thermodynamic corrections to reach a mean statistical uncertainty of 0.05 kcal mol-1. For five organic compounds, the mean unsigned error between direct use of AMOEBA and the AMOEBA/FC dual force field path was only 0.2 kcal mol-1 and not statistically significant. Compared to experimental deposition thermodynamics, the mean unsigned error for AMOEBA/FC (1.4 kcal mol-1) was more than a factor of two smaller than uncorrected OPLS-AA (3.2 kcal mol-1). Overall, the dual force field thermodynamic corrections reduced condensed phase sampling in the expensive force field by a factor of 40, and may prove useful for protein stability or binding thermodynamics in the future.

Free Energies of Redox Half-Reactions from First-Principles Calculations

Quantitative prediction of the energetics of redox half-reactions is still a challenge for modern computational chemistry. Here, we propose a simple scheme for reliable calculations of vertical ionization and attachment energies, as well as of redox potentials of solvated molecules. The approach exploits linear response approximation in the context of explicit solvent simulations with spherical boundary conditions. It is shown that both vertical ionization energies and vertical electron affinities, and, consequently redox potentials, exhibit linear dependence on the inverse radius of the solvation sphere. The explanation of the linear dependence is provided, and an extrapolation scheme is suggested. The proposed approach accounts for the specific short-range interactions within hybrid DFT and effective fragment potential approach as well as for the asymptotic system-size effects. The computed vertical ionization energies and redox potentials are in excellent agreement with the experimental values.

Accurate Prediction of One-Electron Reduction Potentials in Aqueous Solution by Variable-Temperature H-Atom Addition/Abstraction Methodology

A robust and efficient theoretical approach for calculation of the reduction potentials of charged species in aqueous solution is presented. Within this approach, the reduction potential of a charged complex (with a charge |n| ≥ 2) is probed by means of the reduction potential of its neutralized (protonated/deprotonated) cognate, employing one or several H-atom addition/abstraction thermodynamic cycles. This includes a separation of one-electron reduction from protonation/deprotonation through the temperature dependence. The accuracy of the method has been assessed for the set of 15 transition-metal complexes that are considered as highly challenging systems for computational electrochemistry. Unlike the standard computational protocol(s), the presented approach yields results that are in excellent agreement with experimental electrochemical data. Last but not least, the applicability and limitations of the approach are thoroughly discussed.

Copper Oxidation/Reduction in Water and Protein: Studies with DFTB3/MM and VALBOND Molecular Dynamics Simulations

We apply two recently developed computational methods, DFTB3 and VALBOND, to study copper oxidation/reduction processes in solution and protein. The properties of interest include the coordination structure of copper in different oxidation states in water or in a protein (plastocyanin) active site, the reduction potential of the copper ion in different environments, and the environmental response to copper oxidation. The DFTB3/MM and VALBOND simulation results are compared to DFT/MM simulations and experimental results whenever possible. For a copper ion in aqueous solution, DFTB3/MM results are generally close to B3LYP/MM with a medium basis, including both solvation structure and reduction potential for Cu(II); for Cu(I), however, DFTB3/MM finds a two-water coordination, similar to previous Born-Oppenheimer molecular dynamics simulations using BLYP and HSE, whereas B3LYP/MM leads to a tetrahedron coordination. For a tetraammonia copper complex in aqueous solution, VALBOND and DFTB3/MM are consistent in terms of both structural and dynamical properties of solvent near copper for both oxidation states. For copper reduction in plastocyanin, DFTB3/MM simulations capture the key properties of the active site, and the computed reduction potential and reorganization energy are in fair agreement with experiment, especially when the periodic boundary condition is used. Overall, the study supports the value of VALBOND and DFTB3(/MM) for the analysis of fundamental copper redox chemistry in water and protein, and the results also help highlight areas where further improvements in these methods are desirable.

Practically Efficient QM/MM Alchemical Free Energy Simulations: The Orthogonal Space Random Walk Strategy

The difference between free energy changes occurring at two chemical states can be rigorously estimated via alchemical free energy (AFE) simulations. Traditionally, most AFE simulations are carried out under the classical energy potential treatment; then, accuracy and applicability of AFE simulations are limited. In the present work, we integrate a recent second-order generalized ensemble strategy, the orthogonal space random walk (OSRW) method, into the combined quantum mechanical/molecular mechanical (QM/MM) potential based AFE simulation scheme. Thereby, within a commonly affordable simulation length, accurate QM/MM alchemical free energy simulations can be achieved. As revealed by the model study on the equilibrium of a tautomerization process of hydrated 3-hydroxypyrazole and by the model calculations of the redox potentials of two flavin derivatives, lumichrome (LC) and riboflavin (RF) in aqueous solution, the present OSRW-based scheme could be a viable path toward the realization of practically efficient QM/MM AFE simulations.

A Multiscale Treatment of Angeli's Salt Decomposition

Sodium trioxodinitrate's (Na2N2O3, Angeli's salt) unique cardiovascular effects have been associated with its ability to yield HNO upon dissociation under physiological conditions. Due to its potential applications in new therapies for heart failure, the dissociation of Angeli's salt has recently received increased attention. The decomposition mechanism has been previously studied by quantum mechanical methods using a continuum approximation (PCM) for the solvent effects. In this work we use our recently developed interface of the Amber and Gaussian packages via the PUPIL package to study Angeli's salt dissociation in a hybrid QM/MM scheme where the water solvent molecules are treated explicitly with classical mechanics while the solute is treated with full quantum mechanics (UB3LYP/6-31+G(d) and UMP2/6-31+G(d)) level. Multiple steered molecular dynamics was used with the Jarzynski relationship to extract the free energy profile for the process. We obtain 4.8 kcal mol(-1) and 6.4 kcal mol(-1) free energy barriers for the N-N bond breaking for UB3LYP and UMP2, respectively. The geometries and Mulliken charges for reactant, transition state, and products have been characterized through a number of hybrid QM/MM molecular dynamics runs with the N-N distance restrained to representative values of each species. The results highlight the role of individual solvent molecules for the reaction energetics and provide a comparison point against implicit solvation methods.

QM/MM calculations with deMon2k

The density functional code deMon2k employs a fitted density throughout (Auxiliary Density Functional Theory), which offers a great speed advantage without sacrificing necessary accuracy. Powerful Quantum Mechanical/Molecular Mechanical (QM/MM) approaches are reviewed. Following an overview of the basic features of deMon2k that make it efficient while retaining accuracy, three QM/MM implementations are compared and contrasted. In the first, deMon2k is interfaced with the CHARMM MM code (CHARMM-deMon2k); in the second MM is coded directly within the deMon2k software; and in the third the Chemistry in Ruby (Cuby) wrapper is used to drive the calculations. Cuby is also used in the context of constrained-DFT/MM calculations. Each of these implementations is described briefly; pros and cons are discussed and a few recent applications are described briefly. Applications include solvated ions and biomolecules, polyglutamine peptides important in polyQ neurodegenerative diseases, copper monooxygenases and ultra-rapid electron transfer in cryptochromes.

Effect of stacking interactions on the thermodynamics and kinetics of lumiflavin: a study with improved density functionals and density functional tight-binding protocol

The π-π stacking interaction between lumiflavin and a number of π-electron-rich molecules has been studied by density functional theory using several new-generation density functionals. Six known lumiflavin-aromatic adducts were used and the models were evaluated by comparing the geometry and energetics with experimental results. The study found that dispersion-corrected and hybrid functionals with larger (>50%) Hartree-Fock exchanges produced superior results in modeling thermodynamic characteristics of these complexes. The functional producing the best energetics for these model systems was used to study the stacking interactions of lumiflavin with biologically relevant aromatic groups. Additionally, the reduction of flavin-in the presence of both a hydride donor and a nondonor π-electronic system was also studied. Weak interactions were observed in the stacked lumiflavin complexes of benzene, phenol, and indole, mimicking phenyl alanine, tryptophan, and tyrosine side chains, respectively, of an enzyme. The stacked complex of naphthalene and flavin showed little change in flavin's redox potential indicating insignificant effect on the thermodynamics of the hydride transfer reaction. In contrast, the hydride transfer reaction with the hydride donor N-methyl nicotinamide tells a different story, as the transition state was found to be strongly impacted by the stacking interactions. A comparison of performance between the density functional theory (DFT) and the computationally less expensive dispersion-corrected self-consistent density functional tight-binding (SCC-DFTB-D) theory revealed that the latter produces consistent energetics for this hydride transfer reaction and additional DFT-computed perturbative corrections could significantly improve these results.

Importance of polarization effect in the study of metalloproteins: application of polarized protein specific charge scheme in predicting the reduction potential of azurin

Molecular dynamics (MD) simulation is commonly used in the study of protein dynamics, and in recent years, the extension of MD simulation to the study of metalloproteins is gaining much interest. Choice of force field is crucial in MD studies, and the inclusion of metal centers complicates the process of accurately describing the electrostatic environment that surrounds the redox centre. Herein, we would like to explore the importance of including electrostatic contribution from both protein and solvent in the study of metalloproteins. MD simulations with the implementation of thermodynamic integration will be conducted to model the reduction process of azurin from Pseudomonas aeruginosa. Three charge schemes will be used to derive the partial charges of azurin. These charge schemes differ in terms of the amount of immediate environment, respective to copper, considered during charge fitting, which ranges from the inclusion of copper and residues in the first coordination sphere during density functional theory charge fitting to the comprehensive inclusion of protein and solvent effect surrounding the metal centre using polarized protein-specific charge scheme. From the simulations conducted, the relative reduction potential of the mutated azurins respective to that of wild-type azurin (ΔEcal) were calculated and compared with experimental values. The ΔEcal approached experimental value with increasing consideration of environmental effect hence substantiating the importance of polarization effect in the study of metalloproteins. This study also attests the practicality of polarized protein-specific charge as a computational tool capable of incorporating both protein environment and solvent effect into MD simulations.

Multiscale Free Energy Simulations: An Efficient Method for Connecting Classical MD Simulations to QM or QM/MM Free Energies Using Non-Boltzmann Bennett Reweighting Schemes

The reliability of free energy simulations (FES) is limited by two factors: (a) the need for correct sampling and (b) the accuracy of the computational method employed. Classical methods (e.g., force fields) are typically used for FES and present a myriad of challenges, with parametrization being a principle one. On the other hand, parameter-free quantum mechanical (QM) methods tend to be too computationally expensive for adequate sampling. One widely used approach is a combination of methods, where the free energy difference between the two end states is computed by, e.g., molecular mechanics (MM), and the end states are corrected by more accurate methods, such as QM or hybrid QM/MM techniques. Here we report two new approaches that significantly improve the aforementioned scheme; with a focus on how to compute corrections between, e.g., the MM and the more accurate QM calculations. First, a molecular dynamics trajectory that properly samples relevant conformational degrees of freedom is generated. Next, potential energies of each trajectory frame are generated with a QM or QM/MM Hamiltonian. Free energy differences are then calculated based on the QM or QM/MM energies using either a non-Boltzmann Bennett approach (QM-NBB) or non-Boltzmann free energy perturbation (NB-FEP). Both approaches are applied to calculate relative and absolute solvation free energies in explicit and implicit solvent environments. Solvation free energy differences (relative and absolute) between ethane and methanol in explicit solvent are used as the initial test case for QM-NBB. Next, implicit solvent methods are employed in conjunction with both QM-NBB and NB-FEP to compute absolute solvation free energies for 21 compounds. These compounds range from small molecules such as ethane and methanol to fairly large, flexible solutes, such as triacetyl glycerol. Several technical aspects were investigated. Ultimately some best practices are suggested for improving methods that seek to connect MM to QM (or QM/MM) levels of theory in FES.

Charging free energy calculations using the Generalized Solvent Boundary Potential (GSBP) and periodic boundary condition: a comparative analysis using ion solvation and oxidation free energy in proteins

Free energy simulations using a finite sphere boundary condition rather than a periodic boundary condition (PBC) are attractive in the study of very large biomolecular systems. To understand the quantitative impact of various approximations in such simulations, we compare charging free energies in both solution and protein systems calculated in a linear response framework with the Generalized Solvent Boundary Potential (GSBP) and PBC simulations. For simple ions in solution, we find good agreements between GSBP and PBC charging free energies, once the relevant correction terms are taken into consideration. For PBC simulations with the particle-mesh-Ewald for long-range electrostatics, the contribution (ΔG(P-M)) due to the use of a particle rather than molecule based summation scheme in real space is found to be significant, as pointed out by Hünenberger and co-workers. For GSBP, when the inner region is close to be charge neutral, the key correction is the overpolarization of water molecules at the inner/outer dielectric boundary; the magnitude of the correction (ΔG(s-pol)), however, is relatively small. For charging (oxidation) free energy in proteins, the situation is more complex, although good agreement between GSBP and PBC can still be obtained when care is exercised. The smooth dielectric boundary approximation inherent to GSBP tends to make significant errors when the inner region is featured with a high net charge. However, the error can be corrected with Poisson-Boltzmann calculations using snapshots from GSBP simulations in a straightforward and robust manner. Because of the more complex charge and solvent distributions, the magnitudes of ΔG(P-M) and ΔG(s-pol) in protein simulations appear to be different from those derived for solution simulations, leading to uncertainty in directly comparing absolute charging free energies from PBC and GSBP simulations for protein systems. The relative charging/oxidation free energies, however, are robust. With the linear response approximation, for the specific protein system (CueR) studied, the effect of freezing the protein structure in the outer region is found to be small, unless a very small (8 Å) inner region is used; even in the latter case, the result is substantially improved when the nearby metal binding loop is allowed to respond to metal oxidation. The implications of these results to the applicability of GSBP to complex biomolecules and in ab initio QM/MM simulations are discussed.

Fragment-based Quantum Mechanical/Molecular Mechanical Simulations of Thermodynamic and Kinetic Process of the Ru2+-Ru3+ Self-Exchange Electron Transfer

A fragment-based fractional number of electron (FNE) approach, is developed to study entire electron transfer (ET) processes from the electron donor region to the acceptor region in condensed phase. Both regions are described by the density-fragment interaction (DFI) method while FNE as an efficient ET order parameter is applied to simulate the electron transfer process. In association with the QM/MM energy expression, the DFI-FNE method is demonstrated to describe ET processes robustly with the Ru2+-Ru3+ self-exchange ET as a proof-of-concept example. This method allows for systematic calculations of redox free energies, reorganization energies, and electronic couplings, and the absolute ET rate constants within the Marcus regime.

Expanding the view of proton pumping in cytochrome c oxidase through computer simulation

In cytochrome c oxidase (CcO), a redox-driven proton pump, protons are transported by the Grotthuss shuttling via hydrogen-bonded water molecules and protonatable residues. Proton transport through the D-pathway is a complicated process that is highly sensitive to alterations in the amino acids or the solvation structure in the channel, both of which can inhibit proton pumping and enzymatic activity. Simulations of proton transport in the hydrophobic cavity showed a clear redox state dependence. To study the mechanism of proton pumping in CcO, multi-state empirical valence bond (MS-EVB) simulations have been conducted, focusing on the proton transport through the D-pathway and the hydrophobic cavity next to the binuclear center. The hydration structures, transport pathways, effects of residues, and free energy surfaces of proton transport were revealed in these MS-EVB simulations. The mechanistic insight gained from them is herein reviewed and placed in context for future studies.

Reduction potentials and acidity constants of Mn superoxide dismutase calculated by QM/MM free-energy methods

We used two theoretical methods to estimate reduction potentials and acidity constants in Mn superoxide dismutase (MnSOD), namely combined quantum mechanical and molecular mechanics (QM/MM) thermodynamic cycle perturbation (QTCP) and the QM/MM-PBSA approach. In the latter, QM/MM energies are combined with continuum solvation energies calculated by solving the Poisson-Boltzmann equation (PB) or by the generalised Born approach (GB) and non-polar solvation energies calculated from the solvent-exposed surface area. We show that using the QTCP method, we can obtain accurate and precise estimates of the proton-coupled reduction potential for MnSOD, 0.30±0.01 V, which compares favourably with experimental estimates of 0.26-0.40 V. However, the calculated potentials depend strongly on the DFT functional used: The B3LYP functional gives 0.6 V more positive potentials than the PBE functional. The QM/MM-PBSA approach leads to somewhat too high reduction potentials for the coupled reaction and the results depend on the solvation model used. For reactions involving a change in the net charge of the metal site, the corresponding results differ by up to 1.3 V or 24 pK(a) units, rendering the QM/MM-PBSA method useless to determine absolute potentials. However, it may still be useful to estimate relative shifts, although the QTCP method is expected to be more accurate.

Proton storage site in bacteriorhodopsin: new insights from quantum mechanics/molecular mechanics simulations of microscopic pK(a) and infrared spectra

Identifying the group that acts as the proton storage/loading site is a challenging but important problem for understanding the mechanism of proton pumping in biomolecular proton pumps, such as bacteriorhodopsin (bR) and cytochrome c oxidase. Recent experimental studies of bR propelled the idea that the proton storage/release group (PRG) in bR is not an amino acid but a water cluster embedded in the protein. We argue that this idea is at odds with our knowledge of protein electrostatics, since invoking the water cluster as the PRG would require the protein to raise the pK(a) of a hydronium by almost 11 pK(a) units, which is difficult considering known cases of pK(a) shifts in proteins. Our recent quantum mechanics/molecular mechanics (QM/MM) simulations suggested an alternative "intermolecular proton bond" model in which the stored proton is shared between two conserved Glu residues (194 and 204). Here we show that this model leads to microscopic pK(a) values consistent with available experimental data and the functional requirement of a PRG. Extensive QM/MM simulations also show that, independent of a number of technical issues, such as the influence of QM region size, starting X-ray structure, and nuclear quantum effects, the "intermolecular proton bond" model is qualitatively consistent with available spectroscopic data. Potential of mean force calculations show explicitly that the stored proton strongly prefers the pair of Glu residues over the water cluster. The results and analyses help highlight the importance of considering protein electrostatics and provide arguments for why the "intermolecular proton bond" model is likely applicable to the PRG in biomolecular proton pumps in general.

Application of the SCC-DFTB method to neutral and protonated water clusters and bulk water

The self-consistent charge density functional tight-binding (SCC-DFTB) method has been actively employed to study proton-transfer processes in biological systems. Recent studies in the literature employing SCC-DFTB reported that the method favors the Zundel form of the hydrated proton over the Eigen form, both in gas-phase water clusters and in bulk water, in disagreement with both higher-level calculations and experimental data. In this work, we explore the performance of SCC-DFTB for protonated gas-phase water clusters and bulk water (the latter both with and without an excess proton) with a modified O-H repulsive potential reported in our earlier work and with on-site third-order expansion of the DFT energy. Our results show that, with the proper set of published parameters, SCC-DFTB does correctly favor the Eigen form of the hydrated proton as compared to the Zundel form, both in gas-phase clusters and in the bulk; the amphiphilic character of the hydrated proton discussed in the literature has also been observed. The analyses do, however, bring forth remaining limitations in terms of the solvation structure around the hydrated proton as well as the structure of bulk water, which can guide future improvements of the method.

Redox potentials and pKa for benzoquinone from density functional theory based molecular dynamics

The density functional theory based molecular dynamics (DFTMD) method for the computation of redox free energies presented in previous publications and the more recent modification for computation of acidity constants are reviewed. The method uses a half reaction scheme based on reversible insertion/removal of electrons and protons. The proton insertion is assisted by restraining potentials acting as chaperones. The procedure for relating the calculated deprotonation free energies to Brønsted acidities (pK(a)) and the oxidation free energies to electrode potentials with respect to the normal hydrogen electrode is discussed in some detail. The method is validated in an application to the reduction of aqueous 1,4-benzoquinone. The conversion of hydroquinone to quinone can take place via a number of alternative pathways consisting of combinations of acid dissociations, oxidations, or dehydrogenations. The free energy changes of all elementary steps (ten in total) are computed. The accuracy of the calculations is assessed by comparing the energies of different pathways for the same reaction (Hess's law) and by comparison to experiment. This two-sided test enables us to separate the errors related with the restrictions on length and time scales accessible to DFTMD from the errors introduced by the DFT approximation. It is found that the DFT approximation is the main source of error for oxidation free energies.