Drude Force Field Research Articles

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.

Optimization of hydrogen-deuterium exchange ensemble reweighting (HDXer) to aid the biophysical characterization of proteins.

Hydrogen-deuterium exchange mass spectrometry (HDX-MS) experiments provide unique insights into native protein structure and dynamics in solution but have limited spatial resolution. Molecular dynamics (MD) simulations, in contrast, allow for the exploration of protein conformational landscapes with an atomistic level of resolution. However, forcefield inaccuracies, crystallographic biases, and timescale limitations may compromise the validity of the final structural ensemble. Hydrogen Deuterium eXchange Ensemble Reweighting (HDXer) is a maximum entropy approach that was recently developed to (1) quantitatively predict hydrogen-deuterium exchange from a modeled conformational ensemble and (2) refine the candidate ensemble to best fit an experimental HDX-MS dataset. Here, we compare how various computational and experimental parameters impact the performance of HDXer, using extracellular regulated kinase 2 (ERK2) and hen egg white lysozyme (HEWL) as model proteins. Conventional molecular dynamics simulations of ERK2 and HEWL were performed in triplicate using the CHARMM36m, AMBER FF19SB, and Drude forcefields. After predicting the HDX rates of each candidate ensemble, HDXer was used to fit the simulated structures to the corresponding HDX-MS datasets. To assess the influence of simulation length, truncated and strided copies of the full-length simulations were generated. Other factors, including peptide charge state selection, ion mobility mass spectrometry, and independent HDX-MS data acquisitions, were also tested. We find that while some variables (e.g., high vs. low peptide charge states) do not have a significant effect, other parameters (e.g., early vs. late timepoint selection) directly correlate with the disagreement between calculated and experimental deuterium uptake. The results from this benchmarking inform the development of a recommended HDXer workflow and best practices for optimal data acquisition.

Partitioning free energies and induced dipole response of small molecules across a POPC membrane using polarizable molecular dynamics.

The cell membrane is the semi-permeable barrier that regulates flux of materials in and out of the cell. The membrane has a characteristic dielectric gradient due to the amphiphilic nature of its constituent lipids, which have polar head groups and hydrophobic tails. This gradient plays governs many aspects of small-molecule partitioning and membrane protein folding and function. Atomistic details of small-molecule partitioning related to hydration and electronic polarization are still largely unknown. Here, we used polarizable molecular dynamics simulations to understand how the electronic environment influences the process of partitioning. We performed a series of unbiased and umbrella sampling (US) simulations using the Drude polarizable force field (FF) for six different amino-acid sidechain analogs in a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer. Using a polarizable FF such as Drude allows for the investigation of explicit dipole responses as a function of local electric fields along the membrane normal. We assessed molecular dipole moments, which vary as a function of localization within the membrane, and compared outcomes of the Drude FF simulations to those of the nonpolarizable CHARMM36 FF to investigate the differences in the free energy surfaces of partitioning for our different amino-acid analogs. In doing so, we sought to quantify the effects of using the fixed-charge convention of nonpolarizable FFs, which may average out the influence of the membrane dielectric gradient. We show how polarity of small molecules influences their partitioning and describe key differences in the free energy surfaces between polarizable and nonpolarizable systems in these contexts. This work serves as a foundation for future investigations into polarizable simulations of drug partitioning, as well as simulations of membrane proteins.

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.

Deep Neural Network Model to Predict the Electrostatic Parameters in the Polarizable Classical Drude Oscillator Force Field.

The Drude polarizable force field (FF) captures electronic polarization effects via auxiliary Drude particles that are attached to non-hydrogen atoms, distinguishing it from commonly used additive FFs that rely on fixed charges. The Drude FF currently includes parameters for biomolecules such as proteins, nucleic acids, lipids, and carbohydrates and small-molecule representative of those classes of molecules as well as a range of atomic ions. Extension of the Drude FF to novel small druglike molecules is challenging as it requires the assignment of partial charges, atomic polarizabilities, and Thole scaling factors. In the present article, deep neural network (DNN) models are trained on quantum mechanical (QM)-based partial charges and atomic polarizabilities along with Thole scale factors trained to target QM molecular dipole moments and polarizabilities. Training of the DNN model used a collection of 39 421 molecules with molecular weights up to 200 Da and containing H, C, N, O, P, S, F, Cl, Br, or I atoms. The DNN model utilizes bond connectivity, including 1,2, 1,3, 1,4, and 1,5 terms and distances of Drude FF atom types as the feature vector to build the model, allowing it to capture both local and nonlocal effects in the molecules. Novel methods have been developed to determine restrained electrostatic potential (RESP) charges on atoms and external points representing lone pairs and to determine Thole scale factors, which have no QM analogue. A penalty scheme is devised as a performance predictor of the trained model. Validation studies show that these DNN models can precisely predict molecular dipole and polarizabilities of Food and Drug Administration (FDA)-approved drugs compared to reference MP2 calculations. The availability of the DNN model allowing for the rapid estimation of the Drude electrostatic parameters will facilitate its applicability to a wider range of molecular species.

Testing and Optimizing the Drude Polarizable Force Field for Blocked Amino Acids Based on High-Level Quantum-Mechanical Energy Surfaces

The correct reproduction of conformational substates of amino acids was tested for the CHARMM Drude polarizable force field. This was achieved by evaluating the reorganization energies for all low lying energy minima occurring in all 15 neutral blocked amino acids on a quantum-mechanical (QM) energy surface at the MP2/cc-pVDZ level. The results indicate that the bonded parameters of the N-acetyl (ACE) and N-Methylamide (CT3) blocking groups lead to significant discrepancies. A reparametrization of five bond angles significantly improved the agreement with the QM energy surface. The corrected Drude force field exhibits almost the same average reorganization energies relative to the MP2 energy surface as the AM1 and PM3 semi-empirical methods.

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 Ca2+ 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 Ca2+ 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.

CHARMM‐GUI Drude prepper for molecular dynamics simulation using the classical Drude polarizable force field

Explicit treatment of electronic polarizability in empirical force fields (FFs) represents an extension over a traditional additive or pairwise FF and provides a more realistic model of the variations in electronic structure in condensed phase, macromolecular simulations. To facilitate utilization of the polarizable FF based on the classical Drude oscillator model, Drude Prepper has been developed in CHARMM-GUI. Drude Prepper ingests additive CHARMM protein structures file (PSF) and pre-equilibrated coordinates in CHARMM, PDB, or NAMD format, from which the molecular components of the system are identified. These include all residues and patches connecting those residues along with water, ions, and other solute molecules. This information is then used to construct the Drude FF-based PSF using molecular generation capabilities in CHARMM, followed by minimization and equilibration. In addition, inputs are generated for molecular dynamics (MD) simulations using CHARMM, GROMACS, NAMD, and OpenMM. Validation of the Drude Prepper protocol and inputs is performed through conversion and MD simulations of various heterogeneous systems that include proteins, nucleic acids, lipids, polysaccharides, and atomic ions using the aforementioned simulation packages. Stable simulations are obtained in all studied systems, including 5 μs simulation of ubiquitin, verifying the integrity of the generated Drude PSFs. In addition, the ability of the Drude FF to model variations in electronic structure is shown through dipole moment analysis in selected systems. The capabilities and availability of Drude Prepper in CHARMM-GUI is anticipated to greatly facilitate the application of the Drude FF to a range of condensed phase, macromolecular systems.

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.

Benchmark Force Fields for the Molecular Dynamic Simulation of G-Quadruplexes.

G-quadruplexes have drawn widespread attention for serving as a potential anti-cancer target and their application in material science. Molecular dynamics (MD) simulation is the key theoretical tool in the study of GQ’s structure-function relationship. In this article, we systematically benchmarked the five force fields of parmbsc0, parmbsc1, OL15, AMOEBA, and Drude2017 on the MD simulation of G-quadruplex from four aspects: structural stability, central ion channel stability, description of Hoogsteen hydrogen bond network, and description of the main chain dihedral angle. The results show that the overall performance of the Drude force field is the best. Although there may be a certain over-polarization effect, it is still the best choice for the MD simulation of G-quadruplexes.

Ion Binding to the bcl-2 G-Quadruplex from Polarizable Simulations with the Drude Force Field

G-quadruplexes (GQs) are noncanonical nucleic acid structures that can form in guanine-rich sequences and can adopt many diverse conformations depending on their constituent nucleotides. GQs are commonly present in DNA promoter regions, so alterations to their structures can contribute to diseases such as cancer. The B-cell lymphoma 2 (bcl-2) proto-oncogene codes for a protein that inhibits cell apoptosis and overexpression of bcl-2 is associated with cancer progression. Transcription of bcl-2 is controlled by two promoters, one of which contains a GQ-forming sequence. Understanding the stabilizing forces in this structure may lead to novel chemotherapeutic design. There are gaps in our current understanding of under what conditions GQs fold, which is crucial information that can give insight into the equilibrium of GQ conformations. To advance our understanding of the dominant bcl-2 promoter GQ, we performed molecular dynamics (MD) simulations using the Drude-2017 polarizable force field. Our simulation outcomes highlight two distinct changes within the GQ that coincide with ion binding. The first is the recruitment of a bulk K+ ion to the solvent-exposed face of the tetrad stem, which led to changes in base dipole moments and the K+ interaction energies of the guanines within the tetrad core. The second is the emergence of an “electronegative pocket” between the tetrad core and the long loop below the tetrad core. Ion binding in the electronegative pocket induces a backbone conformational change that may be relevant for inducing local conformational changes at K+ concentrations similar to those found in the cell. These results suggest two different locations within the bcl-2 GQ that may be targeted by small molecules to induce conformational changes.

Characterizing the water wire in the Gramicidin channel found by Monte Carlo sampling using continuum electrostatics and in molecular dynamics trajectories with conventional or polarizable force fields

Water molecules play a key role in all biochemical processes. They help define the shape of proteins, and they are reactant or product in many reactions and are released as ligands are bound. They facilitate the transfer of protons through transmembrane proton channel, pump and transporter proteins. Continuum electrostatics (CE) force fields used by program Multiconformation CE (MCCE) capture electrostatic interactions in biomolecules with an implicit solvent, which captures the averaged solvent water equilibrium properties. Hybrid CE methods can use explicit water molecules within the protein surrounded by implicit solvent. These hybrid methods permit the study of explicit hydrogen bond networks within the protein and allow analysis of processes such as proton transfer reactions. Yet hybrid CE methods have not been rigorously tested. Here, we present an explicit treatment of water molecules in the Gramicidin A (gA) channel using MCCE and compare the resulting distributions of water molecules and key hydration features against those obtained with explicit solvent Molecular Dynamics (MD) simulations with the nonpolarizable CHARMM36 and polarizable Drude force fields. CHARMM36 leads to an aligned water wire in the channel characterized by a large absolute net water dipole moment; the MCCE and Drude analysis lead to a small net dipole moment as the water molecules change orientation within the channel. The correct orientation is not as yet known, so these calculations identify an open question.

Characterizing the Water Wire in the Gramicidin Channel Found by Monte Carlo Sampling Using Continuum Electrostatics and in Molecular Dynamics Trajectories with Conventional or Polarizable Force Fields

Water molecules play a key role in all biochemical processes. They help define the shape of proteins, and they are reactant or product in many reactions and are released as ligands are bound. They facilitate the transfer of protons through transmembrane proton channel, pump and transporter proteins. Continuum electrostatics (CE) force fields used by program Multiconformation CE (MCCE) capture electrostatic interactions in biomolecules with an implicit solvent, which captures the averaged solvent water equilibrium properties. Hybrid CE methods can use explicit water molecules within the protein surrounded by implicit solvent. These hybrid methods permit the study of explicit hydrogen bond networks within the protein and allow analysis of processes such as proton transfer reactions. Yet hybrid CE methods have not been rigorously tested. Here, we present an explicit treatment of water molecules in the Gramicidin A (gA) channel using MCCE and compare the resulting distributions of water molecules and key hydration features against those obtained with explicit solvent Molecular Dynamics (MD) simulations with the nonpolarizable CHARMM36 and polarizable Drude force fields. CHARMM36 leads to an aligned water wire in the channel characterized by a large absolute net water dipole moment; the MCCE and Drude analysis lead to a small net dipole moment as the water molecules change orientation within the channel. The correct orientation is not as yet known, so these calculations identify an open question.

Dynamics of the 1:2:1 and 1:6:1 C-Myc G-Quadruplexes with the Drude Polarizable Force Field

G-quadruplexes (GQs) are noncanonical nucleic acid structures that form in guanine-rich sequences, importantly gene promoters and telomeres. The guanines within GQs arrange in square, planar tetrads through Hoogsteen hydrogen bonding. These tetrads stack on top of one another to make a highly ordered structure with an electronegative core. Coordinated cations are required to counteract the core's electronegativity and stabilize the GQ. GQs play important roles in gene expression and genomic stability and are functionally relevant in several human proto-oncogenes. One such proto-oncogene, c-myc, is linked to several human cancers such as Burkitt's lymphoma and multiples types of myeloma due to GQ formation within the promoter region of the gene acting as a gene silencer. As such, the c-myc GQ is a novel drug target for cancer therapy. Two solved GQ structures of interest in the c-myc promoter are the 1:2:1 and 1:6:1 GQs, named for differences in central loop length. Loop regions of GQs are known to affect topology, stability, and overall structure. Moreover, studies have shown nucleolin preferentially binds to the 1:6:1 GQ, prompting our present study on the dynamics of c-myc GQs. To analyze the importance of loop sequence on structural dynamics and gain information on this potential drug target, we performed molecular dynamics simulations utilizing the CHARMM36 non-polarizable and the Drude-2017 polarizable force fields (FF). The Drude FF accounts for electronic degrees of freedom via negatively charged particles placed on every heavy atom in the systems and has provided new insights on GQ dynamics and DNA-ion interactions. Our simulations of the 1:2:1 and 1:6:1 GQs demonstrate that polarization is important to model the stability, ion-ion, and ion-DNA interactions that underly GQ dynamics, emphasizing the importance of polarizable molecular dynamics in studying GQs.

FFParam: Standalone package for CHARMM additive and Drude polarizable force field parametrization of small molecules.

Accurate force-field (FF) parameters are key to reliable prediction of properties obtained from molecular modeling (MM) and molecular dynamics (MD) simulations. With ever-widening applicability of MD simulations, robust parameters need to be generated for a wider range of chemical species. The CHARMM General Force Field program (CGenFF, https://cgenff.umaryland.edu/) is a tool for obtaining initial parameters for a given small molecule based on analogy with the available CGenFF parameters. However, improvement of these parameters is often required and performing their optimization remains tedious and time consuming. In addition, tools for optimization of small molecule parameters in the context of the Drude polarizable FF are not yet available. To overcome these issues, the FFParam package has been designed to facilitate the parametrization process. The package includes a graphical user interface (GUI) created using Qt libraries. FFParam supports Gaussian and Psi4 for performing quantum mechanical calculations and CHARMM and OpenMM for MM calculations. A Monte Carlo simulated annealing (MCSA) algorithm has been implemented for automated fitting of partial atomic charge, atomic polarizabilities and Thole scale parameters. The LSFITPAR program is called for automated fitting of bonded parameters. Accordingly, FFParam provides all the features required for generation and analysis of CHARMM and Drude FF parameters for small molecules. FFParam-GUI includes a text editor, graph plotter, molecular visualization, and text to table converter to meet various requirements of the parametrization process. It is anticipated that FFParam will facilitate wider use of CGenFF as well as promote future use of the Drude polarizable FF.

Structure and Dynamics of the Bcl-2 Promoter G-Quadruplex using the Drude Polarizable Force Field

G-quadruplexes (GQ) are noncanonical nucleic acid structures that form in guanine-rich sequences and exist in a variety of topological forms. The guanine bases in GQ arrange in a square planar configuration and stack to form tetrads via Hoogsteen hydrogen bonding and coordination of ions like K+. These highly ordered secondary structures are implicated in the regulation of gene expression and genome maintenance, and alterations to their structures contribute to diseases such as cancer. The B-cell CLL/lymphoma 2 (bcl-2) oncogene is associated with a variety of cancers and its promoter features a GQ that arises due to two consecutive substitutions of guanine with thymine. Understanding the forces driving GQ folding and stabilization is required for developing novel therapeutics that target these structures. To develop a greater understanding of the atomistic details in the bcl-2 GQ, we performed molecular dynamics simulations using the CHARMM36 additive and Drude-2017 polarizable force fields. In the classical Drude oscillator model, electronic degrees of freedom are modeled as negatively charged particles attached harmonically to all heavy atoms. This approach overcomes limitations of nonpolarizable force fields, which approximate average polarization via assignment of partial charges. We performed simulations of the mutant bcl-2 promoter GQ and created the wild-type sequence by reverting the point mutations. The wild-type sequence was then thermally denatured and allowed to re-fold to determine the structural differences between the mutant bcl-2 GQ and the native, wild-type fold. Our analysis provides insight into the stability of the bcl-2 GQ structure, backbone conformational ensembles in dihedral space, base-ion interaction energetics, and the magnitude and fluctuation of dipole moments. These simulations illustrate that polarization is integral for tetrad stability and suggests that the Drude force field is essential for simulating GQ structures.

Optimized Parameters for the Drude Polarizable Force Field for Small Organic Molecules

Meaningful efforts in computer-aided drug dsign (CADD) need accurate molecular mechanical force fields to quantitatively characterize protein-ligand interactions. The most widely used strategy to ascribe force field parameters for new compounds apply pre-defined tabulated values following knowledge-based rules. Two widely used approaches following this strategy are the General Amber Force Field (GAFF) and the CHARMM General Force Field (CGenFF). While these approaches are extremely useful, an important limitation of using pre-tabulated parameter values is that they may be inadequate in the context of a specific molecule. To resolve this issue, we previously introduced the General Automated Atomic Model Parameterization (GAAMP) to automatically generate the parameters of atomic models of small molecules using the results from ab initio quantum mechanical (QM) calculations as target data [DOI:10.1021/ct4003477]. The GAAMP protocol uses QM data to optimize the bond, valence angle, and dihedral angle internal parameters, and atomic partial charges. However, the Lennard-Jones (LJ) 6-12 parameters cannot easily be refined from QM data and are unchanged from the initial atom types assignments severely limiting the achievable accuracy by these models. Here, we systematically optimize the LJ parameters of the Drude force field to reproduce experimental neat liquid densities and enthalpies of vaporization for a large set of compounds, covering a wide range of chemical functionalities. The optimization is carried out by constructing an objective function and exploiting the analytical gradient information within an optimization framework. The new set of LJ parameters will help to improve the accuracy of the computed physical properties of diverse small organic molecules. These developments will provide the necessary tools to widely apply drude polarizable force field model for drug discovery endeavors.

Improved Modeling of Halogenated Ligand-Protein Interactions Using the Drude Polarizable and CHARMM Additive Empirical Force Fields.

Halogenated ligands can participate in nonbonding interactions with proteins via halogen bond (XB) or halogen-hydrogen bond donor (X-HBD) interactions. In the context of molecular dynamics (MD) simulations, the accuracy of the simulations depends strongly on the force field (FF) used. To ensure good reproduction of XB and X-HBD interactions with proteins, we optimized the previously developed additive CHARMM36/CHARMM General force field (CGenFF) and Drude polarizable force field by including atom pair-specific Lennard-Jones parameters for aromatic halogen-protein interactions. The optimization targeted quantum mechanical interaction energy surfaces with the developed parameters then examined for their ability to reproduce experimental halogen-containing ligand-protein interactions in MD simulations. The calculated halogenated ligand interaction geometries were in good overall agreement with the experimental crystal data for both the polarizable and additive FFs, showing that these models can accurately treat both XB and X-HBD interactions. Analysis of the ligand-protein interactions shows significant contributions of polarizability to binding occurring in the Drude FF, with self-polarization energy making both favorable and unfavorable contributions to binding. Further analysis of the dipole moments from aqueous solution to protein indicates the polarizable FF accounts for subtle changes of the environment of the ligands that can impact binding. The present work demonstrates the utility of the updated additive CHARMM36/CGenFF and polarizable Drude FFs for the study of halogenated ligand-protein interactions in computer-aided drug design.

Comparative Analysis of Protein Hydration from MD simulations with Additive and Polarizable Force Fields

AbstractRecent development of the Drude polarizable (Drude) force field (FF), based on the extension of an induced dipole model, has reached a milestone in the past few years providing a complete set of polarizable parameters for proteins, water, ions, and many lipid types. This FF enables stable simulations up to microseconds, surpassing the capability of other polarizable FFs. The quality of the Drude FF, however, has remained largely untested for modeling the secondary structures of small peptides in explicit solvents compared with classical non‐polarizable FFs. It is critical to benchmark the complex and mutually dependent dynamics of hydrogen‐bond (H‐bond) networks formed by water–water, protein–water, and protein–protein interactions that are expected to have a major impact on the stability of protein structures and their conformational space. Here, a direct comparison is presented between the current Drude FF and the CHARMM‐36 non‐polarizable classical FF for 1) the solvation free energy of mimetics for all amino acid side‐chain equivalents, 2) limited conformational space, 3) protein–water and protein–protein interactions, and 4) the comparative lifetimes of H‐bonds. The impact of counterions on the stabilization of secondary structure in model peptides is additionally discussed and compared between these FFs.