CHARMM Program Research Articles

Exploring the Limits of the Generalized CHARMM and AMBER Force Fields through Predictions of Hydration Free Energy of Small Molecules.

Accurate force field parameters, potential energy functions, and receptor-ligand models are essential for modeling the solvation and binding of drug-like molecules to a receptor. A large and ever-growing chemical space of medicinally relevant scaffolds has also required these factors, especially force field parameters, to be highly transferable. Generalized force fields such as the CHARMM General Force Field (CGenFF) and the generalized AMBER force field (GAFF) have accomplished this feat along with other contemporaneous ones like OPLS. Here, we analyze the limits in the parametrization of drug-like small molecules by CGenFF and GAFF in terms of the various functional groups represented within them. Specifically, we link the presence of specific functional groups to the error in the absolute hydration free energy of over 600 small molecules, predicted by alchemical free energy methods implemented in the CHARMM program. Our investigation reveals that molecules with (i) a nitro group in CGenFF and GAFF are, respectively, over- or undersolubilized in aqueous medium, (ii) amine groups are undersolubilized more so in CGenFF than in GAFF, and (iii) carboxyl groups are more oversolubilized in GAFF than in CGenFF. We present our analyses of the potential factors underlying these trends. We also showcase the use of a machine-learning-based approach combined with the SHapley Additive exPlanations framework to attribute these trends to specific functional groups, which can be easily adopted to explore the limits of other general force fields.

Computation of conformational transitions and activation of Src tyrosine kinase.

A protein molecule can adapt different structural forms depending on solution conditions or ligation state. The stable forms can be defined experimentally, but the transition pathways and metastable states between them cannot. A computational method to explore equilibrium conformational transition pathways is the adaptive bias path optimization (ABPO) algorithm. ABPO, implemented in the CHARMM program, differs from other path methods by using enhanced sampling along a path between two known states without restraining the system to the reduced coordinates of the evolving path. ABPO converges to a physically reasonable path even in the case of a highly rugged energy landscape and is amenable to assessing the power of the reduced variables chosen to explore the transition. ABPO was applied to three cases where a conformational transition is localized to one part of the polypeptide chain. With proper selection of the reduced coordinate set, sampling without restraints converged to a physically reasonable path. We also applied ABPO to a more complex transition of Src tyrosine kinase catalytic domain that undergoes a more complex transition between enzymatically down‑regulated and activated conformations. The results using coarse-grained and all‑atom models identified a switched electrostatic network (SEN) that could confer an entropic advantage to reduce energetic barriers by limiting the search space of the activation loop. The richness in atomic detail provided by the ABPO path reveals how the SEN couples the dynamics of the key αC‑helix and the activation loop, which provided a new rationale for the highly conserved HRD sequence of kinases. Further, the results related to energetics of the rotation of the αC‑helix provide a rationale for the interactions observed in a number of kinase regulatory complexes, some for which the basis for regulation from inspection of the structure along is puzzling.

Two physics-based models for pH-dependent calculations of protein solubility.

When engineering a protein for its biological function, many physicochemical properties are also optimized throughout the engineering process, and the protein's solubility is among the most important properties to consider. Here, we report two novel computational methods to calculate the pH-dependent protein solubility, and to rank the solubility of mutants. The first is an empirical method developed for fast ranking of the solubility of a large number of mutants of a protein. It takes into account electrostatic solvation energy term calculated using Generalized Born approximation, hydrophobic patches, protein charge, and charge asymmetry, as well as the changes of protein stability upon mutation. This method has been tested on over 100 mutations for 17 globular proteins, as well as on 44 variants of five different antibodies. The prediction rate is over 80%. The antibody tests showed a Pearson correlation coefficient, R, with experimental data from .83 to .91. The second method is based on a novel, completely force-field-based approach using CHARMm program modules to calculate the binding energy of the protein to a part of the crystal lattice, generated from X-ray structure. The method predicted with very high accuracy the solubility of Ribonuclease SA and its 3K and 5K mutants as a function of pH without any parameter adjustments of the existing BIOVIA Discovery Studio binding affinity model. Our methods can be used for rapid screening of large numbers of design candidates based on solubility, and to guide the design of solution conditions for antibody formulation.

Polarizable molecular dynamics simulations of ionic liquids: Influence of temperature control.

Ionic liquids are an interesting class of soft matter with viscosities of one or two orders of magnitude higher than that of water. Unfortunately, classical, non-polarizable molecular dynamics (MD) simulations of ionic liquids result in too slow dynamics and demonstrate the need for explicit inclusion of polarizability. The inclusion of polarizability, here via the Drude oscillator model, requires amendments to the employed thermostat, where we consider a dual Nosé-Hoover thermostat, as well as a dual Langevin thermostat. We investigate the effects of the choice of a thermostat and the underlying parameters such as the masses and force constants of the Drude particles on static and dynamic properties of ionic liquids. Here, we show that Langevin thermostats are not suitable for investigating the dynamics of ionic liquids. Since polarizable MD simulations are associated with high computational costs, we employed a self-developed graphics processing unit enhanced code within the MD program CHARMM to keep the overall computational effort reasonable.

IN SILICO AND IN VIVO PHARMACOLOGICAL STUDIES OF CLOZAPINE AND D-AMINO ACID OXIDASE INHIBITOR FOR COGNITIVE ENHANCEMENT

Objective: D-amino acid oxidase inhibitors (DAAOIs) are of particular focus for cognition study. Atypical antipsychotics are known DAAO inhibitors. The present examination was done to check out the binding affinity of atypical antipsychotics by docking toward the DAAO protein; in conclusion, the picked antipsychotic drug was checked for their cognition enhancing activity with scopolamine-induced amnesia.Methods: The crystal structure of DAAO was obtained from Protein Data Bank, the energy minimization was performed with CHARMM program, then active site prediction was made out using Ramachandran plot, and finally, docking examination was finished using Autodock 4.2 tool. For in vivo study, the mice were divided into three groups. Group I - vehicle (Saline) treated, Group II – saline +scopolamine (1 mg/kg, intraperitoneal [i.p]) treated, and Group III - clozapine (20 mg/kg, i.p) + scopolamine (1 mg/kg, i.p).Results: The Autodock examination shows significant binding affinity of - 5.22 for brexpiprazole and least or positive binding affinity of +1 for iloperidone. Clozapine with binding energy of - 2.87 was decided for completing the in vivo cognition study. The in vivo shows up that clozapine (20 mg/kg, i.p) exhibits a change in the impairment of spatial memory.Conclusion: The results recommend that the clozapine produces cognitive enhancement through both DAAOI and antipsychotic action. Clozapine has cognitive improvement potential, favoring its usage in reducing toxic impacts of scopolamine

Classical Molecular Dynamics with Mobile Protons.

An important limitation of standard classical molecular dynamics simulations is the inability to make or break chemical bonds. This restricts severely our ability to study processes that involve even the simplest of chemical reactions, the transfer of a proton. Existing approaches for allowing proton transfer in the context of classical mechanics are rather cumbersome and have not achieved widespread use and routine status. Here we reconsider the combination of molecular dynamics with periodic stochastic proton hops. To ensure computational efficiency, we propose a non-Boltzmann acceptance criterion that is heuristically adjusted to maintain the correct or desirable thermodynamic equilibria between different protonation states and proton transfer rates. Parameters are proposed for hydronium, Asp, Glu, and His. The algorithm is implemented in the program CHARMM and tested on proton diffusion in bulk water and carbon nanotubes and on proton conductance in the gramicidin A channel. Using hopping parameters determined from proton diffusion in bulk water, the model reproduces the enhanced proton diffusivity in carbon nanotubes and gives a reasonable estimate of the proton conductance in gramicidin A.

GAFit: A general-purpose, user-friendly program for fitting potential energy surfaces

We have developed a software package based on a genetic algorithm that fits an analytic function to a given set of data points. The code, called GAFit, was also interfaced with the CHARMM and MOPAC programs in order to facilitate force field parameterizations and fittings of specific reaction parameters (SRP) for semiempirical Hamiltonians. The present tool may be applied to a wide range of fitting problems, though it has been especially designed to significantly reduce the hard work involved in the development of potential energy surfaces for complex systems. For this purpose, it has been equipped with several programs to help the user in the preparation of the input files. We showcase the application of the computational tool to several chemical-relevant problems: force-field parameterization, with emphasis on nonbonded energy terms or intermolecular potentials, derivation of SRP for semiempirical Hamiltonians, and fittings of generic analytical functions. Program summaryProgram title : GAFitProgram Files doi : http://dx.doi.org/10.17632/9gy6bjcwk5.1Licensing provisions : GNU General Public License 3 (GPL)Programming language : Fortran 90, C, Perl and JavaNature of problem : Potential energy surfaces (PESs) have a key role in reaction dynamics and molecular dynamics. Chemical dynamics simulations can be performed by using analytical PESs or by direct dynamics, that is, by computing energies and forces “on-the-fly” by molecular structure calculations. The development of analytical PESs is, quite often, a very complex and tedious task. There are several programs in the literature that may help to develop PESs but, as far as we know, none of them show, at the same time, great generality and flexibility.Solution method : GAFit was designed to help users develop analytical potential energy functions. We made special efforts to write a code that exhibits the three features mentioned above. For this reason, we have also interfaced GAFit with the CHARMM and MOPAC programs. The former is one of the popular packages for molecular dynamics and the latter is a semiempirical quantum chemistry program that may be conveniently used for direct dynamics simulations. Because, in general, the fitting of PESs involves a simultaneous optimization of many nonlinear parameters, we have selected a genetic algorithm as the driver for the parameterizations. As shown in the paper, GAFit may be applied to other fitting problems.

Revisiting sulfur H-bonds in proteins: The example of peroxiredoxin AhpE.

In many established methods, identification of hydrogen bonds (H-bonds) is primarily based on pairwise comparison of distances between atoms. These methods often give rise to systematic errors when sulfur is involved. A more accurate method is the non-covalent interaction index, which determines the strength of the H-bonds based on the associated electron density and its gradient. We applied the NCI index on the active site of a single-cysteine peroxiredoxin. We found a different sulfur hydrogen-bonding network to that typically found by established methods, and we propose a more accurate equation for determining sulfur H-bonds based on geometrical criteria. This new algorithm will be implemented in the next release of the widely-used CHARMM program (version 41b), and will be particularly useful for analyzing water molecule-mediated H-bonds involving different atom types. Furthermore, based on the identification of the weakest sulfur-water H-bond, the location of hydrogen peroxide for the nucleophilic attack by the cysteine sulfur can be predicted. In general, current methods to determine H-bonds will need to be reevaluated, thereby leading to better understanding of the catalytic mechanisms in which sulfur chemistry is involved.

Protein-Protein Docking Using EMAP in CHARMM and Support Vector Machine: Application to Ab/Ag Complexes.

In this work, we have (i) evaluated the ability of the EMAP method implemented in the CHARMM program to generate the correct conformation of Ab/Ag complex structures and (ii) developed a support vector machine (SVM) classifier to detect native conformations among the thousands of refined Ab/Ag configurations using the individual components of the binding free energy based on a thermodynamic cycle as input features in training the SVM. Tests on 24 Ab/Ag complexes from the protein-protein docking benchmark version 3.0 showed that based on CAPRI evaluation criteria, EMAP could generate medium-quality native conformations in each case. Furthermore, the SVM classifier could rank medium/high-quality native conformations mostly in the top six among the thousands of refined Ab/Ag configurations. Thus, Ab-Ag docking can be performed using different levels of protein representations, from grid-based (EMAP) to polar hydrogen (united-atom) to all-atom representation within the same program. The scripts used and the trained SVM are available at the www.charmm.org forum script repository.

Model structure of a full-length anti-IgE complexed with IgE (P3384)

Abstract Modeling antibody/antigen (Ab/Ag) complexes remains a challenging problem and is important for understanding how proteins associate, which in turn can aid drug design. Recently, we have implemented a method (EMAP) for producing the conformations of Ab/Ag complexes in the widely used molecular dynamics program CHARMM. Here, we apply this method to model the structure of IgE in complex with omalizumab, a humanized IgG1 antibody that has been approved in many countries for treating severe, persistent allergic asthma. First, the EMAP method was calibrated by using it to reproduce 24 nonredundant Ab/Ag complexes whose bound and unbound structures have been solved. It could yield near-native conformations for 22 of the 24 complexes, whereas other docking methods such as ZDOCK yield less near-native conformations. Having validated the EMAP method, it was used to dock full-length models of omalizumab and IgE together. The resulting IgE/omalizumab complex structure provides insights into the relative importance of the residues involved in binding. It also provides a plausible explanation for why omalizumab and IgE do not form large immune complexes.

Unfolding Ubiquitin by force: water mediated H-bond destabilization

Using the pull and wait (PNW) simulation protocol at 300 K, we studied the unfolding of a ubiquitin molecule by force. PNW was implemented in the CHARMM program using an integration time step of 1 fs and a uniform dielectric constant of 1. The ubiquitin molecule, initially solvated, was put under mechanical stress, exerting forces from different directions. The rupture of five hydrogen bonds between parallel strands β1 and β5 takes place during the extension from 13 to 15 A, defines a mechanical barrier for unfolding and dominates the point of maximum unfolding force. The simulations described here show that given adequate time, a small applied force can destabilize those five H-bonds relative to the bonds that can be created to water molecules; allowing the formation of stable H-bonds between a single water molecule and the donor and acceptor groups of the interstrand H-bonds. Thus, simulations run with PNW show that the force is not responsible for ripping apart the backbone H-bonds; it merely destabilizes them making them less stable than the H-bonds they can make with water. Additional simulations show that the force necessary to destabilize the H-bonds and allow them to be replaced by H-bonds to water molecules depends strongly on the pulling direction. By using a simulation protocol that allows equilibration at each extension we have been able to observe the details of the events leading to the unfolding of ubiquitin by mechanical force.

Maintain rigid structures in Verlet based Cartesian molecular dynamics simulations

An algorithm is presented to maintain rigid structures in Verlet based cartesian molecular dynamics (MD) simulations. After each unconstrained MD step, the coordinates of selected particles are corrected to maintain rigid structures through an iterative procedure of rotation matrix computation. This algorithm, named as SHAPE and implemented in CHARMM program suite, avoids the calculations of Lagrange multipliers, so that the complexity of computation does not increase with the number of particles in a rigid structure. The implementation of this algorithm does not require significant modification of propagation integrator, and can be plugged into any cartesian based MD integration scheme. A unique feature of the SHAPE method is that it is interchangeable with SHAKE for any object that can be constrained as a rigid structure using multiple SHAKE constraints. Unlike SHAKE, the SHAPE method can be applied to large linear (with three or more centers) and planar (with four or more centers) rigid bodies. Numerical tests with four model systems including two proteins demonstrate that the accuracy and reliability of the SHAPE method are comparable to the SHAKE method, but with much more applicability and efficiency.

Development of a tuned interfacial force field parameter set for the simulation of protein adsorption to silica glass.

Adsorption free energies for eight host-guest peptides (TGTG-X-GTGT, with X = N, D, G, K, F, T, W, and V) on two different silica surfaces [quartz (100) and silica glass] were calculated using umbrella sampling and replica exchange molecular dynamics and compared with experimental values determined by atomic force microscopy. Using the CHARMM force field, adsorption free energies were found to be overestimated (i.e., too strongly adsorbing) by about 5-9 kcal/mol compared to the experimental data for both types of silica surfaces. Peptide adsorption behavior for the silica glass surface was then adjusted using a modified version of the CHARMM program, which we call dual force-field CHARMM, which allows separate sets of nonbonded parameters (i.e., partial charge and Lennard-Jones parameters) to be used to represent intra-phase and inter-phase interactions within a given molecular system. Using this program, interfacial force field (IFF) parameters for the peptide-silica glass systems were corrected to obtain adsorption free energies within about 0.5 kcal/mol of their respective experimental values, while IFF tuning for the quartz (100) surface remains for future work. The tuned IFF parameter set for silica glass will subsequently be used for simulations of protein adsorption behavior on silica glass with greater confidence in the balance between relative adsorption affinities of amino acid residues and the aqueous solution for the silica glass surface.

Glycan reader: Automated sugar identification and simulation preparation for carbohydrates and glycoproteins

Understanding how glycosylation affects protein structure, dynamics, and function is an emerging and challenging problem in biology. As a first step toward glycan modeling in the context of structural glycobiology, we have developed Glycan Reader and integrated it into the CHARMM-GUI, http://www.charmm-gui.org/input/glycan. Glycan Reader greatly simplifies the reading of PDB structure files containing glycans through (i) detection of carbohydrate molecules, (ii) automatic annotation of carbohydrates based on their three-dimensional structures, (iii) recognition of glycosidic linkages between carbohydrates as well as N-/O-glycosidic linkages to proteins, and (iv) generation of inputs for the biomolecular simulation program CHARMM with the proper glycosidic linkage setup. In addition, Glycan Reader is linked to other functional modules in CHARMM-GUI, allowing users to easily generate carbohydrate or glycoprotein molecular simulation systems in solution or membrane environments and visualize the electrostatic potential on glycoprotein surfaces. These tools are useful for studying the impact of glycosylation on protein structure and dynamics.

The distributed diagonal force decomposition method for parallelizing molecular dynamics simulations

Parallelization is an effective way to reduce the computational time needed for molecular dynamics simulations. We describe a new parallelization method, the distributed-diagonal force decomposition method, with which we extend and improve the existing force decomposition methods. Our new method requires less data communication during molecular dynamics simulations than replicated data and current force decomposition methods, increasing the parallel efficiency. It also dynamically load-balances the processors' computational load throughout the simulation. The method is readily implemented in existing molecular dynamics codes and it has been incorporated into the CHARMM program, allowing its immediate use in conjunction with the many molecular dynamics simulation techniques that are already present in the program. We also present the design of the Force Decomposition Machine, a cluster of personal computers and networks that is tailored to running molecular dynamics simulations using the distributed diagonal force decomposition method. The design is expandable and provides various degrees of fault resilience. This approach is easily adaptable to computers with Graphics Processing Units because it is independent of the processor type being used.

A versatile method for systematic conformational searches: application to CheY.

A novel molecular structure prediction method, the Z Method, is described. It provides a versatile platform for the development and use of systematic, grid-based conformational search protocols, in which statistical information (i.e., rotamers) can also be included. The Z Method generates trial structures by applying many changes of the same type to a single starting structure, thereby sampling the conformation space in an unbiased way. The method, implemented in the CHARMM program as the Z Module, is applied here to an illustrative model problem in which rigid, systematic searches are performed in a 36-dimensional conformational space that describes the relative positions of the 10 secondary structural elements of the protein CheY. A polar hydrogen representation with an implicit solvation term (EEF1) is used to evaluate successively larger fragments of the protein generated in a hierarchical build-up procedure. After a final refinement stage, and a total computational time of about two-and-a-half CPU days on AMD Opteron processors, the prediction is within 1.56 Å of the native structure. The errors in the predicted backbone dihedral angles are found to approximately cancel. Monte Carlo and simulated annealing trials on the same or smaller versions of the problem, using the same atomic model and energy terms, are shown to result in less accurate predictions. Although the problem solved here is a limited one, the findings illustrate the utility of systematic searches with atom-based models for macromolecular structure prediction and the importance of unbiased sampling in structure prediction methods.

Automated Umbrella Sampling Simulations for the Calculation of Multidimensional Potentials of Mean Force

The potential of mean force (PMF) describing the conformational changes of biomolecules is a central quantity that determines the function of biomolecular systems. Calculating a multidimensional PMF (free energy landscape) is a time consuming process. Each additional reaction coordinate drastically increases the required simulation time for conformation sampling, making calculations in four or more dimensions practically impossible. However, in most cases, only a small fraction of such multidimensional space is energetically relevant. PMF calculations in high dimensionality could thus be achieved if one could effectively focus the simulation effort on those region of high interest. We have developed a method based on umbrella sampling (US) that determines, using a feedback mechanism, which regions of the multidimensional space is worth exploring. The first aim of such application is to manage the creation and analysis of sampling windows for the calculation of PMFs involving up to four reaction coordinates. While this approach could in some cases be used to find the minimum free energy pathway underlying large conformational changes, it is not as general as other approaches that were specifically developed in that purpose. However, contrary to most other approaches, our method allows for the simultaneous characterization of several pathways, and not only the most probable one. The current implementation consists in a C++ application with a python interface that manages MD simulations performed with the biomolecular simulation program CHARMM. The feedback mechanism and final PMF calculation involves the unbiasing of simulations using the weighted-histogram analysis methods (WHAM).

Study of the Cooperativity of Calcium Binding in Calbindin D9k using 2D Replica-Exchange Umbrella Sampling

Department of Biochemistry and Molecular Biology, The University of Chicago, IL 60637 Argonne National Laboratory, Argonne, IL 60439 Calbindin D9k, a small protein possessing a pair of EF-hands that binds two calcium ions in a cooperative fashion, is chosen as a model system for studying the molecular mechanism of cooperativity.1 Binding is characterized in terms of a potential of mean force (PMF) as a function two variables: the distance r between the ion, and the binding pocket, and the root-mean-square deviation (RMSD) of the conformation of the EF-hand relative to its ion-bound structure. The PMF is calculated using a novel two-dimensional replica-exchange molecular dynamics (MD) umbrella sampling scheme, which is developed and implemented in the program CHARMM to increase the configuration space sampling. Using 2048 replicas on Blue Gene/P, the 2D-PMF/REMD calculation for the binding the second calcium ion converges within 800 ps, with an exchange probability of 30%. In contrast, standard umbrella sampling PMF/MD simulations without replica-exchange did not converge, even after 1.2 ns. The absolute binding free energy of a second ion to a singly occupied calbindin calculated from the 2D-PMF following the statistical mechanical formulation of noncovalent association2 is −9.4 kcal/mol, in excellent agreement with the experimental value3. The 2D-PMF/REMD simulations will be extended to provide important information about the molecular basis of calcium binding cooperativity. Footnotes 1 Marchand S. and Roux B. Protein Struct. Funct. Genet. 1998, 33, 265-284. 2 Woo H. and Roux B., PNAS, 2005, 19, 6825. 3 Linse S. et al. Biochemistry 1991, 30, 154.

Computing ensembles of transitions from stable states: Dynamic importance sampling.

There is an increasing dataset of solved biomolecular structures in more than one conformation and increasing evidence that large-scale conformational change is critical for biomolecular function. In this article, we present our implementation of a dynamic importance sampling (DIMS) algorithm that is directed toward improving our understanding of important intermediate states between experimentally defined starting and ending points. This complements traditional molecular dynamics methods where most of the sampling time is spent in the stable free energy wells defined by these initial and final points. As such, the algorithm creates a candidate set of transitions that provide insights for the much slower and probably most important, functionally relevant degrees of freedom. The method is implemented in the program CHARMM and is tested on six systems of growing size and complexity. These systems, the folding of Protein A and of Protein G, the conformational changes in the calcium sensor S100A6, the glucose-galactose-binding protein, maltodextrin, and lactoferrin, are also compared against other approaches that have been suggested in the literature. The results suggest good sampling on a diverse set of intermediates for all six systems with an ability to control the bias and thus to sample distributions of trajectories for the analysis of intermediate states.

Electron Microscopy and Persistence Length Analysis of Semi-Rigid Smooth Muscle Tropomyosin Strands

The structural mechanics of tropomyosin are essential determinants of its affinity and positioning on F-actin. Thus, tissue-specific differences among tropomyosin isoforms may influence both access of actin-binding proteins along the actin filaments and the cooperativity of actin-myosin interactions. Here, 40 nm long smooth and striated muscle tropomyosin molecules were rotary-shadowed and compared by means of electron microscopy. Electron microscopy shows that striated muscle tropomyosin primarily consists of single molecules or paired molecules linked end-to-end. In contrast, smooth muscle tropomyosin is more a mixture of varying-length chains of end-to-end polymers. Both isoforms are characterized by gradually bending molecular contours that lack obvious signs of kinking. The flexural stiffness of the tropomyosins was quantified and evaluated. The persistence lengths along the shaft of rotary-shadowed smooth and striated muscle tropomyosin molecules are equivalent to each other (approximately 100 nm) and to values obtained from molecular-dynamics simulations of the tropomyosins; however, the persistence length surrounding the end-to-end linkage is almost twofold higher for smooth compared to cardiac muscle tropomyosin. The tendency of smooth muscle tropomyosin to form semi-rigid polymers with continuous and undampened rigidity may compensate for the lack of troponin-based structural support in smooth muscles and ensure positional fidelity on smooth muscle thin filaments.