Simple electrolyte solutions: comparison of DRISM and molecular dynamics results for alkali halide solutions

Unstructured linker regions play a role in the differential splicing activities of paralogous RNA binding proteins PTBP1 and PTBP2

RNA Binding Proteins regulate, in part, alternative pre-mRNA splicing and, in turn, gene expression patterns. Polypyrimidine tract binding proteins PTBP1 and PTBP2 are paralogous RNA binding proteins sharing 74% amino acid sequence identity. Both proteins contain four structured RNA-recognition motifs (RRMs) connected by linker regions and an N-terminal region. Despite their similarities, the paralogs have distinct tissue-specific expression patterns and can regulate discrete sets of target exons. How two highly structurally similar proteins can exert different splicing outcomes is not well understood. Previous studies revealed that PTBP2 is post-translationally phosphorylated in the unstructured N-terminal, Linker 1, and Linker 2 regions that share less sequence identity with PTBP1 signifying a role for these regions in dictating the paralog's distinct splicing activities. To this end, we conducted bioinformatics analysis to determine the evolutionary conservation of RRMs versus linker regions in PTBP1 and PTBP2 across species. To determine the role of PTBP2 unstructured regions in splicing activity, we created hybrid PTBP1-PTBP2 constructs that had counterpart PTBP1 regions swapped to an otherwise PTBP2 protein and assayed on differentially regulated exons. We also conducted molecular dynamics studies to investigate how negative charges introduced by phosphorylation in PTBP2 unstructured regions can alter their physical properties. Collectively, results from our studies reveal an important role for PTBP2 unstructured regions and suggest a role for phosphorylation in the differential splicing activities of the paralogs on certain regulated exons.

Electrocatalytic Mechanisms for an Oxygen Evolution Reaction at a Rhombohedral Boron Monosulfide Electrode/Alkaline Medium Interface

Rhombohedral boron monosulfide (r-BS) with a layer stacking structure is a promising electrocatalyst for an oxygen evolution reaction (OER) within an alkaline solution. We investigated the catalytic mechanisms at the r-BS electrode/alkaline medium interface for an OER using hybrid solvation theory based on the first-principles method combined with classical solution theory. In this study, we elucidate the activities of the OER at the outermost r-BS sheet with and without various surface defects. The Gibbs free energies along the OER path indicate that the boron vacancies at the first and second layers of the r-BS surface (VB1 and VB2) can promote the OER. However, we found that the VB1 is easily occupied by the oxygen atom during the OER, degrading its electrocatalytic performance. In contrast, VB2 is suitable for the active site of the OER due to its structure stability. Next, we applied a bias voltage with the OER potential to the r-BS electrode. The bias voltage incorporates the positive excess surface charge into pristine r-BS and VB2, which can be understood by the relationship between the OER potential and potentials of zero charge at the r-BS electrode. Because the OH- ions are the starting point of the OER, the positively charged surface is kinetically favorable for the electrocatalyst owing to the attractive interaction with the OH- ions. Finally, we qualitatively discuss the flat-band potential at a semiconductor/alkaline solution interface. It suggests that p-type carrier doping could promote the catalytic performance of r-BS. These results explain the previous measurement of the OER performance with the r-BS-based electrode and provide valuable insights into developing a semiconductor electrode/water interface.

Identifying Systematic Force Field Errors Using a 3D-RISM Element Counting Correction

Hydration free energies of small molecules are commonly used as benchmarks for solvation models. However, errors in predicting hydration free energies are partially due to the force fields used and not just the solvation model. To address this, we have used the 3D reference interaction site model (3D-RISM) of molecular solvation and existing benchmark explicit solvent calculations with a simple element count correction (ECC) to identify problems with the non-bond parameters in the general AMBER force field (GAFF). 3D-RISM was used to calculate hydration free energies of all 642 molecules in the FreeSolv database, and a partial molar volume correction (PMVC), ECC, and their combination (PMVECC) were applied to the results. The PMVECC produced a mean unsigned error of 1.01±0.04kcal/mol and root mean squared error of 1.44±0.07kcal/mol, better than the benchmark explicit solvent calculations from FreeSolv, and required less than 15 s of computing time per molecule on a single CPU core. Importantly, parameters for PMVECC showed systematic errors for molecules containing Cl, Br, I, and P. Applying ECC to the explicit solvent hydration free energies found the same systematic errors. The results strongly suggest that some small adjustments to the Lennard-Jones parameters for GAFF will lead to improved hydration free energy calculations for all solvent models.

The Ion-Dipole Correction of the 3DRISM Solvation Model to Accurately Compute Water Distributions around Negatively Charged Biomolecules

The 3D reference interaction site model (3DRISM) provides an efficient grid-based solvation model to compute the structural and thermodynamic properties of biomolecules in aqueous solutions. However, it remains challenging for existing 3DRISM methods to correctly predict water distributions around negatively charged solute molecules. In this paper, we first show that this challenge is mainly due to the orientation of water molecules in the first solvation shell of the negatively charged solute molecules. To properly consider this orientational preference, position-dependent two-body intramolecular correlations of solvent need to be included in the 3DRISM theory, but direct evaluations of these position-dependent two-body intramolecular correlations remain numerically intractable. To address this challenge, we introduce the Ion-Dipole Correction (IDC) to the 3DRISM theory, in which we incorporate the orientation preference of water molecules via an additional solute-solvent interaction term (i.e., the ion-dipole interaction) while keeping the formulism of the 3DRISM equation unchanged. We prove that this newly introduced IDC term is equivalent to an effective direct correlation function which can effectively consider the orientation effect that arises from position dependent two-body correlations. We first quantitatively validate our 3DRISM-IDC theory combined with the PSE3 closure on Cl^-, [ClO]^- (a two-site anion), and [NO₂]^- (a three-site anion). For all three anions, we show that our 3DRISM-IDC theory significantly outperforms the 3DRISM theory in accurately predicting the solvation structures in comparison to MD simulations, including RDFs and 3D water distributions. Furthermore, we have also demonstrated that the 3DRISM-IDC can improve the accuracy of hydration free-energy calculation for Cl^-. We further demonstrate that our 3DRISM-IDC theory yields significant improvements over the 3DRISM theory when applied to compute the solvation structures for various negatively charged solute molecules, including adenosine triphosphate (ATP), a short peptide containing 19 residues, a DNA hairpin containing 24 nucleotides, and a riboswitch RNA molecule with 77 nucleotides. We expect that our 3DRISM-IDC-PSE3 solvation model holds great promise to be widely applied to study solvation properties for nucleic acids and other biomolecules containing negatively charged functional groups.

Structure and thermodynamics of a 2D Lennard-Jones hexagonal fluid

The thermodynamic and structural properties of the 2D hexagonal soft-sites fluid are examined by integral equation theory benchmarked against extensive Monte Carlo simulations. Hexamers are built of six equal Lennard-Jones segments. Site-site integral equation theory is used to compute site-site correlation functions, excess internal energies and isotherms over a wide range of conditions and compared with results obtained from Monte Carlo simulations. Various approaches for computing the pressure are discussed as well. Satisfactory qualitative agreement between theory and simulations is found with details depending on the applied closure relation.

Activity Coefficients and Solubilities of NaCl in Water-Methanol Solutions from Molecular Dynamics Simulations

We obtain activity coefficients and solubilities of NaCl in water-methanol solutions at 298.15 K and 1 bar from molecular dynamics (MD) simulations with the Joung-Cheatham, SPC/E, and OPLS-AA force fields for NaCl, water, and methanol, respectively. The Lorentz-Berthelot combining rules were adopted for the unlike-pair interactions of Na⁺, Cl^-, and the oxygen site in SPC/E water, and geometric combining rules were utilized for the remainder of the cross interactions. We found that the selection of appropriate combining rules is important in obtaining physically realistic solubilities. The solvent compositions studied range from pure water to pure methanol. Several salt concentrations were investigated at each solvent composition, from the lowest concentrations permitted by the system size used up to the experimental solubilities. We first calculated individual ion activity coefficients (IIACs) for Na⁺ and Cl^- from the free energy change due to the gradual insertion of a single cation or anion into the solution, accompanied by a neutralizing background. We obtained the salt solubilities by comparing the chemical potentials in solution with solid NaCl chemical potentials calculated previously using the Einstein crystal method. Mean ionic activity coefficients obtained from the IIACs are in reasonable agreement with experimental data, with deviations increasing for solutions of higher methanol content. Predictions for the salt solubility are in surprisingly good agreement with experimental data, despite well-known challenges in the simultaneous calculation of activity coefficients and solubilities with classical MD simulations. The present study demonstrates that good predictions for these two important phase equilibrium properties can be obtained for mixed-solvent electrolyte solutions using existing nonpolarizable models and further suggests that the previously proposed single ion insertion technique can be extended to complex mixed-solvent solutions as well.

Ion Pairing of the Neurotransmitters Acetylcholine and Glutamate in Aqueous Solutions

Neurotransmitters (NTs) play an important role in neural communication, regulating a variety of functions such as motivation, learning, memory, and muscle contraction. Their intermolecular interactions in biological media are an important factor affecting their biological activity. However, the available information on the features of these interactions is scarce and contradictory, especially, in an estimation of possible ion binding. In this paper, we present the results of a study for two well-known NTs, acetylcholine (ACh) and glutamate (Glu), with relation to the NT-inorganic ion and the NT-NT binding in a water environment. The features of NT pairing are investigated in aqueous AChCl and NaGlu solutions over a wide concentration range using the integral equation method in 1D- and 3D- reference interaction site model (RISM) approaches. The data for ACh are given for its two bioactive TG (trans, gauche) and TT (trans, trans) conformers. As was found, for both NTs, the results indicate the NT-inorganic counterion contact pair to be the predominant associate type in the concentrated solutions. In this case, the counterions occupy the vacated "water" space in the hydration shell of the onium moiety (ACh) or carboxylate groups (Glu). For ACh, the "unfolded" TT conformer demonstrates a slightly greater possibility for counterion pairing in comparison with the "folded" TG conformer. For Glu, the probability of its binding with a counterion is slightly stronger for the "side-chain" carboxylate group than for the "backbone" group. The obtained results also revealed an insignificant probability of Glu^--Glu^- pairing. Namely, the RISM data indicate Glu^--Glu^- binding by NH₃⁺-COO^- interactions. A link between the ion binding of NTs and their biological activity is discussed. This contribution adds new knowledge to our understanding of the interactions between the NTs and their molecular environment, providing further insights into the behavior of these compounds in biological media.

Renormalized site density functional theory for models of ion hydration

The development of accurate statistical mechanics models of molecular liquid systems is a problem of great practical and fundamental importance. Site-density functional theory (SDFT) is one of the promising directions in this area, but its success hinges upon the ability to efficiently reconcile the co-existence of two distinct intra- and inter-molecular interaction regimes in a molecular liquid. The renormalized formulation of SDFT (RSDFT), which we have recently developed, resolves this problem by introducing an additional potential field variable that decouples two interaction scales and maps the molecular liquid problem onto the effective simple liquid mixture. This work provides a critical assessment of RSDFT for the hydrated ion system-a problem that historically has always been one of the most difficult cases for SDFT applications. Using a two-site model of water, we perform a comprehensive analysis of hydrated alkali metal and halogen ions, including both structural and free energy based characteristics. The results indicate that RSDFT provides a significant improvement over conventional three-dimensional reference interaction site model implementations and may prove useful in coarse grained simulations based on two-site solvent models.

Structural and dynamic properties of some aqueous salt solutions

Aqueous salt solutions are utilized and encountered in wide-ranging technological applications and natural settings. Towards improved understanding of the effect of salts on the dynamic properties of such systems, dilute aqueous salt solutions (up to 1 molar concentration) are investigated here, via experiments and molecular simulations. Four salts are considered: sodium chloride, for which published results are readily available for comparison, ammonium acetate, barium acetate and barium nitrate, for which published data are scarce. In the present work, molecular dynamics (MD) simulations are conducted to quantify viscosity and water self-diffusion coefficients, together with rheometry and Pulsed Field Gradient Spin Echo (PFGSE)-NMR experiments for validation. Simulation predictions are consistent with experimental observations in terms of trend and magnitude of salt-specific effects. Combining insights from the approaches considered, an interpretation of the results is proposed whereby the capacity of salts to influence bulk dynamics arises from their molecular interfacial area and strength of interaction with first hydration-shell water molecules. For the concentration range investigated, the interpretation could be useful in formulating aqueous systems for applications including the manufacturing of advanced catalysts.

Individual Ion Activity Coefficients in Aqueous Electrolytes from Explicit-Water Molecular Dynamics Simulations

We compute individual ion activity coefficients (IIACs) in aqueous NaCl, KCl, NaF, and KF solutions from explicit-water molecular dynamics simulations. Free energy changes are obtained from insertion of single ions-accompanied by uniform neutralizing backgrounds-into solution by gradually turning on first Lennard-Jones interactions, followed by Coulombic interactions using Ewald electrostatics. Simulations are performed at multiple system sizes, and all results are extrapolated to the thermodynamic limit, thus eliminating any possible artifacts from the neutralizing backgrounds. Because of controversies associated with measurements of IIACs from electrochemical cells with ion-selective electrodes, the reported experimental data are not widely accepted; thus there remains a knowledge gap with respect to the contributions of individual ions to solution nonidealities. Our results are in good qualitative agreement with these reported measurements, though significantly larger in magnitude. In particular, the relative positioning for the activity coefficients of anions and cations matches the experimental ordering for all four systems. This work establishes a robust thermodynamic framework, without a need to invoke extra hypotheses, that sheds light on the behavior of individual ions and their contributions to nonidealities of aqueous electrolyte solutions.

Failure of molecular dynamics to provide appropriate structures for quantum mechanical description of the aqueous chloride ion charge-transfer-to-solvent ultraviolet spectrum

The lowest band in the charge-transfer-to-solvent ultraviolet absorption spectrum of aqueous chloride ion is studied by experiment and computation. Interestingly, the experiments indicate that at concentrations up to at least 0.25 M, where calculations indicate ion pairing to be significant, there is no notable effect of ionic strength on the spectrum. The experimental spectra are fitted to aid comparison with computations. Classical molecular dynamic simulations are carried out on dilute aqueous Cl-, Na+, and NaCl, producing radial distribution functions in reasonable agreement with experiment and, for NaCl, clear evidence of ion pairing. Clusters are extracted from the simulations for quantum mechanical excited state calculations. Accurate ab initio coupled-cluster benchmark calculations on a small number of representative clusters are carried out and used to identify and validate an efficient protocol based on time-dependent density functional theory. The latter is used to carry out quantum mechanical calculations on thousands of clusters. The resulting computed spectrum is in excellent agreement with experiment for the peak position, with little influence from ion pairing, but is in qualitative disagreement on the width, being only about half as wide. It is concluded that simulation by classical molecular dynamics fails to provide an adequate variety of structures to explain the experimental CTTS spectrum of aqueous Cl-.

Continuum Electrostatic Behavior of a 3D-RISM Theory

The electrostatic response underlying the 3D-RISM theory and its general relationship to models in which the solvent is represented in terms of a dielectric continuum are examined. It is found that the theory provides a coherent picture of solvation, although its behavior is not entirely consistent with the trends that are expected in the limit of a large solute. The electrostatic discrepancy is due to the nature of the isotropic pair additive site-site correlation function associated with the susceptibility response of the uniform fluid. The influence of the discrepancy in the magnitude of the solvation free energy is negligible for a solvent with a large dielectric constant.

Theory of electrolytes including steric, attractive, and hydration interactions

We present a continuum theory of electrolytes composed of a waterlike solvent and univalent ions. First, we start with a density functional F for the coarse-grained solvent, cation, and anion densities, including the Debye-Hückel free energy, the Coulombic interaction, and the direct interactions among these three components. These densities fluctuate obeying the distribution ∝exp(-F/k_BT). Eliminating the solvent density deviation in F, we obtain the effective non-Coulombic interactions among the ions, which consist of the direct ones and the solvent-mediated ones. We then derive general expressions for the ion correlation, the apparent partial volume, and the activity and osmotic coefficients up to linear order in the average salt density n_s. Second, we perform numerical analysis using the Mansoori-Carnahan-Starling-Leland model [J. Chem. Phys. 54, 1523 (1971)] for three-component hardspheres. The effective interactions sensitively depend on the cation and anion sizes due to competition between the steric and hydration effects, which are repulsive between small-large ion pairs and attractive between symmetric pairs. These agree with previous experiments and Collins' rule [Biophys. J. 72, 65 (1997)]. We also give simple approximate expressions for the ionic interaction coefficients valid for any ion sizes.

A molecular reconstruction approach to site-based 3D-RISM and comparison to GIST hydration thermodynamic maps in an enzyme active site

Computed, high-resolution, spatial distributions of solvation energy and entropy can provide detailed information about the role of water in molecular recognition. While grid inhomogeneous solvation theory (GIST) provides rigorous, detailed thermodynamic information from explicit solvent molecular dynamics simulations, recent developments in the 3D reference interaction site model (3D-RISM) theory allow many of the same quantities to be calculated in a fraction of the time. However, 3D-RISM produces atomic-site, rather than molecular, density distributions, which are difficult to extract physical meaning from. To overcome this difficulty, we introduce a method to reconstruct molecular density distributions from atomic-site density distributions. Furthermore, we assess the quality of the resulting solvation thermodynamics density distributions by analyzing the binding site of coagulation Factor Xa with both GIST and 3D-RISM. We find good qualitative agreement between the methods for oxygen and hydrogen densities as well as direct solute-solvent energetic interactions. However, 3D-RISM predicts lower energetic and entropic penalties for moving water from the bulk to the binding site.

Closure for the Ornstein-Zernike equation with pressure and free energy consistency

The Ornstein-Zernike (OZ) integral equation theory is a powerful approach to simple liquids due to its low computational cost and the fact that, when combined with an appropriate closure equation, the theory is thermodynamically complete. However, approximate closures proposed to date exhibit pressure or free energy inconsistencies that produce inaccurate or ambiguous results, limiting the usefulness of the Ornstein-Zernike approach. To address this problem, we combine methods to enforce both pressure and free energy consistency to create a new closure approximation and test it for a single-component Lennard-Jones fluid. The closure is a simple power series in the direct and total correlation functions for which we have derived analytical formulas for the excess Helmholtz free energy and chemical potential. These expressions contain a partial molar volumelike term, similar to excess chemical potential correction terms recently developed. Using our bridge approximation, we have calculated the pressure, Helmholtz free energy, and chemical potential for the Lennard-Jones fluid using the Kirkwood charging, thermodynamic integration techniques, and analytic expressions. These results are compared with those from the hypernetted chain equation and the Verlet-modified closure against Monte Carlo and equations-of-state data for reduced densities of ρ^{*}<1 and temperatures of T^{*}=1.5, 2.74, and 5. Our closure shows consistency among all thermodynamic paths, except for one expression of the Gibbs-Duhem relation, whereas the hypernetted chain equation and the Verlet-modified closure exhibit consistency between only a few relations. Accuracy of the closure is comparable to the Verlet-modified closure and a significant improvement to results obtained from the hypernetted chain equation.

Predicting Site-Binding Modes of Ions and Water to Nucleic Acids Using Molecular Solvation Theory

Site binding of ions and water shapes nucleic acids folding, dynamics, and biological function, complementing the more diffuse, nonspecific "territorial" ion binding. Unlike territorial binding, prediction of site-specific binding to nucleic acids remains an unsolved challenge in computational biophysics. This work presents a new toolset based on the 3D-RISM molecular solvation theory and topological analysis that predicts cation and water site binding to nucleic acids. 3D-RISM is shown to accurately capture alkali cations and water binding to the central channel, transversal loops, and grooves of the Oxytricha nova's telomeres' G-quadruplex ( Oxy-GQ), in agreement with high-resolution crystallographic data. To improve the computed cation occupancy along the Oxy-GQ central channel, it was necessary to refine and validate new cation-oxygen parameters using structural and thermodynamic data available for crown ethers and ion channels. This single set of parameters that describes both localized and delocalized binding to various biological systems is used to gain insight into cation occupancy along the Oxy-GQ channel under various salt conditions. The paper concludes with prospects for extending the method to predict divalent cation binding to nucleic acids. This work advances the forefront of theoretical methods able to provide predictive insight into ion atmosphere effects on nucleic acids function.

Computational characterization of the selective inhibition of human norepinephrine and serotonin transporters by an escitalopram scaffold

Selective inhibition of human norepinephrine and serotonin transporters has been studied by computational approaches. 4 warm spots in hNET and 4 in hSERT were found to exert a pronounced effect on inhibition by the studied ligands.

How Water's Properties Are Encoded in Its Molecular Structure and Energies

How are water's material properties encoded within the structure of the water molecule? This is pertinent to understanding Earth's living systems, its materials, its geochemistry and geophysics, and a broad spectrum of its industrial chemistry. Water has distinctive liquid and solid properties: It is highly cohesive. It has volumetric anomalies-water's solid (ice) floats on its liquid; pressure can melt the solid rather than freezing the liquid; heating can shrink the liquid. It has more solid phases than other materials. Its supercooled liquid has divergent thermodynamic response functions. Its glassy state is neither fragile nor strong. Its component ions-hydroxide and protons-diffuse much faster than other ions. Aqueous solvation of ions or oils entails large entropies and heat capacities. We review how these properties are encoded within water's molecular structure and energies, as understood from theories, simulations, and experiments. Like simpler liquids, water molecules are nearly spherical and interact with each other through van der Waals forces. Unlike simpler liquids, water's orientation-dependent hydrogen bonding leads to open tetrahedral cage-like structuring that contributes to its remarkable volumetric and thermal properties.

Extracting water and ion distributions from solution x-ray scattering experiments

Small-angle X-ray scattering measurements can provide valuable information about the solvent environment around biomolecules, but it can be difficult to extract solvent-specific information from observed intensity profiles. Intensities are proportional to the square of scattering amplitudes, which are complex quantities. Amplitudes in the forward direction are real, and the contribution from a solute of known structure (and from the waters it excludes) can be estimated from theory; hence, the amplitude arising from the solvent environment can be computed by difference. We have found that this "square root subtraction scheme" can be extended to non-zero q values, out to 0.1 Å(-1) for the systems considered here, since the phases arising from the solute and from the water environment are nearly identical in this angle range. This allows us to extract aspects of the water and ion distributions (beyond their total numbers), by combining experimental data for the complete system with calculations for the solutes. We use this approach to test molecular dynamics and integral-equation (3D-RISM (three-dimensional reference interaction site model)) models for solvent structure around myoglobin, lysozyme, and a 25 base-pair duplex DNA. Comparisons can be made both in Fourier space and in terms of the distribution of interatomic distances in real space. Generally, computed solvent distributions arising from the MD simulations fit experimental data better than those from 3D-RISM, even though the total small-angle X-ray scattering patterns are very similar; this illustrates the potential power of this sort of analysis to guide the development of computational models.

Potential of Mean Force Calculations of Solute Permeation Across UT-B and AQP1: A Comparison between Molecular Dynamics and 3D-RISM

Membrane channels facilitate the efficient and selective flux of various solutes across biological membranes. A common approach to investigate the selectivity of a channel has been the calculation of potentials of mean force (PMFs) for solute permeation across the pore. PMFs have been frequently computed from molecular dynamics (MD) simulations, yet the three-dimensional reference interaction site model (3D-RISM) has been suggested as a computationally efficient alternative to MD. Whether the two methods yield comparable PMFs for solute permeation has remained unclear. In this study, we calculated potentials of mean force for water, ammonia, urea, molecular oxygen, and methanol across the urea transporter B (UT-B) and aquaporin-1 (AQP1), using 3D-RISM, as well as using MD simulations and umbrella sampling. To allow direct comparison between the PMFs from 3D-RISM and MD, we ensure that all PMFs refer to a well-defined reference area in the bulk or, equivalently, to a well-defined density of channels in the membrane. For PMFs of water permeation, we found reasonable agreement between the two methods, with differences of ≲3 kJ mol-1. In contrast, we found stark discrepancies for the PMFs for all other solutes. Additional calculations confirm that discrepancies between MD and 3D-RISM are not explained by the choice for the closure relation, the definition the reaction coordinate (center of mass-based versus atomic site-based), details of the molecule force field, or fluctuations of the protein. Comparison of the PMFs suggests that 3D-RISM may underestimate effects from hydrophobic solute-channel interactions, thereby, for instance, missing the urea binding sites in UT-B. Furthermore, we speculate that the orientational averages inherent to 3D-RISM might lead to discrepancies in the narrow channel lumen. These findings suggest that current 3D-RISM solvers provide reasonable estimates for the PMF for water permeation, but that they are not suitable to study the selectivity of membrane channels with respect to uncharged nonwater solutes.

Emergence of Solvent-Separated Na+-Cl- Ion Pair in Salt Water: Photoelectron Spectroscopy and Theoretical Calculations

Solvation of salts in water is a fundamental physical chemical process, but the underlying mechanism remains unclear. We investigated the contact ion pair (CIP) to solvent-separated ion pair (SSIP) transition in NaCl(H₂O)_n clusters with anion photoelectron spectroscopy and ab initio calculations. It is found that the SSIP type of structures show up at n = 2 for NaCl^-(H₂O)_n anions. For neutral NaCl(H₂O)_n, the CIP structures are dominant at n < 9. At n = 9-12, the CIP structures and SSIP structures of NaCl(H₂O)_n are nearly degenerate in energy, coincident to the H₂O:NaCl molar ratio of NaCl saturated solution and implying that the CIP and SSIP structures can coexist in concentrated solutions. These results are useful for understanding the solvation of salts at the molecular level.

A closure relation to molecular theory of solvation for macromolecules

We propose a closure to the integral equations of molecular theory of solvation, particularly suitable for polar and charged macromolecules in electrolyte solution. This includes such systems as oligomeric polyelectrolytes at a finite concentration in aqueous and various non-aqueous solutions, as well as drug-like compounds in solution. The new closure by Kobryn, Gusarov, and Kovalenko (KGK closure) imposes the mean spherical approximation (MSA) almost everywhere in the solvation shell but levels out the density distribution function to zero (with the continuity at joint boundaries) inside the repulsive core and in the spatial regions of strong density depletion emerging due to molecular associative interactions. Similarly to MSA, the KGK closure reduces the problem to a linear equation for the direct correlation function which is predefined analytically on most of the solvation shells and has to be determined numerically on a relatively small (three-dimensional) domain of strong depletion, typically within the repulsive core. The KGK closure leads to the solvation free energy in the form of the Gaussian fluctuation (GF) functional. We first test the performance of the KGK closure coupled to the reference interaction site model (RISM) integral equations on the examples of Lennard-Jones liquids, polar and nonpolar molecular solvents, including water, and aqueous solutions of simple ions. The solvation structure, solvation chemical potential, and compressibility obtained from RISM with the KGK closure favorably compare to the results of the hypernetted chain (HNC) and Kovalenko-Hirata (KH) closures, including their combination with the GF solvation free energy. We then use the KGK closure coupled to RISM to obtain the solvation structure and thermodynamics of oligomeric polyelectrolytes and drug-like compounds at a finite concentration in electrolyte solution, for which no convergence is obtained with other closures. For comparison, we calculate their solvation structure from molecular dynamics (MD) simulations. We further couple the 3D-RISM integral equation with the 3D-version of the KGK closure, and solve it for molecular mixtures as well as oligomeric polyelectrolytes and drug-like molecules in electrolyte solutions.

Designing molecular complexes using free-energy derivatives from liquid-state integral equation theory

Complex formation between molecules in solution is the key process by which molecular interactions are translated into functional systems. These processes are governed by the binding or free energy of association which depends on both direct molecular interactions and the solvation contribution. A design goal frequently addressed in pharmaceutical sciences is the optimization of chemical properties of the complex partners in the sense of minimizing their binding free energy with respect to a change in chemical structure. Here, we demonstrate that liquid-state theory in the form of the solute-solute equation of the reference interaction site model provides all necessary information for such a task with high efficiency. In particular, computing derivatives of the potential of mean force (PMF), which defines the free-energy surface of complex formation, with respect to potential parameters can be viewed as a means to define a direction in chemical space toward better binders. We illustrate the methodology in the benchmark case of alkali ion binding to the crown ether 18-crown-6 in aqueous solution. In order to examine the validity of the underlying solute-solute theory, we first compare PMFs computed by different approaches, including explicit free-energy molecular dynamics simulations as a reference. Predictions of an optimally binding ion radius based on free-energy derivatives are then shown to yield consistent results for different ion parameter sets and to compare well with earlier, orders-of-magnitude more costly explicit simulation results. This proof-of-principle study, therefore, demonstrates the potential of liquid-state theory for molecular design problems.

Small molecule hydration energy and entropy from 3D-RISM

Implicit solvent models offer an attractive way to estimate the effects of a solvent environment on the properties of small or large solutes without the complications of explicit simulations. One common test of accuracy is to compute the free energy of transfer from gas to liquid for a variety of small molecules, since many of these values have been measured. Studies of the temperature dependence of these values (i.e. solvation enthalpies and entropies) can provide additional insights into the performance of implicit solvent models. Here, we show how to compute temperature derivatives of hydration free energies for the 3D-RISM integral equation approach. We have computed hydration free energies of 1123 small drug-like molecules (both neutral and charged). Temperature derivatives were also used to calculate hydration energies and entropies of 74 of these molecules (both neutral and charged) for which experimental data is available. While direct results have rather poor agreement with experiment, we have found that several previously proposed linear hydration free energy correction schemes give good agreement with experiment. These corrections also provide good agreement for hydration energies and entropies though simple extensions are required in some cases.

SAMPL5: 3D-RISM partition coefficient calculations with partial molar volume corrections and solute conformational sampling

Implicit solvent methods for classical molecular modeling are frequently used to provide fast, physics-based hydration free energies of macromolecules. Less commonly considered is the transferability of these methods to other solvents. The Statistical Assessment of Modeling of Proteins and Ligands 5 (SAMPL5) distribution coefficient dataset and the accompanying explicit solvent partition coefficient reference calculations provide a direct test of solvent model transferability. Here we use the 3D reference interaction site model (3D-RISM) statistical-mechanical solvation theory, with a well tested water model and a new united atom cyclohexane model, to calculate partition coefficients for the SAMPL5 dataset. The cyclohexane model performed well in training and testing ([Formula: see text] for amino acid neutral side chain analogues) but only if a parameterized solvation free energy correction was used. In contrast, the same protocol, using single solute conformations, performed poorly on the SAMPL5 dataset, obtaining [Formula: see text] compared to the reference partition coefficients, likely due to the much larger solute sizes. Including solute conformational sampling through molecular dynamics coupled with 3D-RISM (MD/3D-RISM) improved agreement with the reference calculation to [Formula: see text]. Since our initial calculations only considered partition coefficients and not distribution coefficients, solute sampling provided little benefit comparing against experiment, where ionized and tautomer states are more important. Applying a simple [Formula: see text] correction improved agreement with experiment from [Formula: see text] to [Formula: see text], despite a small number of outliers. Better agreement is possible by accounting for tautomers and improving the ionization correction.

Hydration Free Energies of Molecular Ions from Theory and Simulation

We present a theoretical/computational framework for accurate calculation of hydration free energies of ionized molecular species. The method is based on a molecular theory, 3D-RISM, combined with a recently developed pressure correction (PC+). The 3D-RISM/PC+ model can provide ∼3 kcal/mol hydration free energy accuracy for a large variety of ionic compounds, provided that the Galvani potential of water is taken into account. The results are compared with direct atomistic simulations. Several methodological aspects of hydration free energy calculations for charged species are discussed.

Competitive interaction of monovalent cations with DNA from 3D-RISM

The composition of the ion atmosphere surrounding nucleic acids affects their folding, condensation and binding to other molecules. It is thus of fundamental importance to gain predictive insight into the formation of the ion atmosphere and thermodynamic consequences when varying ionic conditions. An early step toward this goal is to benchmark computational models against quantitative experimental measurements. Herein, we test the ability of the three dimensional reference interaction site model (3D-RISM) to reproduce preferential interaction parameters determined from ion counting (IC) experiments for mixed alkali chlorides and dsDNA. Calculations agree well with experiment with slight deviations for salt concentrations >200 mM and capture the observed trend where the extent of cation accumulation around the DNA varies inversely with its ionic size. Ion distributions indicate that the smaller, more competitive cations accumulate to a greater extent near the phosphoryl groups, penetrating deeper into the grooves. In accord with experiment, calculated IC profiles do not vary with sequence, although the predicted ion distributions in the grooves are sequence and ion size dependent. Calculations on other nucleic acid conformations predict that the variation in linear charge density has a minor effect on the extent of cation competition.

Assessment of metal-assisted nucleophile activation in the hepatitis delta virus ribozyme from molecular simulation and 3D-RISM

The hepatitis delta virus ribozyme is an efficient catalyst of RNA 2'-O-transphosphorylation and has emerged as a key experimental system for identifying and characterizing fundamental features of RNA catalysis. Recent structural and biochemical data have led to a proposed mechanistic model whereby an active site Mg(2+) ion facilitates deprotonation of the O2' nucleophile, and a protonated cytosine residue (C75) acts as an acid to donate a proton to the O5' leaving group as noted in a previous study. This model assumes that the active site Mg(2+) ion forms an inner-sphere coordination with the O2' nucleophile and a nonbridging oxygen of the scissile phosphate. These contacts, however, are not fully resolved in the crystal structure, and biochemical data are not able to unambiguously exclude other mechanistic models. In order to explore the feasibility of this model, we exhaustively mapped the free energy surfaces with different active site ion occupancies via quantum mechanical/molecular mechanical (QM/MM) simulations. We further incorporate a three-dimensional reference interaction site model for the solvated ion atmosphere that allows these calculations to consider not only the rate associated with the chemical steps, but also the probability of observing the system in the presumed active state with the Mg(2+) ion bound. The QM/MM results predict that a pathway involving metal-assisted nucleophile activation is feasible based on the rate-controlling transition state barrier departing from the presumed metal-bound active state. However, QM/MM results for a similar pathway in the absence of Mg(2+) are not consistent with experimental data, suggesting that a structural model in which the crystallographically determined Mg(2+) is simply replaced with Na(+) is likely incorrect. It should be emphasized, however, that these results hinge upon the assumption of the validity of the presumed Mg(2+)-bound starting state, which has not yet been definitively verified experimentally, nor explored in depth computationally. Thus, further experimental and theoretical study is needed such that a consensus view of the catalytic mechanism emerges.

Uranyl Solvation by a Three-Dimensional Reference Interaction Site Model

We report an implementation of the three-dimensional reference interaction site model (3D RISM) that in particular addresses the treatment of the long-range Coulomb field of charged species, represented by point charges and/or a distributed charge density. A comparison of 1D and 3D results for atomic ions demonstrates a reasonable accuracy, even for a moderate size of the unit cell and a moderate grid resolution. In an application to uranyl complexes with 4-6 explicit aqua ligands and an implicit bulk solvent modeled by RISM, we show that the 3D technique is not susceptible to the deficiencies of the 1D technique exposed in our previous work [Li, Matveev, Krüger, Rösch, Comp. Theor. Chem. 2015, 1051, 151]. The 3D method eliminates the artificial superposition of explicit aqua ligands and the RISM medium and predicts essentially the same values for uranyl and uranyl-water bond lengths as a state-of-the-art polarizable continuum model. With the first solvation shell treated explicitly, the observables are nearly independent of the order of the closure relationship used when solving the set of integral equations for the various distribution functions. Furthermore, we calculated the activation barrier of water exchange with a hybrid approach that combines the 3D RISM model for the bulk aqueous solvent and a quantum mechanical description (at the level of electronic density functional theory) of uranyl interacting with explicitly represented water molecules. The calculated result agrees very well with experiment and the best theoretical estimates.

Accurate small and wide angle x-ray scattering profiles from atomic models of proteins and nucleic acids

A new method is introduced to compute X-ray solution scattering profiles from atomic models of macromolecules. The three-dimensional version of the Reference Interaction Site Model (RISM) from liquid-state statistical mechanics is employed to compute the solvent distribution around the solute, including both water and ions. X-ray scattering profiles are computed from this distribution together with the solute geometry. We describe an efficient procedure for performing this calculation employing a Lebedev grid for the angular averaging. The intensity profiles (which involve no adjustable parameters) match experiment and molecular dynamics simulations up to wide angle for two proteins (lysozyme and myoglobin) in water, as well as the small-angle profiles for a dozen biomolecules taken from the BioIsis.net database. The RISM model is especially well-suited for studies of nucleic acids in salt solution. Use of fiber-diffraction models for the structure of duplex DNA in solution yields close agreement with the observed scattering profiles in both the small and wide angle scattering (SAXS and WAXS) regimes. In addition, computed profiles of anomalous SAXS signals (for Rb(+) and Sr(2+)) emphasize the ionic contribution to scattering and are in reasonable agreement with experiment. In cases where an absolute calibration of the experimental data at q = 0 is available, one can extract a count of the excess number of waters and ions; computed values depend on the closure that is assumed in the solution of the Ornstein-Zernike equations, with results from the Kovalenko-Hirata closure being closest to experiment for the cases studied here.

3D RISM theory with fast reciprocal-space electrostatics

The calculation of electrostatic solute-solvent interactions in 3D RISM ("three-dimensional reference interaction site model") integral equation theory is recast in a form that allows for a computational treatment analogous to the "particle-mesh Ewald" formalism as used for molecular simulations. In addition, relations that connect 3D RISM correlation functions and interaction potentials with thermodynamic quantities such as the chemical potential and average solute-solvent interaction energy are reformulated in a way that calculations of expensive real-space electrostatic terms on the 3D grid are completely avoided. These methodical enhancements allow for both, a significant speedup particularly for large solute systems and a smoother convergence of predicted thermodynamic quantities with respect to box size, as illustrated for several benchmark systems.

The effect of concentration- and temperature-dependent dielectric constant on the activity coefficient of NaCl electrolyte solutions

Our implicit-solvent model for the estimation of the excess chemical potential (or, equivalently, the activity coefficient) of electrolytes is based on using a dielectric constant that depends on the thermodynamic state, namely, the temperature and concentration of the electrolyte, ε(c, T). As a consequence, the excess chemical potential is split into two terms corresponding to ion-ion (II) and ion-water (IW) interactions. The II term is obtained from computer simulation using the Primitive Model of electrolytes, while the IW term is estimated from the Born treatment. In our previous work [J. Vincze, M. Valiskó, and D. Boda, "The nonmonotonic concentration dependence of the mean activity coefficient of electrolytes is a result of a balance between solvation and ion-ion correlations," J. Chem. Phys. 133, 154507 (2010)], we showed that the nonmonotonic concentration dependence of the activity coefficient can be reproduced qualitatively with this II+IW model without using any adjustable parameter. The Pauling radii were used in the calculation of the II term, while experimental solvation free energies were used in the calculation of the IW term. In this work, we analyze the effect of the parameters (dielectric constant, ionic radii, solvation free energy) on the concentration and temperature dependence of the mean activity coefficient of NaCl. We conclude that the II+IW model can explain the experimental behavior using a concentration-dependent dielectric constant and that we do not need the artificial concept of "solvated ionic radius" assumed by earlier studies.

Multiple time step molecular dynamics in the optimized isokinetic ensemble steered with the molecular theory of solvation: accelerating with advanced extrapolation of effective solvation forces

We develop efficient handling of solvation forces in the multiscale method of multiple time step molecular dynamics (MTS-MD) of a biomolecule steered by the solvation free energy (effective solvation forces) obtained from the 3D-RISM-KH molecular theory of solvation (three-dimensional reference interaction site model complemented with the Kovalenko-Hirata closure approximation). To reduce the computational expenses, we calculate the effective solvation forces acting on the biomolecule by using advanced solvation force extrapolation (ASFE) at inner time steps while converging the 3D-RISM-KH integral equations only at large outer time steps. The idea of ASFE consists in developing a discrete non-Eckart rotational transformation of atomic coordinates that minimizes the distances between the atomic positions of the biomolecule at different time moments. The effective solvation forces for the biomolecule in a current conformation at an inner time step are then extrapolated in the transformed subspace of those at outer time steps by using a modified least square fit approach applied to a relatively small number of the best force-coordinate pairs. The latter are selected from an extended set collecting the effective solvation forces obtained from 3D-RISM-KH at outer time steps over a broad time interval. The MTS-MD integration with effective solvation forces obtained by converging 3D-RISM-KH at outer time steps and applying ASFE at inner time steps is stabilized by employing the optimized isokinetic Nosé-Hoover chain (OIN) ensemble. Compared to the previous extrapolation schemes used in combination with the Langevin thermostat, the ASFE approach substantially improves the accuracy of evaluation of effective solvation forces and in combination with the OIN thermostat enables a dramatic increase of outer time steps. We demonstrate on a fully flexible model of alanine dipeptide in aqueous solution that the MTS-MD/OIN/ASFE/3D-RISM-KH multiscale method of molecular dynamics steered by effective solvation forces allows huge outer time steps up to tens of picoseconds without affecting the equilibrium and conformational properties, and thus provides a 100- to 500-fold effective speedup in comparison to conventional MD with explicit solvent. With the statistical-mechanical 3D-RISM-KH account for effective solvation forces, the method provides efficient sampling of biomolecular processes with slow and/or rare solvation events such as conformational transitions of hydrated alanine dipeptide with the mean life times ranging from 30 ps up to 10 ns for "flip-flop" conformations, and is particularly beneficial for biomolecular systems with exchange and localization of solvent and ions, ligand binding, and molecular recognition.

Acidity in DMSO from the embedded cluster integral equation quantum solvation model

The embedded cluster reference interaction site model (EC-RISM) is applied to the prediction of acidity constants of organic molecules in dimethyl sulfoxide (DMSO) solution. EC-RISM is based on a self-consistent treatment of the solute's electronic structure and the solvent's structure by coupling quantum-chemical calculations with three-dimensional (3D) RISM integral equation theory. We compare available DMSO force fields with reference calculations obtained using the polarizable continuum model (PCM). The results are evaluated statistically using two different approaches to eliminating the proton contribution: a linear regression model and an analysis of pK(a) shifts for compound pairs. Suitable levels of theory for the integral equation methodology are benchmarked. The results are further analyzed and illustrated by visualizing solvent site distribution functions and comparing them with an aqueous environment.