Explicit Solvation Shell Model and Continuum Solvation Models for Solvation Energy and pKa Determination of Amino Acids

Explainable Supervised Machine Learning Model To Predict Solvation Gibbs Energy

Many challenges persist in developing accurate computational models for predicting solvation free energy (ΔGsol). Despite recent developments in Machine Learning (ML) methodologies that outperformed traditional quantum mechanical models, several issues remain concerning explanatory insights for broad chemical predictions with an acceptable speed-accuracy trade-off. To overcome this, we present a novel supervised ML model to predict the ΔGsol for an array of solvent-solute pairs. Using two different ensemble regressor algorithms, we made fast and accurate property predictions using open-source chemical features, encoding complex electronic, structural, and surface area descriptors for every solvent and solute. By integrating molecular properties and chemical interaction features, we have analyzed individual descriptor importance and optimized our model though explanatory information form feature groups. On aqueous and organic solvent databases, ML models revealed the predictive relevance of solutes with increasing polar surface area and decreasing polarizability, yielding better results than state-of-the-art benchmark Neural Network methods (without complex quantum mechanical or molecular dynamic simulations). Both algorithms successfully outperformed previous ΔGsol predictions methods, with a maximum absolute error of 0.22 ± 0.02 kcal mol-1, further validated in an external benchmark database and with solvent hold-out tests. With these explanatory and statistical insights, they allow a thoughtful application of this method for predicting other thermodynamic properties, stressing the relevance of ML modeling for further complex computational chemistry problems.

Preparation of small peptoid-modified zwitterion-exchange/reversed-phase mixed-mode chromatography stationary phases for the separation of arenes and antibiotics

In this study, three small peptoids with different structures, named Sil-peptoids, were developed to improve the separation selectivity of zwitterion-exchange/reversed-phase mixed-mode chromatography stationary phases for multi-component complex drugs. Nonpolar, amphoteric, and alkaline drugs were used as test samples to demonstrate their retention behaviors in reversed-phase, ionic, and mixed-mode interactions. It was observed that different carboxyl anions in the small peptoids of the Sil-peptoids had vast differences in their stereo-selectivity. The stereo-selectivity and the influence of Sil-peptoids on the retention behavior of complex drugs and their interaction mechanism for the drug molecules were effectively evaluated through the combination of chromatographic analysis and molecular modeling. Finally, a mixture of drugs consisting of two polar and six non-polar drugs was used to obtain a separation effect with a resolution >1.5. Two other groups of polar antibiotics were used to verify the separation ability of the Sil-peptoids. The results indicated that the Sil-peptoids could separate multiple substances simultaneously. These novel stationary phases can be applied to the analysis of complex multi-component drugs.

Role of Explicit Solvation in Computational Modeling of Chemical Reactions: Mechanism of Cu(I) Transfer Between Thiolate-Based Chelators in Water

Solvation plays important roles in controlling the thermodynamic and kinetic aspects of chemical reactions. The conventional approaches to treat solvation via electronic structure methods are likely to become inadequate, when the reacting solutes have strong electrostatic and hydrogen bonding interactions with the solvent and undergo significant structural changes during the course of the reaction. In this article, we present evidence of such solvent and structural effects in the computational study of the Cu(I) transfer reaction between thiolate-based chelators dithiobutylamine (DTBA) and dithiotheritol (DTT) in water, inspired from biological copper trafficking phenomena. We propose a general solution to the problem by combining classical molecular dynamics (MD) simulations of the bulk system and static quantum chemistry calculations. The fluctuating solvation shell was estimated from MD, and energetics was assessed by averaging QM energies of a series of molecular clusters constructed from the MD snapshots. Applying this approach, we propose a reaction pathway with estimates of relative intermediate stabilities and barriers, which suggest the overall reaction to be reversible in nature and likely to go through both two and three coordinated intermediates, confirming previous studies of similar protein analogues. An interesting fact that emerged from our study was the strong indication that the rate-determining step is the deprotonation of initial thiol bound Cu(I) complex, without involving any Cu(I)-S bonds. The proposed method will lead to a better treatment of solvations, and these mechanistic insights will aid our understanding of biological copper(I) trafficking.

Determination of Amino Acids' pKa: Importance of Cavity Scaling within Implicit Solvation Models and Choice of DFT Functionals

Protonation states of molecules significantly influence the thermodynamics and kinetics of chemical reactions. This is especially important in biochemical processes, where appropriate protonation states of amino acids control the exo/endoergicity of practically all biochemical cycles. This paper is focused on appraisal of the impact of DFT functionals and PCM solvation models on the accuracy of pKa evaluations for all proteinogenic amino acids. Eight functionals (B3LYP, PBE0, revPBE0, M06-2X, M11, M11-L, TPSSh, and ωB97X-D) and four basis sets are considered, together with four kinds of implicit solvation models when additional attention is paid to a cavity construction. An influence of nonelectrostatic contributions and Wertz's corrections on Gibbs free energy is investigated together with accuracy of provided proton solvation energy. The best model is based on the M06-2X/6-311++G**/D-PCM/UAKS computational level. The fitting procedure is utilized to improve the accuracy of the evaluated models. All of these results are also compared with values obtained from the COSMOtherm program and CCSD(T) calculations. Results for cysteine and histidine are discussed individually, as they can be found in different protonation states at neutral pH.

Solvation energies of the ferrous ion in water and in ammonia at various temperatures

CONTEXT

The solvation of metal ions is crucial to understanding relevant properties in physics, chemistry, or biology. Therefore, we present solvation enthalpies and solvation free energies of the ferrous ion in water and ammonia. Our results agree well with the experimental reports for the hydration free energy and hydration enthalpy. We obtained [Formula: see text] kJ mol[Formula: see text] for the hydration free energy and [Formula: see text] kJ mol[Formula: see text] for the hydration enthalpy of ferrous ion in water at room temperature. At ambient temperature, we obtained [Formula: see text] kJ mol[Formula: see text] as the [Formula: see text] ammoniation free energy and [Formula: see text] kJ mol[Formula: see text] for the ammoniation enthalpy. In addition, the free energy of solvation is deeply affected when the temperature increases. This pattern can be attributed to the rise of entropy when the temperature rises. Besides, the temperature does not affect the ammoniation enthalpies and the hydration enthalpy of the [Formula: see text] ion.

METHOD

All the geometry optimizations are performed at the MP2 methods associated with the 6-31++g(d,p) basis set of Pople. solvated phase structures of [Formula: see text] ion in water or in ammonia are performed using the PCM model. The [Formula: see text] program suite was used to perform all the calculations. The program TEMPO was also used to evaluate the temperature sensitivity of the different obtained geometries.

DOX_BDW: Incorporating Solvation and Desolvation Effects of Cavity Water into Nonfitting Protein-Ligand Binding Affinity Prediction

Accurate prediction of the protein-ligand binding affinity (PLBA) with an affordable cost is one of the ultimate goals in the field of structure-based drug design (SBDD), as well as a great challenge in the computational and theoretical chemistry. Herein, we have systematically addressed the complicated solvation and desolvation effects on the PLBA brought by the difference of the explicit water in the protein cavity before and after ligands bind to the protein-binding site. Based on the new solvation model, a nonfitting method at the first-principles level for the PLBA prediction was developed by taking the bridging and displaced water (BDW) molecules into account simultaneously. The newly developed method, DOX_BDW, was validated against a total of 765 noncovalent and covalent protein-ligand binding pairs, including the CASF2016 core set, Cov_2022 covalent binding testing set, and six testing sets for the hit and lead compound optimization (HLO) simulation. In all of the testing sets, the DOX_BDW method was able to produce PLBA predictions that were strongly correlated with the corresponding experimental data (R = 0.66-0.85). The overall performance of DOX_BDW is better than the current empirical scoring functions that are heavily parameterized. DOX_BDW is particularly outstanding for the covalent binding situation, implying the need for considering an electronic structure in covalent drug design. Furthermore, the method is especially recommended to be used in the HLO scenario of SBDD, where hundreds of similar derivatives need to be screened and refined. The computational cost of DOX_BDW is affordable, and its accuracy is remarkable.

Disulfide radical anion as a super-reductant in biology and photoredox chemistry

Disulfides are involved in a broad range of radical-based synthetic organic and biochemical transformations. In particular, the reduction of a disulfide to the corresponding radical anion, followed by S-S bond cleavage to yield a thiyl radical and a thiolate anion plays critical roles in radical-based photoredox transformations and the disulfide radical anion in conjunction with a proton donor, mediates the enzymatic synthesis of deoxynucleotides from nucleotides within the active site of the enzyme, ribonucleotide reductase (RNR). To gain fundamental thermodynamic insight into these reactions, we have performed experimental measurements to furnish the transfer coefficient from which the standard E0(RSSR/RSSR˙-) reduction potential has been determined for a homologous series of disulfides. The electrochemical potentials are found to be strongly dependent on the structures and electronic properties of the substituents of the disulfides. In the case of cysteine, a standard potential of E0(RSSR/RSSR˙-) = -1.38 V vs. NHE is determined, making the disulfide radical anion of cysteine one of the most reducing cofactors in biology.

Computational Insights into the Regeneration of Ovothiol and Ergothioneine and Their Selenium Analogues by Glutathione

Ovothiol and ergothioneine are powerful antioxidants that readily react with oxidants by forming their respective disulfides. In fact, ovothiol is widely considered one of the most powerful natural antioxidants. However, for these antioxidants to be again involved in reacting with oxidants, they must be regenerated via the reduction of the disulfide bonds. In the present work, the regeneration of the antioxidants ovothiol and ergothioneine and their selenium analogues, by the closed-shell nucleophilic attack of glutathione, was investigated using density functional theory. From the calculated thermodynamic data, the attack of glutathione on OSSO and EYYE (where Y = S and/or Se) will readily occur in solution. Moreover, in comparison to the reference reaction GSH + GSSG → GSSG + GSH, all reactions are expected to be faster. Overall, the results presented herein show that the key antioxidant GSH should readily recycle ovothiol, ovoselenol, ergothioneine, and ergoseloneine from OYYO and EYYE (where Y = S and/or Se).

Modeling Differential Enthalpy of Absorption of CO2 with Piperazine as a Function of Temperature

Temperature-dependent correlations for equilibrium constants (ln K) and heat of absorption (ΔH_abs) of different reactions (i.e., deprotonation, double deprotonation, carbamate formation, protonated carbamate formation, dicarbamate formation) involved in the piperazine (PZ)/CO₂/H₂O system have been calculated using computational chemistry based ln K values input to the Gibbs-Helmholtz equation. This work also presents an extensive study of gaseous phase free energy and enthalpy for different reactions using composite (G3MP2B3, G3MP2, CBS-QB3, and G4MP2) and density functional theory [B3LYP/6-311++G(d,p)] methods. The explicit solvation shell (ESS) model and SM8T solvation free energy coupled with gaseous phase density functional theory calculations give temperature-dependent reaction equilibrium constants for different reactions. Calculated individual and overall reaction equilibrium constants and enthalpies of different reactions involved in CO₂ absorption in piperazine solution are compared against experimental data, where available, in the temperature range 273.15-373 K. Postcombustion CO₂ capture (PCC) is a temperature swing absorption-desorption process. The enthalpy of the solution directly correlates with the steam requirement of the amine regeneration step. Temperature-dependent correlations for ln K and ΔH_abs calculated using computational chemistry tools can help evaluate potential PCC solvents' thermodynamics and cost-efficiency. These correlations can also be employed in thermodynamic models (e.g., e-UNIQUAC, e-NRTL) to better understand postcombustion CO₂ capture solvent chemistry.

Microsecond molecular dynamics studies of cholesterol-mediated myelin sheath degeneration in early Alzheimer's disease

Cholesterol-mediated perturbations of membrane structural integrity are key early events in the molecular pathogenesis of Alzheimer's disease (AD). In AD, protein misfolding (proteopathy) and pro-inflammatory conditions (immunopathy) culminate in neuronal death, a process enabled by altered membrane biophysical properties which render neurons more susceptible to proteopathic and immunopathic cytotoxicities. Since cholesterol is a principal neuronal membrane lipid, normal cholesterol homeostasis is central to membrane health; also, since increased cholesterol composition is especially present in neuronal myelin sheath (i.e. brain "white matter"), recent studies have not surprisingly revealed that white matter atrophy precedes the conventional biomarkers of AD (amyloid plaques, tau tangles). Employing extensive microsecond all-atom molecular dynamics simulations, we investigated biophysical and mechanical properties of myelin sheath membrane as a function of cholesterol mole fraction (χ_CHL). Impaired χ_CHL modulates multiple bilayer properties, including surface area per lipid (APL), chain order, number and mass density profiles, area compressibility and bending moduli, bilayer thickness, lipid tilt angles, H-bonding interactions and tail interdigitation. The increased orientational ordering of both palmitoyl and oleoyl chains in model healthy myelin sheath (HMS) membranes illustrates the condensing effect of cholesterol. With an increase in χ_CHL, number density profiles of water tend to attain bulk water number density more quickly, indicating shrinkage in the interfacial region with increasing χ_CHL. The average tilt value is 11.5° for the C10-C13 angle in cholesterol and 64.2° for the P-N angle in POPC lipids in HMS. These calculations provide a molecular-level understanding of myelin sheath susceptibility to pathology as an early event in the pathogenesis of AD.

Axonal plasma membrane-mediated toxicity of cholesterol in Alzheimer's disease: A microsecond molecular dynamics study

Alzheimer's disease is increasingly being recognized as an immune-mediated disease of brain. Since physiological brain health and brain immune function is dependent upon homeostatic neuronal membrane structure and function, alterations in membrane lipid biochemistry may predispose to disease. Brain is rich in cholesterol, and cholesterol metabolism dysfunction is a known risk factor for AD. Employing extensive microsecond all-atom molecular dynamics simulations, we investigated the properties of model neuronal membranes as a function of cholesterol concentration; phospholipid and phospholipid/cholesterol bilayers were also simulated to compare against available experimental data. Increased cholesterol concentrations compact and stiffen the lipid membrane, reducing permeability while modulating local water densities in the peri-membranous environment. Conversely, lower cholesterol mole fraction yields membranes with increased molecular disorder, enhanced fluidity, higher molecular tilting, and augmented interdigitation between bilayer leaflet lipids. Our findings provide a molecular insight on effect of cholesterol composition on various biochemical processes occurring at neuronal axon plasma membrane. These calculations also endeavor to establish a membrane-based link between cholesterol as an AD risk factor and possible AD pathology.

Understanding Carbamate Formation Reaction Thermochemistry of Amino Acids as Solvents for Postcombustion CO2 Capture

The carbamate stability constant for a data set of 10 amino acids, having potential for being postcombustion CO₂ capture (PCC) solvents, has been calculated using various implicit and explicit solvation shell models. This work also includes an extensive study of gas-phase free energy and enthalpy for the amino acid carbamate formation reaction with the Hartree Fock method, density functional methods [B3LYP/6-311++G(d,p)], and composite methods (G3MP2B3, G3MP2, CBS-QB3, and G4MP2). Ideal PCC solvent properties require finding a profitable tradeoff between various thermodynamic and system optimization parameters. Benchmark gaseous-phase and solution-phase thermodynamic properties given in this work can help in making informed decisions when choosing promising PCC solvents. The temperature dependency of the carbamate stability constant of amino acids is predicted using PCM and SM8T implicit solvation models. PCC is a temperature swing absorption-desorption process, and the high-temperature sensitivity of the ln K_c^AmCOO- value is of vital importance in attaining cost-efficient processes.

Density Functional Calculations for Aqueous Silver Clusters Containing Water and Nitrate Ligands

The incorporation of a nitrate ion into silver-water aqueous clusters has been examined using PBE0 density functional theory with the solvent model density (SMD) solvation model. The Gibbs free energy of solvation and other thermodynamic variables are calculated using the harmonic/rigid rotor/ideal gas model at 298.15 K for aqueous solutes including the effects of solute relaxation in water and with London dispersive forces at the D3 level. Free energies of solvation for Ag⁺ and NO₃^- were found to agree well with experimental values of -118.2 and -60.1 kcal/mol, respectively, calculated using cluster-continuum models with six to eight water molecules and including solute relaxation and London D3 dispersive interactions. An analysis of data of varying cluster size upon calculated free energy is presented. A direct procedure is applied to aqueous clusters such as Ag_n^z(NO₃^-)(H₂O)₅, Ag_n^z(H₂O)₅, and (NO₃^-)(H₂O)₆ n = 1-4; z = 0, +1 in the SMD solvent representation to calculate equilibrium constants for nitrate association with silver clusters in solution that includes fully relaxed solutes. The equilibrium structures of the nitrate-containing clusters involve one or more bonds from nitrate oxygen to positive silver clusters. Water molecules interact with nitrate through H atoms, and overall, the structure represents a silver nitrate cluster with water ligands having similarity to a close ion pair in many aspects. The neutral silver atom is attached to nitrate through H-bonded water molecules. The ratio of nitrate-containing silver clusters to nitrate-free clusters using a calculated equilibrium constant of 0.51 L/mol for Ag⁺ is small in the range of many experiments. Similar values are found for positive silver clusters up to four atoms in size. The resulting procedures were applied to aqueous clusters of Ag_n(NO₃)_m^+(n-m) that have been previously experimentally studied for silver reduction in aqueous solution. A chain-like structure with collinear and bidentate oxygen bonds to silver was found, and the equilibrium constants for clustering were determined. A simplified model calculation for the reduction of Ag(H₂O)₆⁺ clusters in the presence of silver clusters in aqueous media was studied to understand catalytic effects observed in these systems. The reduction potentials vary with silver cluster size indicating a more favorable reduction caused by the presence of larger silver clusters.

Role of Water in the Reaction Mechanism and endo/exo Selectivity of 1,3-Dipolar Cycloadditions Elucidated by Quantum Chemistry and Machine Learning

Asymmetric 1,3-dipolar cycloadditions of azomethine ylides with activated olefins are among the most important and versatile methods for the synthesis of enantioenriched pyrroline and pyrrolidine derivatives. Despite both theoretical and practical importance, the role of water molecules in the reactivity and endo/exo selectivity remains unclear. To explore how water accelerates the reactions and improves the endo/exo selectivity of the cycloadditions of 1,3-dipole phthalazinium-2-dicyanomethanide (1) and two dipolarophiles, an ab initio-quality neural network potential that overcomes the computational bottleneck of explicitly considering water molecules was used. It is demonstrated that not only the nature of both the dipolarophile and the 1,3-dipole, but also the solvent medium, can perturb or even alter the reaction mechanism. An extreme case was found for the reaction of 1,3-dipole 1 with methyl vinyl ketone, in which the reaction mechanism changes from a concerted to a stepwise mode on going from MeCN to H₂ O as solvent, with formation of a zwitterionic intermediate that is a very shallow minimum on the energy surface. Thus, high stereocontrol can still be expected despite the stepwise nature of the mechanism. The results indicate that water can induce global polarization along the reaction coordinate and highlight the role of microsolvation effects and bulk-phase effects in reproducing the experimentally observed aqueous acceleration and enhanced endo/exo selectivity.

The multiple dissociation constants of glutathione disulfide: interpreting experimental pH-titration curves with ab initio MD simulations

The hexapeptide glutathione disulfide (GSSG) has six ionizable groups with six associated dissociation constants. The experimentally measured pH-titration curve, however, does not exhibit the six corresponding equivalence points and bears little resemblance to standard textbook examples of acid-base pH-titration curves. The curve highlights the difficulties in determining multiple pKa values of polyprotic acids - typically proteins and peptides - from experiment. The six pKa values of GSSG can, however, be estimated using Car-Parrinello molecular dynamics (CPMD) simulations in conjunction with metadynamics sampling of the underlying free energy landscape of the dissociation reactions. Ab initio MD simulations were performed on a GSSG molecule solvated by 200 water molecules. Using the estimated pKa values the theoretical titration curve was calculated and found to be in good agreement with experiment. The results clearly highlight how dissociation constants estimated from ab initio MD simulations can facilitate the interpretation of the pH-titration curves of complex chemical and biological systems.

Performance of 3D-RISM-KH in Predicting Hydration Free Energy: Effect of Solute Parameters

The three-dimensional reference interaction site model molecular solvation theory with the Kovalenko-Hirata closure relation has been shown to produce excellent solvation characteristics for a large class of (bio)chemical systems in solution. Correct calculation of hydration free energy is central to successful application of any solvation model. In order to find out the best possible force-field parameters to be used for hydration free energy calculation with the aforementioned theory, we have developed an extended database containing a large number of experimental solvation free energies available in the current literature and used a plethora of theoretical models for assessment. The general Amber force field was found to perform satisfactorily, whereas special care should be taken in solute charge assignment with the universal force field.

Theoretical modeling of pKa's of thiol compounds in aqueous solution

The pK_a's of different kinds of thiols (R-SH) were investigated by using the M06-2X method with a SMD_sSAS model.

Isodesmic reaction for accurate theoretical pKa calculations of amino acids and peptides

Theoretical and quantitative prediction of pKa values at low computational cost is a current challenge in computational chemistry. We report that the isodesmic reaction scheme provides semi-quantitative predictions (i.e. mean absolute errors of 0.5-1.0 pKa unit) for the pKa1 (α-carboxyl), pKa2 (α-amino) and pKa3 (sidechain groups) of a broad set of amino acids and peptides. This method fills the gaps of thermodynamic cycles for the computational pKa calculation of molecules that are unstable in the gas phase or undergo proton transfer reactions or large conformational changes from solution to the gas phase. We also report the key criteria to choose a reference species to make accurate predictions. This method is computationally inexpensive and makes use of standard density functional theory (DFT) and continuum solvent models. It is also conceptually simple and easy to use for researchers not specialized in theoretical chemistry methods.

In vitro study of a glucose attached poly(aryl ether) dendron based gel as a drug carrier for a local anaesthetic

A poly(aryl ether) dendron based gelator as an efficient transdermal drug carrier for the controlled release of prilocaine hydrochloride.

Charge Transfer Modulated Self-Assembly in Poly(aryl ether) Dendron Derivatives with Improved Stability and Transport Characteristics

Alteration of native gelation properties of anthracene and pyrene cored first generation poly(aryl ether) dendrons, G1-An and G1-Py, by introducing a common acceptor, 2,4,7-trinitro-9H-fluoren-9-one (TNF), results in forming charge transfer gels in long chain alcoholic solvents. This strategy leads to significant perturbation of optical and electronic properties within the gel matrix. Consequently, a noticeable increase of their electrical conductivities is observed, making these poly(aryl ether) dendron based gels potential candidates for organic electronics. While the dc-conductivity (σ) value for the native gel from G1-An is 2.8 × 10^-4 S m^-1, the value increased 3 times (σ = 8.7 × 10^-4 S m^-1) for its corresponding charge transfer gel. Further, the dc-conductivity for the native gel self-assembled from G1-Py dramatically enhanced by approximately an order of magnitude from 4.9 × 10^-4 to 1.3 × 10^-3 S m^-1, under the influence of an acceptor. Apart from H-bonding and π···π interactions, charge transfer results in the formation of a robust 3D network of fibers, with improved aspect ratio, providing high thermo-mechanical stability to the gels compared to the native ones. The charge transfer gels self-assembled from G1-An/TNF (1:1) and G1-Py/TNF exhibit a 7.3- and 2.5-fold increase in their yield stress, respectively, compared to their native assemblies. A similar trend follows in the case of their thermal stabilities. This is attributed to the typical bilayer self-assembly of the former which is not present in the case of G1-Py/TNF charge transfer gel. Density functional calculations provide deeper insights accounting for the role of charge transfer interactions in the mode of self-assembly. The 1D potential energy surface for the G1-An/TNF dimer and G1-Py/TNF dimer is found to be 11.8 and 1.9 kcal mol^-1 more stable than their corresponding native gel dimers, G1-An/G1-An and G1-Py/G1-Py, respectively.

Calculating Free Energy Changes in Continuum Solvation Models

We recently showed for a large data set of pKas and reduction potentials that free energies calculated directly within the SMD continuum model compares very well with corresponding thermodynamic cycle calculations in both aqueous and organic solvents [ Phys. Chem. Chem. Phys. 2015 , 17 , 2859 ]. In this paper, we significantly expand the scope of our study to examine the suitability of this approach for calculating general solution phase kinetics and thermodynamics, in conjunction with several commonly used solvation models (SMD-M062X, SMD-HF, CPCM-UAKS, and CPCM-UAHF) for a broad range of systems. This includes cluster-continuum schemes for pKa calculations as well as various neutral, radical, and ionic reactions such as enolization, cycloaddition, hydrogen and chlorine atom transfer, and SN2 and E2 reactions. On the basis of this benchmarking study, we conclude that the accuracies of both approaches are generally very similar-the mean errors for Gibbs free energy changes of neutral and ionic reactions are approximately 5 and 25 kJ mol(-1), respectively. In systems where there are significant structural changes due to solvation, as is the case for certain ionic transition states and amino acids, the direct approach generally afford free energy changes that are in better agreement with experiment.

Ab Initio MD Simulations of the Brønsted Acidity of Glutathione in Aqueous Solutions: Predicting pKa Shifts of the Cysteine Residue

The tripeptide glutathione (GSH) is one of the most abundant peptides and the major repository for nonprotein sulfur in both animal and plant cells. It plays a critical role in intracellular oxidative stress management by the reversible formation of glutathione disulfide with the thiol-disulfide pair acting as a redox buffer. The state of charge of the ionizable groups of GSH can influence the redox couple, and hence the pKa value of the cysteine residue of GSH is critical to its functioning. Here we report ab initio Car-Parrinello molecular dynamics simulations of glutathione solvated by 200 water molecules, all of which are considered in the simulation. We show that the free-energy landscape for the protonation-deprotonation reaction of the cysteine residue of GSH computed using metadynamics sampling provides accurate estimates of the pKa and correctly predicts the shift in the dissociation constant values as compared with the isolated cysteine amino acid.

Ab Initio Molecular Dynamics Simulations of Amino Acids in Aqueous Solutions: Estimating pKa Values from Metadynamics Sampling

Changes in the protonation and deprotonation of amino acid residues in proteins play a key role in many biological processes and pathways. Here, we report calculations of the free-energy profile for the protonation-deprotonation reaction of the 20 canonical α amino acids in aqueous solutions using ab initio Car-Parrinello molecular dynamics simulations coupled with metadynamics sampling. We show here that the calculated change in free energy of the dissociation reaction provides estimates of the multiple pKa values of the amino acids that are in good agreement with experiment. We use the bond-length-dependent number of the protons coordinated to the hydroxyl oxygen of the carboxylic and the amine groups as the collective variables to explore the free-energy profiles of the Bronsted acid-base chemistry of amino acids in aqueous solutions. We ensure that the amino acid undergoing dissociation is solvated by at least three hydrations shells with all water molecules included in the simulations. The method works equally well for amino acids with neutral, acidic and basic side chains and provides estimates of the multiple pKa values with a mean relative error, with respect to experimental results, of 0.2 pKa units.

Calculations of pKa's and redox potentials of nucleobases with explicit waters and polarizable continuum solvation

The SMD implicit solvation model augmented with one and four explicit water molecules was used to calculate pKa's and redox potentials of N-methyl-substituted nucleic acid bases guanine, adenine, cytosine, thymine, and uracil. Calculations were carried out with the B3LYP/6-31+G(d,p) level of theory. The same numbers of water molecules were hydrogen bonded to the neutral, protonated, and deprotonated nucleobases in their unoxidized and oxidized forms. The improvement in pKa1 involving neutrals and cations was modest. By contrast, the improvement in pKa2 involving neutrals and anions was quite significant, reducing the mean absolute error from 4.6 pKa units with no waters, to 2.6 with one water and 1.7 with four waters. For the oxidation of nucleobases, adding explicit waters did little to improve E(X(•),H(+)/XH), possibly because both species in the redox couple are neutral molecules at pH 7.

Predicting pKa in Implicit Solvents: Current Status and Future Directions

Computational prediction of condensed phase acidity is a topic of much interest in the field today. We introduce the methods available for predicting gas phase acidity and pKas in aqueous and non-aqueous solvents including high-level electronic structure methods, empirical linear free energy relationships (LFERs), implicit solvent methods, explicit solvent statistical free energy methods, and hybrid implicit–explicit approaches. The focus of this paper is on implicit solvent methods, and we review recent developments including new electronic structure methods, cluster-continuum schemes for calculating ionic solvation free energies, as well as address issues relating to the choice of proton solvation free energy to use with implicit solvation models, and whether thermodynamic cycles are necessary for the computation of pKas. A comparison of the scope and accuracy of implicit solvent methods with ab initio molecular dynamics free energy methods is also presented. The present status of the theory and future directions are outlined.