Representation of Ion-Protein Interactions Using the Drude Polarizable Force-Field

A non-B DNA binding peptidomimetic channel alters cellular functions

DNA binding transcription factors possess the ability to interact with lipid membranes to construct ion-permeable pathways. Herein, we present a thiazole-based DNA binding peptide mimic TBP2, which forms transmembrane ion channels, impacting cellular ion concentration and consequently stabilizing G-quadruplex DNA structures. TBP2 self-assembles into nanostructures, e.g., vesicles and nanofibers and facilitates the transportation of Na+ and K+ across lipid membranes with high conductance (~0.6 nS). Moreover, TBP2 exhibits increased fluorescence when incorporated into the membrane or in cellular nuclei. Monomeric TBP2 can enter the lipid membrane and localize to the nuclei of cancer cells. The coordinated process of time-dependent membrane or nuclear localization of TBP2, combined with elevated intracellular cation levels and direct G-quadruplex (G4) interaction, synergistically promotes formation and stability of G4 structures, triggering cancer cell death. This study introduces a platform to mimic and control intricate biological functions, leading to the discovery of innovative therapeutic approaches.

FFParam-v2.0: A Comprehensive Tool for CHARMM Additive and Drude Polarizable Force-Field Parameter Optimization and Validation

Developing production quality CHARMM force-field (FF) parameters is a very detailed process involving a variety of calculations, many of which are specific for the molecule of interest. The first version of FFParam was developed as a standalone Python package designed for the optimization of electrostatic and bonded parameters of the CHARMM additive and polarizable Drude FFs by using quantum mechanical (QM) target data. The new version of FFParam has multiple new capabilities for FF parameter optimization and validation, with an emphasis on the ability to use condensed-phase target data in optimization. FFParam-v2 allows optimization of Lennard-Jones (LJ) parameters using potential energy scans of interactions between selected atoms in a molecule and noble gases, viz., He and Ne, and through condensed-phase calculations, from which experimental observables such as heats of vaporization and free energies of solvation may be obtained. This functionality serves as a gold standard for both optimizing parameters and validating the performance of the final parameters. A new bonded parameter optimization algorithm has been introduced to account for simultaneously optimizing multiple molecules sharing parameters. FFParam-v2 also supports the comparison of normal modes and the potential energy distribution of internal coordinates towards each normal mode obtained from QM and molecular mechanics calculations. Such comparison capability is vital to validate the balance among various bonded parameters that contribute to the complex normal modes of molecules. User interaction has been extended beyond the original graphical user interface to include command-line interface capabilities that allow for integration of FFParam in workflows, thereby facilitating the automation of parameter optimization. With these new functionalities, FFParam is a more comprehensive parameter optimization tool for both beginners and advanced users.

Balancing Group I Monatomic Ion-Polar Compound Interactions for Condensed Phase Simulation in the Polarizable Drude Force Field

Molecular dynamics (MD) simulations are a commonly used method for investigating molecular behavior at the atomic level. Achieving reliable MD simulation results necessitates the use of an accurate force field. In the present work, we present a protocol to enhance the quality of group 1 monatomic ions (specifically Li+, Na+, K+, Rb+, and Cs+) with respect to their interactions with common polar model compounds in biomolecules in condensed phases in the context of the Drude polarizable force field. Instead of adjusting preexisting individual parameters for ions, model compounds, and water, we employ atom-pair specific Lennard-Jones (LJ) (known as NBFIX in CHARMM) and through-space Thole dipole screening (NBTHOLE) terms to fine-tune the balance of ion-model compound, ion-water, and model compound-water interactions. This involved establishing a protocol for the optimization of NBFIX and NBTHOLE parameters targeting the difference between molecular mechanical (MM) and quantum mechanical (QM) potential energy scans (PES). It is shown that targeting PES involving complexes that include multiple model compounds and/or ions as trimers and tetramers yields parameters that produce condensed phase properties in agreement with experimental data. Validation of this protocol involved the reproduction of experimental thermodynamic benchmarks, including solvation free energies of ions in methanol and N-methylacetamide, osmotic pressures, ionic conductivities, and diffusion coefficients within the condensed phase. These results show the importance of including more complex ion-model compound complexes beyond dimers in the QM target data to account for many-body effects during parameter fitting. The presented parameters represent a significant refinement of the Drude polarizable force field, which will lead to improved accuracy for modeling ion-biomolecular interactions.

System-Specific Parameter Optimization for Nonpolarizable and Polarizable Force Fields

The accuracy of classical force fields (FFs) has been shown to be limited for the simulation of cation-protein systems despite their importance in understanding the processes of life. Improvements can result from optimizing the parameters of classical FFs or by extending the FF formulation by terms describing charge transfer (CT) and polarization (POL) effects. In this work, we introduce our implementation of the CTPOL model in OpenMM, which extends the classical additive FF formula by adding CT and POL. Furthermore, we present an open-source parametrization tool, called FFAFFURR, that enables the (system-specific) parametrization of OPLS-AA and CTPOL models. The performance of our workflow was evaluated by its ability to reproduce quantum chemistry energies and by molecular dynamics simulations of a zinc-finger protein.

A Fast, Convenient, Polarizable Electrostatic Model for Molecular Dynamics

We present an efficient polarizable electrostatic model, utilizing typed, atom-centered polarizabilities and the fast direct approximation, designed for efficient use in molecular dynamics (MD) simulations. The model provides two convenient approaches for assigning partial charges in the context of atomic polarizabilities. One is a generalization of RESP, called RESP-dPol, and the other, AM1-BCC-dPol, is an adaptation of the widely used AM1-BCC method. Both are designed to accurately replicate gas-phase quantum mechanical electrostatic potentials. Benchmarks of this polarizable electrostatic model against gas-phase dipole moments, molecular polarizabilities, bulk liquid densities, and static dielectric constants of organic liquids show good agreement with the reference values. Of note, the model yields markedly more accurate dielectric constants of organic liquids, relative to a matched nonpolarizable force field. MD simulations with this method, which is currently parametrized for molecules containing elements C, N, O, and H, run only about 3.6-fold slower than fixed charge force fields, while simulations with the self-consistent mutual polarization average 4.5-fold slower. Our results suggest that RESP-dPol and AM1-BCC-dPol afford improved accuracy relative to fixed charge force fields and are good starting points for developing general, affordable, and transferable polarizable force fields. The software implementing these approaches has been designed to utilize the force field fitting frameworks developed and maintained by the Open Force Field Initiative, setting the stage for further exploration of this approach to polarizable force field development.

Polarizable AMOEBA Model for Simulating Mg2+·Protein·Nucleotide Complexes

Molecular mechanics (MM) simulations have the potential to provide detailed insights into the mechanisms of enzymes that utilize nucleotides as cofactors. In most cases, the activities of these enzymes also require the binding of divalent cations to catalytic sites. However, modeling divalent cations in MM simulations has been challenging. The inclusion of explicit polarization was considered promising, but despite improvements over nonpolarizable force fields and despite the inclusion of "Nonbonded-fix (NB-fix)" corrections, errors in interaction energies of divalent cations with proteins remain large. Importantly, the application of these models fails to reproduce the experimental structural data on Mg2+·Protein·ATP complexes. Focusing on these complexes, here we provide a systematic assessment of the polarizable AMOEBA model and recommend critical changes that substantially improve its predictive performance. Our key results are as follows. We first show that our recent revision of the AMOEBA protein model (AMOEBABIO18-HFC), which contains high field corrections (HFCs) to induced dipoles, dramatically improves Mg2+-protein interaction energies, reducing the mean absolute error (MAE) from 17 to 10 kcal/mol. This further supports the general applicability of AMOEBABIO18-HFC. The inclusion of many-body NB-fix corrections further reduces MAE to 6 kcal/mol, which amounts to less than 2% error. The errors are estimated with respect to vdW-inclusive density functional theory that we benchmark against CCSD(T) calculations and experiments. We also present a new model of ATP with revised polarization parameters to better capture its high field response, as well as new vdW and dihedral parameters. The ATP model accurately predicts experimental Mg2+-ATP binding free energy in the aqueous phase and provides new insights into how Mg2+ associates with ATP. Finally, we show that molecular dynamics (MD) simulations of Mg2+·Kinase·ATP complexes carried out with these improvements lead to a better agreement in global and local catalytic site structures between MD and X-ray crystallography.

Water inside the Selectivity Filter of a K+ Ion Channel: Structural Heterogeneity, Picosecond Dynamics, and Hydrogen Bonding

Water inside biological ion channels regulates the key properties of these proteins, such as selectivity, ion conductance, and gating. In this article, we measure the picosecond spectral diffusion of amide I vibrations of an isotope-labeled KcsA potassium channel using two-dimensional infrared (2D IR) spectroscopy. By combining waiting time (100-2000 fs) 2D IR measurements of the KcsA channel including 13C18O isotope-labeled Val76 and Gly77 residues with molecular dynamics simulations, we elucidated the site-specific dynamics of water and K+ ions inside the selectivity filter of KcsA. We observe inhomogeneous 2D line shapes with extremely slow spectral diffusion. Our simulations quantitatively reproduce the experiments and show that water is the only component with any appreciable dynamics, whereas K+ ions and the protein are essentially static on a picosecond timescale. By analyzing simulated and experimental vibrational frequencies, we find that water in the selectivity filter can be oriented to form hydrogen bonds with adjacent or nonadjacent carbonyl groups with the reorientation timescales being three times slower and comparable to that of water molecules in liquid, respectively. Water molecules can reside in the cavity sufficiently far from carbonyls and behave essentially like "free" gas-phase-like water with fast reorientation times. Remarkably, no interconversion between these configurations was observed on a picosecond timescale. These dynamics are in stark contrast with liquid water, which remains highly dynamic even in the presence of ions at high concentrations.

Development of a QM/MM(ABEEM) method for the deprotonation of neutral and cation radicals in the G-tetrad and GGX(8-oxo-G) tetrad

The rapid deprotonation of G˙+ in the DNA strand impedes positive charge (hole) transfer, whereas the slow deprotonation rate of G˙+ in the G-tetrad makes it a more suitable carrier for hole conduction. The QM/MM(ABEEM) combined method, which involves the integration of QM and the ABEEM polarizable force field (ABEEM PFF), was developed to investigate the deprotonation of neutral and cation free radicals in the G-tetrad and GGX(8-oxo-G) tetrad (xanthine and 8-oxoguanine dual substituted G-tetrad). By incorporating valence-state electronegativity piecewise functions χ*(r) and implementing charge local conservation conditions, QM/MM(ABEEM) possesses the advantage of accurately simulating charge transfer and polarization effect during deprotonation. The activation energy calculated by the QM method of X˙ is the lowest among other bases in the GGX(8-oxo-G) tetrad, which is supported by the computation of the average electronegativity calculated by ABEEM PFF. By utilizing QM/MM(ABEEM) with a two-way free energy perturbation method, the deprotonation activation energy of X˙ in the GGX(8-oxo-G) tetrad is determined to be 33.0 ± 2.1 kJ mol-1, while that of G˙+ in the G-tetrad is 20.7 ± 0.6 kJ mol-1, consistent with the experimental measurement of 20 ± 1.0 kJ mol-1. These results manifest that X˙ in the GGX(8-oxo-G) tetrad exhibits a slower deprotonation rate than G˙+ in the G-tetrad, suggesting that the GGX(8-oxo-G) tetrad may serve as a more favorable hole transport carrier. Furthermore, the unequal average electronegativities of bases in the GGX(8-oxo-G) tetrad impede the deprotonation rate. This study provides a potential foundation for investigating the microscopic mechanism of DNA electronic devices.

Water inside the selectivity filter of a K+ ion channel: structural heterogeneity, picosecond dynamics, and hydrogen-bonding

Water inside biological ion channels regulates the key properties of these proteins such as selectivity, ion conductance, and gating. In this Article we measure the picosecond spectral diffusion of amide I vibrations of an isotope labeled KcsA potassium channel using two-dimensional infrared (2D IR) spectroscopy. By combining waiting time (100 - 2000 fs) 2D IR measurements of the KcsA channel including 13C18O isotope labeled Val76 and Gly77 residues with molecular dynamics simulations, we elucidated the site-specific dynamics of water and K+ ions inside the selectivity filter of KcsA. We observe inhomogeneous 2D lineshapes with extremely slow spectral diffusion. Our simulations quantitatively reproduce the experiments and show that water is the only component with any appreciable dynamics, whereas K+ ions and the protein are essentially static on a picosecond timescale. By analyzing simulated and experimental vibrational frequencies, we find that water in the selectivity filter can be oriented to form hydrogen bonds with adjacent, or non-adjacent carbonyl groups with the reorientation timescales being three times slower and comparable to that of water molecules in liquid, respectively. Water molecules can reside in the cavity sufficiently far from carbonyls and behave essentially like "free" gas-phase-like water with fast reorientation times. Remarkably, no interconversion between these configurations were observed on a picosecond timescale. These dynamics are in stark contrast with liquid water that remains highly dynamic even in the presence of ions at high concentrations.

From bonds to interactions: comprehensive molecular characterization via polarizable bond-dipole approach

Accurately characterizing molecular interactions stands as a pivotal requirement for ensuring the reliability of molecular dynamics simulations. In line with our bond-dipole-based interaction model proposed by Gao et al. [X.-C. Gao, Q. Hao and C.-S. Wang, J. Chem. Theory Comput., 2017, 13, 2730-2741.], we have implemented an efficient and concise approach to compute electrostatic potential. This methodology capitalizes on the polarizable nature of chemical bond dipoles, resulting in a model of remarkable simplicity. In this study, we have revised the polarizable bond-dipole-based force field (PBFF) through the meticulous curation of quantum chemical data sets. These data sets encompass a comprehensive collection of 40 000 conformations, including those of water, methylamine, methanol, and N-methylacetamide. Additionally, we incorporate 520 hydrogen-bonded dimers into our data sets. In pursuit of enhanced accuracy in molecular dynamics simulations and a more faithful representation of potential energy landscapes, we undertook the re-optimization of the nonbonded parameters within the PBFF framework. Concurrently, we intricately fine-tuned the bonded parameters. The results of our comprehensive evaluation denote that this newly optimized force field method adeptly and efficiently computes structural characteristics, harmonic frequencies, and interaction energies. Overall, this study provides further validation for the applicability of PBFF in molecular dynamics simulations.

Molecular Modeling of Water-in-Salt Electrolytes: A Comprehensive Analysis of Polarization Effects and Force Field Parameters in Molecular Dynamics Simulations

Accurate modeling of highly concentrated aqueous solutions, such as water-in-salt (WiS) electrolytes in battery applications, requires proper consideration of polarization contributions to atomic interactions. Within the force field molecular dynamics (MD) simulations, the atomic polarization can be accounted for at various levels. Nonpolarizable force fields implicitly account for polarization effects by incorporating them into their van der Waals interaction parameters. They can additionally mimic electron polarization within a mean-field approximation through ionic charge scaling. Alternatively, explicit polarization description methods, such as the Drude oscillator model, can be selectively applied to either a subset of polarizable atoms or all polarizable atoms to enhance simulation accuracy. The trade-off between simulation accuracy and computational efficiency highlights the importance of determining an optimal level of accounting for atomic polarization. In this study, we analyze different approaches to include polarization effects in MD simulations of WiS electrolytes, with an example of a Na-OTF solution. These approaches range from a nonpolarizable to a fully polarizable force field. After careful examination of computational costs, simulation stability, and feasibility of controlling the electrolyte properties, we identify an efficient combination of force fields: the Drude polarizable force field for salt ions and non-polarizable models for water. This cost-effective combination is sufficiently flexible to reproduce a broad range of electrolyte properties, while ensuring simulation stability over a relatively wide range of force field parameters. Furthermore, we conduct a thorough evaluation of the influence of various force field parameters on both the simulation results and technical requirements, with the aim of establishing a general framework for force field optimization and facilitating parametrization of similar systems.

Probing Ion Configurations in the KcsA Selectivity Filter with Single-Isotope Labels and 2D IR Spectroscopy

The potassium ion (K+) configurations of the selectivity filter of the KcsA ion channel protein are investigated with two-dimensional infrared (2D IR) spectroscopy of amide I vibrations. Single 13C-18O isotope labels are used, for the first time, to selectively probe the S1/S2 or S2/S3 binding sites in the selectivity filter. These binding sites have the largest differences in ion occupancy in two competing K+ transport mechanisms: soft-knock and hard-knock. According to the former, water molecules alternate between K+ ions in the selectivity filter while the latter assumes that K+ ions occupy the adjacent sites. Molecular dynamics simulations and computational spectroscopy are employed to interpret experimental 2D IR spectra. We find that in the closed conductive state of the KcsA channel, K+ ions do not occupy adjacent binding sites. The experimental data is consistent with simulated 2D IR spectra of soft-knock ion configurations. In contrast, the simulated spectra for the hard-knock ion configurations do not reproduce the experimental results. 2D IR spectra of the hard-knock mechanism have lower frequencies, homogeneous 2D lineshapes, and multiple peaks. In contrast, ion configurations of the soft-knock model produce 2D IR spectra with a single peak at a higher frequency and inhomogeneous lineshape. We conclude that under equilibrium conditions, in the absence of transmembrane voltage, both water and K+ ions occupy the selectivity filter of the KcsA channel in the closed conductive state. The ion configuration is central to the mechanism of ion transport through potassium channels.

Symmetrized Drude Oscillator Force Fields Improve Numerical Performance of Polarizable Molecular Dynamics

Drude oscillator potentials are a popular and computationally efficient class of polarizable models that represent each polarizable atom as a positively charged Drude core harmonically bound to a negatively charged Drude shell. We show that existing force fields that place all non-Coulomb forces on the Drude core and none on the shell inadvertently couple the dipole to non-Coulombic forces. This introduces errors where interactions with neutral particles can erroneously induce atomic polarization, leading to spurious polarizations in the absence of an electric field, exacerbating violations of equipartition in the employed Carr-Parinello scheme. A suitable symmetrization of the interaction potential that correctly splits the force between the Drude core and shell can correct this shortcoming, improving the stability and numerical performance of Drude oscillator-based simulations. The symmetrization procedure is straightforward and only requires the rescaling of a few force field parameters.

A cooperative knock-on mechanism underpins Ca2+-selective cation permeation in TRPV channels

The selective exchange of ions across cellular membranes is a vital biological process. Ca2+-mediated signaling is implicated in a broad array of physiological processes in cells, while elevated intracellular concentrations of Ca2+ are cytotoxic. Due to the significance of this cation, strict Ca2+ concentration gradients are maintained across the plasma and organelle membranes. Therefore, Ca2+ signaling relies on permeation through selective ion channels that control the flux of Ca2+ ions. A key family of Ca2+-permeable membrane channels is the polymodal signal-detecting transient receptor potential (TRP) ion channels. TRP channels are activated by a wide variety of cues including temperature, small molecules, transmembrane voltage, and mechanical stimuli. While most members of this family permeate a broad range of cations non-selectively, TRPV5 and TRPV6 are unique due to their strong Ca2+ selectivity. Here, we address the question of how some members of the TRPV subfamily show a high degree of Ca2+ selectivity while others conduct a wider spectrum of cations. We present results from all-atom molecular dynamics simulations of ion permeation through two Ca2+-selective and two non-selective TRPV channels. Using a new method to quantify permeation cooperativity based on mutual information, we show that Ca2+-selective TRPV channel permeation occurs by a three-binding site knock-on mechanism, whereas a two-binding site knock-on mechanism is observed in non-selective TRPV channels. Each of the ion binding sites involved displayed greater affinity for Ca2+ over Na+. As such, our results suggest that coupling to an extra binding site in the Ca2+-selective TRPV channels underpins their increased selectivity for Ca2+ over Na+ ions. Furthermore, analysis of all available TRPV channel structures shows that the selectivity filter entrance region is wider for the non-selective TRPV channels, slightly destabilizing ion binding at this site, which is likely to underlie mechanistic decoupling.

Probing ion binding in the selectivity filter of the Cav1.1 channel with molecular dynamics

Cav1.1 is the voltage-gated calcium channel essential for the contraction of skeletal muscles upon membrane potential changes. Structural determination of the Cav1.1 channel opens the avenue toward understanding of the structure-function relationship of voltage-gated calcium channels. Here, we show that there exist two Ca2+-binding sites, termed S1 and S2, within the selectivity filter of Cav1.1 through extensive molecular dynamics simulations on various initial ion arrangement configurations. The formation of both binding sites is associated with the four Glu residues (Glu292/614/1014/1323) that constitute the so-called EEEE locus. At the S1 site near the extracellular side, the Ca2+ ion is coordinated with the negatively charged carboxylic groups of these Glu residues and of the Asp615 residue either in a direct way or via an intermediate water molecule. At the S2 site, Ca2+ binding shows two distinct states: an upper state involving two out of the four Glu residues in the EEEE locus and a lower state involving only one Glu residue. In addition, there exist two recruitment sites for Ca2+ above the entrance of the filter. These findings promote the understanding of mechanism for ion permeation and selectivity in calcium channels.

Conformational dynamics and putative substrate extrusion pathways of the N-glycosylated outer membrane factor CmeC from Campylobacter jejuni

The outer membrane factor CmeC of the efflux machinery CmeABC plays an important role in conferring antibiotic and bile resistance to Campylobacter jejuni. Curiously, the protein is N-glycosylated, with the glycans playing a key role in the effective function of this system. In this work we have employed atomistic equilibrium molecular dynamics simulations of CmeC in a representative model of the C. jejuni outer membrane to characterise the dynamics of the protein and its associated glycans. We show that the glycans are more conformationally labile than had previously been thought. The extracellular loops of CmeC visit the open and closed states freely suggesting the absence of a gating mechanism on this side, while the narrow periplasmic entrance remains tightly closed, regulated via coordination to solvated cations. We identify several cation binding sites on the interior surface of the protein. Additionally, we used steered molecular dynamics simulations to elucidate translocation pathways for a bile acid and a macrolide antibiotic. These, and additional equilibrium simulations suggest that the anionic bile acid utilises multivalent cations to climb the ladder of acidic residues that line the interior surface of the protein.

Inclusion of High-Field Target Data in AMOEBA's Calibration Improves Predictions of Protein-Ion Interactions

The reliability of molecular mechanics simulations to predict effects of ion binding to proteins depends on their ability to simultaneously describe ion-protein, ion-water, and protein-water interactions. Force fields (FFs) to describe protein-water and ion-water interactions have been constructed carefully and have also been refined routinely to improve accuracy. Descriptions for ion-protein interactions have also been refined, although in an a posteriori manner through the use of "nonbonded-fix (NB-fix)" approaches in which parameters from default Lennard-Jones mixing rules are replaced with those optimized against some reference data. However, even after NB-fix corrections, there remains a significant need for improvement. This is also true for polarizable FFs that include self-consistent inducible moments. Our recent studies on the polarizable AMOEBA FF suggested that the problem associated with modeling ion-protein interactions could be alleviated by recalibrating polarization models of cation-coordinating functional groups so that they respond better to the high electric fields present near ions. Here, we present such a recalibration of carbonyls, carboxylates, and hydroxyls in the AMOEBA protein FF and report that it does improve predictions substantially─mean absolute errors in Na⁺-protein and K⁺-protein interaction energies decrease from 8.7 to 5.3 and 9.6 to 6.3 kcal/mol, respectively. Errors are computed with respect to estimates from van der Waals-inclusive density functional theory benchmarked against high-level quantum mechanical calculations and experiments. While recalibration does improve ion-protein interaction energies, they still remain underestimated, suggesting that further improvements can be made in a systematic manner through modifications in classical formalism. Nevertheless, we show that by applying our many-body NB-fix correction to Lennard-Jones components, these errors are further reduced to 2.7 and 2.6 kcal/mol, respectively, for Na⁺ and K⁺ ions. Finally, we show that the recalibrated AMOEBA protein FF retains its intrinsic reliability in predicting protein structure and dynamics in the condensed phase.

Millisecond molecular dynamics simulations of KRas-dimer formation and interfaces

Ras dimers have been proposed as building blocks for initiating the extracellular signal-regulated kinase (ERK)/mitogen-activated protein kinase (MAPK) cellular signaling pathway. To better examine the structure of possible dimer interfaces, the dynamics of Ras dimerization, and its potential signaling consequences, we performed molecular dynamics simulations totaling 1 ms of sampling, using an all-atom model of two full-length, farnesylated, guanosine triphosphate (GTP)-bound, wild-type KRas4b proteins diffusing on 29%POPS (1-palmitoyl-2-oleoyl-sn-glycero-3-phospho-L-serine)-mixed POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) membranes. Our simulations unveil an ensemble of thermodynamically weak KRas dimers spanning multiple conformations. The most stable conformations, having the largest interface areas, involve helix α2 and a hypervariable region (HVR). Among the dimer conformations, we found that the HVR of each KRas has frequent interactions with various parts of the dimer, thus potentially mediating the dimerization. Some dimer configurations have one KRas G-domain elevated above the lipid bilayer surface by residing on top of the other G-domain, thus likely contributing to the recruitment of cytosolic Raf kinases in the context of a stably formed multi-protein complex. We identified a variant of the α4-α5 KRas-dimer interface that is similar to the interfaces obtained with fluorescence resonance energy transfer (FRET) data of HRas on lipid bilayers. Interestingly, we found two arginine fingers, R68 and R149, that directly interact with the beta-phosphate of the GTP bound in KRas, in a manner similar to what is observed in a crystal structure of GAP-HRas complex, which can facilitate the GTP hydrolysis via the arginine finger of GTPase-activating protein (GAP).

Molecular Dynamics Investigation into pH Dependent Metal Binding of the Intrinsically Disordered Worm Jaw Protein, Nvjp-1

Sclerotization of the Nereis virens jaw is mediated by metal binding to the histidine-rich jaw protein, Nvjp-1. Previous studies showed that the mechanical properties of Nvjp-1 hydrogels could be modulated with zinc binding as well as the associated anion. Here, we show that the mechanical properties of Nvjp-1 hydrogels can be modulated by pH and that zinc binding to Nvjp-1 is stable at both acidic and alkaline pH conditions. To probe the mechanism of Zn²⁺ binding to Nvjp-1 at different pH conditions, we utilized all atom molecular dynamics simulations employing a polarizable force field. At low pH conditions, polar residues predominantly interacted with Zn²⁺, with at most two residues interacting with a given zinc ion. Surprisingly, little to no Zn²⁺ binding was observed with the abundant Nvjp-1 acidic residues, which form salt-bridges with the protonated histidines to effectively block their binding to Zn²⁺ ions. As the pH was shifted to alkaline conditions, Zn²⁺ binding residues reconfigured to form additional coordination bonds with histidine, resulting in a reduction in the radius of gyration that correlated with hydrogel sclerotization. Furthermore, acetate ions were shown to facilitate the capture of zinc ions through association with protonated histidines at low pH, freeing acidic residues to interact with Zn²⁺ ions and increasing the number of Zn²⁺ ions that diffuse into the Nvjp-1 interior. Thus, these studies provide valuable molecular insights into how amino acid residues in Nvjp-1 manage metal salt binding and coordination in hydrogels as a function of the pH and ionic environments.

Multiscale molecular modelling: from electronic structure to dynamics of nanosystems and beyond

Important contemporary biological and materials problems often depend on interactions that span orders of magnitude differences in spatial and temporal dimensions. This Tutorial Review attempts to provide an introduction to such fascinating problems through a series of case studies, aimed at beginning researchers, graduate students, postdocs and more senior colleagues who are changing direction to focus on multiscale aspects of their research. The choice of specific examples is highly personal, with examples either chosen from our own work or outstanding multiscale efforts from the literature. I start with various embedding schemes, as exemplified by polarizable continuum models, 3-D RISM, molecular DFT and frozen-density embedding. Next, QM/MM (quantum mechanical/molecular mechanical) techniques are the workhorse of pm-to-nm/ps-to-ns simulations; examples are drawn from enzymes and from nanocatalysis for oil-sands upgrading. Using polarizable force-fields in the QM/MM framework represents a burgeoning subfield; with examples from ion channels and electron dynamics in molecules subject to strong external fields, probing the atto-second dynamics of the electrons with RT-TDDFT (real-time - time-dependent density functional theory) eventually coupled with nuclear motion through the Ehrenfest approximation. This is followed by a section on coarse graining, bridging dimensions from atoms to cells. The penultimate chapter gives a quick overview of multiscale approaches that extend into the meso- and macro-scales, building on atomistic and coarse-grained techniques to enter the world of materials engineering, on the one hand, and cell biology, on the other. A final chapter gives just a glimpse of the burgeoning impact of machine learning on the structure-dynamics front. I aim to capture the excitement of contemporary leading-edge breakthroughs in the description of physico-chemical systems and processes in complex environments, with only enough historical content to provide context and aid the next generation of methodological development. While I aim also for a clear description of the essence of methodological breakthroughs, equations are kept to a minimum and detailed formalism and implementation details are left to the references. My approach is very selective (case studies) rather than exhaustive. I think that these case studies should provide fodder to build as complete a reference tree on multiscale modelling as the reader may wish, through forward and backward citation analysis. I hope that my choices of cases will excite interest in newcomers and help to fuel the growth of multiscale modelling in general.

Harnessing Deep Learning for Optimization of Lennard-Jones Parameters for the Polarizable Classical Drude Oscillator Force Field

The outcomes of computational chemistry and biology research, including drug design, are significantly influenced by the underlying force field (FF) used in molecular simulations. While improved FF accuracy may be achieved via inclusion of explicit treatment of electronic polarization, such an extension must be accompanied by optimization of van der Waals (vdW) interactions, in the context of the Lennard-Jones (LJ) formalism in the present study. This is particularly challenging due to the extensive nature of chemical space combined with the correlated nature of LJ parameters. To address this challenge, a deep learning (DL)-based parametrization framework is developed, allowing for sampling of wide ranges of LJ parameters targeting experimental condensed phase thermodynamic properties. The present work utilizes this framework to develop the LJ parameters for atoms associated with four distinct groups covering 10 different atom types. Final parameter selection was facilitated by quantum mechanical data on rare-gas interactions with the training set molecules. The chosen parameters were then validated through experimental hydration free energies and condensed phase thermodynamic properties of validation set molecules to confirm transferability. The ultimate outcome of utilizing this framework is a set of LJ parameters in the context of the polarizable Drude FF, which demonstrated improvement in the reproduction of both experimental pure solvent and crystal properties and hydration free energies of the molecules compared to the additive CHARMM General FF (CGenFF) including the ability of the Drude FF to accurately reproduce both experimental pure solvent properties and hydration free energies. The study also shows how correlations between difference in the reproduction of condensed phase data between model compounds may be used to direct the selection of new atom types and training set molecules during FF development.

High-Dimensional Parameter Search Method to Determine Force Field Mixing Terms in Molecular Simulations

Molecular dynamics (MD) force fields for lipids and ions are typically developed independently of one another. In simulations consisting of both lipids and ions, lipid-ion interaction energies are estimated using a predefined set of mixing rules for Lennard-Jones (LJ) interactions. This, however, does not guarantee their reliability. In fact, compared to the quantum mechanical reference data, Lorentz-Berthelot mixing rules substantially underestimate the binding energies of Na⁺ ions with small-molecule analogues of lipid headgroups, yielding errors on the order of 80 and 130 kJ/mol, respectively, for methyl acetate and diethyl phosphate. Previously, errors associated with mixing force fields have been reduced using approaches such as "NB-fix" in which LJ interactions are computed using explicit cross terms rather than those from mixing rules. Building on this idea, we derive explicit lipid-ion cross terms that also may implicitly include many-body cooperativity effects. Additionally, to account for the interdependency between cross terms, we optimize all cross terms simultaneously by performing high-dimensional searches using our ParOpt software. The cross terms we obtain reduce the errors due to mixing rules to below 10 kJ/mol. MD simulation of the lipid bilayer conducted using these optimized cross terms resolves the structural discrepancies between our previous simulations and small-angle X-ray and neutron scattering experiments. These results demonstrate that simulations of lipid bilayers with ions that are accurate up to structural data from scattering experiments can be performed without explicit polarization terms. However, it is worth noting that such NB-fix cross terms are not based on any physical principle; a polarizable lipid model would be more realistic and is still desired. Our approach is generic and can be applied to improve the accuracies of simulations employing mixed force fields.

Binding Energy and Free Energy of Calcium Ion to Calmodulin EF-Hands with the Drude Polarizable Force Field

Calcium ions are important messenger molecules in cells, which bind calcium-binding proteins to trigger many biochemical processes. We constructed four model systems, each containing one EF-hand loop of calmodulin with one calcium ion bound, and investigated the binding energy and free energy of Ca²⁺ by the quantum mechanics symmetry-adapted perturbation theory (SAPT) method and the molecular mechanics with the additive CHARMM36m (C36m) and the polarizable Drude force fields (FFs). Our results show that the explicit introduction of polarizability in the Drude not only yields considerably improved agreement with the binding energy calculated from the SAPT method but is also able to capture each component of the binding energies including electrostatic, induction, exchange, and dispersion terms. However, binding free energies computed with the Drude and the C36m FFs both deviated significantly from the experimental measurements. Detailed analysis indicated that one of main reasons might be that the strong interactions between Ca²⁺ and the side chain nitrogen of Asn/Gln in the Drude FF caused the distorted coordination geometries of calcium. Our work illustrated the importance of polarization in modeling ion-protein interactions and the difficulty in generating accurate and balanced FF models to represent the polarization effects.

Global Optimization of the Lennard-Jones Parameters for the Drude Polarizable Force Field

Molecular dynamics (MD) simulations based on atomic models play an important role in the drug-discovery process to screen molecules, estimate binding free energies, and optimize lead compounds in chemical space. Accurate computations of thermodynamic and kinetic properties using MD simulations are highly dependent on the accuracy of the underlying atomic force field. In this context, going beyond the nonpolarizable fixed-charge model by accounting explicitly for induced polarization is highly desirable. The CHARMM polarizable force field based on classical Drude oscillators, in which an auxiliary charged particle is attached via a harmonic spring to its parent nucleus, offers both a computationally convenient and rigorous framework to model explicitly induced electronic polarization in MD simulations. For any molecule of interest, electrostatic partial charges, atomic polarizabilities, and Thole shielding factors, as well as bonded parameters can either be determined from ab initio calculations or ascribed from the knowledge-based library of the CHARMM Generalized force field. While this approach is fairly reliable in general, it is well understood that the overall accuracy of the models with respect to thermodynamic properties such as bulk density, enthalpies, and solvation free energies is particularly sensitive to the nonbonded Lennard-Jones (LJ) parameters. In the present study, we systematically refined the set of LJ parameters for the atom types available in the Drude force field to best match the experimental thermodynamic properties for 416 small druglike organic molecules. To further test the transferability of the optimized parameters, the hydration free energy of 372 molecules was computed. The calculations resulted in a small average error of 0.46 kcal/mol and a Pearson R of 0.9, representing a significant improvement over the additive GAFF force field in our previous study, where an average error of ∼2 kcal/mol was obtained. Such an improvement is consistent with the ability of the polarizable Drude model to more accurately model interactions in different environments. The effort provides a roadmap for the global optimization of force field parameters using experimental data. It is hoped that the present effort will further the application of the Drude polarizable force field in molecular simulations including drug design and discovery.

Polarization Effects in Water-Mediated Selective Cation Transport across a Narrow Transmembrane Channel

Despite the progress in modeling complex molecular systems of ever-increasing complexity, a quantitatively accurate computational treatment of ion permeation through narrow membrane channels remains challenging. An important factor to reach this goal is induced electronic polarization, which is likely to impact the permeation rate of small ions through narrow molecular pores. In this work, we extended the recently developed polarizable force field based on the classical Drude oscillators to assess the role of induced polarization effects on the energetics of sodium and potassium ion transport across the gramicidin A (gA) ion channel. The inclusion of induced polarization lowers barriers present in 1D potential of mean force (PMF) for cation permeation by ∼50% compared to those obtained with the additive force field. Conductance properties calculated with 1D PMFs from Drude simulations are in better agreement with experimental results. Polarization of single-file water molecules and protein atoms forming the narrow pore has a direct impact on the free-energy barriers and cation-specific solid-state NMR chemical shifts. Sensitivity analysis indicates that small changes to water-channel interactions can alter the free energy barrier for ion permeation. These results, illustrating polarization effects present in the complex electrostatic environment of the gA channel, have broad implications for revising proposed mechanisms of ion permeation and selectivity in a variety of ion channels.

Transferable interactions of Li+ and Mg2+ ions in polarizable models

Therapeutic implications of Li⁺, in many cases, stem from its ability to inhibit certain Mg²⁺-dependent enzymes, where it interacts with or substitutes for Mg²⁺. The underlying details of its action are, however, unknown. Molecular simulations can provide insights, but their reliability depends on how well they describe relative interactions of Li⁺ and Mg²⁺ with water and other biochemical groups. Here, we explore, benchmark, and recommend improvements to two simulation approaches: the one that employs an all-atom polarizable molecular mechanics (MM) model and the other that uses a hybrid quantum and MM implementation of the quasi-chemical theory (QCT). The strength of the former is that it describes thermal motions explicitly and that of the latter is that it derives local contributions from electron densities. Reference data are taken from the experiment, and also obtained systematically from CCSD(T) theory, followed by a benchmarked vdW-inclusive density functional theory. We find that the QCT model predicts relative hydration energies and structures in agreement with the experiment and without the need for additional parameterization. This implies that accurate descriptions of local interactions are essential. Consistent with this observation, recalibration of local interactions in the MM model, which reduces errors from 10.0 kcal/mol to 1.4 kcal/mol, also fixes aqueous phase properties. Finally, we show that ion-ligand transferability errors in the MM model can be reduced significantly from 10.3 kcal/mol to 1.2 kcal/mol by correcting the ligand's polarization term and by introducing Lennard-Jones cross-terms. In general, this work sets up systematic approaches to evaluate and improve molecular models of ions binding to proteins.

Improved description of ligand polarization enhances transferability of ion-ligand interactions

The reliability of molecular mechanics (MM) simulations in describing biomolecular ion-driven processes depends on their ability to accurately model interactions of ions simultaneously with water and other biochemical groups. In these models, ion descriptors are calibrated against reference data on ion-water interactions, and it is then assumed that these descriptors will also satisfactorily describe interactions of ions with other biochemical ligands. The comparison against the experiment and high-level quantum mechanical data show that this transferability assumption can break down severely. One approach to improve transferability is to assign cross terms or separate sets of non-bonded descriptors for every distinct pair of ion type and its coordinating ligand. Here, we propose an alternative solution that targets an error-source directly and corrects misrepresented physics. In standard model development, ligand descriptors are never calibrated or benchmarked in the high electric fields present near ions. We demonstrate for a representative MM model that when the polarization descriptors of its ligands are improved to respond to both low and high fields, ligand interactions with ions also improve, and transferability errors reduce substantially. In our case, the overall transferability error reduces from 3.3 kcal/mol to 1.8 kcal/mol. These improvements are observed without compromising on the accuracy of low-field interactions of ligands in gas and condensed phases. Reference data for calibration and performance evaluation are taken from the experiment and also obtained systematically from "gold-standard" CCSD(T) in the complete basis set limit, followed by benchmarked vdW-inclusive density functional theory.

Open and Closed Structures of a Barium-Blocked Potassium Channel

Barium (Ba²⁺) is a classic permeant blocker of potassium (K⁺) channels. The "external lock-in effect" in barium block experiments, whereby the binding of external K⁺ impedes the forward translocation of the blocker, provides a powerful avenue to investigate the selectivity of the binding sites along the pore of potassium channels. Barium block experiments show that the external lock-in site is highly selective for K⁺ over Na⁺. Wild-type KcsA was crystallized in low K⁺ conditions, and the crystals were soaked in solutions containing various concentrations of barium. Structural analysis reveals open and closed gate conformations of the KcsA channel. Anomalous diffraction experiments show that Ba²⁺ primarily binds to the innermost site S4 of the selectivity filter of the open-gate conformation and also the site S2, but no binding is detected with the closed-gate conformation. Alchemical free-energy perturbation calculations indicate that the presence of a Ba²⁺ ion in the selectivity filter boosts the specificity of K⁺ binding relative to Na⁺ in the external sites S0-S2.

Accuracy of Molecular Simulation-Based Predictions of koff Values: A Metadynamics Study

The k_off values of ligands unbinding to proteins are key parameters for drug discovery. Their predictions based on molecular simulation may under- or overestimate experiment in a system- and/or technique-dependent way. Here we use an established method-infrequent metadynamics, based on the AMBER force field-to compute the k_off of the ligand iperoxo (in clinical use) targeting the muscarinic receptor M₂. The ligand charges are calculated by either (i) the Amber standard procedure or (ii) B3LYP-DFT. The calculations using (i) turn out not to provide a reasonable estimation of the transition-state free energy. Those using (ii) differ from experiment by 2 orders of magnitude. On the basis of B3LYP DFT QM/MM simulations, we suggest that the observed discrepancy in (ii) arises, at least in part, from the lack of electronic polarization and/or charge transfer in biomolecular force fields. These issues might be present in other systems, such as DNA-protein complexes.

Investigating dual Ca2+ modulation of the ryanodine receptor 1 by molecular dynamics simulation

The ryanodine receptors (RyR) are essential to calcium signaling in striated muscles. A deep understanding of the complex Ca²⁺ -activation/inhibition mechanism of RyRs requires detailed structural and dynamic information for RyRs in different functional states (eg, with Ca²⁺ bound to activating or inhibitory sites). Recently, high-resolution structures of the RyR isoform 1 (RyR1) were solved by cryo-electron microscopy, revealing the location of a Ca²⁺ binding site for activation. Toward elucidating the Ca²⁺ -modulation mechanism of RyR1, we performed extensive molecular dynamics simulation of the core RyR1 structure in the presence and absence of activating and solvent Ca²⁺ (total simulation time is >5 μs). In the presence of solvent Ca²⁺ , Ca²⁺ binding to the activating site enhanced dynamics of RyR1 with higher inter-subunit flexibility, asymmetric inter-subunit motions, outward domain motions and partial pore dilation, which may prime RyR1 for subsequent channel opening. In contrast, the solvent Ca²⁺ alone reduced dynamics of RyR1 and led to inward domain motions and pore contraction, which may cause inhibition. Combining our simulation with the map of disease mutation sites in RyR1, we constructed a wiring diagram of key domains coupled via specific hydrogen bonds involving the mutation sites, some of which were modulated by Ca²⁺ binding. The structural and dynamic information gained from this study will inform future mutational and functional studies of RyR1 activation and inhibition by Ca²⁺ .

Further Optimization and Validation of the Classical Drude Polarizable Protein Force Field

The CHARMM Drude-2013 polarizable force field (FF) was developed to include the explicit treatment of induced electronic polarizability, resulting in a more accurate description of the electrostatic interactions in molecular dynamics (MD) simulations. While the Drude-2013 protein FF has shown success in improving the folding properties of α-helical peptides and to reproduce experimental observables in simulations up to 1 μs, some limitations were noted regarding the stability of β-sheet structures in simulations longer than 100 ns as well as larger deviations from crystal structures in simulations of a number of proteins compared to the additive CHARMM36 protein FF. The origin of the instability has been identified and appears to be primarily due to overestimated atomic polarizabilities and induced dipole-dipole interactions on the Cβ, Cγ, and Cδ side chain atoms. To resolve this and other issues, a number of aspects of the model were revisited, resulting in Drude-2019 protein FF. Backbone parameters were optimized targeting the conformational properties of the (Ala)5 peptide in solution along with gas phase properties of the alanine dipeptide. Dipeptides that contain N-acetylated and N'-methylamidated termini, excluding Gly, Pro, and Ala, were used as models to optimize the atomic polarizabilities and Thole screening factors on selected Cβ, Cγ, and Cδ carbons by targeting quantum mechanical (QM) dipole moments and molecular polarizabilities. In addition, to obtain better conformational properties, side chain χ1 and χ2 dihedral parameters were optimized targeting QM data for the respective side chain dipeptide conformations as well as Protein Data Bank survey data based on the χ1, χ2 sampling from Hamiltonian replica-exchange MD simulations of (Ala)4-X-(Ala)4 in solution, where X is the amino acid of interest. Further improvements include optimizing nonbonded interactions between charged residues to reproduce QM interaction energies of the charged-protein model compounds and experimental osmotic pressures. Validation of the optimized Drude protein FF includes MD simulations of a collection of peptides and proteins including β-sheet structures, as well as transmembrane ion channels. Results showed that the updated Drude-2019 protein FF yields smaller overall root-mean-square differences of proteins as compared to the additive CHARMM36m and Drude-2013 FFs as well as similar or improved agreement with experimental NMR properties, allowing for long time scale simulation studies of proteins and more complex biomolecular systems in conjunction with the remainder of the Drude polarizable FF.

Impact of electronic polarizability on protein-functional group interactions

Interactions of proteins with functional groups are key to their biological functions, making it essential that they be accurately modeled. To investigate the impact of the inclusion of explicit treatment of electronic polarizability in force fields on protein-functional group interactions, the additive CHARMM and Drude polarizable force field are compared in the context of the Site-Identification by Ligand Competitive Saturation (SILCS) simulation methodology from which functional group interaction patterns with five proteins for which experimental binding affinities of multiple ligands are available, were obtained. The explicit treatment of polarizability produces significant differences in the functional group interactions in the ligand binding sites including overall enhanced binding of functional groups to the proteins. This is associated with variations of the dipole moments of solutes representative of functional groups in the binding sites relative to aqueous solution with higher dipole moments systematically occurring in the latter, though exceptions occur with positively charged methylammonium. Such variation indicates the complex, heterogeneous nature of the electronic environments of ligand binding sites and emphasizes the inherent limitation of fixed charged, additive force fields for modeling ligand-protein interactions. These effects yield more defined orientation of the functional groups in the binding pockets and a small, but systematic improvement in the ability of the SILCS method to predict the binding orientation and relative affinities of ligands to their target proteins. Overall, these results indicate that the physical model associated with the explicit treatment of polarizability along with the presence of lone pairs in a force field leads to changes in the nature of the interactions of functional groups with proteins versus that occurring with additive force fields, suggesting the utility of polarizable force fields in obtaining a more realistic understanding of protein-ligand interactions.

The Ca2+ permeation mechanism of the ryanodine receptor revealed by a multi-site ion model

Ryanodine receptors (RyR) are ion channels responsible for the release of Ca²⁺ from the sarco/endoplasmic reticulum and play a crucial role in the precise control of Ca²⁺ concentration in the cytosol. The detailed permeation mechanism of Ca²⁺ through RyR is still elusive. By using molecular dynamics simulations with a specially designed Ca²⁺ model, we show that multiple Ca²⁺ ions accumulate in the upper selectivity filter of RyR1, but only one Ca²⁺ can occupy and translocate in the narrow pore at a time, assisted by electrostatic repulsion from the Ca²⁺ within the upper selectivity filter. The Ca²⁺ is nearly fully hydrated with the first solvation shell intact during the whole permeation process. These results suggest a remote knock-on permeation mechanism and one-at-a-time occupation pattern for the hydrated Ca²⁺ within the narrow pore, uncovering the basis underlying the high permeability and low selectivity of the RyR channels.

Polarisable force fields: what do they add in biomolecular simulations?

The quality of biomolecular simulations critically depends on the accuracy of the force field used to calculate the potential energy of the molecular configurations. Currently, most simulations employ non-polarisable force fields, which describe electrostatic interactions as the sum of Coulombic interactions between fixed atomic charges. Polarisation of these charge distributions is incorporated only in a mean-field manner. In the past decade, extensive efforts have been devoted to developing simple, efficient, and yet generally applicable polarisable force fields for biomolecular simulations. In this review, we summarise the latest developments in accounting for key biomolecular interactions with polarisable force fields and applications to address challenging biological questions. In the end, we provide an outlook for future development in polarisable force fields.

Toward a Model for Activation of Orai Channel

Store-operated calcium release-activated calcium (CRAC) channels mediate a variety of cellular signaling functions. The CRAC channel pore-forming protein, Orai1, is a hexamer arranged with 3-fold symmetry. Despite its importance in moving Ca²⁺ ions into cells, a detailed mechanistic understanding of Orai1 activation is lacking. Herein, a working model is proposed for the putative open state of Orai from Drosophila melanogaster (dOrai), which involves a “twist-to-open” gating mechanism. The proposed model is supported by energetic, structural, and experimental evidence. Fluorescent imaging demonstrates that each subunit on the intracellular side of the pore is inherently strongly cross-linked, which is important for coupling to STIM1, the pore activator, and graded activation of the Orai1 channel. The proposed model thus paves the way for understanding key aspects of calcium signaling at a molecular level.

•

Mechanical coupling within the calcium channel pore is critical for its activation

•

Molecular modeling could disclose gating mechanism of ion channels at atomic level

•

The predicted open-state structure of the pore was further confirmed by experiments

Interactions of Water and Alkanes: Modifying Additive Force Fields to Account for Polarization Effects

Atomistic biomolecular simulations predominantly utilize additive force fields (FF), where the electrostatic potential is modeled by fixed point charges. Among other consequences, the lack of polarizability in these models undermines the balance of hydrophilic/hydrophobic nonbonded interactions. Simulations of water/alkane systems using the TIP3P water model and CHARMM36 parameters reveal a 1 kcal/mol overestimate of the experimental transfer free energy of water to hexadecane; more recent optimized water models (SPC/E, TIP4P/2005, TIP4P-Ew, TIP3P-FB, TIP4P-FB, OPC, TIP4P-D) overestimate this transfer free energy by approximately 2 kcal/mol. In contrast, the polarizable SWM4-NDP and SWM6 water models reproduce experimental values to within statistical error. As an alternative to explicitly modeling polarizability, this paper develops an efficient automated workflow to optimize pair-specific Lennard-Jones parameters within an additive FF. Water/hexadecane is used as a prototype and the free energy of water transfer to hexadecane as a target. The optimized model yields quantitative agreement with the experimental transfer free energy and improves the water/hexadecane interfacial tension by 6%. Simulations of five different lipid bilayers show a strong increase of water permeabilities compared to the unmodified CHARMM36 lipid FF which consistently improves match with experiment: the order-of-magnitude underestimate for monounsaturated bilayers is rectified and the factor of 2.8-4 underestimate for saturated bilayers is turned into a factor of 1.5-3 overestimate. While agreement with experiment is decreased for the diffusion constant of water in hexadecane, alkane transfer free energies, and the bilayers' area per lipid, the method provides a permeant-specific route to achieve a wide range of heterogeneous observables via rapidly optimized pairwise parameters.

Ion-Hydroxyl Interactions: From High-Level Quantum Benchmarks to Transferable Polarizable Force Fields

Ion descriptors in molecular mechanics models are calibrated against reference data on ion-water interactions. It is then typically assumed that these descriptors will also satisfactorily describe interactions of ions with other functional groups, such as those present in biomolecules. However, several studies now demonstrate that this transferability assumption produces, in many different cases, large errors. Here we address this issue in a representative polarizable model and focus on transferability of cationic interactions from water to a series of alcohols. Both water and alcohols use hydroxyls for ion-coordination, and, therefore, this set of molecules constitutes the simplest possible case of transferability. We obtain gas phase reference data systematically from "gold-standard" quantum Monte Carlo and CCSD(T) methods, followed by benchmarked vdW-corrected DFT. We learn that the original polarizable model yields large gas phase water → alcohol transferability errors - the RMS and maximum errors are 2.3 and 5.1 kcal/mol, respectively. These errors are, nevertheless, systematic in that ion-alcohol interactions are overstabilized, and systematic errors typically imply that some essential physics is either missing or misrepresented. A comprehensive analysis shows that when both low- and high-field responses of ligand dipole polarization are described accurately, then transferability improves significantly - the RMS and maximum errors in the gas phase reduce, respectively, to 0.9 and 2.5 kcal/mol. Additionally, predictions of condensed phase transfer free energies also improve. Nevertheless, within the limits of the extrathermodynamic assumptions necessary to separate experimental estimates of salt dissolution into constituent cationic and anionic contributions, we note that the error in the condensed phase is systematic, which we attribute, at least, partially to the parametrization in long-range electrostatics. Overall, this work demonstrates a rational approach to boosting transferability of ionic interactions that will be applicable broadly to improving other polarizable and nonpolarizable models.

Challenges and advances in atomistic simulations of potassium and sodium ion channel gating and permeation

Ion channels are implicated in many essential physiological events such as electrical signal propagation and cellular communication. The advent of K⁺ and Na⁺ ion channel structure determination has facilitated numerous investigations of molecular determinants of their behaviour. At the same time, rapid development of computer hardware and molecular simulation methodologies has made computational studies of large biological molecules in all-atom representation tractable. The concurrent evolution of experimental structural biology with biomolecular computer modelling has yielded mechanistic details of fundamental processes unavailable through experiments alone, such as ion conduction and ion channel gating. This review is a short survey of the atomistic computational investigations of K⁺ and Na⁺ ion channels, focusing on KcsA and several voltage-gated channels from the K_V and Na_V families, which have garnered many successes and engendered several long-standing controversies regarding the nature of their structure-function relationship. We review the latest advancements and challenges facing the field of molecular modelling and simulation regarding the structural and energetic determinants of ion channel function and their agreement with experimental observations.

Polarizable force field for RNA based on the classical drude oscillator

RNA molecules are highly dynamic and capable of adopting a wide range of complex, folded structures. The factors driving the folding and dynamics of these structures are dependent on a balance of base pairing, hydration, base stacking, ion interactions, and the conformational sampling of the 2'-hydroxyl group in the ribose sugar. The representation of these features is a challenge for empirical force fields used in molecular dynamics simulations. Toward meeting this challenge, the inclusion of explicit electronic polarization is important in accurately modeling RNA structure. In this work, we present a polarizable force field for RNA based on the classical Drude oscillator model, which represents electronic degrees of freedom via negatively charged particles attached to their parent atoms by harmonic springs. Beginning with parametrization against quantum mechanical base stacking interaction energy and conformational energy data, we have extended the Drude-2017 nucleic acid force field to include RNA. The conformational sampling of a range of RNA sequences were used to validate the force field, including canonical A-form RNA duplexes, stem-loops, and complex tertiary folds that bind multiple Mg2+ ions. Overall, the Drude-2017 RNA force field reproduces important properties of these structures, including the conformational sampling of the 2'-hydroxyl and key interactions with Mg2+ ions. © 2018 Wiley Periodicals, Inc.

The Role of Proton Transport in Gating Current in a Voltage Gated Ion Channel, as Shown by Quantum Calculations

Over two-thirds of a century ago, Hodgkin and Huxley proposed the existence of voltage gated ion channels (VGICs) to carry Na⁺ and K⁺ ions across the cell membrane to create the nerve impulse, in response to depolarization of the membrane. The channels have multiple physiological roles, and play a central role in a wide variety of diseases when they malfunction. The first channel structure was found by MacKinnon and coworkers in 1998. Subsequently, the structure of a number of VGICs was determined in the open (ion conducting) state. This type of channel consists of four voltage sensing domains (VSDs), each formed from four transmembrane (TM) segments, plus a pore domain through which ions move. Understanding the gating mechanism (how the channel opens and closes) requires structures. One TM segment (S4) has an arginine in every third position, with one such segment per domain. It is usually assumed that these arginines are all ionized, and in the resting state are held toward the intracellular side of the membrane by voltage across the membrane. They are assumed to move outward (extracellular direction) when released by depolarization of this voltage, producing a capacitive gating current and opening the channel. We suggest alternate interpretations of the evidence that led to these models. Measured gating current is the total charge displacement of all atoms in the VSD; we propose that the prime, but not sole, contributor is proton motion, not displacement of the charges on the arginines of S4. It is known that the VSD can conduct protons. Quantum calculations on the K_v1.2 potassium channel VSD show how; the key is the amphoteric nature of the arginine side chain, which allows it to transfer a proton. This appears to be the first time the arginine side chain has had its amphoteric character considered. We have calculated one such proton transfer in detail: this proton starts from a tyrosine that can ionize, transferring to the NE of the third arginine on S4; that arginine’s NH then transfers a proton to a glutamate. The backbone remains static. A mutation predicted to affect the proton transfer has been qualitatively confirmed experimentally, from the change in the gating current-voltage curve. The total charge displacement in going from a normal closed potential of −70 mV across the membrane to 0 mV (open), is calculated to be approximately consistent with measured values, although the error limits on the calculation require caution in interpretation.

Discrete and Continuum Models for the Salt in Crowded Environments of Suspended Charged Particles

Electrostatic forces greatly affect the overall dynamics and diffusional activities of suspended charged particles in crowded environments. Accordingly, the concentration of counter- or co-ions in a fluid-''the salt"-determines the range, strength, and order of electrostatic interactions between particles. This environment fosters engineering routes for controlling directed assembly of particles at both the micro- and nanoscale. Here, we analyzed two computational modeling schemes that considered salt within suspensions of charged particles, or polyelectrolytes: discrete and continuum. Electrostatic interactions were included through a Green's function formalism, where the confined fundamental solution for Poisson's equation is resolved by the general geometry Ewald-like method. For the discrete model, the salt was considered as regularized point-charges with a specific valence and size, while concentration fields were defined for each ionic species for the continuum model. These considerations were evolved using Brownian dynamics of the suspended charged particles and the discrete salt ions, while a convection-diffusion transport equation, including the Nernst-Planck diffusion mechanism, accounted for the dynamics of the concentration fields. The salt/particle models were considered as suspensions under slit-confinement conditions for creating crowded "macro-ions", where density distributions and radial distribution functions were used to compare and differentiate computational models. Importantly, our analysis shows that disparate length scales or increased system size presented by the salt and suspended particles are best dealt with using concentration fields to model the ions. These findings were then validated by novel simulations of a semipermeable polyelectrolyte membrane, at the mesoscale, from which ionic channels emerged and enable ion conduction.