Polarizable simulations with second order interaction model (POSSIM) force field: Developing parameters for protein side‐chain analogues

Modeling structural and thermodynamic properties of nucleic acids has long been a challenge in the development of force fields. Polarizable force fields are a new generation of potential functions to take charge redistribution and induced dipole into account, and have been proved to be reliable to model small molecules, polypeptides and proteins, but their use on nucleic acids is still rather limited. In this article, the interactions between nucleic acids and a small molecule or ion were modeled by AMOEBAbio09, a modern polarizable force field, and conventional non-polarizable AMBER99sb and CHARMM36 force fields. The resulting intermolecular interaction energies were compared with those calculated by ab initio quantum mechanics methods. Although the test is not sufficient to prove the reliability of the polarizable force field, the results at least validate its capability in modeling energetics of static configurations, which is one basic component in force field parameterization.

Molecular dynamics (MD) simulations have been widely applied to computer-aided drug design (CADD). While MD has been used in a variety of applications such as free energy perturbation and long-time simulations, the accuracy of the results from those methods depends strongly on the force field used. Force fields for small molecules are crucial, as they not only serve as building blocks for developing force fields for larger biomolecules but also act as model compounds that will be transferred to ligands used in CADD. Currently, a wide range of small molecule force fields based on additive or nonpolarizable models have been developed. While these nonpolarizable force fields can produce reasonable estimations of physical properties and have shown success in a variety of systems, there is still room for improvements due to inherent limitations in these models including the lack of an electronic polarization response. For this reason, incorporating polarization effects into the energy function underlying a force field is believed to be an important step forward, giving rise to the development of polarizable force fields. Recent simulations of biological systems have indicated that polarizable force fields are able to provide a better physical representation of intermolecular interactions and, in many cases, better agreement with experimental properties than nonpolarizable, additive force fields. Therefore, this chapter focuses on the development of small molecule force fields with emphasis on polarizable models. It begins with a brief introduction on the importance of small molecule force fields and their evolution from additive to polarizable force fields. Emphasis is placed on the additive CHARMM General Force Field and the polarizable force field based on the classical Drude oscillator. The theory for the Drude polarizable force field and results for small molecules are presented showing their improvements over the additive model. The potential importance of polarization for their application in a wide range of biological systems including CADD is then discussed.

The induced dipole and the classical Drude oscillator represent two major approaches for the explicit inclusion of electronic polarizability into force field-based molecular modeling and simulations. In this work, we explore the equivalency of these two models by comparing condensed phase properties computed using the Drude force field and a multipole and induced dipole (MPID) model. Presented is an approach to map the electrostatic model optimized in the context of the Drude force field onto the MPID model. Condensed phase simulations on water and 15 small model compounds show that without any reparametrization, the MPID model yields properties similar to the Drude force field with both models yielding satisfactory reproduction of a range of experimental values and quantum mechanical data. Our results illustrate that the Drude oscillator model and the point induced dipole model are different representations of essentially the same physical model. However, results indicate the presence of small differences between the use of atomic multipoles and off-center charge sites. Additionally, results on the use of dispersion particle mesh Ewald further support its utility for treating long-range Lennard Jones dispersion contributions in the context of polarizable force fields. The main motivation in demonstrating the transferability of parameters between the Drude and MPID models is that the more than 15 years of development of the Drude polarizable force field can now be used with MPID formalism without the need for dual-thermostat integrators nor self-consistent iterations. This opens up a wide range of new methodological opportunities for polarizable models.

Molecular dynamics simulations, based on molecular mechanics force fields, have been instrumental in increasing our understanding of a range of biomolecular systems. Most of the force fields commonly used for the simulation of biomolecules, however, represent electrostatic properties by a set of fixed partial atomic charges and are referred to as additive force fields. This approach is problematic because it does not explicitly include polarizability, an important component of the electrostatic interaction arising from the response of the molecular dipoles to an external electric field.To overcome this problem, multiple efforts are currently underway to develop force fields including an explicit representation of polarizability: one such effort is the CHARMM Drude polarizable force field. In this chapter, we begin by discussing the development and implementation of the theory associated with the Drude model, as well as a robust scheme for parameter optimization. This is followed by discussion of the way in which these tools are being used to optimize a polarizable force field for the simulation of biomolecules. The primary focus has been the development of parameters for small molecule analogues of functional groups present within biomolecules, but is shifting to the optimization of covalent connections between these small molecules, and testing of the resulting parameters in fully hydrated biomolecules. For the nucleic acids a case study is presented, illustrating that the polarizable force field already reveals atomic-level details not observed with an additive force field.

Condensed-phase computational studies of molecules using molecular mechanics approaches require the use of force fields to describe the energetics of the systems as a function of structure. The advantage of polarizable force fields over nonpolarizable (or additive) models lies in their ability to vary their electronic distribution as a function of the environment. Toward development of a polarizable force field for biological molecules, parameters for a series of sulfur-containing molecules are presented. Parameter optimization was performed to reproduce quantum mechanical and experimental data for gas phase properties including geometries, conformational energies, vibrational spectra, and dipole moments as well as for condensed phase properties such as heats of vaporization, molecular volumes, and free energies of hydration. Compounds in the training set include methanethiol, ethanethiol, propanethiol, ethyl methyl sulfide, and dimethyl disulfide. The molecular volumes and heats of vaporization are in good accordance with experimental values, with the polarizable model performing better than the CHARMM22 nonpolarizable force field. Improvements with the polarizable model were also obtained for molecular dipole moments and in the treatment of intermolecular interactions as a function of orientation, in part due to the presence of lone pairs and anisotropic atomic polarizability on the sulfur atoms. Significant advantage of the polarizable model was reflected in calculation of the dielectric constants, a property that CHARMM22 systematically underestimates. The ability of this polarizable model to accurately describe a range of gas and condensed phase properties paves the way for more accurate simulation studies of sulfur-containing molecules including cysteine and methionine residues in proteins.

Computational approaches are widely used to help discover and develop new drugs, in particular, to understand how these molecules interact with their biomolecular targets. There are well-optimised and validated force field parameters available to describe interactions of common biomolecules, such as proteins, lipids and nucleotides. However, these force fields are not designed to represent heteromolecular ligands such as substrates, inhibitors, co-factors and potential drug molecules. Errors in parameters may result in incorrect predictions of ligand structure, orientation and conformation, which in turn can lead to the failure of computational drug design efforts. With over 25% of all the structures in the Protein Data Bank (PDB) containing ligands and over a million ligand molecules of potential interest in drug design in other databases, ligand validation and parameterisation represents a significant scientific challenge. While several automated parameterisation protocols have been proposed to generate the parameters, none of the current procedures are properly validated. As a response to the high demand for interaction parameters for ligands compatible with the GROMOS force field, a web accessible Automated force field Topology Builder (ATB; http://atb.uq.edu.au/) and Repository was developed. The ATB and Repository is intended to facilitate the development of molecular force fields and generates parameters that can be used in X-ray refinement, structure-based drug design and study of biomolecule-ligand complexes. In this thesis, molecular dynamics (MD) simulations were used to calculate the thermodynamic properties for the validation and refinement of ATB force fields parameters. A fully automated validation protocol of the ATB force field parameters for small organic molecules based on thermodynamic and structural information was developed and incorporated into the ATB. A novel integration and convergence protocol which increases the efficiency of the TI method for free energy calculations was also proposed. The validation has shown good overall agreement between the experimental and computed values and indicated problematic functional groups with parameters to be refined. The generation of parameters for novel molecules that are compatible with a given biomolecular force field can be tedious, time-consuming and error-prone. Here, a novel method to refine parameters in classical force fields in an automated manner that can be extended to the parameterisation of all other atom types is also presented. Single-step perturbation protocols were developed for small halogenated molecules to establish alternative van der Waals parameters which describe the interactions of these molecules with high precision. Parameters were successfully refined against experimental hydration free energies, densities and heats of vaporisation. As a result, alternative van der Waals interactions parameters for chlorine and bromine were proposed. Properties calculated with these proposed parameters closely matched the experiment data. The quality of the ATB force field parameters in validating the X-ray structures and investigating binding of ligands to the protein endothiapepsin was also investigated. Analysis of the MD trajectories showed that the protein structure remained relatively stable after 10 ns of simulations. The results for the validation of the ligands were mixed. While in some cases the simulations reproduced the binding mode and the interactions between the protein and the ligand with high accuracy suggesting the ATB parameters performed well, in other cases the complexes were unstable. This was despite the fact that the electron density for all ligands was well defined and the parameters for the ligand were generated using the same procedure. Endothiapepsin contains a catalytic dyad consisting of two aspartates that are generally assumed to share a proton. We demonstrated that the protonation of both the ligand and the residues with the binding site were critical to the ability to reproduce the crystal complex. To conclude, the work presented in this thesis provided insight into the development and validation of the force field parameters for drug-like molecules and their applications in structure-based drug design.

Development of accurate force field parameters for molecular ions in the context of a polarizable energy function based on the classical Drude oscillator is a crucial step toward an accurate polarizable model for modeling and simulations of biological macromolecules. Toward this goal we have undertaken a hierarchical approach in which force field parameter optimization is initially performed for small molecules for which experimental data exists that serve as building blocks of macromolecular systems. Small molecules representative of the ionic moieties of biological macromolecules include the cationic ammonium and methyl substituted ammonium derivatives, imidazolium, guanidinium and methylguanidinium, and the anionic acetate, phenolate, and alkanethiolates. In the present work, parameters for molecular ions in the context of the Drude polarizable force field are optimized and compared to results from the nonpolarizable additive CHARMM general force field (CGenFF). Electrostatic and Lennard-Jones parameters for the model compounds are developed in the context of the polarizable SWM4-NDP water model, with emphasis on assuring that the hydration free energies are consistent with previously reported parameters for atomic ions. The final parameters are shown to be in good agreement with the selected quantum mechanical (QM) and experimental target data. Analysis of the structure of water around the ions reveals substantial differences between the Drude and additive force fields indicating the important role of polarization in dictating the molecular details of aqueous solvation. The presented parameters represent the foundation for the charged functionalities in future generations of the Drude polarizable force field for biological macromolecules as well as for drug-like molecules.

The reliability of molecular dynamics (MD) simulations in predicting macroscopic properties of complex fluids and soft materials, such as liquid crystals, colloidal suspensions, or polymers, relies on the accuracy of the adopted force field (FF). We present an automated protocol to derive specific and accurate FFs, fully based on ab initio quantum mechanical (QM) data. The integration of the Joyce and Picky procedures, recently proposed by our group to provide an accurate description of simple liquids, is here extended to larger molecules, capable of exhibiting more complex fluid phases. While the standard Joyce protocol is employed to parameterize the intramolecular FF term, a new automated procedure is here proposed to handle the computational cost of the QM calculations required for the parameterization of the intermolecular FF term. The latter is thus obtained by integrating the old Picky procedure with a fragmentation reconstruction method (FRM) that allows for a reliable, yet computationally feasible sampling of the intermolecular energy surface at the QM level. The whole FF parameterization protocol is tested on a benchmark liquid crystal, and the performances of the resulting quantum mechanically derived (QMD) FF were compared with those delivered by a general-purpose, transferable one, and by the third, "hybrid" FF, where only the bonded terms were refined against QM data. Lengthy atomistic MD simulations are carried out with each FF on extended 5CB systems in both isotropic and nematic phases, eventually validating the proposed protocol by comparing the resulting macroscopic properties with other computational models and with experiments. The QMD-FF yields the best performances, reproducing both phases in the correct range of temperatures and well describing their structure, dynamics, and thermodynamic properties, thus providing a clear protocol that may be explored to predict such properties on other complex fluids or soft materials.

This study proposes a Polarizable Transferable Anisotropic Mie (P-TAMie) force field for ethers and mixtures of ethers with n-alkanes and nitrogen. In our approach, polarizability is accounted for using a charge-on-spring model (also referred to as the Drude model or core-shell model) placed on every united-atom site. Parameters for polarizabilities are obtained from ab initio calculations, without the need to introduce empirical scaling factors. The van der Waals parameters and static partial charges are optimized to minimize deviations to experimental vapor pressure and coexisting liquid density data. Vapor-liquid phase equilibria are calculated using Monte Carlo simulations in the grand canonical ensemble, using a MC scheme with multiparticle moves. Correlation results for pure substances and predictions of mixture behavior obtained from the proposed polarizable force field (P-TAMie) are in good agreement with experimental data. We compare the polarizable force field with an analogous nonpolarizable force field, namely, the Transferable Anisotropic Mie (TAMie) force field, where both force fields share the same intramolecular force field and were parametrized using the same objective function. The comparison should thus single out the role of polarizability in predicting phase equilibria. The comparison of P-TAMie with the nonpolarizable TAMie force field, however, shows very similar results of both models. This study does not include water and aqueous mixtures or other pronounced polar substances.

Due to the recent increase in computing power, the molecular modeling community has been focused on improving the accuracy and overall quality of biomolecular simulations. These technological improvements centered on the development of new force fields and simulation packages, allowed for more complex and heavier systems to be tackled, while maintaining, or even improving simulations speed. Some force fields, such as GROMOS, have been parameterized and validated using a reaction field (RF), charge groups and a twin-range cutoff scheme (0.8 and 1.4 nm) to treat long range electrostatics [1]. However, in GROMACS software package, the use of group-based cutoff scheme will be deprecated in future versions. As to properly use the newer and faster versions of this simulation package coupled with GROMOS 54A7 and RF, it is crucial to assess the impact on the system sampling when using a single atomistic cutoff (based on the Verlet method) instead of the twin-range group-based scheme. We reproduced the GROMOS parameterization procedure [1] with both schemes and measured very similar hydration free energy values of small amino acid side chains analogs between both protocols [2]. However, we observed a small, yet significant, difference on the conformational spaces of G1-PAMAM and DMPC when using the atomistic cutoff scheme [2]. The DMPC Al values decrease to a region outside the experimental range, which is not surprising since the force field parameters were optimized for a group-based scheme. Nevertheless, the structural properties for both systems are better converged for the used atomistic cutoff range (1.4-2.0 nm) relative to the group-based cutoff simulation set [2]. The use of a single atomistic cutoff scheme seems a viable approach for MD simulations of biomolecules using G54A7 force field, even if in some cases, new calibration protocols are needed. (1) Oostenbrink, C.; Villa, A.; Mark, A. E.; van Gunsteren, W. F. A biomolecular force field based on the free enthalpy of hydration and solvation: The GROMOS force-field parameter sets 53A5 and 53A6. J. Comput. Chem. 2004, 25, 1656–1676 (2) Silva, T. F. D., Vila-Viçosa, D., Reis, P. B. P. S., Victor, B. L., Diem, M., Oostenbrink, C., and Machuqueiro, M. (2018) "The impact of using single atomistic long range cutoff schemes with the GROMOS 54A7 force field", J. Chem. Theory Comput., 14, 5823-5833

Development of a ReaxFF reactive force field for intrinsic point defects in titanium dioxide

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.

The structural flexibility of DNA plays a key role in many biological processes of DNA, such as protein-DNA interactions, DNA packaging in viruses and nucleosome positioning on genomic DNA. Some experimental techniques have been employed to investigate the structural flexibility of DNA with the combination of elastic models, but these experiments could only provide the macroscopic properties of DNA, and thus, it is still difficult to understand the corresponding microscopic mechanisms. Recently, all-atom molecular dynamics (MD) simulation has emerged as a useful tool to investigate not only the macroscopic properties of DNA, but also the microscopic description of the flexibility of DNA at an atomic level. The most important issue in all-atom MD simulations of DNA is to choose an appropriate force field for simulating DNA. Very recently, a new force field for DNA has been developed based on the last generation force field of Amber bsc0, which was named Amber bsc1. In this work, all-atom MD simulations are employed to study the flexibility of a 30-bp DNA with the force fields of Amber bsc1 and Amber bsc0 in a comparative way. Our aim of the research is to examine the improvement of the new development of force field (Amber bsc1) in the macroscopic and microscopic properties of DNA, in comparison with the corresponding experimental measurements. All the MD simulations are performed with Gromacs 4.6 and lasted with a simulation time of 600 ns. The MD trajectories are analyzed with Curves+ for the last 500 ns, since the system reaches equilibrium approximately after ~100 ns. Our results show that the new force field (Amber bsc1) can lead to the improvements in the macroscopic parameters of DNA flexibility, i.e., stretch modulus S and twist-stretch coupling D become closer to experimental measurements, while bending persistence lengths lp and torsional persistence lengths C from the two force fields (bsc1 and bsc0) are both in good agreement with experimental data. Our microscopic analyses show that the microscopic structure parameters of DNA from the MD simulation with the Amber bsc1 force field are closer to the experimental values than those with the Amber bsc0 force field, except for slide, and the obvious improvements are observed in some microscopic parameters such as twist and inclination. Our further analyses show that the improvements in macroscopic flexibility from the Amber bsc1 force field are tightly related to the microscopic parameters and their fluctuations. This study would be helpful in understanding the performances of Amber bsc1 and bsc0 force fields in the description of DNA flexibility at both macroscopic and microscopic level.

The development of accurate molecular mechanics force fields is a significant challenge that must be addressed for the continued success of molecular simulation. We developed the ForceBalance method to automatically derive accurate force field parameters using flexible combinations of experimental and theoretical reference data. The method is demonstrated in the parametrization of two rigid water models, yielding new parameter sets (TIP3P-FB and TIP4P-FB) that accurately describe many physical properties of water.

We developed a ReaxFF reactive force field for NASICON-type Li1+xAlxTi2-x(PO4)3 (LATP) materials, which is a promising solid-electrolyte that may enable all-solid-state lithium-ion batteries. The force field parameters were optimized based on density functional theory (DFT) data, including equations of state and the heats of formation of ternary metal oxides and metal phosphate crystal phases (e.g., LixTiO2, Al2TiO5, LiAlO2, AlPO4, Li3PO4 and LiTi2(PO4)3 (LTP)), and the energy barriers for Li diffusion in TiO2 and LTP via vacancies and interstitial sites. Using ReaxFF, the structural and the energetic features of LATP were described properly across various compositions - Li occupies more preferentially the interstitial site next to Al than next to Ti. Also, as observed in experimental data, the lattice parameters decrease when Ti is partly substituted by Al because of the smaller size of the Al cation. Using this force field, the diffusion mechanism and the ionic conductivity of Li in LTP and LATP were investigated at T = 300-1100 K. Low ionic conductivity (5.9 × 10-5 S cm-1 at 300 K) was obtained in LTP as previously reported. In LATP at x = 0.2, the ionic conductivity was slightly improved (8.4 × 10-5 S cm-1), but it is still below the experimental value, which is on the order of 10-4 to 10-3 S cm-1 at x = 0.3-0.5. At higher x (higher Al composition), LATP has a configurational diversity due to the Al substitution and the concomitant insertion of Li. By performing a hybrid MC/MD simulation for LATP at x = 0.5, a thermodynamically stable LATP configuration was obtained. The ionic conductivity of this LATP configuration was calculated to be 7.4 × 10-4 S cm-1 at 300 K, which is one order of magnitude higher than the ionic conductivity for LTP and LATP at x = 0.2. This value is in good agreement with our experimental value (2.5 × 10-4 S cm-1 at 300 K) and the literature values. The composition-dependent ionic conductivity of LATP was successfully demonstrated using the ReaxFF reactive force field, verifying the applicability of the LATP force field for the understanding of Li diffusion and the design of highly Li ion conductive solid electrolytes. Furthermore, our results also demonstrate the feasibility of the MC/MD method in modeling LATP configuration, and provide compelling evidence for the solid solution sensitivity on ionic conductivity.