Computation of methodology-independent single-ion solvation properties from molecular simulations. III. Correction terms for the solvation free energies, enthalpies, entropies, heat capacities, volumes, compressibilities, and expansivities of solvated ions

Methodological and force field effects in the molecular dynamics-based prediction of binding free energies of host-guest systems

As a contribution to the understanding and rationalization of methodological and modeling effects in recent host-guest SAMPL challenges, using an alchemical molecular dynamics technique we have examined the impact of force field parameterization and ionic strength in connection with guest charge neutralization on computed dissociation free energies in two typical SAMPL heavily charged macrocyclic hosts encapsulating small protonated amines with disparate binding affinities. We have shown that the methodological treatment for host neutralization, with explicit ions or with the background neutralizing plasma in the context of alchemical calculations under periodic boundary conditions, has a moderate effect on the calculated affinities. On the other hand, we have shown that seemingly small differences in the force field parameterization in highly symmetric hosts can produce systematic effects on the structural features that can have a significant impact on the predicted binding affinities.

Broadening the Scope of Binding Free Energy Calculations Using a Separated Topologies Approach

Binding free energy calculations predict the potency of compounds to protein binding sites in a physically rigorous manner and see broad application in prioritizing the synthesis of novel drug candidates. Relative binding free energy (RBFE) calculations have emerged as an industry-standard approach to achieve highly accurate rank-order predictions of the potency of related compounds; however, this approach requires that the ligands share a common scaffold and a common binding mode, restricting the methods' domain of applicability. This is a critical limitation since complex modifications to the ligands, especially core hopping, are very common in drug design. Absolute binding free energy (ABFE) calculations are an alternate method that can be used for ligands that are not congeneric. However, ABFE suffers from a known problem of long convergence times due to the need to sample additional degrees of freedom within each system, such as sampling rearrangements necessary to open and close the binding site. Here, we report on an alternative method for RBFE, called Separated Topologies (SepTop), which overcomes the issues in both of the aforementioned methods by enabling large scaffold changes between ligands with a convergence time comparable to traditional RBFE. Instead of only mutating atoms that vary between two ligands, this approach performs two absolute free energy calculations at the same time in opposite directions, one for each ligand. Defining the two ligands independently allows the comparison of the binding of diverse ligands without the artificial constraints of identical poses or a suitable atom-atom mapping. This approach also avoids the need to sample the unbound state of the protein, making it more efficient than absolute binding free energy calculations. Here, we introduce an implementation of SepTop. We developed a general and efficient protocol for running SepTop, and we demonstrated the method on four diverse, pharmaceutically relevant systems. We report the performance of the method, as well as our practical insights into the strengths, weaknesses, and challenges of applying this method in an industrial drug design setting. We find that the accuracy of the approach is sufficiently high to rank order ligands with an accuracy comparable to traditional RBFE calculations while maintaining the additional flexibility of SepTop.

The automated optimisation of a coarse-grained force field using free energy data

Atomistic models provide a detailed representation of molecular systems, but are sometimes inadequate for simulations of large systems over long timescales. Coarse-grained models enable accelerated simulations by reducing the number of degrees of freedom, at the cost of reduced accuracy. New optimisation processes to parameterise these models could improve their quality and range of applicability. We present an automated approach for the optimisation of coarse-grained force fields, by reproducing free energy data derived from atomistic molecular simulations. To illustrate the approach, we implemented hydration free energy gradients as a new target for force field optimisation in ForceBalance and applied it successfully to optimise the un-charged side-chains and the protein backbone in the SIRAH protein coarse-grain force field. The optimised parameters closely reproduced hydration free energies of atomistic models and gave improved agreement with experiment.

How phosphorylation of peptides affects their interaction with 14-3-3η domains

Members of the 14-3-3 domain family have important functions as adapter domains. Via an amphipathic groove on their protein surface they typically bind to disordered C-terminals of other proteins. Importantly, binding partners of 14-3-3 domains usually contain a phosphorylated serine or threonine residue at their binding interface and possess one of three different sequence motifs. Binding of the respective unphosphorylated versions of the peptides is typically strongly disfavored. There is a wealth of structural and thermodynamic data available for the phosphorylated forms but not for the unphosphorylated forms as the binding affinities seem to be too weak to be measurable experimentally. Here, we characterized the mechanistic details that govern the preference for the binding of phosphorylated peptides to 14-3-3η domains by means of molecular dynamics (MD) simulations. We found that the phosphate group is ideally coordinated in the binding pocket whereas the respective unphosphorylated side-chain counterpart is not. Thus, the binding preference results from the tight coordination of the phosphorylated residue at the center of the binding interface. Furthermore, MD simulations of 14-3-3η dimers showed a preference for the simultaneous binding of two phosphorylated peptides in agreement with their experimentally observed cooperativity.

Zinc binding alters the conformational dynamics and drives the transport cycle of the cation diffusion facilitator YiiP

YiiP is a secondary transporter that couples Zn2+ transport to the proton motive force. Structural studies of YiiP from prokaryotes and Znt8 from humans have revealed three different Zn2+ sites and a conserved homodimeric architecture. These structures define the inward-facing and outward-facing states that characterize the archetypal alternating access mechanism of transport. To study the effects of Zn2+ binding on the conformational transition, we use cryo-EM together with molecular dynamics simulation to compare structures of YiiP from Shewanella oneidensis in the presence and absence of Zn2+. To enable single-particle cryo-EM, we used a phage-display library to develop a Fab antibody fragment with high affinity for YiiP, thus producing a YiiP/Fab complex. To perform MD simulations, we developed a nonbonded dummy model for Zn2+ and validated its performance with known Zn2+-binding proteins. Using these tools, we find that, in the presence of Zn2+, YiiP adopts an inward-facing conformation consistent with that previously seen in tubular crystals. After removal of Zn2+ with high-affinity chelators, YiiP exhibits enhanced flexibility and adopts a novel conformation that appears to be intermediate between inward-facing and outward-facing states. This conformation involves closure of a hydrophobic gate that has been postulated to control access to the primary transport site. Comparison of several independent cryo-EM maps suggests that the transition from the inward-facing state is controlled by occupancy of a secondary Zn2+ site at the cytoplasmic membrane interface. This work enhances our understanding of individual Zn2+ binding sites and their role in the conformational dynamics that govern the transport cycle.

Developing and Assessing Nonbonded Dummy Models of Magnesium Ion with Different Hydration Free Energy References

A large diversity in the targeted hydration free energies (HFEs) during model parameterization of metal ions was reported in the literature with a difference by dozens of kcal/mol. Here, we developed a series of nonbonded dummy models of the Mg²⁺ ion targeting different HFE references in TIP3P water, followed by assessments of the designed models in the simulations of MgCl₂ solution and biological systems. Together with the comparison of existing models, we conclude that the difference in the targeted HFEs has a limited influence on the model performance, while the usability of these models differs from case to case. The feasibility of reproducing more properties of Mg²⁺ such as diffusion constants and water exchange rates using a nonbonded dummy model is demonstrated. Underestimated activity derivative and osmotic coefficient of MgCl₂ solutions in high concentration reveal a necessity for further optimization of ion-pair interactions. The developed dummy models are applicable to metal coordination with Asp, Glu, and His residues in metalloenzymes, and the performance in predicting monodentate or bidentate binding modes of Asp/Glu residues depends on the complexity of metal centers and the choice of protein force fields. When both the binding modes coexist, the nonbonded dummy models outperform point charge models, probably in need of considering polarization of metal-binding residues by, for instance, charge calibration in classical force fields. This work is valuable for the use and further development of magnesium ion models for simulations of metal-containing systems with good accuracy.

SAMPL7 Host-Guest Challenge Overview: assessing the reliability of polarizable and non-polarizable methods for binding free energy calculations

The SAMPL challenges focus on testing and driving progress of computational methods to help guide pharmaceutical drug discovery. However, assessment of methods for predicting binding affinities is often hampered by computational challenges such as conformational sampling, protonation state uncertainties, variation in test sets selected, and even lack of high quality experimental data. SAMPL blind challenges have thus frequently included a component focusing on host-guest binding, which removes some of these challenges while still focusing on molecular recognition. Here, we report on the results of the SAMPL7 blind prediction challenge for host-guest affinity prediction. In this study, we focused on three different host-guest categories-a familiar deep cavity cavitand series which has been featured in several prior challenges (where we examine binding of a series of guests to two hosts), a new series of cyclodextrin derivatives which are monofunctionalized around the rim to add amino acid-like functionality (where we examine binding of two guests to a series of hosts), and binding of a series of guests to a new acyclic TrimerTrip host which is related to previous cucurbituril hosts. Many predictions used methods based on molecular simulations, and overall success was mixed, though several methods stood out. As in SAMPL6, we find that one strategy for achieving reasonable accuracy here was to make empirical corrections to binding predictions based on previous data for host categories which have been studied well before, though this can be of limited value when new systems are included. Additionally, we found that alchemical free energy methods using the AMOEBA polarizable force field had considerable success for the two host categories in which they participated. The new TrimerTrip system was also found to introduce some sampling problems, because multiple conformations may be relevant to binding and interconvert only slowly. Overall, results in this challenge tentatively suggest that further investigation of polarizable force fields for these challenges may be warranted.

Charge-Changing Perturbations and Path Sampling via Classical Molecular Dynamic Simulations of Simple Guest-Host Systems

Currently, two different methods dominate the field of biomolecular free-energy calculations for the prediction of binding affinities. Pathway methods are frequently used for large ligands that bind on the surface of a host, such as protein-protein complexes. Alchemical methods, on the other hand, are preferably applied for small ligands that bind to deeply buried binding sites. The latter methods are also widely known to be heavily artifacted by the representation of electrostatic energies in periodic simulation boxes, in particular, when net-charge changes are involved. Different methods have been described to deal with these artifacts, including postsimulation correction schemes and instantaneous correction schemes (e.g., co-alchemical perturbation of ions). Here, we use very simple test systems to show that instantaneous correction schemes with no change in the system net charge lower the artifacts but do not eliminate them. Furthermore, we show that free energies from pathway methods suffer from the same artifacts.

Enhancing the promiscuity of a member of the Caspase protease family by rational design

The N-terminal cleavage of fusion tags to restore the native N-terminus of recombinant proteins is a challenging task and up to today, protocols need to be optimized for different proteins individually. Within this work, we present a novel protease that was designed in-silico to yield enhanced promiscuity toward different N-terminal amino acids. Two mutations in the active-site amino acids of human Caspase-2 were determined to increase the recognition of branched amino-acids, which show only poor binding capabilities in the unmutated protease. These mutations were determined by sequential and structural comparisons of Caspase-2 and Caspase-3 and their effect was additionally predicted using free-energy calculations. The two mutants proposed in the in-silico studies were expressed and in-vitro experiments confirmed the simulation results. Both mutants showed not only enhanced activities toward branched amino acids, but also smaller, unbranched amino acids. We believe that the created mutants constitute an important step toward generalized procedures to restore original N-termini of recombinant fusion proteins.

Dielectric constant and density of aqueous alkali halide solutions by molecular dynamics: A force field assessment

The concentration dependence of the dielectric constant and the density of 11 aqueous alkali halide solutions (LiCl, NaCl, KCl, RbCl, CsCl, LiI, NaI, KI, CsI, KF, and CsF) is investigated by molecular simulation. Predictions using eight non-polarizable ion force fields combined with the TIP4P/ε water model are compared to experimental data. The influence of the water model and the temperature on the results for the NaCl brine are also addressed. The TIP4P/ε water model improves the accuracy of dielectric constant predictions compared to the SPC/E water model. The solution density is predicted well by most ion models. Almost all ion force fields qualitatively capture the decline of the dielectric constant with the increase of concentration for all solutions and with the increase of temperature for NaCl brine. However, the sampled dielectric constant is mostly in poor quantitative agreement with experimental data. These results are related to the microscopic solution structure, ion pairing, and ultimately the force field parameters. Ion force fields with excessive contact ion pairing and precipitation below the experimental solubility limit generally yield higher dielectric constant values. An adequate reproduction of the experimental solubility limit should therefore be a prerequisite for further investigations of the dielectric constant of aqueous electrolyte solutions by molecular simulation.

Correcting electrostatic artifacts due to net-charge changes in the calculation of ligand binding free energies

Alchemically derived free energies are artifacted when the perturbed moiety has a nonzero net charge. The source of the artifacts lies in the effective treatment of the electrostatic interactions within and between the perturbed atoms and remaining (partial) charges in the simulated system. To treat the electrostatic interactions effectively, lattice-summation (LS) methods or cutoff schemes in combination with a reaction-field contribution are usually employed. Both methods render the charging component of the calculated free energies sensitive to essential parameters of the system like the cutoff radius or the box side lengths. Here, we discuss the results of three previously published studies of ligand binding. These studies presented estimates of binding free energies that were artifacted due to the charged nature of the ligands. We show that the size of the artifacts can be efficiently calculated and raw simulation data can be corrected. We compare the corrected results with experimental estimates and nonartifacted estimates from path-sampling methods. Although the employed correction scheme involves computationally demanding continuum-electrostatics calculations, we show that the correction estimate can be deduced from a small sample of configurations rather than from the entire ensemble. This observation makes the calculations of correction terms feasible for complex biological systems. To show the general applicability of the proposed procedure, we also present results where the correction scheme was used to correct independent free energies obtained from simulations employing a cutoff scheme or LS electrostatics. In this work, we give practical guidelines on how to apply the appropriate corrections easily.

Optimization of Protein-Ligand Electrostatic Interactions Using an Alchemical Free-Energy Method

We present an explicit solvent alchemical free-energy method for optimizing the partial charges of a ligand to maximize the binding affinity with a receptor. This methodology can be applied to known ligand–protein complexes to determine an optimized set of ligand partial atomic changes. Three protein–ligand complexes have been optimized in this work: FXa, P38, and the androgen receptor. The sets of optimized charges can be used to identify design principles for chemical changes to the ligands which improve the binding affinity for all three systems. In this work, beneficial chemical mutations are generated from these principles and the resulting molecules tested using free-energy perturbation calculations. We show that three quarters of our chemical changes are predicted to improve the binding affinity, with an average improvement for the beneficial mutations of approximately 1 kcal/mol. In the cases where experimental data are available, the agreement between prediction and experiment is also good. The results demonstrate that charge optimization in explicit solvent is a useful tool for predicting beneficial chemical changes such as pyridinations, fluorinations, and oxygen to sulfur mutations.

Structure, Dynamics, and Hydration Free Energy of Carbon Dioxide in Aqueous Solution: A Quantum Mechanical/Molecular Mechanics Molecular Dynamics Thermodynamic Integration (QM/MM MD TI) Simulation Study

The solvation of carbon dioxide in solution represents a key step for the capture and fixation CO₂ in nature, which may be further influenced by the formation of (bi)carbonate species and/or the formation of CO₂ clusters in solution. The latter processes are strongly dependent on the exact environment of the liquid state (e.g., pH value, solvated ions, etc.) and may interfere with the experimental determination of structural, dynamical, and thermodynamic properties. In this work a hybrid quantum mechanical/molecular mechanical (QM/MM) simulation approach at correlated ab initio level of theory resolution-of-identity second-order Møller-Plesset Perturbation Theory (RI-MP2) has been applied in the framework of thermodynamic integration (TI) to study structure, dynamics, and the hydration free energy of a single carbon dioxide molecule in aqueous solution. A detailed analysis of the individual QM/MM potential energy contributions demonstrate that the overall potential remains highly consistent over the entire sampling phase and that no artificial contributions are influencing the determination of the hydration free energy. The latter value of 0.01 ± 0.92 kcal/mol was found in very good agreement with the values of 0.06 and 0.24 kcal/mol obtained via quasi-chemical theory and experimental measurements, respectively. In order to obtain detailed information about the C- and O _C-water interaction, conically restricted regions with respect to the main axis of the CO₂ molecule have been employed in structural analysis. The presented data not only provide detailed information about the hydration properties of CO₂ but act as a critical validation of the simulation technique, which will be beneficial in the study of nonaqueous solvents such as pure and aqueous NH₃ solutions, which have been suggested as potential candidates to capture CO₂ from anthropogenic sources.

Single-Ion Thermodynamics from First Principles: Calculation of the Absolute Hydration Free Energy and Single-Electrode Potential of Aqueous Li+ Using ab Initio Quantum Mechanical/Molecular Mechanical Molecular Dynamics Simulations

A recently proposed thermodynamic integration (TI) approach formulated in the framework of quantum mechanical/molecular mechanical molecular dynamics (QM/MM MD) simulations is applied to study the structure, dynamics, and absolute intrinsic hydration free energy Δ_s G_M⁺,wat^◦ of the Li⁺ ion at a correlated ab initio level of theory. Based on the results, standard values (298.15 K, ideal gas at 1 bar, ideal solute at 1 molal) for the absolute intrinsic hydration free energy [Formula: see text] of the proton, the surface electric potential jump χ_wat^◦ upon entering bulk water, and the absolute single-electrode potential [Formula: see text] of the reference hydrogen electrode are calculated to be -1099.9 ± 4.2 kJ·mol^-1, 0.13 ± 0.08 V, and 4.28 ± 0.04 V, respectively, in excellent agreement with the standard values recommended by Hünenberger and Reif on the basis of an extensive evaluation of the available experimental data (-1100 ± 5 kJ·mol^-1, 0.13 ± 0.10 V, and 4.28 ± 0.13 V). The simulation results for Li⁺ are also compared to those for Na⁺ and K⁺ from a previous study in terms of relative hydration free energies ΔΔ_s G_M⁺,wat^◦ and relative electrode potentials [Formula: see text]. The calculated values are found to agree extremely well with the experimental differences in standard conventional hydration free energies ΔΔ_s G_M⁺,wat^• and redox potentials [Formula: see text]. The level of agreement between simulation and experiment, which is quantitative within error bars, underlines the substantial accuracy improvement achieved by applying a highly demanding QM/MM approach at the resolution-of-identity second-order Møller-Plesset perturbation (RIMP2) level over calculations relying on purely molecular mechanical or density functional theory (DFT) descriptions. A detailed analysis of the structural and dynamical properties of the Li⁺ hydrate indicates that a correct description of the solvation structure and dynamics is achieved as well at this level of theory. Consideration of the QM/MM potential-energy components also shows that the partitioning into QM and MM zones does not induce any significant energetic artifact for the system considered.

Binding affinities of the farnesoid X receptor in the D3R Grand Challenge 2 estimated by free-energy perturbation and docking

We have studied the binding of 102 ligands to the farnesoid X receptor within the D3R Grand Challenge 2016 blind-prediction competition. First, we employed docking with five different docking software and scoring functions. The selected docked poses gave an average root-mean-squared deviation of 4.2 Å. Consensus scoring gave decent results with a Kendall's τ of 0.26 ± 0.06 and a Spearman's ρ of 0.41 ± 0.08. For a subset of 33 ligands, we calculated relative binding free energies with free-energy perturbation. Five transformations between the ligands involved a change of the net charge and we implemented and benchmarked a semi-analytic correction (Rocklin et al., J Chem Phys 139:184103, 2013) for artifacts caused by the periodic boundary conditions and Ewald summation. The results gave a mean absolute deviation of 7.5 kJ/mol compared to the experimental estimates and a correlation coefficient of R 2 = 0.1. These results were among the four best in this competition out of 22 submissions. The charge corrections were significant (7-8 kJ/mol) and always improved the results. By employing 23 intermediate states in the free-energy perturbation, there was a proper overlap between all states and the precision was 0.1-0.7 kJ/mol. However, thermodynamic cycles indicate that the sampling was insufficient in some of the perturbations.

Extending the Nonbonded Cationic Dummy Model to Account for Ion-Induced Dipole Interactions

Modeling metalloproteins often requires classical molecular dynamics (MD) simulations in order to capture their relevant motions, which in turn necessitates reliable descriptions of the metal centers involved. One of the most successful approaches to date is provided by the "cationic dummy model", where the positive charge of the metal ion is transferred toward dummy particles that are bonded to the central metal ion in a predefined coordination geometry. While this approach allows for ligand exchange, and captures the correct electrostatics as demonstrated for different divalent metal ions, current dummy models neglect ion-induced dipole interactions. In the present work, we resolve this weakness by taking advantage of the recently introduced 12-6-4 type Lennard-Jones potential to include ion-induced dipole interactions. We revise our previous dummy model for Mg²⁺ and demonstrate that the resulting model can simultaneously reproduce the experimental solvation free energy and metal-ligand distances without the need for artificial restraints or bonds. As ion-induced dipole interactions become particularly important for highly charged metal ions, we develop dummy models for the biologically relevant ions Al³⁺, Fe³⁺, and Cr³⁺. Finally, the effectiveness of our new models is demonstrated in MD simulations of several diverse (and highly challenging to simulate) metalloproteins.

Energetics and mechanism of anion permeation across formate-nitrite transporters

Formate-nitrite transporters (FNTs) facilitate the translocation of monovalent polyatomic anions, such as formate and nitrite, across biological membranes. FNTs are widely distributed among pathogenic bacteria and eukaryotic parasites, but they lack human homologues, making them attractive drug targets. The mechanisms and energetics involved in anion permeation across the FNTs have remained largely unclear. Both, channel and transporter mode of function have been proposed, with strong indication of proton coupling to the permeation process. We combine molecular dynamics simulations, quantum mechanical calculations, and pK _a calculations, to compute the energetics of the complete permeation cycle of an FNT. We find that anions as such, are not able to traverse the FNT pore. Instead, anion binding into the pore is energetically coupled to protonation of a centrally located histidine. In turn, the histidine can protonate the permeating anion, thereby enabling its release. Such mechanism can accommodate the functional diversity among the FNTs, as it may facilitate both, export and import of substrates, with or without proton co-transport. The mechanism excludes proton leakage via the Grotthuss mechanism, and it rationalises the selectivity for weak acids.

Polarizable and Non-Polarizable Force Field Representations of Ferric Cation and Validations

The AMOEBA polarizable force field of ferric ion was optimized and applied to study the hydration of ferric ion and its complexation with porphine in the aqueous phase. The nonpolarizable force field was also optimized for comparison. The AMOEBA force field was found to give a more accurate hydration free energy than the nonpolarizable force field with respect to experimental data, and correctly predict the most stable electronic state of hydrated Fe³⁺, which is the sextet state, and of the Fe^(III)-Por complex, which is the quartet state, consistent with the literature that was carried out using the DFT method. The explicit inclusion of charge transfer between Fe³⁺ and ligand was found to be important in order to obtain a precise picture of polarization energy and van der Waals energy, which otherwise deviate from the corresponding energy components derived from ab initio calculations. The successful application of the AMOEBA force field in the characterization of aquo Fe^(III)-Por complexes suggests that its use may be extended to the study of the dynamics of metalloenzyme containing highly charged metal ions in the condensed phase with reliable treatment of the interactions between metal atom and protein.

Is the Solution Activity Derivative Sufficient to Parametrize Ion-Ion Interactions? Ions for TIP5P Water

Biomolecular processes involve hydrated ions, and thus molecular simulations of such processes require accurate force-field parameters for these ions. In the best force-fields, both ion-water and anion-cation interactions are explicitly parametrized. First, the ion Lennard-Jones parameters are optimized to reproduce, for example, single ion solvation free energies; then ion-pair interactions are adjusted to match experimental activity or activity derivatives. Here, we apply this approach to derive optimized parameters for concentrated NaCl, KCl, MgCl₂, and CaCl₂ salt solutions, to be used with the TIP5P water model. These parameters are of interest because of a number of desirable properties of the TIP5P water model, especially for the simulation of carbohydrates. The results show, that this state of the art approach is insufficient, because the activity derivative often reaches a plateau near the target experimental value, for a wide range of parameter values. The plateau emerges from the interconversion between different types of ion pairs, so parameters leading to equally good agreement with the target solution activity or activity derivative yield very different solution structures. To resolve this indetermination, a second target property, such as the experimentally determined ion-ion coordination number, is required to uniquely determine anion-cation interactions. Simulations show that combining activity derivatives and coordination number as experimental target properties to parametrize ion-ion interactions, is a powerful method for reliable ion-water force field parametrization, and gives insight into the concentration of contact or solvent shared ion pairs in a wide range of salt concentrations. For the alkali and halide ions Li⁺, Rb⁺, Cs⁺, F^-, Br^-, and I^-, we present ion-water parameters appropriate at infinite dilution only.

Finite-size effect on the charging free energy of protein in explicit solvent

The finite-size effect in periodic system is examined for the charging free energy of protein in explicit solvent over a variety of charged states. The key to the finite-size correction is the self-energy, which is defined as the interaction energy of the solute with its own periodic images and the neutralizing background. By employing the thermodynamic-integration method with systematically varied sizes of the unit cell of molecular dynamics (MD) simulations, we show for ubiquitin that the self-energy corrects the finite-size effect on the charging free energy within 1 kcal/mol at total charges of -5e, -1e, neutral, and +1e and within 5 kcal/mol even for a highly charged state with +8e. We then sought the additional correction from the solvation effect using the numerical solution to the Poisson equation of the protein with implicit solvent. This correction reduces the cell-size dependence of the charging free energy at +8e to 3 kcal/mol and is well expressed as the self-energy divided by the dielectric constant of solvent water.

Structural, dynamical, and transport properties of the hydrated halides: How do At(-) bulk properties compare with those of the other halides, from F(-) to I(-)?

The properties of halides from the lightest, fluoride (F(-)), to the heaviest, astatide (At(-)), have been studied in water using a polarizable force-field approach based on molecular dynamics (MD) simulations at the 10 ns scale. The selected force-field explicitly treats the cooperativity within the halide-water hydrogen bond networks. The force-field parameters have been adjusted to ab initio data on anion/water clusters computed at the relativistic Möller-Plesset second-order perturbation theory level of theory. The anion static polarizabilities of the two heaviest halides, I(-) and At(-), were computed in the gas phase using large and diffuse atomic basis sets, and taking into account both electron correlation and spin-orbit coupling within a four-component framework. Our MD simulation results show the solvation properties of I(-) and At(-) in aqueous phase to be very close. For instance, their first hydration shells are structured and encompass 9.2 and 9.1 water molecules at about 3.70 ± 0.05 Å, respectively. These values have to be compared to the F(-), Cl(-), and Br(-) ones, i.e., 6.3, 8.4, and 9.0 water molecules at 2.74, 3.38, and 3.55 Å, respectively. Moreover our computations predict the solvation free energy of At(-) in liquid water at ambient conditions to be 68 kcal mol(-1), a value also close the I(-) one, about 70 kcal mol(-1). In all, our simulation results for I(-) are in excellent agreement with the latest neutron- and X-ray diffraction studies. Those for the At(-) ion are predictive, as no theoretical or experimental data are available to date.

Origin of Asymmetric Solvation Effects for Ions in Water and Organic Solvents Investigated Using Molecular Dynamics Simulations: The Swain Acity-Basity Scale Revisited

The asymmetric solvation of ions can be defined as the tendency of a solvent to preferentially solvate anions over cations or cations over anions, at identical ionic charge magnitudes and effective sizes. Taking water as a reference, these effects are quantified experimentally for many solvents by the relative acity (A) and basity (B) parameters of the Swain scale. The goal of the present study is to investigate the asymmetric solvation of ions using molecular dynamics simulations, and to connect the results to this empirical scale. To this purpose, the charging free energies of alkali and halide ions, and of their hypothetical oppositely charged counterparts, are calculated in a variety of solvents. In a first set of calculations, artificial solvent models are considered that present either a charge or a shape asymmetry at the molecular level. The solvation asymmetry, probed by the difference in charging free energy between the two oppositely charged ions, is found to encompass a term quadratic in the ion charge, related to the different solvation structures around the anion and cation, and a term linear in the ion charge, related to the solvation structure around the uncharged ion-sized cavity. For these simple solvent models, the two terms are systematically counteracting each other, and it is argued that only the quadratic term should be retained when comparing the results of simulations involving physical solvents to experimental data. In a second set of calculations, 16 physical solvents are considered. The theoretical estimates for the acity A are found to correlate very well with the Swain parameters, whereas the correlation for B is very poor. Based on this observation, the Swain scale is reformulated into a new scale involving an asymmetry parameter Σ, positive for acitic solvents and negative for basitic ones, and a polarity parameter Π. This revised scale has the same predictive power as the original scale, but it characterizes asymmetry in an absolute sense, the atomistic simulations playing the role of an extra-thermodynamic assumption, and is optimally compatible with the simulation results. Considering the 55 solvents in the Swain set, it is observed that a moderate basity (Σ between -0.9 and -0.3, related to electronic polarization) represents the baseline for most solvents, while a highly variable acity (Σ between 0.0 and 3.0, related to hydrogen-bond donor capacity modulated by inductive effects) represents a landmark of protic solvents.

A classical force field for tetrahedral oxyanions developed using hydration properties: The examples of pertechnetate (TcO4(-)) and sulfate (SO4(2-))

Radioactive pertechnetate, (99)TcO4(-), is one of the most problematic ionic species in the context of the clean up and storage of nuclear waste. Molecular simulations can be used to understand the behavior of TcO4(-) in dilute aqueous solutions, providing reliable potentials are available. This work outlines the development of a new potential model for TcO4(-) and competing SO4(2-), optimized using their hydration properties, such as the Gibbs hydration free energy (calculated using Bennett's acceptance ratio method). The findings show that the TcO4(-) oxyanion has a very low hydration free energy (-202 kJ mol(-1)) compared to other anions (Cl(-), I(-), SO4(2-)) leading to fast water exchange dynamics and explaining its observed high mobility in the aqueous environment. Its hydrated structure, investigated using ion-water radial distribution functions, shows that it is unique amongst the other anions in that it does not possess well-defined hydration shells. Since contaminants and ubiquitous species in the aqueous environment are often present as tetrahedral oxyanions, it is proposed that the approach could easily be extended to a whole host of other species.

Nuclear quantum effects in water exchange around lithium and fluoride ions

We employ classical and ring polymer molecular dynamics simulations to study the effect of nuclear quantum fluctuations on the structure and the water exchange dynamics of aqueous solutions of lithium and fluoride ions. While we obtain reasonably good agreement with experimental data for solutions of lithium by augmenting the Coulombic interactions between the ion and the water molecules with a standard Lennard-Jones ion-oxygen potential, the same is not true for solutions of fluoride, for which we find that a potential with a softer repulsive wall gives much better agreement. A small degree of destabilization of the first hydration shell is found in quantum simulations of both ions when compared with classical simulations, with the shell becoming less sharply defined and the mean residence time of the water molecules in the shell decreasing. In line with these modest differences, we find that the mechanisms of the exchange processes are unaffected by quantization, so a classical description of these reactions gives qualitatively correct and quantitatively reasonable results. We also find that the quantum effects in solutions of lithium are larger than in solutions of fluoride. This is partly due to the stronger interaction of lithium with water molecules, partly due to the lighter mass of lithium and partly due to competing quantum effects in the hydration of fluoride, which are absent in the hydration of lithium.

Estimating Translational and Orientational Entropies Using the k-Nearest Neighbors Algorithm

Inhomogeneous fluid solvation theory (IFST) and free energy perturbation (FEP) calculations were performed for a set of 20 solutes to compute the hydration free energies. We identify the weakness of histogram methods in computing the IFST hydration entropy by showing that previously employed histogram methods overestimate the translational and orientational entropies and thus underestimate their contribution to the free energy by a significant amount. Conversely, we demonstrate the accuracy of the k-nearest neighbors (KNN) algorithm in computing these translational and orientational entropies. Implementing the KNN algorithm within the IFST framework produces a powerful method that can be used to calculate free-energy changes for large perturbations. We introduce a new KNN approach to compute the total solute-water entropy with six degrees of freedom, as well as the translational and orientational contributions. However, results suggest that both the solute-water and water-water entropy terms are significant and must be included. When they are combined, the IFST and FEP hydration free energies are highly correlated, with an R(2) of 0.999 and a mean unsigned difference of 0.9 kcal/mol. IFST predictions are also highly correlated with experimental hydration free energies, with an R(2) of 0.997 and a mean unsigned error of 1.2 kcal/mol. In summary, the KNN algorithm is shown to yield accurate estimates of the combined translational-orientational entropy and the novel approach of combining distance metrics that is developed here could be extended to provide a powerful method for entropy estimation in numerous contexts.

Toward the correction of effective electrostatic forces in explicit-solvent molecular dynamics simulations: restraints on solvent-generated electrostatic potential and solvent polarization

Despite considerable advances in computing power, atomistic simulations under nonperiodic boundary conditions, with Coulombic electrostatic interactions and in systems large enough to reduce finite-size associated errors in thermodynamic quantities to within the thermal energy, are still not affordable. As a result, periodic boundary conditions, systems of microscopic size and effective electrostatic interaction functions are frequently resorted to. Ensuing artifacts in thermodynamic quantities are nowadays routinely corrected a posteriori, but the underlying configurational sampling still descends from spurious forces. The present study addresses this problem through the introduction of on-the-fly corrections to the physical forces during an atomistic molecular dynamics simulation. Two different approaches are suggested, where the force corrections are derived from special potential energy terms. In the first approach, the solvent-generated electrostatic potential sampled at a given atom site is restrained to a target value involving corrections for electrostatic artifacts. In the second approach, the long-range regime of the solvent polarization around a given atom site is restrained to the Born polarization, i.e., the solvent polarization corresponding to the ideal situation of a macroscopic system under nonperiodic boundary conditions and governed by Coulombic electrostatic interactions. The restraints are applied to the explicit-water simulation of a hydrated sodium ion, and the effect of the restraints on the structural and energetic properties of the solvent is illustrated. Furthermore, by means of the calculation of the charging free energy of a hydrated sodium ion, it is shown how the electrostatic potential restraint translates into the on-the-fly consideration of the corresponding free-energy correction terms. It is discussed how the restraints can be generalized to situations involving several solute particles. Although the present study considers a very simple system only, it is an important step toward the on-the-fly elimination of finite-size and approximate-electrostatic artifacts during atomistic molecular dynamics simulations.

Cluster-continuum quasichemical theory calculation of the lithium ion solvation in water, acetonitrile and dimethyl sulfoxide: an absolute single-ion solvation free energy scale

A solvation free energy scale excluding the net electrostatic potential inside the solute cavity is presented.

Hydrothermal properties of the COS/D2 water model: a polarizable charge-on-spring water model, at elevated temperatures and pressures

Molecular simulations have been conducted to assess the pVT properties and static permittivity of the charge-on-spring polarizable water model COS/D2 at hydrothermal conditions from 300 to 450 °C and bulk densities of 0.001 to 1.0 g cm⁻³.

Absolute binding-free energies between standard RNA/DNA nucleobases and amino-acid sidechain analogs in different environments

Despite the great importance of nucleic acid-protein interactions in the cell, our understanding of their physico-chemical basis remains incomplete. In order to address this challenge, we have for the first time determined potentials of mean force and the associated absolute binding free energies between all standard RNA/DNA nucleobases and amino-acid sidechain analogs in high- and low-dielectric environments using molecular dynamics simulations and umbrella sampling. A comparison against a limited set of available experimental values for analogous systems attests to the quality of the computational approach and the force field used. Overall, our analysis provides a microscopic picture behind nucleobase/sidechain interaction preferences and creates a unified framework for understanding and sculpting nucleic acid-protein interactions in different contexts. Here, we use this framework to demonstrate a strong relationship between nucleobase density profiles of mRNAs and nucleobase affinity profiles of their cognate proteins and critically analyze a recent hypothesis that the two may be capable of direct, complementary interactions.

Practical Aspects of Free-Energy Calculations: A Review

Free-energy calculations in the framework of classical molecular dynamics simulations are nowadays used in a wide range of research areas including solvation thermodynamics, molecular recognition, and protein folding. The basic components of a free-energy calculation, that is, a suitable model Hamiltonian, a sampling protocol, and an estimator for the free energy, are independent of the specific application. However, the attention that one has to pay to these components depends considerably on the specific application. Here, we review six different areas of application and discuss the relative importance of the three main components to provide the reader with an organigram and to make nonexperts aware of the many pitfalls present in free energy calculations.

Net charge changes in the calculation of relative ligand-binding free energies via classical atomistic molecular dynamics simulation

The calculation of binding free energies of charged species to a target molecule is a frequently encountered problem in molecular dynamics studies of (bio-)chemical thermodynamics. Many important endogenous receptor-binding molecules, enzyme substrates, or drug molecules have a nonzero net charge. Absolute binding free energies, as well as binding free energies relative to another molecule with a different net charge will be affected by artifacts due to the used effective electrostatic interaction function and associated parameters (e.g., size of the computational box). In the present study, charging contributions to binding free energies of small oligoatomic ions to a series of model host cavities functionalized with different chemical groups are calculated with classical atomistic molecular dynamics simulation. Electrostatic interactions are treated using a lattice-summation scheme or a cutoff-truncation scheme with Barker-Watts reaction-field correction, and the simulations are conducted in boxes of different edge lengths. It is illustrated that the charging free energies of the guest molecules in water and in the host strongly depend on the applied methodology and that neglect of correction terms for the artifacts introduced by the finite size of the simulated system and the use of an effective electrostatic interaction function considerably impairs the thermodynamic interpretation of guest-host interactions. Application of correction terms for the various artifacts yields consistent results for the charging contribution to binding free energies and is thus a prerequisite for the valid interpretation or prediction of experimental data via molecular dynamics simulation. Analysis and correction of electrostatic artifacts according to the scheme proposed in the present study should therefore be considered an integral part of careful free-energy calculation studies if changes in the net charge are involved.

Calculating the binding free energies of charged species based on explicit-solvent simulations employing lattice-sum methods: an accurate correction scheme for electrostatic finite-size effects

The calculation of a protein-ligand binding free energy based on molecular dynamics (MD) simulations generally relies on a thermodynamic cycle in which the ligand is alchemically inserted into the system, both in the solvated protein and free in solution. The corresponding ligand-insertion free energies are typically calculated in nanoscale computational boxes simulated under periodic boundary conditions and considering electrostatic interactions defined by a periodic lattice-sum. This is distinct from the ideal bulk situation of a system of macroscopic size simulated under non-periodic boundary conditions with Coulombic electrostatic interactions. This discrepancy results in finite-size effects, which affect primarily the charging component of the insertion free energy, are dependent on the box size, and can be large when the ligand bears a net charge, especially if the protein is charged as well. This article investigates finite-size effects on calculated charging free energies using as a test case the binding of the ligand 2-amino-5-methylthiazole (net charge +1 e) to a mutant form of yeast cytochrome c peroxidase in water. Considering different charge isoforms of the protein (net charges -5, 0, +3, or +9 e), either in the absence or the presence of neutralizing counter-ions, and sizes of the cubic computational box (edges ranging from 7.42 to 11.02 nm), the potentially large magnitude of finite-size effects on the raw charging free energies (up to 17.1 kJ mol(-1)) is demonstrated. Two correction schemes are then proposed to eliminate these effects, a numerical and an analytical one. Both schemes are based on a continuum-electrostatics analysis and require performing Poisson-Boltzmann (PB) calculations on the protein-ligand system. While the numerical scheme requires PB calculations under both non-periodic and periodic boundary conditions, the latter at the box size considered in the MD simulations, the analytical scheme only requires three non-periodic PB calculations for a given system, its dependence on the box size being analytical. The latter scheme also provides insight into the physical origin of the finite-size effects. These two schemes also encompass a correction for discrete solvent effects that persists even in the limit of infinite box sizes. Application of either scheme essentially eliminates the size dependence of the corrected charging free energies (maximal deviation of 1.5 kJ mol(-1)). Because it is simple to apply, the analytical correction scheme offers a general solution to the problem of finite-size effects in free-energy calculations involving charged solutes, as encountered in calculations concerning, e.g., protein-ligand binding, biomolecular association, residue mutation, pKa and redox potential estimation, substrate transformation, solvation, and solvent-solvent partitioning.