Simulating the Chelate Effect

Despite the rich history of experimental studies focusing on the thermochemistry and kinetics associated with the chelate effect, molecular-level computational studies on the chelate ring opening/ring closure are scarce. The challenge lies in an accurate description of both the metal ion and its aqueous environment. Herein, we demonstrate that an optimized 12-6-4 Lennard-Jones (LJ) model can capture the thermodynamics and provide detailed structural and mechanistic insights into the formation of ethylenediamine (en) complexes with metal ions. The water molecules in the first solvation shell of the metal ion are found to facilitate the chelate ring formation. The optimized parameters further simulate the formation of bis and tris(en) complexes representing the wide applicability of the model to simulate coordination chemistry and self-assembly processes.

GPU-Accelerated Molecular Dynamics and Free Energy Methods in Amber18: Performance Enhancements and New Features

We report progress in graphics processing unit (GPU)-accelerated molecular dynamics and free energy methods in Amber18. Of particular interest is the development of alchemical free energy algorithms, including free energy perturbation and thermodynamic integration methods with support for nonlinear soft-core potential and parameter interpolation transformation pathways. These methods can be used in conjunction with enhanced sampling techniques such as replica exchange, constant-pH molecular dynamics, and new 12-6-4 potentials for metal ions. Additional performance enhancements have been made that enable appreciable speed-up on GPUs relative to the previous software release.

Anharmonic vibrations of N-H in Cl(-)(N-methylacetamide)1(H2O)(0-2)Ar2 cluster ions. Combined IRPD experiments and BOMD simulations

Infrared Predissociation (IRPD) spectra of Cl(-)(NMA)1(H2O)0-2Ar2 combined with Born-Oppenheimer Molecular Dynamics (BOMD) IR spectra have been acquired, providing the structure and dynamics of these systems. We show that the chloride ion is bound to the hydrogen of the amide N-H group, forming a strong ionic hydrogen bond, weakening the N-H stretch, and shifting it to lower frequency. The presence of water molecules enhances the ionic hydrogen bond by binding to the amide carbonyl oxygen of NMA and shifts the N-H stretch further to lower frequency. The BOMD IR spectra can recapture all, but about 100 cm(-1), of the 600 to 700 cm(-1) shifts due to the strong N-H stretch anharmonicities observed in experiments. This residual error was found to be due to the lack of zero point energy in the classical treatment of motion in the BOMD method.

Exploring the ion selectivity properties of a large number of simplified binding site models

The ability to discriminate between different cations efficiently is essential for the proper physiological functioning of many membrane transport proteins. One obvious mechanism of ion selectivity is when a binding site is structurally constrained by the protein architecture and its geometry is precisely adapted to fit an ion of a given size. This mechanism is not effective in the case of flexible protein binding sites that are able to deform structurally or to adapt to a bound ion. In this study, the concept of nontrivial ion selectivity arising in a highly flexible protein binding site conceptually represented as a microdroplet of ligands confined to a small volume is explored. The environment imposed by the spatial confinement is a critical feature of the reduced models. A large number of reduced binding site models (1077) comprising typical ion-coordinating ligands (carbonyl, hydroxyl, carboxylate, water) are constructed and characterized for Na(+)/K(+) and Ca(2+)/Ba(2+) size selectivity using free energy perturbation molecular dynamics simulations. Free energies are highly correlated with the sum of ion-ligand and ligand-ligand mean interactions, but the relative balance of those two contributions is different for K(+)-selective and Na(+)-selective binding sites. The analysis indicates that both the number and the type of ligands are important factors in ion selectivity.

Ion selectivity of the KcsA channel: a perspective from multi-ion free energy landscapes

Potassium (K(+)) channels are specialized membrane proteins that are able to facilitate and regulate the conduction of K(+) through cell membranes. Comprising five specific cation binding sites (S(0)-S(4)) formed by the backbone carbonyl groups of conserved residues common to all K(+) channels, the narrow selectivity filter allows fast conduction of K(+) while being highly selective for K(+) over Na(+). To extend our knowledge of the microscopic mechanism underlying selectivity in K(+) channels, we characterize the free energy landscapes governing the entry and translocation of a Na(+) or a K(+) from the extracellular side into the selectivity filter of KcsA. The entry process of an extracellular ion is examined in the presence of two additional K(+) in the pore, and the three-ion potential of mean force is computed using extensive all-atom umbrella sampling molecular dynamics simulations. A comparison of the potentials of mean force yields a number of important results. First, the free energy minima corresponding to configurations with extracellular K(+) or Na(+) in binding site S(0) or S(1) are similar in depth, suggesting that the thermodynamic selectivity governed by the free energy minima for those two binding sites is insignificant. Second, the free energy barriers between stable multi-ion configurations are generally higher for Na(+) than for K(+), implying that the kinetics of ion conduction is slower when a Na(+) enters the pore. Third, the region corresponding to binding site S(2) near the center of the narrow pore emerges as the most selective for K(+) over Na(+). In particular, while there is a stable minimum for K(+) in site S(2), Na(+) faces a steep free energy increase with no local free energy well in this region. Lastly, analysis shows that selectivity is not correlated with the overall coordination number of the ion entering the pore, but is predominantly affected by changes in the type of coordinating ligands (carbonyls versus water molecules). These results further highlight the importance of the central region near binding site S(2) in the selectivity filter of K(+) channels.

Molecular modeling and dynamics studies with explicit inclusion of electronic polarizability. Theory and applications

A current emphasis in empirical force fields is on the development of potential functions that explicitly treat electronic polarizability. In the present article, the commonly used methodologies for modelling electronic polarization are presented along with an overview of selected application studies. Models presented include induced point-dipoles, classical Drude oscillators, and fluctuating charge methods. The theoretical background of each method is followed by an introduction to extended Langrangian integrators required for computationally tractable molecular dynamics simulations using polarizable force fields. The remainder of the review focuses on application studies using these methods. Emphasis is placed on water models, for which numerous examples exist, with a more thorough discussion presented on the recently published models associated with the Drude-based CHARMM and the AMOEBA force fields. The utility of polarizable models for the study of ion solvation is then presented followed by an overview of studies of small molecules (e.g. CCl(4), alkanes, etc) and macromolecule (proteins, nucleic acids and lipid bilayers) application studies. The review is written with the goal of providing a general overview of the current status of the field and to facilitate future application and developments.

Interfacing Q-Chem and CHARMM to perform QM/MM reaction path calculations

A hybrid quantum mechanical/molecular mechanical (QM/MM) potential energy function with Hartree-Fock, density functional theory (DFT), and post-HF (RIMP2, MP2, CCSD) capability has been implemented in the CHARMM and Q-Chem software packages. In addition, we have modified CHARMM and Q-Chem to take advantage of the newly introduced replica path and the nudged elastic band methods, which are powerful techniques for studying reaction pathways in a highly parallel (i.e., parallel/parallel) fashion, with each pathway point being distributed to a different node of a large cluster. To test our implementation, a series of systems were studied and comparisons were made to both full QM calculations and previous QM/MM studies and experiments. For instance, the differences between HF, DFT, MP2, and CCSD QM/MM calculations of H2O...H2O, H2O...Na+, and H2O...Cl- complexes have been explored. Furthermore, the recently implemented polarizable Drude water model was used to make comparisons to the popular TIP3P and TIP4P water models for doing QM/MM calculations. We have also computed the energetic profile of the chorismate mutase catalyzed Claisen rearrangement at various QM/MM levels of theory and have compared the results with previous studies. Our best estimate for the activation energy is 8.20 kcal/mol and for the reaction energy is -23.1 kcal/mol, both calculated at the MP2/6-31+G(d)//MP2/6-31+G(d)/C22 level of theory.

Absolute hydration free energy scale for alkali and halide ions established from simulations with a polarizable force field

A polarizable potential function for the hydration of alkali and halide ions is developed on the basis of the recent SWM4-DP water model [Lamoureux, G.; MacKerell, A. D., Jr.; Roux, B. J. Chem. Phys. 2003, 119, 5185]. Induced polarization is incorporated using classical Drude oscillators that are treated as auxiliary dynamical degrees of freedom. The ions are represented as polarizable Lennard-Jones centers, whose parameters are optimized to reproduce the binding energies of gas-phase monohydrates and the hydration free energies in the bulk liquid. Systematic exploration of the parameters shows that the monohydrate binding energies can be consistent with a unique hydration free energy scale if the computed hydration free energies incorporate the contribution from the air/water interfacial electrostatic potential (-540 mV for SWM4-DP). The final model, which can satisfyingly reproduce both gas and bulk-phase properties, corresponds to an absolute scale in which the intrinsic hydration free energy of the proton is -247 kcal/mol.

A theoretical study of aqueous solvation of K comparing ab initio, polarizable, and fixed-charge models

The hydration of K(+) is studied using a hierarchy of theoretical approaches, including ab initio Born-Oppenheimer molecular dynamics and Car-Parrinello molecular dynamics, a polarizable force field model based on classical Drude oscillators, and a nonpolarizable fixed-charge potential based on the TIP3P water model. While models based more directly on quantum mechanics offer the possibility to account for complex electronic effects, polarizable and fixed-charges force fields allow for simulations of large systems and the calculation of thermodynamic observables with relatively modest computational costs. A particular emphasis is placed on investigating the sensitivity of the polarizable model to reproduce key aspects of aqueous K(+), such as the coordination structure, the bulk hydration free energy, and the self diffusion of K(+). It is generally found that, while the simple functional form of the polarizable Drude model imposes some restrictions on the range of properties that can simultaneously be fitted, the resulting hydration structure for aqueous K(+) agrees well with experiment and with more sophisticated computational models. A counterintuitive result, seen in Car-Parrinello molecular dynamics and in simulations with the Drude polarizable force field, is that the average induced molecular dipole of the water molecules within the first hydration shell around K(+) is slightly smaller than the corresponding value in the bulk. In final analysis, the perspective of K(+) hydration emerging from the various computational models is broadly consistent with experimental data, though at a finer level there remain a number of issues that should be resolved to further our ability in modeling ion hydration accurately.

Polarization effects and charge transfer in the KcsA potassium channel

The electronic structure of the selectivity filter of KcsA K(+) channel is investigated by density functional theory (DFT/BLYP) and QM/MM methods. The quantum part includes the selectivity filter, which is polarized by the electrostatic field of the environment treated with the Amber force field. The details of the electronic structure were investigated using the maximally localized Wannier function centers of charge and Bader's atoms in molecules charge analysis. Our results show that the channel backbone carbonyl groups are largely polarized and that there is a sizeable charge transfer from the backbone to the cations. These effects are expected to be important for an accurate description of the carbonyl groups and the ion-ion electrostatic repulsion, which have been proposed to play a central role for the selectivity mechanism of the channel [S.Y. Noskov, S. Berneche, B. Roux, Control of ion selectivity in potassium channels by electrostatic and dynamic properties of carbonyl ligands. Nature 431 (2004) 830-834].

Changes in phosphatidylcholine headgroup tilt and water order induced by monovalent salts: molecular dynamics simulations

The association between monovalent salts and neutral lipid bilayers is known to influence global bilayer structural properties such as headgroup conformational fluctuations and the dipole potential. The local influence of the ions, however, has been unknown due to limited structural resolution of experimental methods. Molecular dynamics simulations are used here to elucidate local structural rearrangements upon association of a series of monovalent Na(+) salts to a palmitoyl-oleoyl-phosphatidylcholine bilayer. We observe association of all ion types in the interfacial region. Larger anions, which are meant to rationalize data regarding a Hofmeister series of anions, bind more deeply within the bilayer than either Cl(-) or Na(+). Although the simulations are able to reproduce experimentally measured quantities, the analysis is focused on local properties currently invisible to experiments, which may be critical to biological systems. As such, for all ion types, including Cl(-), we show local ion-induced perturbations to headgroup tilt, the extent and direction of which is sensitive to ion charge and size. Additionally, we report salt-induced ordering of the water well beyond the interfacial region, which may be significant in terms of hydration repulsion between stacked bilayers.

Dynamical behavior of one-dimensional water molecule chains in zeolites: Nanosecond time-scale molecular dynamics simulations of bikitaite

Nanosecond scale molecular dynamics simulations of the behavior of the one-dimensional water molecule chains adsorbed in the parallel nanochannels of bikitaite, a rare lithium containing zeolite, were performed at different temperatures and for the fully and partially hydrated material. New empirical potential functions have been developed for representing lithium-water interactions. The structure and the vibrational spectrum of bikitaite were in agreement both with experimental data and Car-Parrinello molecular dynamics results. Classical molecular dynamics simulations were extended to the nanosecond time scale in order to study the flip motion of water molecules around the hydrogen bonds connecting adjacent molecules in the chains, which has been observed by NMR experiments, and the dehydration mechanism at high temperature. Computed relaxation times of the flip motion follow the Arrhenius behavior found experimentally, but the activation energy of the simulated system is slightly underestimated. Based on the results of the simulations, it may be suggested that the dehydration proceeds by a defect-driven stepwise diffusion. The diffusive mechanism appears as a single-file motion: the molecules never pass one another, even at temperatures as high as about 1000 K, nor can they switch between different channels. However, the mean square displacement (MSD) of the molecules, computed with respect to the center of mass of the simulated system, shows an irregular trend from which the single-file diffusion cannot be clearly evidenced. If the MSDs are evaluated with respect to the center of mass of the molecules hosted in each channel, the expected dependence on the square root of time finally appears.

A QM-MM interface between CHARMM and TURBOMOLE: implementation and application to systems in bulk phase and biologically active systems

The implementation of a hybrid QM-MM approach combining ab initio and density functional methods of TURBOMOLE with the molecular mechanics program package CHARMM is described. An interface has been created to allow data exchange between the two applications. With this method the efficient multiprocessor capabilities of TURBOMOLE can be utilized with CHARMM running as a single processor application. Therefore, features of nonparallel running code in CHARMM like the TRAVEL module for locating saddle points or VIBRAN for the calculation of second derivatives can be exploited by running the CPU intensive QM calculations in parallel. To test the methodology, several small systems are studied with both Hartree-Fock and density functional methods and varying QM-MM boundaries. Also, the computationally efficient RI-J method has been examined for use in QM-MM applications. A B(12) cofactor containing cobalt has been studied, to examine systems with a large QM region and transition metals. All tested methods perform satisfactory in comparison with pure quantum calculations. Additionally, algorithms for the characterization of saddle points have been tested for their potential use in QM-MM problems. The TRAVEL module of CHARMM has been applied to the Menshutkin reaction in the condensed phase, and a saddle point was located. This saddle point was verified by calculation of a steepest descent path connecting educt, transition state, and product, and by calculation of vibrational modes.

Understanding the Hofmeister effect in interactions between chaotropic anions and lipid bilayers: molecular dynamics simulations

A set of all-atom molecular dynamics simulations have been performed to better understand critical phenomena regarding a Hofmeister series of anions and lipid bilayers. The simulations isolate the effect of anion size and show clear differences in the interactions with the dipolar phoshpatidylcholine headgroup. Cl- anions penetrate into the headgroup region of the bilayer, but the simulations confirm theories which predict that larger anions penetrate more deeply, into a more heterogeneous and hydrophobic molecular region. That anion size leads to such differences in partitioning in the bilayer provides atomic-level support to hypotheses inspired by several experimental studies. The ability of larger anions to bury deep within the bilayer is correlated with a less well-structured hydration shell, shedding of which upon penetration incurs a smaller penalty for the larger anions than for Cl-.

Modeling permeation energetics in the KcsA potassium channel

The thermodynamics of cation permeation through the KcsA K(+) channel selectivity filter is studied from the perspective of a physically transparent semimicroscopic model using Monte Carlo free energy integration. The computational approach chosen permits dissection of the separate contributions to ionic stabilization arising from different parts of the channel (selectivity filter carbonyls, single-file water, cavity water, reaction field of bulk water, inner helices, ionizable residues). All features play important roles; their relative significance varies with the ion's position in the filter. The cavity appears to act as an electrostatic buffer, shielding filter ions from structural changes in the inner pore. The model exhibits K(+) vs. Na(+) selectivity, and roughly isoenergetic profiles for K(+) and Rb(+), and discriminates against Cs(+), all in agreement with experimental data. It also indicates that Ba(2+) and Na(+) compete effectively with permeant ions at a site near the boundary between the filter and the cavity, in the vicinity of the barium blocker site.

Potassium permeation through the KcsA channel: a density functional study

We present a theoretical study on structural and electronic aspects of K+ permeation through the binding sites of the KcsA channel's selectivity filter. Density functional calculations are carried out on models taken from selected snapshots of a molecular dynamics simulation recently reported [FEBS Lett. 477 (2000) 37]. During the translocation process from one binding site to the other, the coordination number of the permeating K+ ion turns out to decrease and K+ ion polarizes significantly its ligands, backbone carbonyl groups and a water molecule. K+-induced polarization increases significantly at the transition state (TS) between the two binding sites. These findings suggest that polarization effects play a significant role in the microscopic mechanisms regulating potassium permeation.

Computational studies of the gramicidin channel

Ion channels are highly specific membrane-spanning protein structures which serve to facilitate the passage of selected ions across the lipid barrier. In the past decade, molecular dynamics simulations based on atomic models and realistic microscopic interactions with explicit solvent and membrane lipids have been used to gain insight into the function of these complex systems. These calculations have considerably expanded our view of ion permeation at the microscopic level. This Account will mainly focus on computational studies of the gramicidin A channel, one of the simplest and best characterized molecular pore.

Experimental validation of Gaussian-3 lithium cation affinities of amides: implications for the gas-phase lithium cation basicity scale

Using a refined Gaussian-3 (G3) protocol, the highest level of ab initio calculations reported so far, we have established the Li+ cation binding enthalpy (affinity) at 0 K (in kJ mol-1) for formamide (195.7), N-methylformamide (209.2), N,N'-dimethylformamide (220.0), acetamide (211.7), N-methylacetamide (222.5), and N,N'-dimethylacetamide (230.1), with an estimated maximum uncertainty of +/-8 kJ mol-1. With these six theoretical lithium cation binding affinities as reference values, the absolute Li+ affinities of imidazole and dimethoxyethane were determined by the extended kinetic method, and by adopting the statistical data treatment protocol recently proposed by Armentrout. The Li+ affinities obtained for these two ligands are in good agreement (within 6 kJ mol-1) with recent values determined by the threshold collision-induced dissociation method, and consistent with the Li+ basicity values first reported by Taft and co-workers in 1990. Our study confirms that the previously suggested, and recently implemented, downward revision of Taft's original basicity scale by 10.9 kJ mol-1 is justified for ligands with revised basicities less than 151 kJ mol-1. However, for selected ligands with Li+ basicities greater than 151 kJ mol-1, including some of the six amides studied in this work, the reported discrepancy between theoretical and experimental estimates in the revised Li+ basicity scale of Burk et al. is likely to arise from experimental uncertainties.

Molecular dynamics of the KcsA K(+) channel in a bilayer membrane

Molecular dynamics (MD) simulations of an atomic model of the KcsA K(+) channel embedded in an explicit dipalmitoylphosphatidylcholine (DPPC) phospholipid bilayer solvated by a 150 mM KCl aqueous salt solution are performed and analyzed. The model includes the KcsA K(+) channel, based on the recent crystallographic structure of, Science. 280:69-77), 112 DPPC, K(+) and Cl(-) ions, as well as over 6500 water molecules for a total of more than 40,000 atoms. Three K(+) ions are explicitly included in the pore. Two are positioned in the selectivity filter on the extracellular side and one in the large water-filled cavity. Different starting configurations of the ions and water molecules in the selectivity filter are considered, and MD trajectories are generated for more than 4 ns. The conformation of KcsA is very stable in all of the trajectories, with a global backbone root mean square (RMS) deviation of less than 1.9 A with respect to the crystallographic structure. The RMS atomic fluctuations of the residues surrounding the selectivity filter on the extracellular side of the channel are significantly lower than those on the intracellular side. The motion of the residues with aromatic side chains surrounding the selectivity filter (Trp(67), Trp(68), Tyr(78), and Tyr(82)) is anisotropic with the smallest RMS fluctuations in the direction parallel to the membrane plane. A concerted dynamic transition of the three K(+) ions in the pore is observed, during which the K(+) ion located initially in the cavity moves into the narrow part of the selectivity filter, while the other two K(+) ions move toward the extracellular side. A single water molecule is stabilized between each pair of ions during the transition, suggesting that each K(+) cation translocating through the narrow pore is accompanied by exactly one water molecule, in accord with streaming potential measurements (, Biophys. J. 55:367-371). The displacement of the ions is coupled with the structural fluctuations of Val(76) and Gly(77), in the selectivity filter, as well as the side chains of Glu(71), Asp(80), and Arg(89), near the extracellular side. Thus the mechanical response of the channel structure at distances as large as 10-20 A from the ions in the selectivity filter appears to play an important role in the concerted transition.

The binding site of sodium in the gramicidin A channel

The available information concerning the structure and location of the main binding site for sodium in the gramicidin A channel is reviewed and discussed. Results from molecular dynamics simulations using an atomic model of the channel embedded in a lipid bilayer are compared with experimental observations. The combined information from experiment and simulation suggests that the main binding sites for sodium are near the channel's mouth, approximately 9.2 A from the centre of the dimer channel, although the motion along the axis could be as large as 1 to 2 A. In the binding site, the sodium ion is lying off axis, making contact with two carbonyl oxygens and two single-file water molecules. The main channel ligand is provided by the carbonyl group of the Leu10-Trp11 peptide linkage, which exhibits the largest deflection from the ion-free channel structure.

Ionic interactions in multiply occupied channels

A significant number of physiologically important ion channels function via multi-ion mechanisms where repulsion between ions at slightly separated locations is believed to be critical for permeation. We apply the semi-microscopic Monte Carlo approach and analyse how multiple occupancy affects permeation energetics and ion-water-peptide correlations. We consider double occupancy in idealized models of two systems: gramicidin A and the KcsA K+ channel. We focus on the excess repulsion energy due to ion-water and ion-peptide correlations (repulsion energy adjusted for direct ion-ion interaction). Gramicidin, where multiple occupancy is marginally important functionally, is ideal for correlating structure and ion interactions. Pair occupancy is stabilized by interaction with bulk solvent, destabilized by interaction with both the channel water and, as binding sites are far apart, the peptide backbone. In the KcsA K+ channel, double occupancy is promoted by the uneven spacing and the large ion-water separations in the selectivity filter. The carbonyls forming the binding cavities are equally important for pair stabilization. Due to the binding pocket's design, net ionic repulsion is approximately 25-30% of what it would be in a gramicidin-like structure with the same interionic spacing.