On the Implementation of Friedman Boundary Conditions in Liquid Water Simulations

Abstract

We have applied the image approximation to the reaction field as suggested by H.L. Friedman [Mol. Phys., 29, 1533 (1975)] by investigating appropriate cavity sizes and system parameters for use in molecular simulations. The energy of and the structure around a central simple point charge (SPC) water molecule in a dielectric cavity was found to be in good agreement with the properties of a liquid sample. To confine the water molecules within the cavity, we introduced a short-range repulsion between a real charge and its image as the Lennard-Jones repulsive potential between oxygen atoms of the SPC potential. For a system of 65 water molecules a cavity radius of 10.45 Å is appropriate; this radius is altered to 12.00 Å for a cavity surrounding 113 molecules. The effect of the boundary is restricted to the outer-most water layer which is in contact with the dielectric continuum.     

Adsorption of Eu3+ and Am3+ ion towards hard donor-based diglycolamic acid-functionalised carbon nanotubes: density functional theory guided experimental verification

Structure, bonding, energetic and thermodynamic parameters of trivalent Eu³⁺ and Am³⁺ with diglycolamic acid-functionalised carbon nanotubes (CNT–DGA) in gas and solvent phases are reported in order to understand their complexation and extraction behaviour. The calculation was performed with generalised gradient approximated BP86 density functional and hybrid B3LYP functional using SVP/TZVP basis set. The free energy of extraction, ΔG_ext, of Eu³⁺ and Am³⁺ was computed using the Born–Haber thermodynamic cycle in conjunction with conductor-like screening model solvation approach. The calculated free energy of extraction was found to be exergonic using an explicit cluster water model for hydrating the ions and it was found to be higher for Eu³⁺ ion over Am³⁺ ion. From the estimated distribution constant using synthesised CNT–DGA, it is indisputably established that the Eu³⁺ ion is preferentially extracted over Am³⁺ ion and hence confirms the acceptance of the explicit cluster model for ion solvation free energy and thermodynamic cycle for the evaluation of free energy of extraction in solution phase. This is perhaps the first of its kind that a combined theoretical and experimental study has been performed on the free energy of extraction of Eu³⁺ and Am³⁺ using CNT–DGA.     

Effective water model for Monte Carlo simulations of proteins

We present an effective theory for water. Our goal is to formulate on accurate model for the effects of solvation on protein dynamics, without incurring the huge computational cost and the slow temporal evolution typical of molecular dynamics simulations of liquids. We replace the individual water molecules in an all-atom potential with a local dielectric density field, with self interactions given by the Landau-Ginzburg free energy and external interactions by Lennard-Jones forces at the surface of the protein atoms. We explore conformational space with finite temperature Monte Carlo dynamics, using parallel Langevin and Fourier acceleration algorithms well suited to data-parallel computer architectures such as the Connection Machine. To establish the validity of our approximations, we compare our electrostatic contribution to the solvalion energy with the results of Lim, Bashford, and Karplus using a conventional static continuum dielectric cavity model, and the non electrostatic contributions with estimates of hydrophohic surface free energy. Our model can also accommodate ionic charges and temperature fluctuations, We propose future investigations extending our effective theory of solvation to include explicit orientational entropy and hydroxen-bonding terms. © 1995 John Wiley & Sons, Inc.     

Monte Carlo simulation of the water in a channel with charges.

A Monte Carlo simulation of water in a channel with charges suggests the existence of water in immobile, high density, essentially glasslike form near the charges. The channel model has a conical section with an opening through which water molecules can pass, at the narrow end of the cone, and a cylindrical section at the other end. When the charges are placed near the narrow section of the model, the "glass" effectively blocks the channel; with the charges removed, the channel opens. The effect can be determined from the rate of passage of the water molecules through the pore, from the average orientation of the water molecule, and from distortion of the distribution of molecules. In the simulations carried out to date, no external ions have been considered. In addition to the energy, the Helmholtz free energy has been calculated.     

The role of glutathione in cadmium ion detoxification: coordination modes and binding properties--a density functional study

In this investigation the reduced form of glutathione molecule (GSH) was considered as a model in the coordination chemistry of Cd²⁺ and group thiol-containing peptides. Three different forms of GSH, corresponding to the prevalent ones in gas-phase and in aqueous solution, were taken into account in the metallation process. The obtained complexes were characterized at the density functional B3LYP level with the purpose to give better insight in the chelation mechanism of GSH with heavy metal ions in living organisms. Solvent effects, whether with explicit water molecules or with polarizable continuum model (PCM), were considered on the most stable cadmium-complexes obtained by every GSH charged species examined.     

Hydration energy landscape of the active site cavity in cytochrome P450cam

Hydration of protein cavities influences protein stability, dynamics, and function. Protein active sites usually contain water molecules that, upon ligand binding, are either displaced into bulk solvent or retained to mediate protein–ligand interactions. The contribution of water molecules to ligand binding must be accounted for to compute accurate values of binding affinities. This requires estimation of the extent of hydration of the binding site. However, it is often difficult to identify the water molecules involved in the binding process when ligands bind on the surface of a protein. Cytochrome P450cam is, therefore, an ideal model system because its substrate binds in a buried active site, displacing partially disordered solvent, and the protein is well characterized experimentally. We calculated the free energy differences for having five to eight water molecules in the active site cavity of the unliganded enzyme from molecular dynamics simulations by thermodynamic integration employing a three-stage perturbation scheme. The computed free energy differences between the hydration states are small (within 12 kJ mol⁻¹) but distinct. Consistent with the crystallographic determination and studies employing hydrostatic pressure, we calculated that, although ten water molecules could in principle occupy the volume of the active site, occupation by five to six water molecules is thermodynamically most favorable. Proteins 32:381–396, 1998. © 1998 Wiley-Liss, Inc.     

Electronic structure calculation study of metal complexes with a phytosiderophore mugineic acid

Structural and thermodynamic properties of biologically important metal-mugineic acid complexes have been studied from the theoretical side in order to understand the metal-chelating mechanism of phytosiderophore mugineic acid at an atomic level. Density-functional theory methods combined with the polarizable continuum model (PCM) have been employed to obtain free energies of complex formation and redox potentials for metal-mugineic acid complexes in solution. It has been found that the free energies of complex formation calculated at the B3LYP/PCM level of theory are in moderate agreement with available experimental results. The inclusion of explicit water molecules interacting with the carboxylic groups in deprotonated mugineic acid through strong hydrogen-bonds is found to further improve the calculated free energies of complex formation.     

Contributions of water transfer energy to protein‐ligand association and dissociation barriers: Watermap analysis of a series of p38α MAP kinase inhibitors

In our previous work, we proposed that desolvation and resolvation of the binding sites of proteins can serve as the slowest steps during ligand association and dissociation, respectively, and tested this hypothesis on two protein‐ligand systems with known binding kinetics behavior. In the present work, we test this hypothesis on another kinetically‐determined protein‐ligand system—that of p38α and eight Type II BIRB 796 inhibitor analogs. The k_on values among the inhibitor analogs are narrowly distributed (10⁴ ≤ k_on ≤ 10⁵ M^?1 s^?1), suggesting a common rate‐determining step, whereas the k_off values are widely distributed (10^?1 ≤ k_off ≤ 10^?6 s^?1), suggesting a spectrum of rate‐determining steps. We calculated the solvation properties of the DFG‐out protein conformation using an explicit solvent molecular dynamics simulation and thermodynamic analysis method implemented in WaterMap to predict the enthalpic and entropic costs of water transfer to and from bulk solvent incurred upon association and dissociation of each inhibitor. The results suggest that the rate‐determining step for association consists of the transfer of a common set of enthalpically favorable solvating water molecules from the binding site to bulk solvent. The rate‐determining step for inhibitor dissociation consists of the transfer of water from bulk solvent to specific binding site positions that are unfavorably solvated in the apo protein, and evacuated during ligand association. Different sets of unfavorable solvation are evacuated by each ligand, and the observed dissociation barriers are qualitatively consistent with the calculated solvation free energies of those sets.     

Temperature and pressure denaturation of chignolin: Folding and unfolding simulation by multibaric-multithermal molecular dynamics method

A multibaric‐multithermal molecular dynamics (MD) simulation of a 10‐residue protein, chignolin, was performed. All‐atom model with the Amber parm99SB force field was used for the protein and the TIP3P model was used for the explicit water molecules. This MD simulation covered wide ranges of temperature between 260 and 560 K and pressure between 0.1 and 600 MPa and sampled many conformations without getting trapped in local‐minimum free‐energy states. Folding events to the native β‐hairpin structure occurred five times and unfolding events were observed four times. As the temperature and/or pressure increases, fraction of folded chignolin decreases. The partial molar enthalpy change ΔH and partial molar volume change ΔV of unfolding were calculated as ΔH = 24.1 ± 4.9 kJ/mol and ΔV = ?5.6 ± 1.5 cm³/mol, respectively. These values agree well with recent experimental results. Illustrating typical local‐minimum free‐energy conformations, folding and unfolding pathways were revealed. When chignolin unfolds from the β‐hairpin structure, only the C terminus or both C and N termini open first. It may undergo an α‐helix or 3₁₀‐helix structure and finally unfolds to the extended structure. Difference of the mechanism between temperature denaturation and pressure denaturation is also discussed. Temperature denaturation is caused by making the protein transferred to a higher entropy state and making it move around more with larger space. The reason for pressure denaturation is that water molecules approach the hydrophobic residues, which are not well hydrated at the folded state, and some hydrophobic contacts are broken. Proteins 2012;. © 2012 Wiley Periodicals, Inc.     

Molecular Dynamics Free Energy Calculations to Assess the Possibility of Water Existence in Protein Nonpolar Cavities

Are protein nonpolar cavities filled with water molecules? Although many experimental and theoretical investigations have been done, particularly for the nonpolar cavity of IL-1β, the results are still conflicting. To study this problem from the thermodynamic point of view, we calculated hydration free energies of four protein nonpolar cavities by means of the molecular dynamics thermodynamic integration method. In addition to the IL-1β cavity (69 ų), we selected the three largest nonpolar cavities of AvrPphB (81 ų), Trp repressor (87 ų), and hemoglobin (108 ų) from the structural database, in view of the simulation result from another study that showed larger nonpolar cavities are more likely to be hydrated. The calculations were performed with flexible and rigid protein models. The calculated free energy changes were all positive; hydration of the nonpolar cavities was energetically unfavorable for all four cases. Because hydration of smaller cavities should happen more rarely, we conclude that existing protein nonpolar cavities are not likely to be hydrated. Although a possibility remains for much larger nonpolar cavities, such cases are not found experimentally. We present a hypothesis to explain this: hydrated nonpolar cavities are quite unstable and the conformation could not be maintained.     

Molecular Parameterisation and the Theoretical Calculation of Electrode Potentials

The difference in reduction potentials between ortho and para-benzoquinones has been calculated. The employs gas phase ab initio and semi-empirical computations in combination with free energy perturbation theory applied to gas and solution phase Monte Carlo simulations. The effects on calculated results of altering solute electrostatic parameterisation in solution phase simulations is examined. Atom centred charges derived from the molecular electrostatic potentials, MEPs, from optimised ab initio wavefunctions and charges generated by consideration of hydrogen bonded complexes are considered. Parameterisation of hydroxyl torsions in hydroquinone molecules is treated in a physically realistic manner. The coupled torsional system of the ortho-hydrobenzoquinone molecule is described by a potential energy surface calculated using gas phase AM1 semi-empirical computations rather than the simple torsional energy functions frequently employed in such calculations. Calculated differences in electrode potentials show that the electrostatic interactions of quinone and hydroquinone molecules in aqueous solution are not well described by atom centred charges derived from ab initio calculated MEPs. Moreover, results in good agreement with the experimental reduction potential difference can be obtained by employing high level ab initio calculations and solution phase electrostatic parameters developed by consideration of hydrogen bonded complexes.     

Hybrid MC/QC simulations of water-assisted proton transfer in nucleosides. Guanosine and its analog acyclovir

To provide an in-depth insight into the molecular basis of spontaneous tautomerism in DNA and RNA base pairs, a hybrid Monte Carlo (MC)–quantum chemical (QC) methodology is implemented to map two-dimensional potential energy surfaces along the reaction coordinates of solvent-assisted proton transfer processes in guanosine and its analog acyclovir in aqueous solution. The solvent effects were simulated by explicit inclusion of water molecules that model the relevant part of the first hydration shell around the solute. The position of these water molecules was estimated by carrying out a classical Metropolis Monte Carlo simulation of dilute water solutions of the guanosine (Gs) and acyclovir (ACV) and subsequently analyzing solute–solvent intermolecular interactions in the statistically-independent MC-generated configurations. The solvent-assisted proton transfer processes were further investigated using two different ab initio MP2 quantum chemical approaches. In the first one, potential energy surfaces of the ‘bare’ finite solute–solvent clusters containing Gs/ACV and four water molecules (MP2/6-31+G(d,p) level) were explored, while within the second approach, these clusters were embedded in ‘bulk’ solvent treated as polarizable continuum (C-PCM/MP2/6-31+G(d,p) level of theory). It was found that in the gas phase and in water solution, the most stable tautomer for guanosine and acyclovir is the 1H-2-amino-6-oxo form followed by the 2-amino-6-(sZ)-hydroxy form. The energy barriers of the water-assisted proton transfer reaction in guanosine and in acyclovir are found to be very similar – 11.74 kcal mol^?1 for guanosine and 11.16 kcal mol^?1 for acyclovir, and the respective rate constants (k = 1.5?×?10¹ s^?1, guanosine and k = 4.09?×?10¹ s^?1, acyclovir), are sufficiently large to generate the 2-amino-6-(sZ)-hydroxy tautomer. The analysis of the reaction profiles in both compounds shows that the proton transfer processes occur through the asynchronous concerted mechanism.     

1H-NMR spin-spin relaxation time assessment of the reflection coefficient of the Triticale seed cell wall–plasmalemma barrier to mannitol

The ¹H-NMR spin-spin relaxation time (T₂) in Triticale seeds swelling in external osmotica, polyethylene glycol 8000 or mannitol can identify both bound and free water. At the same water content, the free water spin-spin relaxation time increases for seeds imbibed with the mannitol solution, demonstrating inadequate water potential adjustment. The exchange rate of free/bound water molecules is apparently influenced by the driving force for water flow. The reciprocal lifetime of free water molecules, as a measure of water flow through the main cell barrier, was obtained. From a model of the seed as a resistance–capacitor network for water flow, a method was derived for calculating the reflection coefficient σ as a lifetime ratio of the free water molecules in seeds imbibed with two different osmotica (one penetrating across the main cell barrier and one not penetrating) at the same water potential. The ¹H-NMR method and the classical method based on volume rate changes yielded reflection coefficients for mannitol for the cell wall–plasmalemma barrier of 0.78 ± 0.08 and 0.68 ± 0.06, respectively.     

Cadmium(II) sorption from water samples by powdered marble wastes

Abstract

A series of batch adsorption experiments were carried out, with the aim of removing cadmium ions from aqueous solutions and water samples using powdered marble wastes (PMW) as an effective inorganic sorbent. PMW is inexpensive, widespread, and may be considered as environmental problem. The main parameters (i.e. solution pH, sorbent and cadmium concentrations, stirring time, and temperature) influencing the sorption process were investigated. The results obtained for sorption of cadmium ions onto PMW are well described by the Freundlich and Langmuir models. The Dubinin-Radushkevick (D–R) isotherm model was applied to describe the nature of the adsorption of the metal ion; it was found that the adsorption process was chemical in nature. The thermodynamic parameters were also calculated from the Gibbs free energy change (ΔG°), enthalpy (AH°) and entropy (ΔS°). These parameters indicated that the adsorption process of cadmium(II) ions on PMW was spontaneous and endothermic in nature. Under the optimum experimental conditions employed the removal of ca ~100% of Cd²⁺ ions was attained. The procedure was successfully applied to removal of the cadmium ions from aqueous and various natural water samples. The adsorption mechanism is discussed.     

On the calculation of binding free energies using continuum methods: Application to MHC class I protein-peptide interactions

This paper describes a methodology to calculate the binding free energy (ΔG) of a protein-ligand complex using a continuum model of the solvent. A formal thermodynamic cycle is used to decompose the binding free energy into electrostatic and non-electrostatic contributions. In this cycle, the reactants are discharged in water, associated as purely nonpolar entities, and the final complex is then recharged. The total electrostatic free energies of the protein, the ligand, and the complex in water are calculated with the finite difference Poisson-Boltzmann (FDPB) method. The nonpolar (hydrophobic) binding free energy is calculated using a free energy-surface area relationship, with a single alkane/water surface tension coefficient (γ_aw). The loss in backbone and side-chain configurational entropy upon binding is estimated and added to the electrostatic and the nonpolar components of ΔG. The methodology is applied to the binding of the murine MHC class I protein H-2K^b with three distinct peptides, and to the human MHC class I protein HLA-A2 in complex with five different peptides. Despite significant differences in the amino acid sequences of the different peptides, the experimental binding free energy differences (ΔΔG_exp) are quite small (<0.3 and <2.7 kcal/mol for the H-2K^b and HLA-A2 complexes, respectively). For each protein, the calculations are successful in reproducing a fairly small range of values for ΔΔG_calc (<4.4 and <5.2 kcal/mol, respectively) although the relative peptide binding affinities of H-2K^b and HLA-A2 are not reproduced. For all protein-peptide complexes that were treated, it was found that electrostatic interactions oppose binding whereas nonpolar interactions drive complex formation. The two types of interactions appear to be correlated in that larger nonpolar contributions to binding are generally opposed by increased electrostatic contributions favoring dissociation. The factors that drive the binding of peptides to MHC proteins are discussed in light of our results.     

Quantum calculations on water in the KcsA channel cavity with permeant and non-permeant ions

Different ions in the pore of the KcsA channel behave differently, and we relate this to their solvation. We show that the selectivity is dependent, in part, on the solvation in the cavity (sometimes referred to as the vestibule, it is the region containing water molecules between the intracellular gate and the selectivity filter at the extracellular end of the pore). We have shown earlier that potassium is more dependent at the upper end of the cavity region on solvation by the threonines there, while sodium ion has more water molecules as ligands. In addition, sodium ion is placed asymmetrically, while potassium is nearly exactly symmetric with respect to the four-fold symmetry of the channel. We have now extended these calculations to rubidium and cesium ions, and find that rubidium solvation resembles that of potassium (and both are permeant ions), while cesium resembles sodium (and both are non-permeant), in terms of the geometry of up to eight hydrating, and four non-hydrating, water molecules. In each case, a maximum of 12 water molecules are relevant to the calculation. The placement of the water molecules in the two cases is essentially the same as found from the electron density in the X-ray structure of Zhou and MacKinnon. For Na⁺ and K⁺, we show that energy decreases from bulk to the cavity to the lowest position in the selectivity filter (accurate energy could not be calculated for the heavier ions). A separate calculation shows that fixing the Na⁺ ion at the position of the K⁺ minimum, followed by re-optimization produced a significantly modified system, not something that could be produced by thermal fluctuations. Moving the K⁺ into the Na⁺ position in the upper cavity led to a small increase in energy, ≈ 3 k_BT, but was accompanied by large shifts in the positions of hydrating waters, which would create a major kinetic barrier. Therefore, thermal fluctuations could not invalidate the conclusions of the main calculations.     

Theoretical evaluation of adiabatic and vertical electron affinity of some radiosensitizers in solution using FEP, ab initio and DFT methods

The biological activity of radiosensitizers is associated to their electron affinity (EA), which can be divided in two main processes: vertical and adiabatic. In this work, we calculated the EAs of nitrofurans and nitroimidazoles (Fig. 2) using Hartree–Fock (HF) and density functional theory (DFT) methods and evaluated solvent effects (water and carbon tetrachloride) on EAs. For water, we combined the polarized continuum model (PCM) and free energy perturbation (FEP) (finite difference thermodynamic integration, FDTI) methods. For carbon tetrachloride, we used the FDTI method. The values of adiabatic EA obtained are in agreement with experimental data (deviations of 0.013 eV). The vertical EAs were calculated according to Cederbaum's outer valence Green function (OVGF) method. This methodology, which relies on theoretical aspects of free energy calculations on charged molecules in solution, was used to select potential selective radiosensitizers from recently reported compounds and could be helpful in the rational design of new and more selective bioreductive anticancer drugs.     

Transesterification thio effects of phosphate diesters: free energy barriers and kinetic and equilibrium isotope effects from density-functional theory

Primary and secondary kinetic and equilibrium isotope effects are calculated with density-functional methods for the in-line dianionic methanolysis of the native (unsubstituted) and thio-substituted cyclic phosphates. These reactions represent reverse reaction models for RNA transesterification under alkaline conditions. The effect of solvent is treated with explicit (single and double) water molecules and self-consistently with an implicit (continuum) solvation model. Singly substituted reactions at the nonbridging O(P1) position and bridging O(2)('), O(3)('), and O(5)(') positions and a doubly substituted reaction at the O(P1) and O(P2) positions were considered. Aqueous free energy barriers are calculated, and the structures and bond orders of the rate-controlling transition states are characterized. The results are consistent with available experimental data and provide useful information for the interpretation of measured isotope and thio effects used to probe mechanism in phosphoryl transfer reactions catalyzed by enzymes and ribozymes.