Simulation of Liquid and Supercritical Hydrogen Sulfide and of Alkali Ions in the Pure and Aqueous Liquid

Corrections in the CHARMM36 Parametrization of Chloride Interactions with Proteins, Lipids, and Alkali Cations, and Extension to Other Halide Anions

The nonpolarizable CHARMM force field is one of the most widely used energy functions for all-atom biomolecular simulations. Chloride is the only halide ion included in the latest version, CHARMM36m, and is used widely in simulation studies, often as an electrolyte ion but also as the biological substrate of transport proteins and enzymes. Here, we find that existing parameters systematically underestimate the interaction of Cl^- with proteins and lipids. Accordingly, when examined in solution, little to no Cl^-association can be observed with most components of the protein, including backbone, polar side chains and aromatic rings. The strength of the interaction with cationic side chains and with alkali ions is also incongruent with experimental measurements, specifically osmotic coefficients of concentrated solutions. Consistent with these findings, a 4-μs trajectory of the Cl^--specific transport protein CLC-ec1 shows irreversible Cl^- dissociation from the so-called S_cen binding site, even in a 150 mM NaCl buffer. To correct for these deficiencies, we formulate a series of pair-specific Lennard-Jones parameters that override those resulting from the conventional Lorentz-Berthelot combination rules. These parameters, referred to as NBFIX, are systematically calibrated against available experimental data as well as ab initio geometry optimizations and energy evaluations, for a wide set of binary and ternary Cl^- complexes with protein and lipid analogs and alkali cations. Analogously, we also formulate parameter sets for the other three biological halide ions, namely, fluoride, bromide, and iodide. The resulting parameters are used to calculate the potential of mean force defining the interaction of each anion and each of the protein and lipid analogues in bulk water, revealing association free energies in the range of -0.3 to -3.3 kcal/mol, with the F^- complexes being the least stable. The NBFIX corrections also preserve the Cl^- occupancy of CLC-ec1 in a second 4-μs trajectory. We posit that these optimized molecular-mechanics models provide a more realistic foundation for all-atom simulation studies of processes entailing changes in hydration, recognition, or transport of halide anions.

Computational analysis of vibrational frequencies and rovibrational spectroscopic constants of hydrogen sulfide dimer using MP2 and CCSD(T)

Previous studies have shown that the weakly bonded H₂S dimer demands high level quantum chemical calculations to reproduce experimental values. We investigated the hydrogen bonding of H₂S dimer using MP2 and CCSD(T) levels of theory in combination with aug-cc-pV(D,T,Q)Z basis sets. More precisely, the binding energies, potential energy curves, rovibrational spectroscopic constants, decomposition lifetime, and normal vibrational frequencies were calculated. In addition, we introduced the local mode analysis of Konkoli-Cremer to quantify the hydrogen bonding in the H₂S dimer as well as providing for the first time the comprehensive decomposition of normal vibrational modes into local modes contributions, and a decomposition lifetime based on rate constant. The local mode force constant of the H₂S dimer hydrogen bond is smaller than that of the water dimer, in accordance with the weaker hydrogen bonding in the H₂S dimer.

New Molecular-Mechanics Model for Simulations of Hydrogen Fluoride in Chemistry and Biology

Hydrogen fluoride (HF) is the most polar diatomic molecule and one of the simplest molecules capable of hydrogen-bonding. HF deviates from ideality both in the gas phase and in solution and is thus of great interest from a fundamental standpoint. Pure and aqueous HF solutions are broadly used in chemical and industrial processes, despite their high toxicity. HF is a stable species also in some biological conditions, because it does not readily dissociate in water unlike other hydrogen halides; yet, little is known about how HF interacts with biomolecules. Here, we set out to develop a molecular-mechanics model to enable computer simulations of HF in chemical and biological applications. This model is based on a comprehensive high-level ab initio quantum chemical investigation of the structure and energetics of the HF monomer and dimer; (HF)_n clusters, for n = 3-7; various clusters of HF and H₂O; and complexes of HF with analogs of all 20 amino acids and of several commonly occurring lipids, both neutral and ionized. This systematic analysis explains the unique properties of this molecule: for example, that interacting HF molecules favor nonlinear geometries despite being diatomic and that HF is a strong H-bond donor but a poor acceptor. The ab initio data also enables us to calibrate a three-site molecular-mechanics model, with which we investigate the structure and thermodynamic properties of gaseous, liquid, and supercritical HF in a wide range of temperatures and pressures; the solvation structure of HF in water and of H₂O in liquid HF; and the free diffusion of HF across a lipid bilayer, a key process underlying the high cytotoxicity of HF. Despite its inherent simplifications, the model presented significantly improves upon previous efforts to capture the properties of pure and aqueous HF fluids by molecular-mechanics methods and to our knowledge constitutes the first parameter set calibrated for biomolecular simulations.

Cooperative nature of the sulfur centered hydrogen bond: investigation of (H2S)n (n = 2-4) clusters using an affordable yet accurate level of theory

Existing studies have shown that appreciably high level quantum chemical calculations are required to reproduce experimental energetic and geometric features of a H₂S dimer. This condition severely restricts any practical possibility of obtaining reliable results for higher order H₂S clusters. We have shown here that the binding energies calculated at the CCSD(T)/CBS level with counterpoise corrected geometries calculated at the MP2/aug-cc-pV(Q+d)Z level of theory excellently match with the experimental results for the H₂S dimer. Subsequently, the above mentioned levels of theory were used for trimers and tetramers. (H₂S)_n (n = 2-4) clusters were found to show cooperative strengthening of S-HS hydrogen bonds, which is clearly evident from the evolution of binding energies and hydrogen bond lengths, with increasing cluster size. Localized molecular orbital energy decomposition analyses have been carried out to understand how the contributions of various energy components modulate with the size of the clusters and what are their relative contributions towards the overall stabilization of the clusters. Natural bond orbital and atoms in molecules analyses were also carried out in order to look into the evolution of the electronic charge transfer and electron density topology with cluster size.

Computational insight into hydrogen persulfide and a new additive model for chemical and biological simulations

S-Sulfhydration of cysteine to the Cys-SSH persulfide is an oxidative post-translational modification that plays an important regulatory role in many physiological systems. Though hydrogen persulfide (H2S2) has recently been established as a signaling and cellular sulfhydration reagent, the chemistry and chemical biology of persulfides remain poorly explored. We first report an extensive high-level ab initio quantum chemical investigation of (H2S2)n, (H2S2)m·H2O, and (H2O)m·H2S2 clusters (n = 1-3 and m = 1, 2) and of H2S2 complexes with 19 compounds that model the side chains of naturally-occurring amino acids. The high polarizability of S necessitates the use of large, very diffuse, basis sets for proper description of H2S2 and its complexes. H2S2 possesses a skewed equilibrium geometry, with nonpolar trans and more polar cis conformers 6 and 8 kcal mol-1 higher in energy, respectively; the skewed conformation is preserved in all neutral and cationic complexes while a cis geometry prevails in some anionic complexes. H2S2 is found to be a better H-bond donor and a poorer acceptor than H2S, and that in complexes with H2O, alcohols and amines, H2S2 is a better H-bond donor. Radical delocalization on both S atoms stabilizes the perthiyl (HSS˙) over the thiyl (HS˙) radical and results in a ∼20 kcal mol-1 lower S-H homolytic bond dissociation in H2S2, making it a potential antioxidant. A simple additive model is optimized for H2S2 and used together with the TIP3P model and the CHARMM36 all-atom force field (FF) to investigate the structure and thermodynamic properties of liquid H2S2 and the solubility of H2S2 in water, and to model H2S2-protein interactions (for which new FF parameters are further developed). Very weak H-bonding characterizes liquid H2S2 and it is found immiscible in liquid water with a trend in H-bonding strengths between H2S2 and H2O in the order O-HO ≫ S-HO > O-HS. This work does not only provide a thorough description of the structure and energetics of H2S2 and its various complexes, but also yields a reliable FF for investigating H2S2 in chemistry and biology.

A Simple Additive Potential Model for Simulating Hydrogen Peroxide in Chemical and Biological Systems

Hydrogen peroxide (H₂O₂) has numerous industrial, environmental, medical, cosmetic, and biological applications. Given its importance, we provide a simple model as an alternative to experiment for studying the properties of pure liquid H₂O₂ and its concentrated aqueous solutions, which are hazardous, and for understanding the biological roles of H₂O₂ at the molecular level. A four-site additive model is calibrated for H₂O₂ based on the ab initio and experimental properties of the gaseous monomer and the density and heat of vaporization of liquid H₂O₂ at 0 °C. Our model together with the TIP3P water model reproduce the ab initio binding energies of (H₂O₂) _m, H₂O₂· nH₂O, and nH₂O₂·H₂O clusters ( m = 2, 3 and n = 1, 2) calculated at the MP2 level using the 6-311++G(d,p) or the 6-311++G(3df,3pd) basis set. It yields structure, the self-diffusion coefficient, heat capacity, and densities at temperatures up to 200 °C of the pure liquid in good agreement with experiment. The model correctly predicts the hydration free energy of H₂O₂ and reproduces the experimental density of aqueous H₂O₂ solutions at 0-96 °C. Investigation of the solvation of H₂O₂ and H₂O in aqueous H₂O₂ solutions reveals that, as in the gas phase, H₂O₂ is a better H-bond donor but poorer acceptor than H₂O and the bonding stability follows the order O_p-H_p···O_w > O_w-H_w···O_w ≥ O_p-H_p···O_p > O_w-H_w···O_p. Stronger H-bonding in H₂O/H₂O₂ mixtures than in the pure liquids is consistent with exothermic heats of mixing and explains why the observed density and vapor pressure of the aqueous solutions are higher and lower, respectively, than expected from ideal mixing. Results also show that H₂O₂ adopts a skewed equilibrium geometry in gas and liquid phases but more polar cis and nonpolar trans conformations also are accessible and will stabilize H₂O₂ in environments of different polarity. In sum, our simple model presents a reliable tool for simulating H₂O₂ in chemistry and biology.

Modeling Protein S-Aromatic Motifs Reveals Their Structural and Redox Flexibility

S-aromatic motifs are important noncovalent forces for protein stability and function but remain poorly understood. Hence, we performed quantum calculations at the MP2(full)/6-311++G(d,p) level on complexes between Cys (H₂S, MeSH) and Met (Me₂S) models with models of Phe (benzene, toluene), Trp (indole, 3-methylindole), Tyr (phenol, 4-methylphenol), and His (imidazole, 4-methylimidazole). The most stable gas-phase conformers exhibit binding energies of -2 to -6 kcal/mol, and the S atom lies perpendicular to the ring plane. This reveals preferential interaction with the ring π-system, except in the imidazoles where S binds edge-on to an N atom. Complexation tunes the gas-phase vertical ionization potentials of the ligands over as much as 1 eV, and strong σ- or π-type H-bonding supports charge transfer to the H-bond donor, rendering it more oxidizable. When the S atom acts as an H-bond acceptor (N/O-H_ar···S), calibration of the CHARMM36 force field (by optimizing pair-specific Lennard-Jones parameters) is required. Implementing the optimized parameters in molecular dynamics simulations in bulk water, we find stable S-aromatic complexes with binding free energies of -0.6 to -1.1 kcal/mol at ligand separations up to 8 Å. The aqueous S-aromatics exhibit flexible binding conformations, but edge-on conformers are less stable in water. Reflecting this, only 0.3 to 10% of the S-indole, S-phenol, and S-imidazole structures are stabilized by N/O-H_ar···S or S-H···O_ar/N_ar σ-type H-bonding. The wide range of energies and geometries found for S-aromatic interactions and their tunable redox properties expose the versatility and variability of the S-aromatic motif in proteins and allow us to predict a number of their reported properties.

Solubility Parameters of Permanent Gases

The solubility parameters,δH(Tb), of nonreactive permanent gases at their boiling pointsTb(<290 K) are calculated from individually discussed values of their molar enthalpies of vaporization and densities obtained from the literature. These values are tabulated and where available the coefficients of the temperature dependence expressionδH(T)are also tabulated. The trends noted in theδH(Tb)values are dealt with and the values are compared with those reported in the literature and derived from the solubilities of the gases in various solvents. TheδH(Tb)values are shown to correlate linearly with the depths of the potential wells (attractive interaction energies,ε/kB) for binary collisions of the gaseous molecules and with the surface tensions,σ(Tb), of the liquefied gases.