CHARMM: The biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983. © 2009 Wiley Periodicals, Inc.J Comput Chem, 2009.Peer Reviewedhttp://deepblue.lib.umich.edu/bitstream/2027.42/63074/1/21287_ftp.pd

Molecular Dynamics Simulations using Advanced Sampling and Polarizable Force Fields

Molecular dynamics (MD) simulations were carried out for aqueous dipeptides, water over self-assembled monolayer (SAM) surfaces, and the nicotinic acetylcholine receptor (nAChR) ion channel. The main goal is to use advanced methods to increase the accuracy of molecular dynamics simulations while seeking solutions to problems relevant to chemistry, biophysics and materials science. In addition, activation energies of several cyclodimerization reactions were studied quantum mechanically. The simulations of the aqueous dipeptides and SAM surfaces involve modeling and detailed analysis of interfacial water, which is of interest to a range of fields from biology to materials science. For example, water has a central role in biology and medicine since biomolecules cannot function without water. Both sets of simulations were performed using both polarizable and nonpolarizable force fields. These systems were used as a test ground to assess the effects of explicit incorporation of polarizability and also to determine whether the models can adequately reproduce the experimental data, in particular, the aggregation data of aqueous dipeptides and contact angles of water over SAMs of different chemical character. Since the systems are well-characterized and relatively simple, they provide excellent models to test polarizable force fields to increase the accuracy of molecular dynamics simulations. Polarizable water was depolarized around dipeptide solutes and also at the interface with different SAM surfaces, reflecting its ability to adapt to heterogeneous electrostatic environments. Although the water shows more realistic structure and dynamics in the polarizable simulations, the peptide aggregation behavior agrees less well with the experiment. In this case, neither model successfully reproduces the experimental degree of aggregation. In the case of SAM surfaces, both sets of simulations produce fairly similar results. More studies are suggested to further test and improve the polarizable force fields. The third system studied is the modeling of wild-type and mutant nAChR ion channel proteins. Adaptive biasing force method was used to achieve improved sampling, and subsequently increase the efficiency and accuracy of MD simulations. The nAChR channels are involved in a number of cognitive and brain functions including learning and memory. Dysfunction in these receptors are associated in a variety of neuronal diseases including epilepsy, schizophrenia and Alzheimer\u27s Disease. The present study models the wild-type and two physiologically-relevant mutant structures to assess the effects of mutations on ion translocation energetics and the geometry of the channel. Open channel (conducting, active) structures were obtained from the available closed channel structure. One of the mutants was found to increase the energetic barrier for ion translocation, while the other one decreased the barrier. The ion channel structures were analyzed in detail to understand the structural changes that took place during the channel opening. The channel opening was found to be mediated by large-scale helix motions rather than small-scale side chain motions. Aside from the MD simulations, the final project involves quantum mechanical simulations, which are often needed in parametrization of molecular dynamics force fields. Density functional theory (DFT) calculations were employed to calculate the activation energies of three cyclodimerization reactions of trifluorovinyl ether monomers. The results agree with and further explain the experimentally observed reactivity in these types of reactions

Amyloid Fibril Nucleation In Reverse Micelles

The 40-residue amyloid beta protein (Abeta) is the unstructured cleavage product of a common membrane protein that is produced in large quantities, but normally cleared from the brain before it exerts any apparent toxicity. Under some conditions, however, it undergoes a conformational change and aggregates into fibrils. These fibrils then coalesce into amyloid plaques, which are the pathognomonic brain lesions of Alzheimer‘s disease. The plaques are centers of active oxidative stress and neuronal death, so the conditions under which fibrils form is of high interest. When Abeta is encapsulated in a reverse micelle, its infrared spectrum indicates that it spontaneously adopts a fibril-like structure, which is remarkable because only one Abeta strand is present in each reverse micelle. That observation suggests that some aspect of the reverse micelle environment such as crowding, dehydration, proximity to a membrane, or high ionic strength may induce Abeta to nucleate amyloid fibril formation. Therefore, an understanding of the factors that induce Abeta to adopt fibril-like structure in reverse micelles may reveal what causes amyloid fibrils to form in Alzheimer\u27s disease. Molecular dynamics simulations of Abeta in reverse micelles have been performed to identify and understand these factors. Results indicate that Abeta side chains penetrate the reverse micelle surface, anchoring the peptide in the membrane. Other interactions between peptide and membrane stabilize intrachain hydrogen bond formation and secondary structure. These interactions may be important factors in the formation of amyloid fibrils and the pathogenesis of Alzheimer‘s disease

Methods Development and Force Field Evaluation for Molecular Simulations of Iinteractions Between Structured Peptides and Functionalized Material Surfaces

The process of protein adsorption to material surfaces is highly complex and it is one of the most fundamental concepts upon which progress in the field of bioengineering is based. The strategic design of material surfaces for optimal utility in specific biological environments is absolutely dependent upon a thorough understanding of the mechanisms underlying protein adsorption, yet there is still a very limited understanding of these mechanisms. The primary reason for this lack of understanding is that protein adsorption is a dynamic process which occurs at the atomic and macromolecular scale, where experimental analyses provide a view that is static and too coarse to elucidate the stepwise processes behind this critical biochemical phenomenon. In recent years, continual improvements in speed and efficiency of computational hardware and simulation techniques have enabled the use of molecular simulation for studying systems of the size necessary for examining the mechanistic details of protein adsorption (tens to hundreds of thousands of atoms). Of the various forms of molecular simulation, all-atom empirical force field molecular dynamics (MD) simulation has shown the greatest potential for exploring the nature of protein adsorption because it offers a dynamic view of nanosecond-scale processes with atomistic detail. However, a shortcoming of the application of MD in studying protein adsorption is that the most widely used MD force fields (i.e., equations and parameter sets used for calculating structural and energetic properties) have been designed and validated for simulations of solvated molecular systems in the absence of solid surfaces. To address this shortcoming of an otherwise extremely powerful research tool, an initial evaluation of the applicability of existing MD force fields to model systems of structured peptides interacting with functionalized material surfaces is warranted. The work presented here encompasses that initial evaluation of force fields. Numerous detailed analyses of water, ions, and peptides were completed in order to provide the most accurate and comprehensive examination of simulated peptide adsorption available. As a result of this work, simulation methods for these unique systems were tested and determined to be appropriate for accurately representing experimental results. Also, a comparative evaluation of force field performance identified the force field that most consistently reflects experimental findings

ATOMISTIC SIMULATIONS OF INTRINSICALLY DISORDERED PROTEIN FOLDING AND DYNAMICS

Intrinsically disordered proteins (IDPs) are crucial in biology and human diseases, necessitating a comprehensive understanding of their structure, dynamics, and interactions. Atomistic simulations have emerged as a key tool for unraveling the molecular intricacies and establishing mechanistic insights into how these proteins facilitate diverse biological functions. However, achieving accurate simulations requires both an appropriate protein force field capable of describing the energy landscape of functionally relevant IDP conformations and sufficient conformational sampling to capture the free energy landscape of IDP dynamics. These factors are fundamental in comprehending potential IDP structures, dynamics, and interactions. I first conducted explicit solvent simulations to assess the performance of two state-of-the-art protein force fields, namely CHARMM36m and a99SB-disp, in capturing the stability of small protein-protein interactions. To evaluate their accuracy, I selected a set of 46 amino acid backbone and side chain pairs with representative configurations and computed the free energy profiles of their interactions. The results demonstrated that CHARMM36m consistently predicted stronger protein-protein interactions compared to a99SB-disp. Notably, the most significant overestimation in CHARMM36m occurred in charged pairs involving Arg and Glu side chains, with an overestimation of up to 2.9 kcal/mol. Through free energy decomposition analysis, I determined that these overestimations were primarily driven by protein-water electrostatic interactions rather than van der Waals (vdW) interactions. Consequently, these findings suggest that careful rebalancing of electrostatic interactions should be considered in the further optimization of protein force fields. In order to enhance the conformational sampling of IDPs, I developed an integrated approach that combines an improved implicit solvent model called Generalized Born with molecular volume and solvent accessible surface area (GBMV2/SA) with a multiscale enhanced sampling (MSES) technique. To make this approach more efficient, I implemented it as a standalone OpenMM plugin on Graphics Processing Units (GPUs). The results demonstrated that the GPU-GBMV2/SA model achieved numerical equivalence to the original CPU-GBMV2/SA models, while providing a remarkable ~60x speedup on a single NVIDIA TITAN X (Pascal) graphics card for molecular dynamic simulations of both folded and unstructured proteins. This significant acceleration greatly facilitated the application of the approach in biomolecular simulations. In addition, I conducted an evaluation of the reliability of GBMV2/SA models in simulating both folded and unfolded proteins. The results revealed that the GBMV2/SA model accurately describes small proteins, but its applicability is limited when it comes to larger proteins such as KID and p53-TAD proteins. This limitation can be attributed to the absence of long-range solute-solvent dispersion interactions in the model. To address this issue, I introduced a comprehensive treatment of nonpolar solvation free energy called GBMV2/NP model. Unfortunately, the GBMV2/NP model exhibited a destabilizing effect on well-folded proteins, particularly larger ones, due to an inaccurate representation of the repulsive solvent accessible surface area (SASA) model caused by the utilization of unphysical van der Waals volume. This observation highlights the need for further improvements in accurately describing the nonpolar term in the model

Interaction of DNA with groove binding ligands

Synthetic molecules that target the major groove in a sequence-selective way are a major goal in molecular medicine. Recently a major step has been taken toward achieving this goal: a novel cylinder has been developed that binds strongly into the major groove of DNA. Experimental techniques have provided some information regarding the binding strength and preferred binding sites of the cylinder on DNA. From all the experimental data it is clear that the parent cylinder binds in the major groove and is able to induce dramatic conformational changes in the DNA; these are unprecedented effects with synthetic DNA binders. However, gaining molecular level information in such a macromolecular system is challenging. Molecular dynamics (MD) simulations can provide information at the molecular level that is complementary to experiment and therefore are an ideal way to get a better understanding of this system. In this work we present the results of various MD simulations designed to probe the DNA-cylinder system. We have studied the effect of using CHARMM22 and CHARMM27 as the force-field for the simulations. Results showed that uncomplexed DNA simulated with CHARMM22 was less stable in the B-form than the comparable strand of DNA simulated with CHARMM27. Investigations into the effects of the cylinders charge and shape are also reported. Multi-nanosecond simulations were performed using two related synthetic cylinders, one with two Fe(II) metal centers and the other with two Cu(I) centers, and DNA. Finally the role of DNA within the system was investigated by performing a series of simulations of the cylinders with d(ATATATATATAT)2, d(CGCGCGCGCGCG)2 and d(CGCGCATATACG). Simulations with these DNA strands has only produced one system (CCu2+ with d(ATATATATATAT)) where the cylinder causes a conformational change in the DNA

Metal Cations in Protein Force Fields: From Data Set Creation and Benchmarks to Polarizable Force Field Implementation and Adjustment

Metal cations are essential to life. About one-third of all proteins require metal cofactors to accurately fold or to function. Computer simulations using empirical parameters and classical molecular mechanics models (force fields) are the standard tool to investigate proteins’ structural dynamics and functions in silico. Despite many successes, the accuracy of force fields is limited when cations are involved. The focus of this thesis is the development of tools and strategies to create system-specific force field parameters to accurately describe cation-protein interactions. The accuracy of a force field mainly relies on (i) the parameters derived from increasingly large quantum chemistry or experimental data and (ii) the physics behind the energy formula. The first part of this thesis presents a large and comprehensive quantum chemistry data set on a consistent computational footing that can be used for force field parameterization and benchmarking. The data set covers dipeptides of the 20 proteinogenic amino acids with different possible side chain protonation states, 3 divalent cations (Ca2+, Mg2+, and Ba2+), and a wide relative energy range. Crucial properties related to force field development, such as partial charges, interaction energies, etc., are also provided. To make the data available, the data set was uploaded to the NOMAD repository and its data structure was formalized in an ontology. Besides a proper data basis for parameterization, the physics covered by the terms of the additive force field formulation model impacts its applicability. The second part of this thesis benchmarks three popular non-polarizable force fields and the polarizable Drude model against a quantum chemistry data set. After some adjustments, the Drude model was found to reproduce the reference interaction energy substantially better than the non-polarizable force fields, which showed the importance of explicitly addressing polarization effects. Tweaking of the Drude model involved Boltzmann-weighted fitting to optimize Thole factors and Lennard-Jones parameters. The obtained parameters were validated by (i) their ability to reproduce reference interaction energies and (ii) molecular dynamics simulations of the N-lobe of calmodulin. This work facilitates the improvement of polarizable force fields for cation-protein interactions by quantum chemistry-driven parameterization combined with molecular dynamics simulations in the condensed phase. While the Drude model exhibits its potential simulating cation-protein interactions, it lacks description of charge transfer effects, which are significant between cation and protein. The CTPOL model extends the classical force field formulation by charge transfer (CT) and polarization (POL). Since the CTPOL model is not readily available in any of the popular molecular-dynamics packages, it was implemented in OpenMM. Furthermore, an open-source parameterization tool, called FFAFFURR, was implemented that enables the (system specific) parameterization of OPLS-AA and CTPOL models. Following the method established in the previous part, the performance of FFAFFURR was evaluated by its ability to reproduce quantum chemistry energies and molecular dynamics simulations of the zinc finger protein. In conclusion, this thesis steps towards the development of next-generation force fields to accurately describe cation-protein interactions by providing (i) reference data, (ii) a force field model that includes charge transfer and polarization, and (iii) a freely-available parameterization tool.Metallkationen sind für das Leben unerlässlich. Etwa ein Drittel aller Proteine benötigen Metall-Cofaktoren, um sich korrekt zu falten oder zu funktionieren. Computersimulationen unter Verwendung empirischer Parameter und klassischer Molekülmechanik-Modelle (Kraftfelder) sind ein Standardwerkzeug zur Untersuchung der strukturellen Dynamik und Funktionen von Proteinen in silico. Trotz vieler Erfolge ist die Genauigkeit der Kraftfelder begrenzt, wenn Kationen beteiligt sind. Der Schwerpunkt dieser Arbeit liegt auf der Entwicklung von Werkzeugen und Strategien zur Erstellung systemspezifischer Kraftfeldparameter zur genaueren Beschreibung von Kationen-Protein-Wechselwirkungen. Die Genauigkeit eines Kraftfelds hängt hauptsächlich von (i) den Parametern ab, die aus immer größeren quantenchemischen oder experimentellen Daten abgeleitet werden, und (ii) der Physik hinter der Kraftfeld-Formel. Im ersten Teil dieser Arbeit wird ein großer und umfassender quantenchemischer Datensatz auf einer konsistenten rechnerischen Grundlage vorgestellt, der für die Parametrisierung und das Benchmarking von Kraftfeldern verwendet werden kann. Der Datensatz umfasst Dipeptide der 20 proteinogenen Aminosäuren mit verschiedenen möglichen Seitenketten-Protonierungszuständen, 3 zweiwertige Kationen (Ca2+, Mg2+ und Ba2+) und einen breiten relativen Energiebereich. Wichtige Eigenschaften für die Entwicklung von Kraftfeldern, wie Wechselwirkungsenergien, Partialladungen usw., werden ebenfalls bereitgestellt. Um die Daten verfügbar zu machen, wurde der Datensatz in das NOMAD-Repository hochgeladen und seine Datenstruktur wurde in einer Ontologie formalisiert. Neben einer geeigneten Datenbasis für die Parametrisierung beeinflusst die Physik, die von den Termen des additiven Kraftfeld-Modells abgedeckt wird, dessen Anwendbarkeit. Der zweite Teil dieser Arbeit vergleicht drei populäre nichtpolarisierbare Kraftfelder und das polarisierbare Drude-Modell mit einem Datensatz aus der Quantenchemie. Nach einigen Anpassungen stellte sich heraus, dass das Drude-Modell die Referenzwechselwirkungsenergie wesentlich besser reproduziert als die nichtpolarisierbaren Kraftfelder, was zeigt, wie wichtig es ist, Polarisationseffekte explizit zu berücksichtigen. Die Anpassung des Drude-Modells umfasste eine Boltzmann-gewichtete Optimierung der Thole-Faktoren und Lennard-Jones-Parameter. Die erhaltenen Parameter wurden validiert durch (i) ihre Fähigkeit, Referenzwechselwirkungsenergien zu reproduzieren und (ii) Molekulardynamik-Simulationen des Calmodulin-N-Lobe. Diese Arbeit demonstriert die Verbesserung polarisierbarer Kraftfelder für Kationen-Protein-Wechselwirkungen durch quantenchemisch gesteuerte Parametrisierung in Kombination mit Molekulardynamiksimulationen in der kondensierten Phase. Während das Drude-Modell sein Potenzial bei der Simulation von Kation - Protein - Wechselwirkungen zeigt, fehlt ihm die Beschreibung von Ladungstransfereffekten, die zwischen Kation und Protein von Bedeutung sind. Das CTPOL-Modell erweitert die klassische Kraftfeldformulierung um den Ladungstransfer (CT) und die Polarisation (POL). Da das CTPOL-Modell in keinem der gängigen Molekulardynamik-Pakete verfügbar ist, wurde es in OpenMM implementiert. Außerdem wurde ein Open-Source-Parametrisierungswerkzeug namens FFAFFURR implementiert, welches die (systemspezifische) Parametrisierung von OPLS-AA und CTPOL-Modellen ermöglicht. In Anlehnung an die im vorangegangenen Teil etablierte Methode wurde die Leistung von FFAFFURR anhand seiner Fähigkeit, quantenchemische Energien und Molekulardynamiksimulationen des Zinkfingerproteins zu reproduzieren, bewertet. Zusammenfassend lässt sich sagen, dass diese Arbeit einen Schritt in Richtung der Entwicklung von Kraftfeldern der nächsten Generation zur genauen Beschreibung von Kationen-Protein-Wechselwirkungen darstellt, indem sie (i) Referenzdaten, (ii) ein Kraftfeldmodell, das Ladungstransfer und Polarisation einschließt, und (iii) ein frei verfügbares Parametrisierungswerkzeug bereitstellt

The severity of osteogenesis imperfecta: A comparison to the relative free energy differences of collagen model peptides

Molecular dynamics simulations were carried out to calculate free energy differences between the folded and unfolded states of wild type and mutant collagen model peptides. The calculated stability of the collagen models was compared with the severity of osteogenesis imperfecta. Free energy differences of Gly → Xaa (Xaa: Ser, Cys, Glu, and Asp) mutations between the wild type and the mutants at position 15 of the model peptide were 3.8, 4.2, 5.6, and 8.8 kcal/mol, respectively. The corresponding free energy differences of a second Gly mutation at the same position in different chains were, on average, 1.3, 1.5, 2.9, and 5.4 kcal/mol, respectively. Free energy simulations were also performed to estimate the relative stability between an oxidized form and a reduced form of the mutants containing two Cys residues, which indicated that the mutant of the collagen-like peptide containing an intramolecular disulfide bond was more stable than the mutant containing one Cys residue but less stable than the wild type. The calculated free energy differences between an oxidized and a reduced form of the mutants containing two Cys residues are 0.8 and 2.6 kcal/mol for the disulfide bonds between Chains A and B and between Chains A and C, respectively. © 2010 Wiley Periodicals, Inc. Biopolymers 95: 182–193, 2011.Peer Reviewedhttp://deepblue.lib.umich.edu/bitstream/2027.42/78503/1/21552_ftp.pd

MOLECULAR DYNAMICS SIMULATIONS OF BIOLOGICAL MACROMOLECULES: APPLICATIONS TO STRUCTURAL VACCINOLOGY AND PEPTIDE DESIGN

This thesis work is splitted into two parts. The first one is about a computational method for epitope predictions on antigenic proteins, while the second one is related to the characterization of folding/unfolding processes of small natural polypeptides. Starting with the first topic, an increasing number of functional studies of proteins have shown that sequence and structural similarities alone may not be sufficient for reliable prediction of their interaction properties. This is particularly true for proteins recognizing specific antibodies, where the prediction of antibody-binding sites, called epitopes, has proven challenging. The antibody-binding properties of an antigen depend on its structure and related dynamics. Aiming to predict the antibody-binding regions of a protein, we investigate a new approach based on the integrated analysis of the dynamical and energetic properties of antigens, to identify nonoptimized, low-intensity energetic interaction networks in the protein structure isolated in solution. The method is based on the idea that recognition sites may correspond to localized regions with low-intensity energetic couplings with the rest of the protein, which allows them to undergo conformational changes, to be recognized by a binding partner, and to tolerate mutations with minimal energetic expense. Upon analyzing the results on isolated proteins and benchmarking against antibody complexes, it is found that the method successfully identifies binding sites located on the protein surface that are accessible to putative binding partners. The combination of dynamics and energetics can thus discriminate between epitopes and other substructures based only on physical properties. A public web server (BEPPE) has been implemented with MLCE method in order to make it available to the scientific community. Changing topic to folding/unfolding, the analysis of the folding mechanism in peptides adopting well defined secondary structure is fundamental to understand protein folding. Herein, we describe the thermal unfolding of two 15-mer polypeptides (called QK and QK-L10A) homologue to the vascular endothelial growth factor binding region. In particular, on the basis of the temperature dependencies, we characterize the molecules through the combination of spectroscopic (CD and NMR) and computational analyses (MD) highlighting their folding/unfolding steps and how these structures can be used in peptide design

Structural basis for the reactivity of Nitrophorin 7 with diatomic ligands and its interaction with membranes

The present thesis reports exclusively about the main research thread concerning Rhodnius prolixus Nitrophorin 7, a peculiar member of a family of nitric oxide trasporting proteins. This particular protein displays a distinctive behaviour, when compared with the other, better known members of its family, in terms of NO transport and release mechanisms and of interactions with phospholipidic membranes. The computational studies, aimed at clarifying aspects such as the structure of the system of cavities and tunnels inside the protein, the pH-triggered conformational transition, the membrane interaction and its influence on the NO stocking and release, were conducted in collaboration with prof. F. Javier Luque’s group (Universitat de Barcelona) and dr. Markus Knipp’s group at the Max Planck Institut für chemische Energiekonversion (Mülheim a. d. Ruhr)