CHARMM: the biomolecular simulation program

Base-flipping mechanism in postmismatch recognition by MutS

DNA mismatch recognition and repair is vital for preserving the fidelity of the genome. Conserved across prokaryotes and eukaryotes, MutS is the primary protein that is responsible for recognizing a variety of DNA mismatches. From molecular dynamics simulations of the Escherichia coli MutS-DNA complex, we describe significant conformational dynamics in the DNA surrounding a G·T mismatch that involves weakening of the basepair hydrogen bonding in the basepair adjacent to the mismatch and, in one simulation, complete base opening via the major groove. The energetics of base flipping was further examined with Hamiltonian replica exchange free energy calculations revealing a stable flipped-out state with an initial barrier of ~2 kcal/mol. Furthermore, we observe changes in the local DNA structure as well as in the MutS structure that appear to be correlated with base flipping. Our results suggest a role of base flipping as part of the repair initiation mechanism most likely leading to sliding-clamp formation.

Molecular dynamics simulations and in silico peptide ligand screening of the Elk-1 ETS domain

BACKGROUND

The Elk-1 transcription factor is a member of a group of proteins called ternary complex factors, which serve as a paradigm for gene regulation in response to extracellular signals. Its deregulation has been linked to multiple human diseases including the development of tumours. The work herein aims to inform the design of potential peptidomimetic compounds that can inhibit the formation of the Elk-1 dimer, which is key to Elk-1 stability. We have conducted molecular dynamics simulations of the Elk-1 ETS domain followed by virtual screening.

RESULTS

We show the ETS dimerisation site undergoes conformational reorganisation at the α1β1 loop. Through exhaustive screening of di- and tri-peptide libraries against a collection of ETS domain conformations representing the dynamics of the loop, we identified a series of potential binders for the Elk-1 dimer interface. The di-peptides showed no particular preference toward the binding site; however, the tri-peptides made specific interactions with residues: Glu17, Gln18 and Arg49 that are pivotal to the dimer interface.

CONCLUSIONS

We have shown molecular dynamics simulations can be combined with virtual peptide screening to obtain an exhaustive docking protocol that incorporates dynamic fluctuations in a receptor. Based on our findings, we suggest experimental binding studies to be performed on the 12 SILE ranked tri-peptides as possible compounds for the design of inhibitors of Elk-1 dimerisation. It would also be reasonable to consider the score-ranked tri-peptides as a comparative test to establish whether peptide size is a determinant factor of binding to the ETS domain.

Molecular dynamics simulations of the Cx26 hemichannel: evaluation of structural models with Brownian dynamics

The recently published crystal structure of the Cx26 gap junction channel provides a unique opportunity for elucidation of the structure of the conductive connexin pore and the molecular determinants of its ion permeation properties (conductance, current-voltage [I-V] relations, and charge selectivity). However, the crystal structure was incomplete, most notably lacking the coordinates of the N-terminal methionine residue, which resides within the pore, and also lacking two cytosolic domains. To allow computational studies for comparison with the known channel properties, we completed the structure. Grand canonical Monte Carlo Brownian dynamics (GCMC/BD) simulations of the completed and the published Cx26 hemichannel crystal structure indicate that the pore is too narrow to permit significant ion flux. The GCMC/BD simulations predict marked inward current rectification and almost perfect anion selectivity, both inconsistent with known channel properties. The completed structure was refined by all-atom molecular dynamics (MD) simulations (220 ns total) in an explicit solvent and POPC membrane system. These MD simulations produced an equilibrated structure with a larger minimal pore diameter, which decreased the height of the permeation barrier formed by the N terminus. GCMC/BD simulations of the MD-equilibrated structure yielded more appropriate single-channel conductance and less anion/cation selectivity. However, the simulations much more closely matched experimentally determined I-V relations when the charge effects of specific co- and posttranslational modifications of Cx26 previously identified by mass spectrometry were incorporated. We conclude that the average equilibrated structure obtained after MD simulations more closely represents the open Cx26 hemichannel structure than does the crystal structure, and that co- and posttranslational modifications of Cx26 hemichannels are likely to play an important physiological role by defining the conductance and ion selectivity of Cx26 channels. Furthermore, the simulations and data suggest that experimentally observed heterogeneity in Cx26 I-V relations can be accounted for by variation in co- and posttranslational modifications.

Protein-protein binding affinity prediction on a diverse set of structures

MOTIVATION

Accurate binding free energy functions for protein-protein interactions are imperative for a wide range of purposes. Their construction is predicated upon ascertaining the factors that influence binding and their relative importance. A recent benchmark of binding affinities has allowed, for the first time, the evaluation and construction of binding free energy models using a diverse set of complexes, and a systematic assessment of our ability to model the energetics of conformational changes.

RESULTS

We construct a large set of molecular descriptors using commonly available tools, introducing the use of energetic factors associated with conformational changes and disorder to order transitions, as well as features calculated on structural ensembles. The descriptors are used to train and test a binding free energy model using a consensus of four machine learning algorithms, whose performance constitutes a significant improvement over the other state of the art empirical free energy functions tested. The internal workings of the learners show how the descriptors are used, illuminating the determinants of protein-protein binding.

AVAILABILITY

The molecular descriptor set and descriptor values for all complexes are available in the Supplementary Material. A web server for the learners and coordinates for the bound and unbound structures can be accessed from the website: http://bmm.cancerresearchuk.org/~Affinity.

CONTACT

paul.bates@cancer.org.uk.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

MATCH: an atom-typing toolset for molecular mechanics force fields

We introduce a toolset of program libraries collectively titled multipurpose atom-typer for CHARMM (MATCH) for the automated assignment of atom types and force field parameters for molecular mechanics simulation of organic molecules. The toolset includes utilities for the conversion of multiple chemical structure file formats into a molecular graph. A general chemical pattern-matching engine using this graph has been implemented whereby assignment of molecular mechanics atom types, charges, and force field parameters are achieved by comparison against a customizable list of chemical fragments. While initially designed to complement the CHARMM simulation package and force fields by generating the necessary input topology and atom-type data files, MATCH can be expanded to any force field and program, and has core functionality that makes it extendable to other applications such as fragment-based property prediction. In this work, we demonstrate the accurate construction of atomic parameters of molecules within each force field included in CHARMM36 through exhaustive cross validation studies illustrating that bond charge increment rules derived from one force field can be transferred to another. In addition, using leave-one-out substitution it is shown that it is also possible to substitute missing intra and intermolecular parameters with ones included in a force field to complete the parameterization of novel molecules. Finally, to demonstrate the robustness of MATCH and the coverage of chemical space offered by the recent CHARMM general force field (Vanommeslaeghe, et al., J Comput Chem 2010, 31, 671), one million molecules from the PubChem database of small molecules are typed, parameterized, and minimized.

Design of a modified mouse protein with ligand binding properties of its human analog by molecular dynamics simulations: the case of C3 inhibition by compstatin

The peptide compstatin and its derivatives inhibit the complement-component protein C3 in primate mammals and are potential therapeutic agents against the unregulated activation of complement in humans, but are inactive against C3 from lower mammals. Recent molecular dynamics (MD) simulations showed that the most potent compstatin analog comprised entirely of natural amino acids (W4A9) had a smaller affinity for rat C3, due to reproducible changes in the rat protein structure with respect to the human protein, which eliminated or weakened specific protein-ligand interactions seen in the human C3:W4A9 complex. Here, we study by MD simulations three W4A9 complexes with the mouse C3 protein, and two "transgenic" mouse derivatives, containing a small number (6-9) of human C3 substitutions. The mouse complex experiences the conformational changes and affinity reduction of the rat complex. In the "transgenic" complexes, the conformation remains closer to that of the human complex, the protein-ligand interactions are improved, and the affinity for compstatin becomes "human-like." The present work creates new avenues for a compstatin-sensitive animal model. A similar strategy, involving the comparison of a series of complexes by MD simulations, could be used to design "transgenic" sequences in other systems.

Influence of the membrane dipole potential on peptide binding to lipid bilayers

The implicit membrane model IMM1 is extended to include the membrane dipole potential and applied to molecular dynamics simulations of the helical peptides alamethicin, WALP23, influenza hemagglutinin fusion peptide, HIV fusion peptide, magainin, and the pre-sequence of cytochrome c oxidase subunit IV (p25). The results show that the orientation of the peptides in the membrane can be influenced by the dipole potential. The binding affinity of all peptides except for the hemagglutinin fusion peptide decreases upon increase of the dipole potential. The changes in both orientation and binding affinity are explained by the interaction of the dipole potential with the helix backbone dipole and ionic side-chains. In general, peptides that tend to insert the N-terminus in the membrane and/or have positively charged side chains will lose binding affinity upon increase of the dipole potential.

Molecular dynamics approaches estimate the binding energy of HIV-1 integrase inhibitors and correlate with in vitro activity

The design of novel integrase (IN) inhibitors has been aided by recent crystal structures revealing the binding mode of these compounds with a full-length prototype foamy virus (PFV) IN and synthetic viral DNA ends. Earlier docking studies relied on incomplete structures and did not include the contribution of the viral DNA to inhibitor binding. Using the structure of PFV IN as the starting point, we generated a model of the corresponding HIV-1 complex and developed a molecular dynamics (MD)-based approach that correlates with the in vitro activities of novel compounds. Four well-characterized compounds (raltegravir, elvitegravir, MK-0536, and dolutegravir) were used as a training set, and the data for their in vitro activity against the Y143R, N155H, and G140S/Q148H mutants were used in addition to the wild-type (WT) IN data. Three additional compounds were docked into the IN-DNA complex model and subjected to MD simulations. All three gave interaction potentials within 1 standard deviation of values estimated from the training set, and the most active compound was identified. Additional MD analysis of the raltegravir- and dolutegravir-bound complexes gave internal and interaction energy values that closely match the experimental binding energy of a compound related to raltegravir that has similar activity. These approaches can be used to gain a deeper understanding of the interactions of the inhibitors with the HIV-1 intasome and to identify promising scaffolds for novel integrase inhibitors, in particular, compounds that retain activity against a range of drug-resistant mutants, making it possible to streamline synthesis and testing.

Resistance of Akt kinases to dephosphorylation through ATP-dependent conformational plasticity

Phosphorylation of a threonine residue (T308 in Akt1) in the activation loop of Akt kinases is a prerequisite for deregulated Akt activity frequently observed in neoplasia. Akt phosphorylation in vivo is balanced by the opposite activities of kinases and phosphatases. Here we describe that targeting Akt kinase to the cell membrane markedly reduced sensitivity of phosphorylated Akt to dephosphorylation by protein phosphatase 2A. This effect was amplified by occupancy of the ATP binding pocket by either ATP or ATP-competitive inhibitors. Mutational analysis revealed that R273 in Akt1 and the corresponding R274 in Akt2 are essential for shielding T308 in the activation loop against dephosphorylation. Thus, occupancy of the nucleotide binding pocket of Akt kinases enables intramolecular interactions that restrict phosphatase access and sustain Akt phosphorylation. This mechanism provides an explanation for the "paradoxical" Akt hyperphosphorylation induced by ATP-competitive inhibitor, A-443654. The lack of phosphatase resistance further contributes insight into the mechanism by which the human Akt2 R274H missense mutation may cause autosomal-dominant diabetes mellitus.

Specific interactions of sodium salts with alanine dipeptide and tetrapeptide in water: insights from molecular dynamics

We examine computationally the dipeptide and tetrapeptide of alanine in pure water and solutions of sodium chloride (NaCl) and iodide (NaI), with salt concentrations up to 3 M. Enhanced sampling of the configuration space is achieved by the replica exchange method. In agreement with other works, we observe preferential sodium interactions with the peptide carbonyl groups, which are enhanced in the NaI solutions due to the increased affinity of the less hydrophilic iodide anion for the peptide methyl side-chains and terminal blocking groups. These interactions have been associated with a decrease in the helicities of more complex peptides. In our simulations, both salts have a small effect on the dipeptide, but consistently stabilize the intramolecular hydrogen-bonding interactions and "α-helical" conformations of the tetrapeptide. This behavior, and an analysis of the intermolecular interaction energies show that ion-peptide interactions, or changes in the peptide hydration due to salts, are not sufficient determining factors of the peptide conformational preferences. Additional simulations suggest that the observed stabilizing effect is not due to the employed force-field, and that it is maintained in short peptides but is reversed in longer peptides. Thus, the peptide conformational preferences are determined by an interplay of energetic and entropic factors, arising from the peptide sequence and length and the composition of the solution.

NaCl interactions with phosphatidylcholine bilayers do not alter membrane structure but induce long-range ordering of ions and water

It is generally accepted that ions interact directly with lipids in biological membranes. Decades of biophysical studies on pure lipid bilayer systems have shown that only certain types of ions, most significantly large anions and multivalent cations, can fundamentally alter the structure and dynamics of lipid bilayers. It has long been accepted that at physiological concentrations NaCl ions do not alter the physical behavior or structure of bilayers composed solely of zwitterionic phosphatidylcholine (PC) lipids. Recent X-ray scattering experiments have reaffirmed this dogma, showing that below 1 M concentration, NaCl does not significantly alter bilayer structure. However, despite this history, there is an ongoing controversy within the molecular dynamics (MD) simulation community regarding NaCl/PC interactions. In particular, the CHARMM and GROMOS force fields show dramatically different behavior, including the effect on bilayer structure, surface potential, and the ability to form stable, coordinated ion-lipid complexes. Here, using long-timescale, constant-pressure simulations under the newest version of the CHARMM force field, we find that Na⁺ and Cl⁻ associate with PC head groups in a POPC bilayer with approximately equal, though weak, affinity, and that the salt has a negligible effect on bilayer structure, consistent with earlier CHARMM results and more recent X-ray data. The results suggest that interpretation of simulations where lipids interact with charged groups of any sort, including charged proteins, must be carefully scrutinized.

Comparison of two adaptive temperature-based replica exchange methods applied to a sharp phase transition of protein unfolding-folding

Temperature-based replica exchange (T-ReX) enhances sampling of molecular dynamics simulations by autonomously heating and cooling simulation clients via a Metropolis exchange criterion. A pathological case for T-ReX can occur when a change in state (e.g., folding to unfolding of a protein) has a large energetic difference over a short temperature interval leading to insufficient exchanges amongst replica clients near the transition temperature. One solution is to allow the temperature set to dynamically adapt in the temperature space, thereby enriching the population of clients near the transition temperature. In this work, we evaluated two approaches for adapting the temperature set: a method that equalizes exchange rates over all neighbor temperature pairs and a method that attempts to induce clients to visit all temperatures (dubbed "current maximization") by positioning many clients at or near the transition temperature. As a test case, we simulated the 57-residue SH3 domain of alpha-spectrin. Exchange rate equalization yielded the same unfolding-folding transition temperature as fixed-temperature ReX with much smoother convergence of this value. Surprisingly, the current maximization method yielded a significantly lower transition temperature, in close agreement with experimental observation, likely due to more extensive sampling of the transition state.

High performance screening, structural and molecular dynamics analysis to identify H1 inhibitors from TCM Database@Taiwan

New-type oseltamivir-resistant H1N1 influenza viruses have been a major threat to human health since the 2009 flu pandemic. To resolve the drug resistance issue, we aimed to identify a new type of inhibitors against H1 from traditional Chinese medicine (TCM) by employing the world's largest TCM database () for virtual screening and molecular dynamics (MD). From the virtual screening results, sodium (+)-isolaricireinol-2 alpha-sulfate, sodium 3,4-dihydroxy-5-methoxybenzoic acid methyl ester-4-sulfate, sodium (E)-7-hydroxy-1,7-bis(4-hydroxyphenyl)hept-5-ene-3S-sulfonate, and 3-methoxytyramine-betaxanthin were identified as potential drug-like compounds. MD simulation of the binding poses with the key residues Asp103 and Glu83, as well as other binding site residues, identified higher numbers of hydrogen bonds than N-Acetyl-D-Glucosamine (NAG), the natural ligand of the esterase domain in H1. Ionic bonds, salt bridges, and electrostatic energy also contribute to binding stability. Key binding residues include Lys71, Glu83, Asp103, and Arg238. Structural moieties promoting H-bond or salt bridge formations at these locations greatly contribute to a stable ligand-protein complex. An available sodium atom for ionic interactions with Asp103 can further stabilize the ligands. Based on virtual screening, MD simulation, and interaction energy evaluation, TCM candidates demonstrate good potential as novel H1 inhibitors. In addition, the identified stabilizing features can provide insights for designing highly stable H1 inhibitors.

Glu-286 rotation and water wire reorientation are unlikely the gating elements for proton pumping in cytochrome C oxidase

One of the key unresolved issues regarding proton pumping in cytochrome c oxidase (CcO) is the identity of the gating element that prevents the backflow of protons. In this study, we analyze two popular proposals for this element: isomerization of the key branching residue (Glu-286) and (re)orientation of water molecules in the hydrophobic cavity. Using a multifaceted set of computational analyses that involve CcO embedded in either an implicit or explicit treatment of lipid membrane, we show that neither Glu-286 nor active-site water likely constitutes the gating element. Detailed energetic and structural analyses of the simulation results indicate that the gating-relevant properties of these structural motifs observed in previous work are likely a result of the simplified computational models employed in those studies.

Unraveling the structural complexity in a single-stranded RNA tail: implications for efficient ligand binding in the prequeuosine riboswitch

Single-stranded RNAs (ssRNAs) are ubiquitous RNA elements that serve diverse functional roles. Much of our understanding of ssRNA conformational behavior is limited to structures in which ssRNA directly engages in tertiary interactions or is recognized by proteins. Little is known about the structural and dynamic behavior of free ssRNAs at atomic resolution. Here, we report the collaborative application of nuclear magnetic resonance (NMR) and replica exchange molecular dynamics (REMD) simulations to characterize the 12 nt ssRNA tail derived from the prequeuosine riboswitch. NMR carbon spin relaxation data and residual dipolar coupling measurements reveal a flexible yet stacked core adopting an A-form-like conformation, with the level of order decreasing toward the terminal ends. An A-to-C mutation within the polyadenine tract alters the observed dynamics consistent with the introduction of a dynamic kink. Pre-ordering of the tail may increase the efficacy of ligand binding above that achieved by a random-coil ssRNA. The REMD simulations recapitulate important trends in the NMR data, but suggest more internal motions than inferred from the NMR analysis. Our study unmasks a previously unappreciated level of complexity in ssRNA, which we believe will also serve as an excellent model system for testing and developing computational force fields.

Delineation of folding pathways of a β-sheet miniprotein

Several methods have been developed in the past few years for the analysis of molecular dynamics simulations of biological (macro)molecules whose complexity is difficult to capture by simple projections of the free-energy surface onto one or two geometric variables. The locally scaled diffusion map (LSDMap) method is a nonlinear dimensionality reduction technique for describing the dynamics of complex systems in terms of a few collective coordinates. Here, we compare LSDMap to two previously developed approaches for the characterization of the configurational landscape associated with the folding dynamics of a three-stranded antiparallel β-sheet peptide, termed Beta3s. The analysis is aided by an improved procedure for extracting pathways from the equilibrium transition network, which enables calculation of pathway-specific cut-based free energy profiles. We find that the results from LSDMap are consistent with analysis based on transition networks and allow a coherent interpretation of metastable states and folding pathways in terms of different time scales of transitions between minima on the free energy projections.

Molecular dynamics trajectory compression with a coarse-grained model

Molecular dynamics trajectories are very data-intensive thereby limiting sharing and archival of such data. One possible solution is compression of trajectory data. Here, trajectory compression based on conversion to the coarse-grained model PRIMO is proposed. The compressed data is about one third of the original data and fast decompression is possible with an analytical reconstruction procedure from PRIMO to all-atom representations. This protocol largely preserves structural features and to a more limited extent also energetic features of the original trajectory.

pH replica-exchange method based on discrete protonation states

We propose a new algorithm for obtaining proton titration curves of ionizable residues. The algorithm is a pH replica-exchange method (PHREM), which is based on the constant pH algorithm of Mongan et al. (J Comput Chem 2004;25:2038-2048). In the original replica-exchange method, simulations of different replicas are performed at different temperatures, and the temperatures are exchanged between the replicas. In our PHREM, simulations of different replicas are performed at different pH values, and the pHs are exchanged between the replicas. The PHREM was applied to a blocked amino acid and to two protein systems (snake cardiotoxin and turkey ovomucoid third domain), in conjunction with a generalized Born implicit solvent. The performance and accuracy of this algorithm and the original constant pH method (PHMD) were compared. For a single set of simulations at different pHs, the use of PHREM yields more accurate Hill coefficients of titratable residues. By performing multiple sets of constant pH simulations started with different initial states, the accuracy of predicted pK(a) values and Hill coefficients obtained with PHREM and PHMD methods becomes comparable. However, the PHREM algorithm exhibits better samplings of the protonation states of titratable residues and less scatter of the titration points and thus better precision of measured pK(a) values and Hill coefficients. In addition, PHREM exhibits faster convergence of individual simulations than the original constant pH algorithm.

Surfactant effects on amyloid aggregation kinetics

There is strong experimental evidence of the influence of surfactants (e.g., fatty acids) on the kinetics of amyloid fibril formation. However, the structures of mixed assemblies and interactions between surfactants and fibril-forming peptides are still not clear. Here, coarse-grained simulations are employed to study the aggregation kinetics of amyloidogenic peptides in the presence of amphiphilic lipids. The simulations show that the lower the fibril formation propensity of the peptides, the higher the influence of the surfactants on the peptide self-assembly kinetics. In particular, the lag phase of weakly aggregating peptides increases because of the formation of mixed oligomers, which are promoted by hydrophobic interactions and favorable entropy of mixing. A transient peak in the number of surfactants attached to the growing fibril is observed before reaching the mature fibril in some of the simulations. This peak originates from transient fibrillar defects consisting of exposed hydrophobic patches on the fibril surface, which provide a possible explanation for the temporary maximum of fluorescence observed sometimes in kinetic traces of the binding of small-molecule dyes to amyloid fibrils.

A multilayer evaluation approach for protein structure prediction and model quality assessment

Protein tertiary structures are essential for studying functions of proteins at molecular level. An indispensable approach for protein structure solution is computational prediction. Most protein structure prediction methods generate candidate models first and select the best candidates by model quality assessment (QA). In many cases, good models can be produced, but the QA tools fail to select the best ones from the candidate model pool. Because of incomplete understanding of protein folding, each QA method only reflects partial facets of a structure model and thus has limited discerning power with no one consistently outperforming others. In this article, we developed a set of new QA methods, including two QA methods for evaluating target/template alignments, a molecular dynamics (MD)-based QA method, and three consensus QA methods with selected references to reveal new facets of protein structures complementary to the existing methods. Moreover, the underlying relationship among different QA methods were analyzed and then integrated into a multilayer evaluation approach to guide the model generation and model selection in prediction. All methods are integrated and implemented into an innovative and improved prediction system hereafter referred to as MUFOLD. In CASP8 and CASP9, MUFOLD has demonstrated the proof of the principles in terms of both QA discerning power and structure prediction accuracy.

Identification of potent EGFR inhibitors from TCM Database@Taiwan

Overexpression of epidermal growth factor receptor (EGFR) has been associated with cancer. Targeted inhibition of the EGFR pathway has been shown to limit proliferation of cancerous cells. Hence, we employed Traditional Chinese Medicine Database (TCM Database@Taiwan) (http://tcm.cmu.edu.tw) to identify potential EGFR inhibitor. Multiple Linear Regression (MLR), Support Vector Machine (SVM), Comparative Molecular Field Analysis (CoMFA), and Comparative Molecular Similarities Indices Analysis (CoMSIA) models were generated using a training set of EGFR ligands of known inhibitory activities. The top four TCM candidates based on DockScore were 2-O-caffeoyl tartaric acid, Emitine, Rosmaricine, and 2-O-feruloyl tartaric acid, and all had higher binding affinities than the control Iressa®. The TCM candidates had interactions with Asp855, Lys716, and Lys728, all which are residues of the protein kinase binding site. Validated MLR (r² = 0.7858) and SVM (r² = 0.8754) models predicted good bioactivity for the TCM candidates. In addition, the TCM candidates contoured well to the 3D-Quantitative Structure-Activity Relationship (3D-QSAR) map derived from the CoMFA (q² = 0.721, r² = 0.986) and CoMSIA (q² = 0.662, r² = 0.988) models. The steric field, hydrophobic field, and H-bond of the 3D-QSAR map were well matched by each TCM candidate. Molecular docking indicated that all TCM candidates formed H-bonds within the EGFR protein kinase domain. Based on the different structures, H-bonds were formed at either Asp855 or Lys716/Lys728. The compounds remained stable throughout molecular dynamics (MD) simulation. Based on the results of this study, 2-O-caffeoyl tartaric acid, Emitine, Rosmaricine, and 2-O-feruloyl tartaric acid are suggested to be potential EGFR inhibitors.

CHARMM additive all-atom force field for carbohydrate derivatives and its utility in polysaccharide and carbohydrate-protein modeling

Monosaccharide derivatives such as xylose, fucose, N-acetylglucosamine (GlcNAc), N-acetylgalactosamine (GlaNAc), glucuronic acid, iduronic acid, and N-acetylneuraminic acid (Neu5Ac) are important components of eukaryotic glycans. The present work details development of force-field parameters for these monosaccharides and their covalent connections to proteins via O-linkages to serine or threonine sidechains and via N-linkages to asparagine sidechains. The force field development protocol was designed to explicitly yield parameters that are compatible with the existing CHARMM additive force field for proteins, nucleic acids, lipids, carbohydrates, and small molecules. Therefore, when combined with previously developed parameters for pyranose and furanose monosaccharides, for glycosidic linkages between monosaccharides, and for proteins, the present set of parameters enables the molecular simulation of a wide variety of biologically-important molecules such as complex carbohydrates and glycoproteins. Parametrization included fitting to quantum mechanical (QM) geometries and conformational energies of model compounds, as well as to QM pair interaction energies and distances of model compounds with water. Parameters were validated in the context of crystals of relevant monosaccharides, as well NMR and/or x-ray crystallographic data on larger systems including oligomeric hyaluronan, sialyl Lewis X, O- and N-linked glycopeptides, and a lectin:sucrose complex. As the validated parameters are an extension of the CHARMM all-atom additive biomolecular force field, they further broaden the types of heterogeneous systems accessible with a consistently-developed force-field model.

Microscopic Symmetry Imposed by Rotational Symmetry Boundary Conditions in Molecular Dynamics Simulation

A large number of viral capsids, as well as other macromolecular assemblies, have icosahedral structure or structures with other rotational symmetries. This symmetry can be exploited during molecular dynamics (MD) to model in effect the full viral capsid using only a subset of primary atoms plus copies of image atoms generated from rotational symmetry boundary conditions (RSBC). A pure rotational symmetry operation results in both primary and image atoms at short range, and within nonbonded interaction distance of each other, so that nonbonded interactions can not be specified by the minimum image convention and explicit treatment of image atoms is required. As such an unavoidable consequence of RSBC is that the enumeration of nonbonded interactions in regions surrounding certain rotational axes must include both a primary atom and its copied image atom, thereby imposing microscopic symmetry for some forces. We examined the possibility of artifacts arising from this imposed microscopic symmetry of RSBC using two simulation systems: a water shell and human rhinovirus 14 (HRV14) capsid with explicit water. The primary unit was a pentamer of the icosahedron, which has the advantage of direct comparison of icosahedrally equivalent spatial regions, for example regions near a 2-fold symmetry axis with imposed symmetry and a 2-fold axis without imposed symmetry. Analysis of structural and dynamic properties of water molecules and protein atoms found similar behavior near symmetry axes with imposed symmetry and where the minimum image convention fails compared with that in other regions in the simulation system, even though an excluded volume effect was detected for water molecules near the axes with imposed symmetry. These results validate the use of RSBC for icosahedral viral capsids or other rotationally symmetric systems.

Self-templating 2D supramolecular networks: a new avenue to reach control over a bilayer formation

One of the greatest challenges in 2D self-assembly at interfaces is the ability to grow spatially controlled supramolecular motifs in the third dimension, exploiting the surface as a template. In this manuscript a concentration-dependent study by scanning tunneling microscopy at the solid-liquid interface, corroborated by Molecular Dynamics (MD) simulations, reveals the controlled generation of mono- or bilayer self-assembled Kagomé networks based on a fully planar tetracarboxylic acid derivative. By programming the backbone of the molecular building blocks, we present a strategy to gain spatial control over the adlayer structure by conferring self-templating capacity to the 2D self-assembled network.

Predicting pKa values with continuous constant pH molecular dynamics

Knowledge of pK(a) values is important for understanding structure and function relationships in proteins. Over the past two decades, theoretical methods for pK(a) calculations have been mainly based on macroscopic models, in which the protein is considered as a low-dielectric cavity embedded in a high-dielectric continuum. In recent years, constant pH molecular dynamics methods have been developed based on a microscopic description of the protein. We describe here the methodology of continuous constant pH molecular dynamics (CPHMD), which has emerged as one of the most robust and accurate tools for predicting protein pK(a)s and for the study of pH-modulated conformational dynamics. We illustrate the utility of CPHMD by the calculation of pK(a)s for surface residues in ribonuclease A, buried residues in staphylococcal nuclease, and titratable groups in the intrinsically flexible protein α-lactalbumin. We will compare the CPHMD results with experimental data as well as calculations from PB-based and empirical methods. These examples demonstrate the accuracy and robustness of the CPHMD method and its ability to capture the correlation between ionization equilibria and conformational dynamics as well as the local dielectric response to structural rearrangement. Finally, we discuss future improvement of the CPHMD method.

MDpocket: open-source cavity detection and characterization on molecular dynamics trajectories

MOTIVATION

A variety of pocket detection algorithms are now freely or commercially available to the scientific community for the analysis of static protein structures. However, since proteins are dynamic entities, enhancing the capabilities of these programs for the straightforward detection and characterization of cavities taking into account protein conformational ensembles should be valuable for capturing the plasticity of pockets, and therefore allow gaining insight into structure-function relationships.

RESULTS

This article describes a new method, called MDpocket, providing a fast, free and open-source tool for tracking small molecule binding sites and gas migration pathways on molecular dynamics (MDs) trajectories or other conformational ensembles. MDpocket is based on the fpocket cavity detection algorithm and a valuable contribution to existing analysis tools. The capabilities of MDpocket are illustrated for three relevant cases: (i) the detection of transient subpockets using an ensemble of crystal structures of HSP90; (ii) the detection of known xenon binding sites and migration pathways in myoglobin; and (iii) the identification of suitable pockets for molecular docking in P38 Map kinase.

AVAILABILITY

MDpocket is free and open-source software and can be downloaded at http://fpocket.sourceforge.net.

CONTACT

pschmidtke@ub.edu

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

The Class D beta-lactamase family: residues governing the maintenance and diversity of function

Class D β-lactamases, a major source of bacterial resistance to β-lactam antibiotic therapies, represent a distinct subset of the β-lactamase superfamily. They share a serine hydrolase mechanism with Classes A/C vs. Class B. Further understanding of their sequence-structure-function relationships would benefit efforts to design a new generation of antibiotics as well as to predict evolutionary mechanisms in response to such therapies. Here we describe analyses based on our high-resolution multiple sequence alignment and phylogenetic tree of ∼80 Class D β-lactamases that leverage several 3D structures of these enzymes. We observe several sequence clusters on the phylogenetic tree, some that are species specific while others include several species from α-, β- and γ-proteobacteria. Residues characteristic of a specific cluster were identified and shown to be located just outside the active site, possibly modulating the function of the catalytic residues to facilitate reactions with specific types of β-lactams. Most significant was the discovery of a likely disulfide bond in a large group composed of α-, β- and γ-proteobacteria that would contribute to enzyme stability and hence bacterial viability under antibiotic assault. A network of co-evolving residues was identified which suggested the importance of maintaining a surface for binding a highly conserved Phe69.

Surveying implicit solvent models for estimating small molecule absolute hydration free energies

Implicit solvent models are powerful tools in accounting for the aqueous environment at a fraction of the computational expense of explicit solvent representations. Here, we compare the ability of common implicit solvent models (TC, OBC, OBC2, GBMV, GBMV2, GBSW, GBSW/MS, GBSW/MS2 and FACTS) to reproduce experimental absolute hydration free energies for a series of 499 small neutral molecules that are modeled using AMBER/GAFF parameters and AM1-BCC charges. Given optimized surface tension coefficients for scaling the surface area term in the nonpolar contribution, most implicit solvent models demonstrate reasonable agreement with extensive explicit solvent simulations (average difference 1.0-1.7 kcal/mol and R(2)=0.81-0.91) and with experimental hydration free energies (average unsigned errors=1.1-1.4 kcal/mol and R(2)=0.66-0.81). Chemical classes of compounds are identified that need further optimization of their ligand force field parameters and others that require improvement in the physical parameters of the implicit solvent models themselves. More sophisticated nonpolar models are also likely necessary to more effectively represent the underlying physics of solvation and take the quality of hydration free energies estimated from implicit solvent models to the next level.

Simulation of solution phase electron transfer in a compact donor-acceptor dyad

Charge separation (CS) and charge recombination (CR) rates in photosynthetic architectures are difficult to control, yet their ratio can make or break photon-to-current conversion efficiencies. A rational design approach to the enhancement of CS over CR requires a mechanistic understanding of the underlying electron-transfer (ET) process, including the role of the environment. Toward this goal, we introduce a QM/MM protocol for ET simulations and use it to characterize CR in the formanilide-anthraquinone dyad (FAAQ). Our simulations predict fast recombination of the charge-transfer excited state, in agreement with recent experiments. The computed electronic couplings show an electronic state dependence and are weaker in solution than in the gas phase. We explore the role of cis-trans isomerization on the CR kinetics, and we find strong correlation between the vertical energy gaps of the full simulations and a collective solvent polarization coordinate. Our approach relies on constrained density functional theory to obtain accurate diabatic electronic states on the fly for molecular dynamics simulations, while orientational and electronic polarization of the solvent is captured by a polarizable force field based on a Drude oscillator model. The method offers a unified approach to the characterization of driving forces, reorganization energies, electronic couplings, and nonlinear solvent effects in light-harvesting systems.

Multiple functions of aromatic-carbohydrate interactions in a processive cellulase examined with molecular simulation

Proteins employ aromatic residues for carbohydrate binding in a wide range of biological functions. Glycoside hydrolases, which are ubiquitous in nature, typically exhibit tunnels, clefts, or pockets lined with aromatic residues for processing carbohydrates. Mutation of these aromatic residues often results in significant activity differences on insoluble and soluble substrates. However, the thermodynamic basis and molecular level role of these aromatic residues remain unknown. Here, we calculate the relative ligand binding free energy by mutating tryptophans in the Trichoderma reesei family 6 cellulase (Cel6A) to alanine. Removal of aromatic residues near the catalytic site has little impact on the ligand binding free energy, suggesting that aromatic residues immediately upstream of the active site are not directly involved in binding, but play a role in the glucopyranose ring distortion necessary for catalysis. Removal of aromatic residues at the entrance and exit of the Cel6A tunnel, however, dramatically impacts the binding affinity, suggesting that these residues play a role in chain acquisition and product stabilization, respectively. The roles suggested from differences in binding affinity are confirmed by molecular dynamics and normal mode analysis. Surprisingly, our results illustrate that aromatic-carbohydrate interactions vary dramatically depending on the position in the enzyme tunnel. As aromatic-carbohydrate interactions are present in all carbohydrate-active enzymes, these results have implications for understanding protein structure-function relationships in carbohydrate metabolism and recognition, carbon turnover in nature, and protein engineering strategies for biomass utilization. Generally, these results suggest that nature employs aromatic-carbohydrate interactions with a wide range of binding affinities for diverse functions.

Influence of solvent and intramolecular hydrogen bonding on the conformational properties of o-linked glycopeptides

A detailed investigation of the conformational properties of all the biologically relevant O-glycosidic linkages using the Hamiltonian replica exchange (HREX) simulation methodology and the recently developed CHARMM carbohydrate force field parameters is presented. Fourteen biologically relevant O-linkages between the five sugars N-acetylgalactosamine (GalNAc), N-acetylglucosamine (GlcNAc), D-glucose (Glc), D-mannose (Man), and L-fucose (Fuc) and the amino acids serine and threonine were studied. The force field was tested by comparing the simulation results of the model glycopeptides to various NMR (3)J coupling constants, NOE distances, and data from molecular dynamics with time-averaged restraints (tar-MD). The results show the force field to be in overall agreement with experimental and previous tar-MD simulations, although some small limitations are identified. An in-depth hydrogen bond and bridging water analysis revealed an interplay of hydrogen bonding and bridge water interactions influencing the geometry of the underlying peptide backbone, with the O-linkages favoring extended β-sheet and polyproline type II (PPII) conformations over the compact α(R)-helical conformation. The newly developed parameters were also able to identify hydrogen bonding and water mediated interactions between O-linked sugars and proteins. These results indicate that the newly developed parameters in tandem with HREX conformational sampling provide the means to study glycoproteins in the absence of targeted NMR restraint data.

An antiparallel actin dimer is associated with the endocytic pathway in mammalian cells

The dynamic rearrangement of the actin cytoskeleton plays a key role in several cellular processes such as cell motility, endocytosis, RNA processing and chromatin organization. However, the supramolecular actin structures involved in the different processes remain largely unknown. One of the less studied forms of actin is the lower dimer (LD). This unconventional arrangement of two actin molecules in an antiparallel orientation can be detected by chemical crosslinking at the onset of polymerization in vitro. Moreover, evidence for a transient incorporation of LD into growing filaments and its ability to inhibit nucleation of F-actin filament assembly implicate that the LD pathway contributes to supramolecular actin patterning. However, a clear link from this actin species to a specific cellular function has not yet been established. We have developed an antibody that selectively binds to LD configurations in supramolecular actin structures assembled in vitro. This antibody allowed us to unveil the LD in different mammalian cells. In particular, we show an association of the antiparallel actin arrangement with the endocytic compartment at the cellular and ultrastructural level. Taken together, our results strongly support a functional role of LD in the patterning of supramolecular actin assemblies in mammalian cells.

A Kirkwood-Buff force field for the aromatic amino acids

In a continuation of our efforts to develop a united atom non-polarizable protein force field based upon the solution theory of Kirkwood and Buff i.e., the Kirkwood-Buff Force Field (KBFF) approach, we present KBFF models for the side chains of phenylalanine, tyrosine, tryptophan, and histidine, including both tautomers of neutral histidine and doubly-protonated histidine. The force fields were specifically designed to reproduce the thermodynamic properties of mixtures over the full composition range in an attempt to provide an improved description of intermolecular interactions. The models were developed by careful parameterization of the solution phase partial charges to reproduce the experimental Kirkwood-Buff integrals for mixtures of solutes representative of the amino acid sidechains in solution. The KBFF parameters and simulated thermodynamic and structural properties are presented for the following eleven binary mixtures: benzene + methanol, benzene + toluene, toluene + methanol, toluene + phenol, toluene + p-cresol, pyrrole + methanol, indole + methanol, pyridine + methanol, pyridine + water, histidine + water, and histidine hydrochloride + water. It is argued that the present approach and models provide a reasonable description of intermolecular interactions which ensures that the required balance between solute-solute, solute-solvent, and solvent-solvent distributions is obtained.

Hybrid Quantum and Classical Simulations of the Dihydrofolate Reductase Catalyzed Hydride Transfer Reaction on an Accurate Semi-Empirical Potential Energy Surface

Dihydrofolate reductase (DHFR) catalyzes the reduction of 7,8-dihydrofolate by nicotinamide adenine dinucleotide phosphate hydride (NADPH) to form 5,6,7,8-tetrahydrofolate and oxidized nicotinamide. DHFR is a small, flexible, monomeric protein with no metals or SS bonds and serves as one of the enzymes commonly used to examine basic aspects in enzymology. In the current work, we present extensive benchmark calculations for several model reactions in the gas phase that are relevant to the DHFR catalyzed hydride transfer. To this end, we employ G4MP2 and CBS-QB3 ab initio calculations as well as numerous density functional theory methods. Using these results, we develop two specific reaction parameter (SRP) Hamiltonians based on the semiempirical AM1 method. The first generation SRP Hamiltonian does not account for dispersion, while the second generation SRP accounts for dispersion implicitly via the AM1 core-repulsion functions. These SRP semiempirical Hamiltonians are subsequently used in hybrid quantum mechanics/molecular mechanics simulations of the DHFR catalyzed reaction. Finally, kinetic isotope effects are computed using a mass-perturbation-based path-integral approach.

An exact expression to calculate the derivatives of position-dependent observables in molecular simulations with flexible constraints

In this work, we introduce an algorithm to compute the derivatives of physical observables along the constrained subspace when flexible constraints are imposed on the system (i.e., constraints in which the constrained coordinates are fixed to configuration-dependent values). The presented scheme is exact, it does not contain any tunable parameter, and it only requires the calculation and inversion of a sub-block of the Hessian matrix of second derivatives of the function through which the constraints are defined. We also present a practical application to the case in which the sought observables are the Euclidean coordinates of complex molecular systems, and the function whose minimization defines the flexible constraints is the potential energy. Finally, and in order to validate the method, which, as far as we are aware, is the first of its kind in the literature, we compare it to the natural and straightforward finite-differences approach in a toy system and in three molecules of biological relevance: methanol, N-methyl-acetamide and a tri-glycine peptide.

Protonation of key acidic residues is critical for the K⁺-selectivity of the Na/K pump

The sodium-potassium (Na/K) pump is a P-type ATPase that generates Na⁺ and K⁺ concentration gradients across the cell membrane. For each ATP molecule, the pump extrudes three Na⁺ and imports two K⁺ by alternating between outward- and inward-facing conformations that preferentially bind K⁺ or Na⁺, respectively. Remarkably, the selective K⁺ and Na⁺ binding sites share several residues, and how the pump is able to achieve the selectivity required for the functional cycle is unclear. Here, free energy perturbation molecular dynamics (FEP/MD) simulations based on the crystal structures of the Na/K pump in a K⁺-loaded state (E2·P_i) reveal that protonation of the high-field acidic side-chains involved in the binding sites is critical to achieve the proper K⁺ selectivity. This prediction is tested with electrophysiological experiments showing that the selectivity of the E2P state for K⁺ over Na⁺ is affected by extracellular pH.

Structural variants of IFNα preferentially promote antiviral functions

IFNα, a cytokine with multiple functions in innate and adaptive immunity and a potent inhibitor of HIV, exerts antiviral activity, in part, by enhancing apolipoprotein B mRNA-editing enzyme-catalytic polypeptide-like 3 (APOBEC3) family members. Although IFNα therapy is associated with reduced viral burden, this cytokine also mediates immune dysfunction and toxicities. Through detailed mapping of IFNα receptor binding sites, we generated IFNα hybrids and mutants and determined that structural changes in the C-helix alter the ability of IFN to limit retroviral activity. Selective IFNα constructs differentially block HIV replication and their directional magnitude of inhibition correlates with APOBEC3 levels. Importantly, certain mutants exhibited reduced toxicity as reflected by induced indoleamine 2,3-dioxygenase (IDO), suggesting discreet and shared intracellular signaling pathways. Defining IFN structure and function relative to APOBEC and other antiviral genes may enable design of novel IFN-related molecules preserving beneficial antiviral roles while minimizing negative effects.

Domain cooperativity in multidomain proteins: what can we learn from molecular alignment in anisotropic media?

Many proteins have modular design with multiple globular domains connected via flexible linkers. As a simple model of such system, we study a tandem construct consisting of two identical SH3 domains and a variable-length Gly/Ser linker. When the linker is short, this construct represents a dumbbell-shaped molecule with limited amount of domain-domain mobility. Due to its elongated shape, this molecule efficiently aligns in steric alignment media. As the length of the linker increases, the two domains become effectively uncoupled and begin to behave as independent entities. Consequently, their degree of alignment drops, approaching that found in the (near-spherical) isolated SH3 domains. To model the dependence of alignment parameters on the length of the interdomain linker, we have generated in silico a series of conformational ensembles representing SH3 tandems with different linker length. These ensembles were subsequently used as input for alignment prediction software PALES. The predicted alignment tensors were compared with the results of experimental measurements using a series of tandem-SH3 samples in PEG/hexanol alignment media. This comparison broadly confirmed the expected trends. At the same time, it has been found that the isolated SH3 domain aligns much stronger than expected. This finding can be attributed to complex morphology of the PEG/hexanol media and/or to weak site-specific interactions between the protein and the media. In the latter case, there are strong indications that electrostatic interactions may play a role. The fact that PEG/hexanol does not behave as a simple steric media should serve as a caution for studies that use PALES as a quantitative prediction tool (especially for disordered proteins). Further progress in this area depends on our ability to accurately model the anisotropic media and its site-specific interactions with protein molecules. Once this ability is improved, it should be possible to use the alignment parameters as a measure of domain-domain cooperativity, thus identifying the situations where two domains transiently interact with each other or become coupled through a partially structured linker.

Scoring and lessons learned with the CSAR benchmark using an improved iterative knowledge-based scoring function

Based on a statistical mechanics-based iterative method, we have extracted a set of distance-dependent, all-atom pairwise potentials for protein-ligand interactions from the crystal structures of 1300 protein-ligand complexes. The iterative method circumvents the long-standing reference state problem in knowledge-based scoring functions. The resulted scoring function, referred to as ITScore 2.0, has been tested with the CSAR (Community Structure-Activity Resource, 2009 release) benchmark of 345 diverse protein-ligand complexes. ITScore 2.0 achieved a Pearson correlation of R(2) = 0.54 in binding affinity prediction. A comparative analysis has been done on the scoring performances of ITScore 2.0, the van der Waals (VDW) scoring function, the VDW with heavy atoms only, and the force field (FF) scoring function of DOCK which consists of a VDW term and an electrostatic term. The results reveal several important factors that affect the scoring performances, which could be helpful for the improvement of scoring functions.