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

Shadow Molecular Dynamics for a Charge-Potential Equilibration Model

We introduce a shadow molecular dynamics (MD) approach based on the Atom-Condensed Kohn-Sham second-order (ACKS2) charge-potential equilibration model. In contrast to regular flexible charge models, the ACKS2 model includes both flexible atomic charges and potential fluctuation parameters that allow for physically correct charge fragmentation and improved scaling of the polarizability. Our shadow MD scheme is based on an approximation of the ACKS2's flexible charge-potential energy function, in combination with extended Lagrangian Born-Oppenheimer MD. Utilizing this shadow charge-potential equilibration approach mitigates the costly overhead and stability problems associated with finding well-converged iterative solutions to the charges and potential fluctuations of the ACKS2 model in an MD simulation. Our work provides a robust and versatile framework for efficient, high-fidelity MD simulations of diverse physical phenomena and applications.

Shadow Molecular Dynamics with a Machine Learned Flexible Charge Potential

We present an extended Lagrangian shadow molecular dynamics scheme with an interatomic Born-Oppenheimer potential determined by the relaxed atomic charges of a second-order charge equilibration model. To parametrize the charge equilibration model, we use machine learning with neural networks to determine the environment-dependent electronegativities and chemical hardness parameters for each atom, in addition to the charge-independent energy and force terms. The approximate shadow molecular dynamics potential in combination with the extended Lagrangian formulation improves the numerical stability and reduces the number of Coulomb potential calculations required to evaluate accurate conservative forces. We demonstrate efficient and accurate simulations with excellent long-term stability of the molecular dynamics trajectories. The significance of choosing fixed or environment-dependent electronegativities and chemical hardness parameters is evaluated. Finally, we compute the infrared spectrum of molecules via the dipole autocorrelation function and compare to experiments to highlight the accuracy of the shadow molecular dynamics scheme with a machine learned flexible charge potential.

Electronic Polarizability Tunes the Function of the Human Bestrophin 1 Cl- Channel

Mechanisms of anion permeation within ion channels and nanopores remain poorly understood. Recent cryo-electron microscopy structures of the human bestrophin 1 Cl- channel (hBest1) provide an opportunity to evaluate ion interactions predicted by molecular dynamics (MD) simulations against experimental observations. Here, we implement the fully polarizable force field AMOEBA in MD simulations on different conformations of hBest1. This force field models multipole moments up to the quadrupole. Using this approach, we model key biophysical properties of the channel that can only be simulated when electronic polarization is included in the molecular models and show that Cl- permeation through the neck of the pore is achieved through hydrophobic solvation concomitant with partial ion dehydration. Furthermore, we demonstrate how such polarizable simulations can help determine the identity of ion-like densities within high-resolution cryo-EM structures and demonstrate that neglecting polarization places Cl- at positions that do not correspond to their experimentally resolved location. Overall, our results demonstrate the importance of including electronic polarization in realistic and physically accurate models of biological systems, especially channels and pores that selectively permeate anions.

Grand canonical Monte Carlo and deep learning assisted enhanced sampling to characterize the distribution of Mg2+ and influence of the Drude polarizable force field on the stability of folded states of the twister ribozyme

Molecular dynamics simulations are crucial for understanding the structural and dynamical behavior of biomolecular systems, including the impact of their environment. However, there is a gap between the time scale of these simulations and that of real-world experiments. To address this problem, various enhanced simulation methods have been developed. In addition, there has been a significant advancement of the force fields used for simulations associated with the explicit treatment of electronic polarizability. In this study, we apply oscillating chemical potential grand canonical Monte Carlo and machine learning methods to determine reaction coordinates combined with metadynamics simulations to explore the role of Mg2+ distribution and electronic polarizability in the context of the classical Drude oscillator polarizable force field on the stability of the twister ribozyme. The introduction of electronic polarizability along with the details of the distribution of Mg2+ significantly stabilizes the simulations with respect to sampling the crystallographic conformation. The introduction of electronic polarizability leads to increased stability over that obtained with the additive CHARMM36 FF reported in a previous study, allowing for a distribution of a wider range of ions to stabilize twister. Specific interactions contributing to stabilization are identified, including both those observed in the crystal structures and additional experimentally unobserved interactions. Interactions of Mg2+ with the bases are indicated to make important contributions to stabilization. Notably, the presence of specific interactions between the Mg2+ ions and bases or the non-bridging phosphate oxygens (NBPOs) leads to enhanced dipole moments of all three moieties. Mg2+-NBPO interactions led to enhanced dipoles of the phosphates but, interestingly, not in all the participating ions. The present results further indicate the importance of electronic polarizability in stabilizing RNA in molecular simulations and the complicated nature of the relationship of Mg2+-RNA interactions with the polarization response of the bases and phosphates.

Parameterization of a Fluctuating Charge Model for Complexes Containing 3d Transition Metals

Metalloproteins widely exist in biology, playing pivotal roles in diverse life processes. Meanwhile, molecular dynamics (MD) simulations based on classical force fields has emerged as an important tool in scientific research. Partial charges are critical parameters within classical force fields and usually derived from quantum mechanical (QM) calculations. However, QM calculations are often time-consuming and prone to basis set dependence. Alternatively, fluctuating charge (FQ) models offer another avenue for partial charge derivation, which has significant speed advantages and can be used for large-scale screening. Building upon our previous work, which introduced an FQ model for zinc-containing complexes, herein we extend this model to include additional 3d transition metals which are important to the life sciences, namely chromium, manganese, iron, cobalt, and nickel. Employing CM5 charges as target for parametrization, our FQ model accurately reproduces partial charges for 3d metal complexes featuring biologically relevant ligands. Furthermore, by using atomic charges derived by our FQ model, MD simulations have been performed. These charges showed excellent performance in simulating proteomic metal sites housing multiple metal ions, specifically, a metalloprotein containing an iron-sulfur cluster and another containing a dimanganese metal site, showcasing comparable performance to those of RESP charges. We anticipate that our study can accelerate the parametrization of atomic charges for metalloproteins featuring 3d transition metals, thereby facilitating simulations of relevant systems.

Getting the intermolecular forces correct: introducing the ASTA strategy for a water model

Having a force field for water providing good bulk properties is paramount for modern studies of most biological systems. Some of the most common three-site force fields are TIP3, SPC/ε or OPC3, providing a decent range of bulk properties. That does not mean though, that they have realistic inter-atomic forces. These force fields have been parameterized with a top-down approach, meaning, by fitting the force field parameters to the experimental bulk properties. This approach has been the governing strategy also for many variants of four- and more-site models. We test a bottom-up approach, in which the force field is parameterized by optimizing the non-bonded inter-atomic forces. Our philosophy is that correct inter-atomic forces lead to correct geometrical and dynamical properties. The first system we try to optimize with the accurately system tailored atomic (ASTA) approach is water, but we aim to eventually probe other systems in the future as well. We applied our ASTA strategy to find a good set of parameters providing accurate bulk properties for the simple three-site force field forms, and also for AMOEBA, a more detailed and polarizable force field. Even though our bottom-up approach did not provide satisfactory results for the simple three-site force fields (with fixed charges), for the case of the AMOEBA force field it led to a modification of the original strategy, giving very good intra- and inter-molecular forces, as compared to accurate quantum chemically calculated reference forces. At the same time, important bulk properties, in this study restricted to the density and diffusion, were accurately reproduced with respect to the experimental values.

Electronic Polarizability Tunes the Function of the Human Bestrophin 1 Cl- Channel

Mechanisms of anion permeation within ion channels and nanopores remain poorly understood. Recent cryo-electron microscopy structures of the human bestrophin 1 Cl- channel (hBest1) provide an opportunity to evaluate ion interactions predicted by molecular dynamics (MD) simulations against experimental observations. Here, we implement the fully polarizable forcefield AMOEBA in MD simulations on different conformations of hBest1. This forcefield models multipole moments up to the quadrupole; therefore, it captures induced dipole and anion-π interactions. We show that key biophysical properties of the channel can only be simulated when electronic polarization is included in the molecular models and that Cl- permeation through the neck of the pore is achieved through hydrophobic solvation concomitant with partial ion dehydration. Furthermore, we demonstrate how such polarizable simulations can help determine the identity of ion-like densities within high-resolution cryo-EM structures and that neglecting polarization places Cl- at positions that do not correspond with their experimentally resolved location. Overall, our results demonstrate the importance of including electronic polarization in realistic and physically accurate models of biological systems, especially channels and pores that selectively permeate anions.

Condensed Matter Systems Exposed to Radiation: Multiscale Theory, Simulations, and Experiment

This roadmap reviews the new, highly interdisciplinary research field studying the behavior of condensed matter systems exposed to radiation. The Review highlights several recent advances in the field and provides a roadmap for the development of the field over the next decade. Condensed matter systems exposed to radiation can be inorganic, organic, or biological, finite or infinite, composed of different molecular species or materials, exist in different phases, and operate under different thermodynamic conditions. Many of the key phenomena related to the behavior of irradiated systems are very similar and can be understood based on the same fundamental theoretical principles and computational approaches. The multiscale nature of such phenomena requires the quantitative description of the radiation-induced effects occurring at different spatial and temporal scales, ranging from the atomic to the macroscopic, and the interlinks between such descriptions. The multiscale nature of the effects and the similarity of their manifestation in systems of different origins necessarily bring together different disciplines, such as physics, chemistry, biology, materials science, nanoscience, and biomedical research, demonstrating the numerous interlinks and commonalities between them. This research field is highly relevant to many novel and emerging technologies and medical applications.

Efficient Sampling of Cavity Hydration in Proteins with Nonequilibrium Grand Canonical Monte Carlo and Polarizable Force Fields

Prediction of the hydration levels of protein cavities and active sites is important to both mechanistic analysis and ligand design. Due to the unique microscopic environment of these buried water molecules, a polarizable model is expected to be crucial for an accurate treatment of protein internal hydration in simulations. Here we adapt a nonequilibrium candidate Monte Carlo approach for conducting grand canonical Monte Carlo simulations with the Drude polarizable force field. The GPU implementation enables the efficient sampling of internal cavity hydration levels in biomolecular systems. We also develop an enhanced sampling approach referred to as B-walking, which satisfies detailed balance and readily combines with grand canonical integration to efficiently calculate quantitative binding free energies of water to protein cavities. Applications of these developments are illustrated in a solvent box and the polar ligand binding site in trypsin. Our simulation results show that including electronic polarization leads to a modest but clear improvement in the description of water position and occupancy compared to the crystal structure. The B-walking approach enhances the range of water sampling in different chemical potential windows and thus improves the accuracy of water binding free energy calculations.

Development of a Fluctuating Charge Model for Zinc-Containing Metalloproteins

Metalloproteins widely exist in biology and play important roles in various processes. To accurately simulate metalloprotein systems, modeling polarization and charge transfer effects is vital. The fluctuating charge (FQ) model can efficiently generate atomic charges and simulate the charge transfer effect; it has been developed for a wide range of applications, but few models have been specifically tailored for metalloproteins. In this study, we present a fluctuating charge model specifically for zinc-containing metalloproteins based on the extended charge equilibration (EQeq) scheme. Our model was parametrized to reproduce CM5 charges instead of RESP/CHELPG charges because the former is less dependent on the conformation or basis set, does not suffer from unphysical charges for buried atoms, and is still being able to well reproduce the molecular dipoles. During our study, we found that adding the Pauling-bond-order-like term (referred to as the "+C term" in a previous study) between the zinc ion and ligating atoms significantly improves the model's performance. Although our model was trained for four-coordinated zinc sites only, our results indicated it can well describe the atomic charges in diverse zinc sites. Morever, our model was used to generate partial charges for the metal sites in three different zinc-containing metalloproteins (with four-, five-, and six-coordinated metal sites, respectively). These charges exhibited performance comparable to that of the RESP charges in molecular dynamics (MD) simulations. Additional tests indicated our model could also well reproduce the CM5 charges when geometric changes were involved. Those results indicate that our model can efficiently calculate the atomic charges for metal sites and well simulate the charge transfer effect, which marks an important step toward developing versatile polarizable force fields for metalloproteins.

The effect of monomer polarizability on the stability and salt partitioning in model coacervates

Coacervation of charged polymer chains has been a topic of major interest in both polymer and biological sciences, as it is a subset of a phenomenon called liquid-liquid phase separation (LLPS). In this process the polymer-rich phase separates from the polymer-lean supernatant while still maintaining its liquid-like properties. LLPS has been shown to play a crucial role in cellular homeostasis by driving the formation of membraneless organelles. It also has the potential to be harnessed to aid in novel therapeutical applications. Recent studies have demonstrated that there is no one simple mechanism which drives LLPS, which is instead a result of the combined effect of electrostatic, dipolar, hydrophobic, and other weak interactions. Using coarse-grained polymer simulations we investigate the relatively unexplored effects of monomer polarizability and spatially varying dielectric constant on LLPS propensity, and these factors affect the properties of the resulting condensates. In order to produce spatial variations in the dielectric constant, all our simulations include explicit solvent and counterions. We demonstrate that polarizability has only a minor effect on the bulk behaviour of the condensates but plays a major role when ion partitioning and microstructure are considered. We observe that the major contribution comes from the nature of the neutral blocks as endowing them with an induced dipole changes their character from hydrophobic to hydrophilic. We hypothesize that the results of this work can aid in guiding future studies concerned with LLPS by providing a general framework and by highlighting important factors which influence LLPS.

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.

Influence of electronic polarization on the binding of anions to a chloride-pumping rhodopsin

The functional properties of some biological ion channels and membrane transport proteins are proposed to exploit anion-hydrophobic interactions. Here, we investigate a chloride-pumping rhodopsin as an example of a membrane protein known to contain a defined anion binding site composed predominantly of hydrophobic residues. Using molecular dynamics simulations, we explore Cl- binding to this hydrophobic site and compare the dynamics arising when electronic polarization is neglected (CHARMM36 [c36] fixed-charge force field), included implicitly (via the prosECCo force field), or included explicitly (through the polarizable force field, AMOEBA). Free energy landscapes of Cl- moving out of the binding site and into bulk solution demonstrate that the inclusion of polarization results in stronger ion binding and a second metastable binding site in chloride-pumping rhodopsin. Simulations focused on this hydrophobic binding site also indicate longer binding durations and closer ion proximity when polarization is included. Furthermore, simulations reveal that Cl- within this binding site interacts with an adjacent loop to facilitate rebinding events that are not observed when polarization is neglected. These results demonstrate how the inclusion of polarization can influence the behavior of anions within protein binding sites and can yield results comparable with more accurate and computationally demanding methods.

Quantum chemical studies on hydrogen bonds in helical secondary structures

We present a brief review of our recent computational studies of hydrogen bonds (H-bonds) in helical secondary structures of proteins, α-helix and 3₁₀-helix, using a Negative Fragmentation Approach with density functional theory. We found that the depolarized electronic structures of the carbonyl oxygen of the ith residue and the amide hydrogen of the (i + 4)th residue cause weaker H-bond in an α-helix than in an isolated H-bond. Our calculations showed that the H-bond energies in the 3₁₀-helix were also weaker than those of the isolated H-bonds. In the 3₁₀-helices, the adjacent N-H group at the (i + 1)th residue was closer to the C=O group of the H-bond pair than the adjacent C=O group in the 3₁₀-helices, whereas the adjacent C=O group at the (i + 1)th residue was close to the H-bond acceptor in α-helices. Therefore, the destabilization of the H-bond is attributed to the depolarization caused by the adjacent residue of the helical backbone connecting the H-bond donor and acceptor. The differences in the change in electron density revealed that such depolarizations were caused by the local electronic interactions in their neighborhood inside the helical structure and redistributed the electron density. We also present the improvements in the force field of classical molecular simulation, based on our findings.

Supplementary Information

The online version contains supplementary material available at 10.1007/s12551-022-01034-5.

Transferable Potential Function for Flexible H2O Molecules Based on the Single-Center Multipole Expansion

A potential function is presented for describing a system of flexible H₂O molecules based on the single-center multipole expansion (SCME) of the electrostatic interaction. The model, referred to as SCME/f, includes the variation of the molecular quadrupole moment as well as the dipole moment with changes in bond length and angle so as to reproduce results of high-level electronic structure calculations. The multipole expansion also includes fixed octupole and hexadecapole moments, as well as anisotropic dipole-dipole, dipole-quadrupole, and quadrupole-quadrupole polarizability tensors. The model contains five adjustable parameters related to the repulsive interaction and damping functions in the electrostatic and dispersion interactions. Their values are adjusted to reproduce the lowest energy isomers of small clusters, (H₂O)_n with n = 2-6, as well as measured properties of the ice Ih crystal. Subsequent calculations of the energy difference between the various isomer configurations of the clusters show that SCME/f gives good agreement with results of electronic structure calculations and represents a significant improvement over the previously presented rigid SCME potential function. Analysis of the vibrational frequencies of the clusters and structural properties of ice Ih crystal show the importance of accurately describing the variation of the quadrupole moment with molecular structures.

Fast Polarizable Water Model for Atomistic Simulations

Simulating water accurately has been a major challenge in atomistic simulations for decades. Inclusion of electronic polarizability effects holds considerable promise, yet existing approaches suffer from significant computational overheads compared to the widely used nonpolarizable water models. We have developed a globally optimal polarizable water model, OPC3-pol, that explicitly accounts for electronic polarizability with minimal impact on the computational efficiency. OPC3-pol reproduces five key bulk water properties at room temperature with an average relative error of 0.6%. In atomistic simulations, OPC3-pol's computational efficiency is in between that of 3- and 4-point nonpolarizable models; the model supports increased (4 fs) integration time step. OPC3-pol is tested in simulations of globular protein ubiquitin and a B-DNA dodecamer with several AMBER force fields, ff99SB, ff14SB, ff19SB, and OL15, demonstrating structure stability close to reference on multi-microsecond time scale. Simulation of an intrinsically disordered amyloid β-peptide yields an ensemble with the radius of gyration of a random coil. The proposed water model can be trivially adopted by any package that supports standard nonpolarizable force fields and water models; its intended use is in long classical atomistic simulations where water polarization effects are expected to be important.

Classical Exchange Polarization: An Anisotropic Variable Polarizability Model

Polarizability, or the tendency of the electron distribution to distort under an electric field, often depends on the local chemical environment. For example, the polarizability of a chloride ion is larger in gas phase compared to a chloride ion solvated in water. This effect is due to the restriction the Pauli exclusion principle places on the allowed electron states. Because no two electrons can occupy the same state, when a highly polarizable atom comes in close contact with other atoms or molecules, the space of allowed states can dramatically decrease. This constraint suggests that an accurate molecular mechanics polarizability model should depend on the radial distance between neighboring atoms. This paper introduces a variable polarizability model within the framework of the HIPPO (Hydrogen-like Intermolecular Polarizable Potential) force field, by damping the polarizability as a function of the orbital overlap of two atoms. This effectively captures the quantum mechanical exchange polarization effects, without explicit utilization of antisymmetrized wave functions. We show that the variable polarizability model remarkably improves the two-body polarization energies and three-body energies of ion-ion and ion-water systems. Under this model, no manual tuning of atomic polarizabilities for monatomic ions is required; the gas-phase polarizability can be used because an appropriate damping function is able to correct the polarizability at short range.

Extension of the CHARMM Classical Drude Polarizable Force Field to N- and O-Linked Glycopeptides and Glycoproteins

Molecular dynamic simulations are an effective tool to study complex molecular systems and are contingent upon the availability of an accurate and reliable molecular mechanics force field. The Drude polarizable force field, which allows for the explicit treatment of electronic polarization in a computationally efficient fashion, has been shown to reproduce experimental properties that were difficult or impossible to reproduce with the CHARMM additive force field, including peptide folding cooperativity, RNA hairpin structures, and DNA base flipping. Glycoproteins are essential components of glycoconjugate vaccines, antibodies, and many pharmaceutically important molecules, and an accurate polarizable force field that includes compatibility between the protein and carbohydrate aspect of the force field is essential to study these types of systems. In this work, we present an extension of the Drude polarizable force field to glycoproteins, including both N- and O-linked species. Parameter optimization focused on the dihedral terms using a reweighting protocol targeting NMR solution J-coupling data for model glycopeptides. Validation of the model include eight model glycopeptides and four glycoproteins with multiple N- and O-linked glycosylations. The new glycoprotein carbohydrate force field can be used in conjunction with the remainder of Drude polarizable force field through a variety of MD simulation programs including GROMACS, OPENMM, NAMD, and CHARMM and may be accessed through the Drude Prepper module in the CHARMM-GUI.

Depolarizing Effects in Hydrogen Bond Energy in 310-Helices Revealed by Quantum Chemical Analysis

Hydrogen-bond (H-bond) energies in 3₁₀-helices of short alanine peptides were systematically examined by precise DFT calculations with the negative fragmentation approach (NFA), a modified method based on the molecular tailoring approach. The contribution of each H-bond was evaluated in detail from the 3₁₀-helical conformation of total energies (whole helical model, WH_3-10 model), and the results were compared with the property of H-bond in α-helix from our previous study. The H-bond energies of the WH_3-10 model exhibited tendencies different from those exhibited by the α-helix in that they depended on the helical position of the relevant H-bond pair. H-bond pairs adjacent to the terminal H-bond pairs were observed to be strongly destabilized. The analysis of electronic structures indicated that structural characteristics cause the destabilization of the H-bond in 3₁₀-helices. We also found that the longer the helix length, the more stable the H-bond in the terminal pairs of the WH_3-10 model, suggesting the action of H-bond cooperativity.

Quantitative molecular simulations

All-atom simulations can provide molecular-level insights into the dynamics of gas-phase, condensed-phase and surface processes. One important requirement is a sufficiently realistic and detailed description of the underlying intermolecular interactions. The present perspective provides an overview of the present status of quantitative atomistic simulations from colleagues' and our own efforts for gas- and solution-phase processes and for the dynamics on surfaces. Particular attention is paid to direct comparison with experiment. An outlook discusses present challenges and future extensions to bring such dynamics simulations even closer to reality.

Electronic Polarization Is Essential for the Stabilization and Dynamics of Buried Ion Pairs in Staphylococcal Nuclease Mutants

Buried charged residues play important roles in the modulation of protein stabilities and conformational dynamics and make crucial contributions to protein functions. Considering the generally nonpolar nature of protein interior, a key question concerns the contribution of electronic polarization to the stabilization and properties of buried charges. We answer this question by conducting free energy simulations using the latest polarizable CHARMM force field based on Drude oscillators for a series of Staphylococcal nuclease mutants that involve a buried Glu-Lys pair in different titration states and orientations. While a nonpolarizable model suggests that the ionized form of the buried Glu-Lys pair is more than 40 kcal/mol less stable than the charge-neutral form, the two titration states are comparable in stability when electronic polarization is included explicitly, a result better reconcilable with available experimental data. Analysis of free energy components suggests that additional stabilization of the ionized Glu-Lys pair has contributions from both the enhanced salt-bridge strength and stronger interaction between the ion-pair and surrounding protein residues and penetrated water. Despite the stronger direct interaction between Glu and Lys, the ion-pair exhibits considerably larger and faster structural fluctuations when polarization is included, due to compensation of interactions in the cavity. Collectively, observations from this work provide compelling evidence that electronic polarization is essential to the stability, hydration, dynamics, and therefore function of buried charges in proteins. Therefore, our study advocates for the explicit consideration of electronic polarization for mechanistic and engineering studies that implicate buried charged residues, such as enzymes and ion transporters.

Modelling of the dynamic polarizability of macromolecules for single-molecule optical biosensing

The structural dynamics of macromolecules is important for most microbiological processes, from protein folding to the origins of neurodegenerative disorders. Noninvasive measurements of these dynamics are highly challenging. Recently, optical sensors have been shown to allow noninvasive time-resolved measurements of the dynamic polarizability of single-molecules. Here we introduce a method to efficiently predict the dynamic polarizability from the atomic configuration of a given macromolecule. This provides a means to connect the measured dynamic polarizability to the underlying structure of the molecule, and therefore to connect temporal measurements to structural dynamics. To illustrate the methodology we calculate the change in polarizability as a function of time based on conformations extracted from molecular dynamics simulations and using different conformations of motor proteins solved crystalographically. This allows us to quantify the magnitude of the changes in polarizablity due to thermal and functional motions.

Using atomic charges to model molecular polarization

We review different models for introducing electric polarization in force fields, with special focus on methods where polarization is modelled at the atomic charge level. While electric polarization has been included in several force fields, the common approach has been to focus on atomic dipole polarizability. Several approaches allow modelling electric polarization by using charge-flow between charge sites instead, but this has been less exploited, despite that atomic charges and charge-flow is expected to be more important than atomic dipoles and dipole polarizability. A number of challenges are required to be solved for charge-flow models to be incorporated into polarizable force fields, for example how to parameterize the models and how to make them computational efficient.

Singular value decomposition analysis of the electron density changes occurring upon electrostatic polarization of water

In-depth elucidation of how molecules are electrically polarized would be one key factor for understanding the properties of those molecules under various thermodynamic and/or spatial conditions. Here this problem is tackled for the case of hydrogen-bonded water by conducting singular value decomposition of the electron density changes that occur upon electrostatic polarization. It is shown that all those electron density changes are approximately described as linear combinations of ten orthonormal basis “vectors”. One main component is the interatomic charge transfer through each OH bond, while some others are characterized as the atomic dipolar polarizations, meaning that both of these components are important for the electrostatic polarization of water. The interaction parameters that reasonably well reproduce the induced dipole moments are derived, which indicate the extent of mixing of the two components in electrostatic polarization.

The main features of the electron density changes that occur upon electrostatic polarization of water are elucidated by conducting singular value decomposition analysis of those changes.

Does Explicit Polarizability Improve Simulations of Phase Behavior of Ionic Liquids?

Molecular dynamics simulations are performed for a test set of 20 aprotic ionic liquids to investigate whether including an explicit polarizability model in the force field leads to higher accuracy and reliability of the calculated phase behavior properties, especially the enthalpy of fusion. A classical nonpolarizable all-atom optimized potentials for liquid simulations (OPLS) force-field model developed by Canongia Lopes and Pádua (CL&P) serves as a reference level of theory. Polarizability is included either in the form of Drude oscillators, resulting in the CL&P-D models, or in the framework of the atomic multipole optimized energetics for biomolecular application (AMOEBA) force field with polarizable atomic sites. Benchmarking of the calculated fusion enthalpy values against the experimental data reveals that overall the nonpolarizable CL&P model and polarizable CL&P-D models perform similarly with average deviations of about 30%. However, fusion enthalpies from the CL&P-D models exhibit a stronger correlation with their experimental counterparts. The least successful predictions are interestingly obtained from AMOEBA (deviation ca. 60%), which may indicate that a reparametrization of this force-field model is needed to achieve improved predictions of the fusion enthalpy. In general, all FF models tend to underestimate the fusion enthalpies. In addition, quantum chemical calculations are used to compute the electronic cohesive energies of the crystalline phases of the ionic liquids and of the interaction energies within the ion pair. Significant positive correlations are found between the fusion enthalpy and the cohesive energies. The character of the present anions predetermines the magnitude of individual mechanistic components of the interaction energy and related enthalpic and cohesive properties.

Development and application of quantum mechanics/molecular mechanics methods with advanced polarizable potentials

Quantum mechanics/molecular mechanics (QM/MM) simulations are a popular approach to study various features of large systems. A common application of QM/MM calculations is in the investigation of reaction mechanisms in condensed-phase and biological systems. The combination of QM and MM methods to represent a system gives rise to several challenges that need to be addressed. The increase in computational speed has allowed the expanded use of more complicated and accurate methods for both QM and MM simulations. Here, we review some approaches that address several common challenges encountered in QM/MM simulations with advanced polarizable potentials, from methods to account for boundary across covalent bonds and long-range effects, to polarization and advanced embedding potentials.

How Hierarchical Interactions Make Membraneless Organelles Tick Like Clockwork

Biomolecular condensates appear throughout the cell, serving many different biochemical functions. We argue that condensate functionality is optimized when the interactions driving condensation vary widely in affinity. Strong interactions provide structural specificity needed to encode functional properties but carry the risk of kinetic arrest, while weak interactions allow the system to remain dynamic but do not restrict the conformational ensemble enough to sustain specific functional features. To support our opinion, we describe illustrative examples of the interplay of strong and weak interactions that are found in the nucleolus, SPOP/DAXX condensates, polySUMO/polySIM condensates, chromatin, and stress granules. The common feature of these systems is a hierarchical assembly motif in which weak, transient interactions condense structurally defined functional units.

Quantum Simulations of Hydrogen Bonding Effects in Glycerol Carbonate Electrolyte Solutions

The need for environmentally friendly nonaqueous solvents in electrochemistry and other fields has motivated recent research into the molecular-level solvation structure, thermodynamics, and dynamics of candidate organic liquids. In this paper, we present the results of quantum density functional theory simulations of glycerol carbonate (GC), a molecule that has been proposed as a solvent for green industrial chemistry, nonaqueous alternatives for biocatalytic reactions, and liquid media in energy storage devices. We investigate the structure and dynamics of both the pure GC liquid and electrolyte solutions containing KF and KCl ion pairs. These simulations reveal the importance of hydrogen bonding that controls the structural and dynamic behavior of the pure liquid and ion association in the electrolyte solutions. The results illustrate the difficulties associated with classical modeling of complex organic solvents. The simulations lead to a better understanding of the underlying mechanisms behind the previously observed peculiar ion-specific behavior in GC electrolyte solutions.

Extended Lagrangian Born-Oppenheimer molecular dynamics for orbital-free density-functional theory and polarizable charge equilibration models

Extended Lagrangian Born-Oppenheimer molecular dynamics (XL-BOMD) [A. M. N. Niklasson, Phys. Rev. Lett. 100, 123004 (2008)] is formulated for orbital-free Hohenberg-Kohn density-functional theory and for charge equilibration and polarizable force-field models that can be derived from the same orbital-free framework. The purpose is to introduce the most recent features of orbital-based XL-BOMD to molecular dynamics simulations based on charge equilibration and polarizable force-field models. These features include a metric tensor generalization of the extended harmonic potential, preconditioners, and the ability to use only a single Coulomb summation to determine the fully equilibrated charges and the interatomic forces in each time step for the shadow Born-Oppenheimer potential energy surface. The orbital-free formulation has a charge-dependent, short-range energy term that is separate from long-range Coulomb interactions. This enables local parameterizations of the short-range energy term, while the long-range electrostatic interactions can be treated separately. The theory is illustrated for molecular dynamics simulations of an atomistic system described by a charge equilibration model with periodic boundary conditions. The system of linear equations that determines the equilibrated charges and the forces is diagonal, and only a single Ewald summation is needed in each time step. The simulations exhibit the same features in accuracy, convergence, and stability as are expected from orbital-based XL-BOMD.

Molecular Dynamics Simulations Based on Polarizable Models Show that Ion Permeation Interconverts between Different Mechanisms as a Function of Membrane Thickness

Different mechanisms have been proposed to explain the permeation of charged compounds through lipid membranes. Overall, it is expected that an ion-induced defect permeation mechanism, where substantial membrane deformations accompany ion movement, should be dominant in thin membranes but that a solubility-diffusion mechanism, where ions partition into the membrane core with large associated dehydration energy costs, becomes dominant in thicker membranes. However, while this physical picture is intuitively reasonable, capturing the interconversion between these two permeation mechanisms in molecular dynamics (MD) simulations based on atomic models is challenging. In particular, simulations relying on nonpolarizable force fields are artificially unfavorable to the solubility-diffusion mechanism, as induced polarization of the nonpolar hydrocarbon is ignored, causing overestimated free energy costs for charged molecules to enter into this region of the membrane. In this study, all-atom MD simulations based on nonpolarizable and polarizable force fields are used to quantitatively characterize the permeation process for the arginine side chain analog methyl-guanidinium through bilayer membranes of mono-unsaturated phosphatidylcholine lipids with and without cholesterol, resulting in thicknesses spanning from ∼24 to ∼42 Å. With simulations based on a nonpolarizable force field, ion translocation can take place solely through an ion-induced defect mechanism, with free energy barriers increasing linearly from 14 to 40 kcal/mol, depending on the thickness. However, with simulations based on a polarizable force field, ion translocation is predominantly dominated by an ion-induced defect mechanism in thin membranes, which progressively converts to a solubility-diffusion mechanism as the membranes get thicker. The transition between the two mechanisms occurs at a thickness of ∼29 Å, with lipid tails of 22 or more carbon atoms. This situation appears to represent the upper limit for ion-induced defect permeation within the current polarizable models. Beyond this thickness, it becomes energetically preferable for the ion to dehydrate and partition into the membrane core-a phenomenon that cannot be captured using the nonpolarizable models. Induced electronic polarizability therefore leads not just to a shift in permeation energetics but to an interconversion between two strikingly different physical mechanisms. The result highlights the importance of induced polarizability in modeling lipid membranes.

Statistical mechanics of polarizable force fields based on classical Drude oscillators with dynamical propagation by the dual-thermostat extended Lagrangian

Polarizable force fields based on classical Drude oscillators offer a practical and computationally efficient avenue to carry out molecular dynamics (MD) simulations of large biomolecular systems. To treat the polarizable electronic degrees of freedom, the Drude model introduces a virtual charged particle that is attached to its parent nucleus via a harmonic spring. Traditionally, the need to relax the electronic degrees of freedom for each fixed set of nuclear coordinates is achieved by performing an iterative self-consistent field (SCF) calculation to satisfy a selected tolerance. This is a computationally demanding procedure that can increase the computational cost of MD simulations by nearly one order of magnitude. To avoid the costly SCF procedure, a small mass is assigned to the Drude particles, which are then propagated as dynamic variables during the simulations via a dual-thermostat extended Lagrangian algorithm. To help clarify the significance of the dual-thermostat extended Lagrangian propagation in the context of the polarizable force field based on classical Drude oscillators, the statistical mechanics of a dual-temperature canonical ensemble is formulated. The conditions for dynamically maintaining the dual-temperature properties in the case of the classical Drude oscillator are analyzed using the generalized Langevin equation.

Balanced polarizable Drude force field parameters for molecular anions: phosphates, sulfates, sulfamates, and oxides

Polarizable force fields are emerging as a more accurate alternative to additive force fields in terms of modeling and simulations of a variety of chemicals including biomolecules. Explicit treatment of induced polarization in charged species such as phosphates and sulfates offers the potential for achieving an improved atomistic understanding of the physical forces driving their interactions with their environments. To help achieve this, in this study we present balanced Drude polarizable force field parameters for molecular ions including phosphates, sulfates, sulfamates, and oxides. Better balance was primarily achieved in the relative values of minimum interaction energies and distances of the anionic model compounds with water at the Drude and quantum mechanical (QM) model chemistries. Parametrization involved reoptimizing available parameters as well as extending the force field to new molecules with the goal of achieving self-consistency with respect to the Lennard-Jones and electrostatic parameters targeting QM and experimental hydration free energies. The resulting force field parameters achieve consistent treatment across the studied anions, facilitating more balanced simulations of biomolecules and small organic molecules in the context of the classical Drude polarizable force field. Graphical abstract.

Induced Polarization in Molecular Dynamics Simulations of the 5-HT3 Receptor Channel

Ion channel proteins form water-filled nanoscale pores within lipid bilayers, and their properties are dependent on the complex behavior of water in a nanoconfined environment. Using a simplified model of the pore of the 5-HT3 receptor (5HT3R) which restrains the backbone structure to that of the parent channel protein from which it is derived, we compare additive with polarizable models in describing the behavior of water in nanopores. Molecular dynamics simulations were performed with four conformations of the channel: two closed state structures, an intermediate state, and an open state, each embedded in a phosphatidylcholine bilayer. Water density profiles revealed that for all water models, the closed and intermediate states exhibited strong dewetting within the central hydrophobic gate region of the pore. However, the open state conformation exhibited varying degrees of hydration, ranging from partial wetting for the TIP4P/2005 water model to complete wetting for the polarizable AMOEBA14 model. Water dipole moments calculated using polarizable force fields also revealed that water molecules remaining within dewetted sections of the pore resemble gas phase water. Free energy profiles for Na+ and for Cl- ions within the open state pore revealed more rugged energy landscapes using polarizable force fields, and the hydration number profiles of these ions were also sensitive to induced polarization resulting in a substantive reduction of the number of waters within the first hydration shell of Cl- while it permeates the pore. These results demonstrate that induced polarization can influence the complex behavior of water and ions within nanoscale pores and provides important new insights into their chemical properties.

Developing an effective polarizable bond method for small molecules with application to optimized molecular docking

Electrostatic interaction plays an essential role in protein-ligand binding. Due to the polarization effect, electrostatic interactions are largely impacted by their local environments. However, traditional force fields use fixed point charge-charge interactions to describe electrostatic interactions but is unable to include the polarization effect. The lack of the polarization effect in the force field representation can result in substantial error in biomolecular studies, such as molecular dynamics and molecular docking. Docking programs usually employ traditional force fields to estimate the binding energy between a ligand and a protein for pose selection or scoring. The intermolecular interaction energy mainly consists of van der Waals and electrostatic interaction in the force field representation. In the current study, we developed an Effective Polarizable Bond (EPB) method for small organic molecules and applied this EPB method to optimize protein-ligand docking in computational tests for a variety of protein-ligand systems. We tested the method on a set of 38 cocrystallized structures taken from the Protein Data Bank (PDB) and found that the maximum error was reduced from 7.98 Å to 2.03 Å when using EPB Dock, providing strong evidence that the use of EPB charges is important. We found that our optimized docking approach with EPB charges could improve the docking performance, sometimes dramatically, and the maximum error was reduced from 12.88 Å to 1.57 Å in Optimized Docking (in the case of 1fqx). The average RMSD decreased from 2.83 Å to 1.85 Å. Further investigations showed that the use of the EBP method could enhance intermolecular hydrogen bonding, which is a major contributing factor to improved docking performance. Developed tools for the calculation of the polarized ligand charge from a protein-ligand complex structure with the EPB method are freely available on GitHub (https://github.com/Xundrug/EPB).

Interplay between packing, dimer interaction energy and morphology in a series of tricyclic imide crystals

Crystal morphology is a very important feature in many industrial applications. Tricyclic imides, derivatives of 10-oxa-4-azatricyclo[5.2.1.0^2,6]dec-8-ene-3,5-dione with differing small hydrophobic groups (Me, Et), were studied and grouped based on Etter's rule. Using experimental X-ray studies, dimer energy calculations, framework analysis and periodic DFT-D calculations, it is shown that knowledge of the hydrogen-bond pattern can be used to determine the final crystal shape. Molecules forming a ring hydrogen-bond motif crystallize as plate crystals with the {100} facet as the slowest growing, whereas those molecules forming an infinite hydrogen-bond motif in the crystal structure crystallize as needles with the {101} facet having the largest surface area.

Pairwise-additive and polarizable atomistic force fields for molecular dynamics simulations of proteins

Protein force fields have been undergoing continual development since the first complete parameter sets were introduced nearly four decades ago. The functional forms that underlie these models have many common elements for the treatment of bonded and nonbonded forces, which are reviewed here. The most widely used force fields to date use a fixed-charge convention in which electronic polarization effects are treated via a mean-field approximation during partial charge assignment. Despite success in modeling folded proteins over many years, the fixed-charge assumption has limitations that cannot necessarily be overcome within their potential energy equations. To overcome these limitations, several force fields have recently been derived that explicitly treat electronic polarization effects with straightforward extensions of the potential energy functions used by nonpolarizable force fields. Here, we review the history of the most popular nonpolarizable force fields (AMBER, CHARMM, OPLS, and GROMOS) as well as studies that have validated them and applied them to studies of protein folding and misfolding. Building upon these force fields are more recent polarizable interaction potentials, including fluctuating charge models, POSSIM, AMOEBA, and the classical Drude oscillator. These force fields differ in their implementations but all attempt to model electronic polarization in a computationally tractable manner. Despite their recent emergence in the field of protein folding, several studies have already applied these polarizable models to challenging problems in this domain, including the role of polarization in folding free energies and sequence-specific effects on the stability of α-helical structures.

Polarizable Embedding with a Transferable H2O Potential Function I: Formulation and Tests on Dimer

The incorporation of mutual polarization in multiscale simulations where different regions of the system are treated at different level of theory is important in studies of, for example, electronic excitations and charge transfer processes. We present here an energy functional for describing a quantum mechanics/molecular mechanics (QM/MM) scheme that includes reciprocal polarization between the two subsystems. The inclusion of polarization alleviates shortcomings inherent in electrostatic embedding QM/MM models based on point-charge force fields. A density functional theory (DFT) description of the QM subsystem is coupled to a single center multipole expansion (SCME) description of H₂O molecules in the MM subsystem that includes anisotropic dipole and quadrupole polarizability as well as static multipoles up to and including the hexadecapole. The energy functional and the coupling scheme is general and can be extended to arbitrary order in terms of both the static and induced moments. Tests of the energy surface for the H₂O dimer show that the QM/MM results lie in between the pure DFT and pure SCME values. The consistency of the many-body contributions to the energy and analytical forces is demonstrated for an H₂O pentamer.

Accurate Biomolecular Simulations Account for Electronic Polarization

In this perspective, we discuss where and how accounting for electronic many-body polarization affects the accuracy of classical molecular dynamics simulations of biomolecules. While the effects of electronic polarization are highly pronounced for molecules with an opposite total charge, they are also non-negligible for interactions with overall neutral molecules. For instance, neglecting these effects in important biomolecules like amino acids and phospholipids affects the structure of proteins and membranes having a large impact on interpreting experimental data as well as building coarse grained models. With the combined advances in theory, algorithms and computational power it is currently realistic to perform simulations with explicit polarizable dipoles on systems with relevant sizes and complexity. Alternatively, the effects of electronic polarization can also be included at zero additional computational cost compared to standard fixed-charge force fields using the electronic continuum correction, as was recently demonstrated for several classes of biomolecules.

A force consistent method for electrostatic energy calculation in fluctuating charge model

A practical approach to include the polarization effect in a molecular force field is the fluctuating charge method in which atomic charges vary as the configuration of the molecular system changes. However, the use of the Coulomb formula to evaluate energy in a fluctuating charge method is theoretically inconsistent with the forces given by the fluctuating method. In this work, we propose a force-consistent method to correctly calculate electrostatic energies of molecular systems using a fluctuating charge model (Effective Polarizable Bond or EPB). In this protocol, the electrostatic energy is obtained by numerical interaction of the atomic forces along the MD trajectory, rather than using the default Coulomb formula in the EPB model. Test study on the benchmark Barnase-Barstar protein-protein interaction system demonstrates that although the total electrostatic energy of the system shows little deviation due to the averaging effect, specific residue-residue electrostatic interaction energy is affected and the level of the effect depends on the charges of the interacting residues with charged residues showing pronounced differences in calculated energies between using the current protocol and the standard Coulomb formula. It is recommended that the proposed numerical interaction method should be preferred in the calculation of electrostatic energy in fluctuating charge models used in molecular dynamics simulations.

Hydrogen bond donors and acceptors are generally depolarized in α-helices as revealed by a molecular tailoring approach

Hydrogen-bond (H-bond) interaction energies in α-helices of short alanine peptides were systematically examined by precise density functional theory calculations, followed by a molecular tailoring approach. The contribution of each H-bond interaction in α-helices was estimated in detail from the entire conformation energies, and the results were compared with those in the minimal H-bond models, in which only H-bond donors and acceptors exist with the capping methyl groups. The former interaction energies were always significantly weaker than the latter energies, when the same geometries of the H-bond donors and acceptors were applied. The chemical origin of this phenomenon was investigated by analyzing the differences among the electronic structures of the local peptide backbones of the α-helices and those of the minimal H-bond models. Consequently, we found that the reduced H-bond energy originated from the depolarizations of both the H-bond donor and acceptor groups, due to the repulsive interactions with the neighboring polar peptide groups in the α-helix backbone. The classical force fields provide similar H-bond energies to those in the minimal H-bond models, which ignore the current depolarization effect, and thus they overestimate the actual H-bond energies in α-helices. © 2019 The Authors. Journal of Computational Chemistry published by Wiley Periodicals, Inc.

Computational and Experimental Studies of Inhibitor Design for Aldolase A

Glycolytic enzyme fructose-bisphosphate aldolase A is an emerging therapeutic target in cancer. Recently, we have solved the crystal structure of murine aldolase in complex with naphthalene-2,6-diyl bisphosphate (ND1) that served as a template of the design of bisphosphate-based inhibitors. In this work, a series of ND1 analogues containing difluoromethylene (-CF₂), methylene (-CH₂), or aldehyde substitutions were designed. All designed compounds were studied using molecular dynamics (MD) simulations with the AMOEBA force field. Both energetics and structural analyses have been done to understand the calculated binding free energies. The average distances between ligand and protein atoms for ND1 were very similar to those for the ND1 crystal structure, which indicates that our MD simulation is sampling the correct conformation well. CF₂ insertion lowers the binding free energy by 10-15 kcal/mol, while CF₂ substitution slightly increases the binding free energy, which matches the experimental measurement. In addition, we found that NDB with two CF₂ insertions, the strongest binder, is entropically driven, while others including NDA with one CF₂ insertion are all enthalpically driven. This work provides insights into the mechanisms underlying protein-phosphate binding and enhances the capability of applying computational and theoretical frameworks to model, predict, and design diagnostic strategies targeting cancer.

Modeling Charge Flux by Interpolating Atomic Partial Charges

Although the charge flux effect or the geometric dependence of the atomic partial charges have been known for a long time, how it can be effectively handled is not yet established. Here, we present a charge interpolation scheme as a new general tool for representing the charge flux in an analytically well-defined manner. By applying it to the anionic GFP chromophore with the diabatically represented atomic charges, we show that the charge interpolation provides a substantial improvement on the accuracy of the geometry-dependent changes in the molecular dipole moments in the gas phase. We also test the scheme toward describing the electrostatic term in the solvation energy in the aqueous environment and observe that it is also improved but that the extent of the improvement is somewhat limited. We show that the remaining errors can be largely corrected by introducing atomic polarizabilities. Overall, our results show that charge interpolation is an amenable approach for describing the charge flux effect and that its description in the condensed phase should be accompanied by proper treatments of polarization effects.