Fast docking using the CHARMM force field with EADock DSS

Tenebrio molitor antifreeze protein (TmAFP) was simulated with growing ice-water interfaces at a realistic melting temperature using TIP4P/Ice water model. To test compatibility of protein force fields (FFs) with TIP4P/Ice water, CHARMM, AMBER, and OPLS FFs were applied. CHARMM and AMBER FFs predict more β-sheet structure and lower diffusivity of TmAFP at the ice-water interface than does OPLS FF, indicating that β-sheet structure is important for the TmAFP-interface binding and antifreeze activity. In particular, CHARMM FF more clearly distinguishes the strengths of hydrogen bonds in the ice-binding and non-ice-binding sites of TmAFP than do other FFs, in agreement with experiments, implying that CHARMM FF can be a reasonable choice to simulate proteins with TIP4P/Ice water. Simulations of mutated TmAFPs show that for the same density of Thr residues, continuous arrangement of Thr with the distance of 0.4~0.6 nm induces the higher extent of antifreeze activity than does intermittent arrangement of Thr with larger distances. These findings suggest the choice of CHARMM FF for AFP-TIP4P/Ice simulations and help explain the relationship between Thr-residue arrangement and antifreeze activity.

Structural Equilibrium of DNA Represented with Different Force Fields

Physics-based force fields are the backbone of molecular dynamics simulations. In recent years, significant progress has been made in the assessment and improvement of commonly used force fields for describing conformational dynamics of folded proteins. However, the accuracy for the unfolded states remains unclear. The latter is however important for detailed studies of protein folding pathways, conformational transitions involving unfolded states, and dynamics of intrinsically disordered proteins. In this work, we compare the three commonly used force fields, AMBER ff99SB-ILDN, CHARMM22/CMAP, and CHARMM36, for modeling the natively unfolded fragment peptides, NTL9(1-22) and NTL9(6-17), using explicit-solvent replica-exchange molecular dynamics simulations. All three simulations show that NTL9(6-17) is completely unstructured, while NTL9(1-22) transiently samples various β-hairpin states, reminiscent of the first β-hairpin in the structure of the intact NTL9 protein. The radius of gyration of the two peptides is force field independent but likely underestimated due to the current deficiency of additive force fields. Compared to the CHARMM force fields, ff99SB-ILDN gives slightly higher β-sheet propensity and more native-like residual structures for NTL9(1-22), which may be attributed to its known β preference. Surprisingly, only two sequence-local pairs of charged residues make appreciable ionic contacts in the simulations of NTL9(1-22), which are sampled slightly more by the CHARMM force fields. Taken together, these data suggest that the current CHARMM and AMBER force fields are globally in agreement in modeling the unfolded states corresponding to β-sheet in the folded structure, while differing in details such as the native-likeness of the residual structures and interactions.

The drug discovery process has been deeply transformed recently by the use of computational ligand-based or structure-based methods, helping the lead compounds identification and optimization, and finally the delivery of new drug candidates more quickly and at lower cost. Structure-based computational methods for drug discovery mainly involve ligand-protein docking and rapid binding free energy estimation, both of which require force field parameterization for many drug candidates. Here, we present a fast force field generation tool, called SwissParam, able to generate, for arbitrary small organic molecule, topologies, and parameters based on the Merck molecular force field, but in a functional form that is compatible with the CHARMM force field. Output files can be used with CHARMM or GROMACS. The topologies and parameters generated by SwissParam are used by the docking software EADock2 and EADock DSS to describe the small molecules to be docked, whereas the protein is described by the CHARMM force field, and allow them to reach success rates ranging from 56 to 78%. We have also developed a rapid binding free energy estimation approach, using SwissParam for ligands and CHARMM22/27 for proteins, which requires only a short minimization to reproduce the experimental binding free energy of 214 ligand-protein complexes involving 62 different proteins, with a standard error of 2.0 kcal mol(-1), and a correlation coefficient of 0.74. Together, these results demonstrate the relevance of using SwissParam topologies and parameters to describe small organic molecules in computer-aided drug design applications, together with a CHARMM22/27 description of the target protein. SwissParam is available free of charge for academic users at www.swissparam.ch.

Here we show the newest implementation of Flexible Artificial Intelligence Docking (FlexAID) allowing its scoring function to consider the conformational entropy of ligand and biomolecules complexes. The higher accuracy of FlexAID1 on complex cases, the addition of novel features, i.e. the conformational entropy, its accessibility and its easy-to-use graphical user interface place FlexAID in an interesting position to tackle biologically and pharmacologically relevant situations currently ignored by other methods. FlexAID 1 is available as a command-line pre-compiled executable (available at http://bcb.med.usherbrooke.ca/flexaid for Windows, macOS & Linux) or through the NRGsuite, a PyMOL integrated user interface allowing the user to use FlexAID in an intuitive manner with real time visualization. Both the NRGsuite 2 and FlexAID 1 are distributed as open-source software. Introduction A major goal of molecular docking is to predict the experimentally observed binding mode between a biomolecule, i.e. a polymer of amino or nucleic acids, and a ligand— e.g. small molecules, peptides or nucleic acids. This computational method is used to study the structure of the molecular interactions involved in cell’s essential biological functions. Actual molecular docking methods are developed to evaluate the molecular interactions of a single conformation, a single pose at a time, and they are trained to estimate the enthalpic fraction of the binding free energy. Consequently, most current molecular docking methods fail to efficiently model the entropic contributions, especially those of conformational nature, who are fundamental in molecular recognition events. Our research group develops FlexAID 1 , an accessible and competitive ligand and biomolecule molecular docking software whose focus is on molecular flexibility. Here we introduce FlexAID’s newest feature that allows its scoring function to estimate the conformational entropy by redefining the static binding mode usually predicted in molecular docking into a dynamic collection of similar poses evaluated altogether. Methods We implemented the new scoring function in FlexAID, allowing its genetic algorithm to select less favourable, but frequently observed binding modes, with a probability following a Boltzmann distribution. This implementation allows FlexAID to consider conformational entropy of the complexes during the molecular docking simulation. The core of the implementation of the conformational entropy in FlexAID resides within an unsupervised and density-based molecular classification algorithm charged to group similar poses together into binding modes, thus redefining a binding mode as a dynamic collection of poses scored altogether as it an be seen in Figure 1. Results & Conclusions We present the impact of FlexAID’s newest feature on its accuracy in binding mode prediction using three increasingly complex scenarios: the Astex Diverse Set 4 , the Astex Non Native Set 5 and the HAP2 6 dataset. We show that FlexAID outperforms other open-source molecular docking methods when molecular flexibility is crucial. Furthermore, FlexAID now outputs multiple conformations per binding mode, a novelty that allows the user to visualize the dynamics of the complex studied. We believe that its higher accuracy in complex scenarios, the addition of novel features, e.g. the conformational entropy, its accessibility and its easy-to-use graphical user interface, the NRGsuite 2 , place FlexAID in an interesting position to tackle biologically and pharmacologically relevant situations currently ignored by other molecular docking methods.

In this work, we present atomic parameters of perfluoroalkanes for use within the CHARMM force field. Perfluorinated alkanes represent a special class of molecules. On the one hand, they are considerably more hydrophobic than lipids, but on the other hand, they are not lipophilic either. Instead, they represent an independent class of philicity, enabling a whole portfolio of applications within both materials science and biochemistry. We performed a thorough parametrization of all bonded and nonbonded parameters with a particular focus on van der Waals parameters. Here, the general framework of the CHARMM and CGenFF force fields has been followed. The van der Waals parameters have been fitted to experimental densities over a wide range of temperatures and pressures. This newly parametrized class of molecules will open the gate for a variety of simulations of biologically relevant systems within the CHARMM force field. A particular perspective for the present work is the influence of polyphilic transmembrane molecules on membrane properties, aggregation phenomena, and transmembrane channels.

Dissolution study of active pharmaceutical ingredients using molecular dynamics simulations with classical force fields

From my PhD research, under Amber99-Chen-Garcia force field, umbrella sampling and MD simulation helped us found that the 2’-5’ phosphodiester RNA backbone linkage modification on in vitro selected neomycin RNA aptamer can transfer it to ‘amikacin RNA aptamer’ by changing the conformation and backbone rigidity of the RNA binding pocket, which further influenced the binding mode and binding pathway of amikacin. Furthermore, under the same force field, metadynamics was found more powerful for simulating the binding pathway, binding mechanism and predicting binding modes of small molecule ligand recognizing its RNA target by filling up all energy basins on free energy surface of the binding pathway with extra Gaussian potential, which was well-trained on five small molecule/RNA binding pocket systems. This metadynamics prediction method training study also discovered a theory of ‘conformational-selection model of small molecule binding with different degree of induced-fit effect on RNA target based on the rigidity of small molecule ligand’ for the binding model of small molecule ligand/RNA binding pocket recognition system. Furthermore, this well-trained metadynamics simulation prediction method was firstly applied to help us found that 2’-5’ linkage modifications on FMN RNA aptamer narrowed the width of the frontside opening of RNA binding pocket, at the same time, the backside opening became wider, which changed the binding pathway and binding mode of FMN to ‘upside down’. Secondly, this well-trained metadynamics simulation prediction method helped us refining the NMR binding mode of small molecule splicing modifier SMN-C5/RNA duplex and further predicted a new mechanism of splicing modifier SMN-C5 promoting SMN2 exon 7 inclusion for treating spinal muscular atrophy (SMA). In sum, Amber99-Chen-Garcia force field showed its power and potential of capturing RNA-small molecule ligand interactions correctly. Umbrella sampling showed its power and potential of simulating unbinding events and calculating ΔG bind/unbind for RNA-small molecule ligand complex correctly. With a good force field of Amber99-Chen-Garcia and good CVs selection, metadynamics showed its power and potential of predicting binding modes, binding pathway and binding mechanism of RNA-small molecule ligand recognition correctly.

Part of early stage drug discovery involves determining how molecules may bind to the target protein. Through understanding where and how molecules bind, chemists can begin to build ideas on how to design improvements to increase binding affinities. In this retrospective study, we compare how computational approaches like docking, molecular dynamics (MD) simulations, and a non-equilibrium candidate Monte Carlo (NCMC) based method (NCMC+MD) perform in predicting binding modes for a set of 12 fragment-like molecules which bind to soluble epoxide hydrolase. We evaluate each method's effectiveness in identifying the dominant binding mode and finding any additional binding modes (if any). Then, we compare our predicted binding modes to experimentally obtained X-ray crystal structures.We dock each of the 12 small molecules into the apo-protein crystal structure and then run simulations up to 1 microsecond each. Small and fragment-like molecules likely have smaller energy barriers separating different binding modes by virtue of relatively fewer and weaker interactions relative to drug-like molecules, and thus likely undergo more rapid binding mode transitions. We expect, thus, to see more rapid transitions betweeen binding modes in our study. Following this, we build Markov State Models (MSM) to define our stable ligand binding modes. We investigate if adequate sampling of ligand binding modes and transitions between them can occur at the microsecond timescale using traditional MD or a hybrid NCMC+MD simulation approach. Our findings suggest that even with small fragment-like molecules, we fail to sample all the crystallographic binding modes using microsecond MD simulations, but using NCMC+MD we have better success in sampling the crystal structure while obtaining the correct populations.

Part of early stage drug discovery involves determining how molecules may bind to the target protein. Through understanding where and how molecules bind, chemists can begin to build ideas on how to design improvements to increase binding affinities. In this retrospective study, we compare how computational approaches like docking, molecular dynamics (MD) simulations, and a non-equilibrium candidate Monte Carlo (NCMC)-based method (NCMC + MD) perform in predicting binding modes for a set of 12 fragment-like molecules, which bind to soluble epoxide hydrolase. We evaluate each method's effectiveness in identifying the dominant binding mode and finding additional binding modes (if any). Then, we compare our predicted binding modes to experimentally obtained X-ray crystal structures. We dock each of the 12 small molecules into the apo-protein crystal structure and then run simulations up to 1 μs each. Small and fragment-like molecules likely have smaller energy barriers separating different binding modes by virtue of relatively fewer and weaker interactions relative to drug-like molecules and thus likely undergo more rapid binding mode transitions. We expect, thus, to see more rapid transitions between binding modes in our study. Following this, we build Markov State Models to define our stable ligand binding modes. We investigate if adequate sampling of ligand binding modes and transitions between them can occur at the microsecond timescale using traditional MD or a hybrid NCMC+MD simulation approach. Our findings suggest that even with small fragment-like molecules, we fail to sample all the crystallographic binding modes using microsecond MD simulations, but using NCMC+MD, we have better success in sampling the crystal structure while obtaining the correct populations.

New generation of docking programs: Supercomputer validation of force fields and quantum-chemical methods for docking

Studies on light-harvesting (LH) systems have attracted much attention after the finding of long-lived quantum coherences in the exciton dynamics of the Fenna-Matthews-Olson (FMO) complex. In this complex, excitation energy transfer occurs between the bacteriochlorophyll a (BChl a) pigments. Two quantum mechanics/molecular mechanics (QM/MM) studies, each with a different force-field and quantum chemistry approach, reported different excitation energy distributions for the FMO complex. To understand the reasons for these differences in the predicted excitation energies, we have carried out a comparative study between the simulations using the CHARMM and AMBER force field and the Zerner intermediate neglect of differential orbital (ZINDO)/S and time-dependent density functional theory (TDDFT) quantum chemistry methods. The calculations using the CHARMM force field together with ZINDO/S or TDDFT always show a wider spread in the energy distribution compared to those using the AMBER force field. High- or low-energy tails in these energy distributions result in larger values for the spectral density at low frequencies. A detailed study on individual BChl a molecules in solution shows that without the environment, the density of states is the same for both force field sets. Including the environmental point charges, however, the excitation energy distribution gets broader and, depending on the applied methods, also asymmetric. The excitation energy distribution predicted using TDDFT together with the AMBER force field shows a symmetric, Gaussian-like distribution.

Some of the main challenges faced in drug discovery are pocket flexibility and binding mode prediction. In this work, we explored the aromatic cage flexibility of the histone methyllysine reader protein Spindlin1 and its impact on binding mode prediction by means of in silico approaches. We first investigated the Spindlin1 aromatic cage plasticity by analyzing the available crystal structures and through molecular dynamic simulations. Then we assessed the ability of rigid docking and flexible docking to rightly reproduce the binding mode of a known ligand into Spindlin1, as an example of a reader protein displaying flexibility in the binding pocket. The ability of induced fit docking was further probed to test if the right ligand binding mode could be obtained through flexible docking regardless of the initial protein conformation. Finally, the stability of generated docking poses was verified by molecular dynamic simulations. Accurate binding mode prediction was obtained showing that the herein reported approach is a highly promising combination of in silico methods able to rightly predict the binding mode of small molecule ligands in flexible binding pockets, such as those observed in some reader proteins.

A systematic molecular dynamics simulation study of temperature dependent bilayer structural properties

For the engineering and process design of chemical and pharmaceutical plants, the knowledge of thermophysical properties is essential. Here, glass transition temperature (Tg), curves of heat capacity (Cp), isotropic thermal expansion (∝p), and isothermal compressibility (βT) are computed for amorphous/paracrystalline (Am‐Par) structures of cellulose over a wide range of temperature (380–680 K) using molecular dynamics with the CHARMM36 (C36) force field (FF). The fluctuation method under the NPT ensemble is used to calculate Cp, ∝p, and βT, whereas Tg is computed by monitoring specific volume versus temperature. Here, the fluctuation method is used with a quantum mechanical correction term for the calculation of Cp. Results of Cp, ∝p, and βT values at 298 K using extrapolation from these curves are also obtained. The thermophysical properties values from the simulations are compared with experimental data for cellulose with different degree of crystallinity and with those obtained by prominent FFs suggested for cellulose, such as GLYCAM06 and COMPASS. The findings reveal that ∝p, βT, and Tg are somewhat better reproduced than Cp with C36 over the studied temperature range. From this study, it is inferred that, for accurate modeling of heat capacity of pure Am‐Par celluloses with large fragments of glucose, the C36 FF needs re‐parameterization.