CHARMM: the biomolecular simulation program

Theoretical studies of the base pair fidelity of selenium-modified DNA

The introduction of selenium into DNA in the place of oxygen provides a unique opportunity for studying the fidelity of DNA replication, as well as providing advantages in the growth of DNA crystals and the greater resolution of their structures. However, the atomic mechanisms of the relative stability and base pair recognition of the selenium-modified DNA are poorly understood. In the present study, quantum mechanics calculations were performed on base pairings, base stacking, and base-water interactions for both unmodified thymine and thymine with the 2-exo-oxygen replaced with selenium, and the results were used to develop and validate CHARMM force field parameters for the 2-Se-thymine. Subsequently, molecular dynamics simulations and free-energy perturbation calculations were performed on 11-base DNA sequences containing native thymine and the 2-Se-thymine. The calculated relative free-energy values are in good agreement with experimentally determined relative stability, where the 2-Se-thymine offers similar stability to T-A in cognate DNA, while it dramatically destabilizes the DNA containing the T-G mismatch base pair when 2-Se-thymine is incorporated. Thus, 2-Se-thymine largely increases the native T-A base pair fidelity by discouraging the T-G wobble pair. Insights into the high pairing specificity and the relative stability of selenium-modified DNA were obtained based on detailed structural and energetic analysis of molecular dynamics trajectories. Our studies move one step further toward an understanding of the high base pair fidelity and thermodynamic properties of Se-DNA in solution and in protein-DNA complexes.

Importance of domain closure for the autoactivation of ERK2

Extracellular signal-regulated kinases 1 and 2 (ERK1 and -2, respectively) play a critical role in regulating cell division and have been implicated in cancer. In addition to activation by MAPK/ERK kinases 1 and 2 (MEK1 and -2, respectively), certain mutants of ERK2 can be activated by autophosphorylation. To identify the mechanism of autoactivation, we have performed a series of molecular dynamics simulations of ERK1 and -2 in various stages of activation as well as the constitutively active Q103A, I84A, L73P, and R65S ERK2 mutants. Our simulations indicate the importance of domain closure for autoactivation and activity regulation, with that event occurring prior to folding of the activation lip and of loop L16. Results indicate that the second phosphorylation event, that of T183, disrupts hydrogen bonding involving D334, thereby allowing the kinase to lock into the active conformation. On the basis of the simulations, three predictions were made. G83A was suggested to impede activation; K162M was suggested to perturb the interface between the N- and C-domains leading to activation, and Q64C was hypothesized to stop folding of loop L16, thereby perturbing the homodimerization interface. Functional analysis of the mutants validated the predictions concerning the G83A and Q64C mutants. The K162M mutant did not autoactivate as predicted, however, which may be due to the location of the residue on the protein surface near the ED substrate docking domain.

Motion of a disordered polypeptide chain as studied by paramagnetic relaxation enhancements, 15N relaxation, and molecular dynamics simulations: how fast is segmental diffusion in denatured ubiquitin?

Molecular dynamics (MD) simulations have been widely used to analyze dynamic conformational equilibria of folded proteins, especially in relation to NMR observables. However, this approach found little use in the studies of disordered proteins, where the sampling of vast conformational space presents a serious problem. In this paper, we demonstrate that the latest advances in computation technology make it possible to overcome this limitation. The experimentally validated (calibrated) MD models allow for new insights into structure/dynamics of disordered proteins. As a test system, we have chosen denatured ubiquitin in solution with 8 M urea at pH 2. High-temperature MD simulations in implicit solvent have been carried out for the wild-type ubiquitin as well as MTSL-tagged Q2C, D32C, and R74C mutants. To recalibrate the MD data (500 K) in relation to the experimental conditions (278 K, 8 M urea), the time axes of the MD trajectories were rescaled. The scaling factor was adjusted such as to maximize the agreement between the simulated and experimental (15)N relaxation rates. The resulting effective length of the trajectories, 311 μs, ensures good convergence properties of the MD model. The constructed MD model was validated against the array of experimental data, including additional (15)N relaxation parameters, multiple sets of paramagnetic relaxation enhancements (PREs), and the radius of gyration. In each case, a near-quantitative agreement has been obtained, suggesting that the model is successful. Of note, the MD-based approach rigorously predicts the quantities that are inherently dynamic, i.e., dependent on the motional correlation times. This cannot be accomplished, other than in empirical fashion, on the basis of static structural models (conformational ensembles). The MD model was further used to investigate the relative translational motion of the MTSL label and the individual H(N) atoms. The derived segmental diffusion coefficients proved to be nearly uniform along the peptide chain, averaging to D = 0.49-0.55 × 10(-6) cm(2)/s. This result was verified by direct analysis of the experimental PRE data using the recently proposed Ullman-Podkorytov model. In this model, MTSL and H(N) moieties are treated as two tethered spheres undergoing mutual diffusion in a harmonic potential. The fitting of the experimental data involving D as a single adjustable parameter leads to D = 0.45 × 10(-6) cm(2)/s, in good agreement with the MD-based analyses. This result can be compared with the range of estimates obtained from the resonance energy transfer experiments, D = 0.2-6.0 × 10(-6) cm(2)/s.

Chemical synthesis, crystal structure and enzymatic evaluation of a dinucleotide spore photoproduct analogue containing a formacetal linker

Spore photoproduct (SP) is the exclusive DNA photodamage product found in bacterial endospores. Its photoformation and repair by a metalloenzyme spore photoproduct lyase (SPL) composes the unique SP biochemistry. Despite the fact that the SP was discovered almost 50 years ago, its crystal structure is still unknown and the lack of structural information greatly hinders the study of SP biochemistry. Employing a formacetal linker and organic synthesis, we successfully prepared a dinucleotide SP isostere 5R-CH(2) SP, which contains a neutral CH(2) moiety between the two thymine residues instead of a phosphate. The neutral linker dramatically facilitates the crystallization process, allowing us to obtain the crystal structure for this intriguing thymine dimer half a century after its discovery. Further ROESY spectroscopic, DFT computational, and enzymatic studies of this 5R-CH(2) SP compound prove that it possesses similar properties with the 5R-SP species, suggesting that the revealed structure truly reflects that of SP generated in Nature.

Molecular mechanism of a green-shifted, pH-dependent red fluorescent protein mKate variant

Fluorescent proteins that can switch between distinct colors have contributed significantly to modern biomedical imaging technologies and molecular cell biology. Here we report the identification and biochemical analysis of a green-shifted red fluorescent protein variant GmKate, produced by the introduction of two mutations into mKate. Although the mutations decrease the overall brightness of the protein, GmKate is subject to pH-dependent, reversible green-to-red color conversion. At physiological pH, GmKate absorbs blue light (445 nm) and emits green fluorescence (525 nm). At pH above 9.0, GmKate absorbs 598 nm light and emits 646 nm, far-red fluorescence, similar to its sequence homolog mNeptune. Based on optical spectra and crystal structures of GmKate in its green and red states, the reversible color transition is attributed to the different protonation states of the cis-chromophore, an interpretation that was confirmed by quantum chemical calculations. Crystal structures reveal potential hydrogen bond networks around the chromophore that may facilitate the protonation switch, and indicate a molecular basis for the unusual bathochromic shift observed at high pH. This study provides mechanistic insights into the color tuning of mKate variants, which may aid the development of green-to-red color-convertible fluorescent sensors, and suggests GmKate as a prototype of genetically encoded pH sensors for biological studies.

Multi-Site λ-dynamics for simulated Structure-Activity Relationship studies

Multi-Site λ-dynamics (MSλD) is a new free energy simulation method that is based on λ-dynamics. It has been developed to enable multiple substituents at multiple sites on a common ligand core to be modeled simultaneously and their free energies assessed. The efficacy of MSλD for estimating relative hydration free energies and relative binding affinties is demonstrated using three test systems. Model compounds representing multiple identical benzene, dihydroxybenzene and dimethoxybenzene molecules show total combined MSλD trajectory lengths of ~1.5 ns are sufficient to reliably achieve relative hydration free energy estimates within 0.2 kcal/mol and are less sensitive to the number of trajectories that are used to generate these estimates for hybrid ligands that contain up to ten substituents modeled at a single site or five substituents modeled at each of two sites. Relative hydration free energies among six benzene derivatives calculated from MSλD simulations are in very good agreement with those from alchemical free energy simulations (with average unsigned differences of 0.23 kcal/mol and R(2)=0.991) and experiment (with average unsigned errors of 1.8 kcal/mol and R(2)=0.959). Estimates of the relative binding affinities among 14 inhibitors of HIV-1 reverse transcriptase obtained from MSλD simulations are in reasonable agreement with those from traditional free energy simulations and experiment (average unsigned errors of 0.9 kcal/mol and R(2)=0.402). For the same level of accuracy and precision MSλD simulations are achieved ~20-50 times faster than traditional free energy simulations and thus with reliable force field parameters can be used effectively to screen tens to hundreds of compounds in structure-based drug design applications.

Mammalian testis-determining factor SRY and the enigma of inherited human sex reversal: frustrated induced fit in a bent protein-DNA complex

Mammalian testis-determining factor SRY contains a high mobility group box, a conserved eukaryotic motif of DNA bending. Mutations in SRY cause XY gonadal dysgenesis and somatic sex reversal. Although such mutations usually arise de novo in spermatogenesis, some are inherited and so specify male development in one genetic background (the father) but not another (the daughter). Here, we describe the biophysical properties of a representative inherited mutation, V60L, within the minor wing of the L-shaped domain (box position 5). Although the stability and DNA binding properties of the mutant domain are similar to those of wild type, studies of SRY-induced DNA bending by subnanosecond time-resolved fluorescence resonance energy transfer (FRET) revealed enhanced conformational fluctuations leading to long range variation in bend angle. (1)H NMR studies of the variant protein-DNA complex demonstrated only local perturbations near the mutation site. Because the minor wing of SRY folds on DNA binding, the inherited mutation presumably hinders induced fit. Stopped-flow FRET studies indicated that such frustrated packing leads to accelerated dissociation of the bent complex. Studies of SRY-directed transcriptional regulation in an embryonic gonadal cell line demonstrated partial activation of downstream target Sox9. Our results have demonstrated a nonlocal coupling between DNA-directed protein folding and protein-directed DNA bending. Perturbation of this coupling is associated with a genetic switch poised at the threshold of activity.

Mutant alcohol dehydrogenase leads to improved ethanol tolerance in Clostridium thermocellum

Clostridium thermocellum is a thermophilic, obligately anaerobic, gram-positive bacterium that is a candidate microorganism for converting cellulosic biomass into ethanol through consolidated bioprocessing. Ethanol intolerance is an important metric in terms of process economics, and tolerance has often been described as a complex and likely multigenic trait for which complex gene interactions come into play. Here, we resequence the genome of an ethanol-tolerant mutant, show that the tolerant phenotype is primarily due to a mutated bifunctional acetaldehyde-CoA/alcohol dehydrogenase gene (adhE), hypothesize based on structural analysis that cofactor specificity may be affected, and confirm this hypothesis using enzyme assays. Biochemical assays confirm a complete loss of NADH-dependent activity with concomitant acquisition of NADPH-dependent activity, which likely affects electron flow in the mutant. The simplicity of the genetic basis for the ethanol-tolerant phenotype observed here informs rational engineering of mutant microbial strains for cellulosic ethanol production.

Influence of water-protein hydrogen bonding on the stability of Trp-cage miniprotein. A comparison between the TIP3P and TIP4P-Ew water models

We report extensive replica exchange molecular dynamics (REMD) simulations on the folding/unfolding equilibrium of Trp-cage miniprotein using the Amber ff99SB all atom forcefield and TIP3P and TIP4P-Ew explicit water solvent models. REMD simulation-lengths in the 500 ns to the microsecond regime per replica are required to adequately sample the folding/unfolding equilibrium. We observe that this equilibrium is significantly affected by the choice of the water model. Compared with experimental data, simulations using the TIP3P solvent describe the stability of the Trp-cage quite realistically, providing a melting point which is just a few Kelvins above the experimental transition temperature of 317 K. The TIP4P-Ew model shifts the equilibrium towards the unfolded state and lowers the free energy of unfolding by about 3 kJ mol(-1) at 280 K, demonstrating the need to fine-tune the protein-forcefield depending on the chosen water model. We report evidence that the main difference between the two water models is mostly due to the different solvation of polar groups of the peptide. The unfolded state of the Trp-cage is stabilized by an increasing number of hydrogen bonds, destabilizing the α-helical part of the molecule and opening the R-D salt bridge. By reweighting the strength of solvent-peptide hydrogen bonds by adding a hydrogen bond square well potential, we can fully recover the effect of the different water models and estimate the shift in population as due to a difference in hydrogen bond-strength of about 0.4 kJ mol(-1) per hydrogen bond.

MEMO associated with an ErbB2 receptor phosphopeptide reveals a new phosphotyrosine motif

Tyrosine phosphorylations are essential in signal transduction. Recently, a new type of phosphotyrosine binding protein, MEMO (Mediator of ErbB2-driven cell motility), has been reported to bind specifically to an ErbB2-derived phosphorylated peptide encompassing Tyr-1227 (PYD). Structural and functional analyses of variants of this peptide revealed the minimum sequence required for MEMO recognition. Using a docking approach we have generated a structural model for MEMO/PYD complex and compare this new phosphotyrosine motif to SH2 and PTB phosphotyrosine motives.

An ensemble dynamics approach to decipher solid-state NMR observables of membrane proteins

Solid-state NMR (SSNMR) is an invaluable tool for determining orientations of membrane proteins and peptides in lipid bilayers. Such orientational descriptions provide essential information about membrane protein functions. However, when a semi-static single conformer model is used to interpret various SSNMR observables, important dynamics information can be missing, and, sometimes, even orientational information can be misinterpreted. In addition, over the last decade, molecular dynamics (MD) simulation and semi-static SSNMR interpretation have shown certain levels of discrepancies in terms of transmembrane helix orientation and dynamics. Dynamic fitting models have recently been proposed to resolve these discrepancies by taking into account transmembrane helix whole body motions using additional parameters. As an alternative approach, we have developed SSNMR ensemble dynamics (SSNMR-ED) using multiple conformer models, which generates an ensemble of structures that satisfies the experimental observables without any fitting parameters. In this review, various computational methods for determining transmembrane helix orientations are discussed, and the distributions of VpuTM (from HIV-1) and WALP23 (a synthetic peptide) orientations from SSNMR-ED simulations are compared with those from MD simulations and semi-static/dynamic fitting models. Such comparisons illustrate that SSNMR-ED can be used as a general means to extract both membrane protein structure and dynamics from the SSNMR measurements. This article is part of a Special Issue entitled: Membrane protein structure and function.

A model of calcium activation of the cardiac thin filament

The cardiac thin filament regulates actomyosin interactions through calcium-dependent alterations in the dynamics of cardiac troponin and tropomyosin. Over the past several decades, many details of the structure and function of the cardiac thin filament and its components have been elucidated. We propose a dynamic, complete model of the thin filament that encompasses known structures of cardiac troponin, tropomyosin, and actin and show that it is able to capture key experimental findings. By performing molecular dynamics simulations under two conditions, one with calcium bound and the other without calcium bound to site II of cardiac troponin C (cTnC), we found that subtle changes in structure and protein contacts within cardiac troponin resulted in sweeping changes throughout the complex that alter tropomyosin (Tm) dynamics and cardiac troponin--actin interactions. Significant calcium-dependent changes in dynamics occur throughout the cardiac troponin complex, resulting from the combination of the following: structural changes in the N-lobe of cTnC at and adjacent to sites I and II and the link between them; secondary structural changes of the cardiac troponin I (cTnI) switch peptide, of the mobile domain, and in the vicinity of residue 25 of the N-terminus; secondary structural changes in the cardiac troponin T (cTnT) linker and Tm-binding regions; and small changes in cTnC-cTnI and cTnT-Tm contacts. As a result of these changes, we observe large changes in the dynamics of the following regions: the N-lobe of cTnC, the mobile domain of cTnI, the I-T arm, the cTnT linker, and overlapping Tm. Our model demonstrates a comprehensive mechanism for calcium activation of the cardiac thin filament consistent with previous, independent experimental findings. This model provides a valuable tool for research into the normal physiology of cardiac myofilaments and a template for studying cardiac thin filament mutations that cause human cardiomyopathies.

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.

Selection of herbal therapeutics against deltatoxin mediated Clostridial infections

Clostridium perfringens (a versatile pathogenic bacterium) secretes enterotoxins (the deltatoxin, virulent factor) and causes food borne gastroenteritis and gasgangrene. The organism was isolated and characterized from improperly cooked meat and poultry samples. The isolated organism showed multiple drug resistance indicating that the treatment is challenging. Hence, there is need for improved therapeutic agents. The rational design of improved therapeutics requires the crystal structure for the toxin. However, the structure for the toxin is not yet available in its native form. Thus, we modeled the toxin structure using α- hemolysin of Staphylococcus aureus (PDB: 3M4D chain A) as template. The docking of the toxin with the herbal extract curcumin (1,7-bis(4-hydroxy-3- methoxyphenyl)hepta-1,6-diene-3,5-dione) showed a binding energy of -8.6 Kcal/mol, in comparison to the known antibiotic Linezolid with binding energy of -6.1 Kcal/mol. This data finds application in the design and development of novel compounds against the deltatoxin from Clostridium perfringens.

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.

A missense mutation in CLIC2 associated with intellectual disability is predicted by in silico modeling to affect protein stability and dynamics

Large-scale next generation resequencing of X chromosome genes identified a missense mutation in the CLIC2 gene on Xq28 in a male with X-linked intellectual disability (XLID) and not found in healthy individuals. At the same time, numerous nsSNPs (nonsynonomous SNP) have been reported in the CLIC2 gene in healthy individuals indicating that the CLIC2 protein can tolerate amino acid substitutions and be fully functional. To test the possibility that p.H101Q is a disease-causing mutation, we performed in silico simulations to calculate the effects of the p.H101Q mutation on CLIC2 stability, dynamics, and ionization states while comparing the effects obtained for presumably harmless nsSNPs. It was found that p.H101Q, in contrast with other nsSNPs, (a) lessens the flexibility of the joint loop which is important for the normal function of CLIC2, (b) makes the overall 3D structure of CLIC2 more stable and thus reduces the possibility of the large conformational change expected to occur when CLIC2 moves from a soluble to membrane form, and (c) removes the positively charged residue, H101, which may be important for the membrane association of CLIC2. The results of in silico modeling, in conjunction with the polymorphism analysis, suggest that p.H101Q may be a disease-causing mutation, the first one suggested in the CLIC family.

Structural and biophysical properties of a synthetic channel-forming peptide: designing a clinically relevant anion selective pore

The design, synthesis, modeling and in vitro testing of channel-forming peptides derived from the cys-loop superfamily of ligand-gated ion channels are part of an ongoing research focus. Over 300 different sequences have been prepared based on the M2 transmembrane segment of the spinal cord glycine receptor α-subunit. A number of these sequences are water-soluble monomers that readily insert into biological membranes where they undergo supramolecular assembly, yielding channels with a range of selectivities and conductances. Selection of a sequence for further modifications to yield an optimal lead compound came down to a few key biophysical properties: low solution concentrations that yield channel activity, greater ensemble conductance, and enhanced ion selectivity. The sequence NK(4)-M2GlyR T19R, S22W (KKKKPARVGLGITTVLTMRTQW) addressed these criteria. The structure of this peptide has been analyzed by solution NMR as a monomer in detergent micelles, simulated as five-helix bundles in a membrane environment, modified by cysteine-scanning and studied for insertion efficiency in liposomes of selected lipid compositions. Taken together, these results define the structural and key biophysical properties of this sequence in a membrane. This model provides an initial scaffold from which rational substitutions can be proposed and tested to modulate anion selectivity. This article is part of a Special Issue entitled: Protein Folding in Membranes.

MDAnalysis: a toolkit for the analysis of molecular dynamics simulations

MDAnalysis is an object-oriented library for structural and temporal analysis of molecular dynamics (MD) simulation trajectories and individual protein structures. It is written in the Python language with some performance-critical code in C. It uses the powerful NumPy package to expose trajectory data as fast and efficient NumPy arrays. It has been tested on systems of millions of particles. Many common file formats of simulation packages including CHARMM, Gromacs, Amber, and NAMD and the Protein Data Bank format can be read and written. Atoms can be selected with a syntax similar to CHARMM's powerful selection commands. MDAnalysis enables both novice and experienced programmers to rapidly write their own analytical tools and access data stored in trajectories in an easily accessible manner that facilitates interactive explorative analysis. MDAnalysis has been tested on and works for most Unix-based platforms such as Linux and Mac OS X. It is freely available under the GNU General Public License from http://mdanalysis.googlecode.com.

Integrating diffusion maps with umbrella sampling: application to alanine dipeptide

Nonlinear dimensionality reduction techniques can be applied to molecular simulation trajectories to systematically extract a small number of variables with which to parametrize the important dynamical motions of the system. For molecular systems exhibiting free energy barriers exceeding a few k(B)T, inadequate sampling of the barrier regions between stable or metastable basins can lead to a poor global characterization of the free energy landscape. We present an adaptation of a nonlinear dimensionality reduction technique known as the diffusion map that extends its applicability to biased umbrella sampling simulation trajectories in which restraining potentials are employed to drive the system into high free energy regions and improve sampling of phase space. We then propose a bootstrapped approach to iteratively discover good low-dimensional parametrizations by interleaving successive rounds of umbrella sampling and diffusion mapping, and we illustrate the technique through a study of alanine dipeptide in explicit solvent.

Toward canonical ensemble distribution from self-guided Langevin dynamics simulation

This work derives a quantitative description of the conformational distribution in self-guided Langevin dynamics (SGLD) simulations. SGLD simulations employ guiding forces calculated from local average momentums to enhance low-frequency motion. This enhancement in low-frequency motion dramatically accelerates conformational search efficiency, but also induces certain perturbations in conformational distribution. Through the local averaging, we separate properties of molecular systems into low-frequency and high-frequency portions. The guiding force effect on the conformational distribution is quantitatively described using these low-frequency and high-frequency properties. This quantitative relation provides a way to convert between a canonical ensemble and a self-guided ensemble. Using example systems, we demonstrated how to utilize the relation to obtain canonical ensemble properties and conformational distributions from SGLD simulations. This development makes SGLD not only an efficient approach for conformational searching, but also an accurate means for conformational sampling.

Force-induced unfolding simulations of the human Notch1 negative regulatory region: possible roles of the heterodimerization domain in mechanosensing

Notch receptors are core components of the Notch signaling pathway and play a central role in cell fate decisions during development as well as tissue homeostasis. Upon ligand binding, Notch is sequentially cleaved at the S2 site by an ADAM protease and at the S3 site by the γ-secretase complex. Recent X-ray structures of the negative regulatory region (NRR) of the Notch receptor reveal an auto-inhibited fold where three protective Lin12/Notch repeats (LNR) of the NRR shield the S2 cleavage site housed in the heterodimerization (HD) domain. One of the models explaining how ligand binding drives the NRR conformation from a protease-resistant state to a protease-sensitive one invokes a mechanical force exerted on the NRR upon ligand endocytosis. Here, we combined physics-based atomistic simulations and topology-based coarse-grained modeling to investigate the intrinsic and force-induced folding and unfolding mechanisms of the human Notch1 NRR. The simulations support that external force applied to the termini of the NRR disengages the LNR modules from the heterodimerization (HD) domain in a well-defined, largely sequential manner. Importantly, the mechanical force can further drive local unfolding of the HD domain in a functionally relevant fashion that would provide full proteolytic access to the S2 site prior to heterodimer disassociation. We further analyzed local structural features, intrinsic folding free energy surfaces, and correlated motions of the HD domain. The results are consistent with a model in which the HD domain possesses inherent mechanosensing characteristics that could be utilized during Notch activation. This potential role of the HD domain in ligand-dependent Notch activation may have implications for understanding normal and aberrant Notch signaling.

Lagrangian formulation with dissipation of Born-Oppenheimer molecular dynamics using the density-functional tight-binding method

An important element determining the time requirements of Born-Oppenheimer molecular dynamics (BOMD) is the convergence rate of the self-consistent solution of Roothaan equations (SCF). We show here that improved convergence and dynamics stability can be achieved by use of a Lagrangian formalism of BOMD with dissipation (DXL-BOMD). In the DXL-BOMD algorithm, an auxiliary electronic variable (e.g., the electron density or Fock matrix) is propagated and a dissipative force is added in the propagation to maintain the stability of the dynamics. Implementation of the approach in the self-consistent charge density functional tight-binding method makes possible simulations that are several hundred picoseconds in lengths, in contrast to earlier DFT-based BOMD calculations, which have been limited to tens of picoseconds or less. The increase in the simulation time results in a more meaningful evaluation of the DXL-BOMD method. A comparison is made of the number of iterations (and time) required for convergence of the SCF with DXL-BOMD and a standard method (starting with a zero charge guess for all atoms at each step), which gives accurate propagation with reasonable SCF convergence criteria. From tests using NVE simulations of C(2)F(4) and 20 neutral amino acid molecules in the gas phase, it is found that DXL-BOMD can improve SCF convergence by up to a factor of two over the standard method. Corresponding results are obtained in simulations of 32 water molecules in a periodic box. Linear response theory is used to analyze the relationship between the energy drift and the correlation of geometry propagation errors.

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.

G117C MelB, a mutant melibiose permease with a changed conformational equilibrium

Replacement of the glycine at position 117 by a cysteine in the melibiose permease creates an interesting phenotype: while the mutant transporter shows still transport activity comparable to the wild type its pre steady-state kinetic properties are drastically altered. The transient charge displacements after substrate concentration jumps are strongly reduced and the fluorescence changes disappear. Together with its maintained transport activity this indicates that substrate translocation in G117C melibiose permease is not impaired but that the initial conformation of the mutant transporter differs from that of the wild type permease. A kinetic model for the G117C melibiose permease based on a rapid dynamic equilibrium of the substrate free transporter is proposed. Implications of the kinetic model for the transport mechanism of the wild type permease are discussed.

Molecular dynamics simulation and binding energy calculation for estimation of oligonucleotide duplex thermostability in RNA-based therapeutics

For oligonucleotide-based therapeutics, a thorough understanding of the thermodynamic properties of duplex formation is critical to developing stable and potent drugs. For unmodified small interfering RNA (siRNA), DNA antisense oligonucleotide (AON) and locked nucleic acid (LNA), DNA/LNA modified oligonucleotides, nearest neighbor (NN) methods can be effectively used to quickly and accurately predict duplex thermodynamic properties such as melting point. Unfortunately, for chemically modified olignonucleotides, there has been no accurate prediction method available. Here we describe the potential of estimating melting temperature (T(m)) for nonstandard oligonucleotides by using the correlation of the experimental T(m) with the calculated duplex binding energy (BE) for oligonucleotides of a given length. This method has been automated into a standardized molecular dynamics (MD) protocol through Pipeline Pilot (PP) using the CHARMm component in Discovery Studio (DS). Results will be presented showing the correlation of the predicted data with experiment for both standard and chemically modified siRNA and AON.

Synergistic folding of two intrinsically disordered proteins: searching for conformational selection

Intrinsically disordered proteins (IDPs) lack stable structures under physiological conditions but often fold into stable structures upon specific binding. These coupled binding and folding processes underlie the organization of cellular regulatory networks, and a mechanistic understanding is thus of fundamental importance. Here, we investigated the synergistic folding of two IDPs, namely, the NCBD domain of transcription coactivator CBP and the p160 steroid receptor coactivator ACTR, using a topology-based model that was carefully calibrated to balance intrinsic folding propensities and intermolecular interactions. As one of the most structured IDPs, NCBD is a plausible candidate that interacts through conformational selection-like mechanisms, where binding is mainly initiated by pre-existing folded-like conformations. Indeed, the simulations demonstrate that, even though binding and folding of both NCBD and ACTR is highly cooperative on the baseline level, the tertiary folding of NCBD is best described by the "extended conformational selection" model that involves multiple stages of selection and induced folding. The simulations further predict that the NCBD/ACTR recognition is mainly initiated by forming a mini folded core that includes the second and third helices of NCBD and ACTR. These predictions are fully consistent with independent physics-based atomistic simulations as well as a recent experimental mapping of the H/D exchange protection factors. The current work thus adds to the limited number of existing mechanistic studies of coupled binding and folding of IDPs, and provides a first direct demonstration of how conformational selection might contribute to efficient recognition of IDPs. Interestingly, even for highly structured IDPs like NCBD, the recognition is initiated by the more disordered C-terminal segment and with substantial contribution from induced folding. Together with existing studies of IDP interaction mechanisms, this argues that induced folding is likely prevalent in IDP-protein interaction, and emphasizes the importance of understanding how IDPs manage to fold efficiently upon (nonspecific) binding. Success of the current study also further supports the notion that, with careful calibration, topology-based models can be effective tools for mechanistic study of IDP interaction and regulation, especially when combined with physics-based atomistic simulations and experiments.

Impact of 2'-hydroxyl sampling on the conformational properties of RNA: update of the CHARMM all-atom additive force field for RNA

Here, we present an update of the CHARMM27 all-atom additive force field for nucleic acids that improves the treatment of RNA molecules. The original CHARMM27 force field parameters exhibit enhanced Watson-Crick base pair opening which is not consistent with experiment, whereas analysis of molecular dynamics (MD) simulations show the 2'-hydroxyl moiety to almost exclusively sample the O3' orientation. Quantum mechanical (QM) studies of RNA related model compounds indicate the energy minimum associated with the O3' orientation to be too favorable, consistent with the MD results. Optimization of the dihedral parameters dictating the energy of the 2'-hydroxyl proton targeting the QM data yielded several parameter sets, which sample both the base and O3' orientations of the 2'-hydroxyl to varying degrees. Selection of the final dihedral parameters was based on reproduction of hydration behavior as related to a survey of crystallographic data and better agreement with experimental NMR J-coupling values. Application of the model, designated CHARMM36, to a collection of canonical and noncanonical RNA molecules reveals overall improved agreement with a range of experimental observables as compared to CHARMM27. The results also indicate the sensitivity of the conformational heterogeneity of RNA to the orientation of the 2'-hydroxyl moiety and support a model whereby the 2'-hydroxyl can enhance the probability of conformational transitions in RNA.

Interconversion of functional motions between mesophilic and thermophilic adenylate kinases

Dynamic properties are functionally important in many proteins, including the enzyme adenylate kinase (AK), for which the open/closed transition limits the rate of catalytic turnover. Here, we compare our previously published coarse-grained (double-well Gō) simulation of mesophilic AK from E. coli (AKmeso) to simulations of thermophilic AK from Aquifex aeolicus (AKthermo). In AKthermo, as with AKmeso, the LID domain prefers to close before the NMP domain in the presence of ligand, but LID rigid-body flexibility in the open (O) ensemble decreases significantly. Backbone foldedness in O and/or transition state (TS) ensembles increases significantly relative to AKmeso in some interdomain backbone hinges and within LID. In contact space, the TS of AKthermo has fewer contacts at the CORE-LID interface but a stronger contact network surrounding the CORE-NMP interface than the TS of AKmeso. A “heated” simulation of AKthermo at 375K slightly increases LID rigid-body flexibility in accordance with the “corresponding states” hypothesis. Furthermore, while computational mutation of 7 prolines in AKthermo to their AKmeso counterparts produces similar small perturbations, mutation of these sites, especially positions 8 and 155, to glycine is required to achieve LID rigid-body flexibility and hinge flexibilities comparable to AKmeso. Mutating the 7 sites to proline in AKmeso reduces some hinges' flexibilities, especially hinge 2, but does not reduce LID rigid-body flexibility, suggesting that these two types of motion are decoupled in AKmeso. In conclusion, our results suggest that hinge flexibility and global functional motions alike are correlated with but not exclusively determined by the hinge residues. This mutational framework can inform the rational design of functionally important flexibility and allostery in other proteins toward engineering novel biochemical pathways.

Dynamic properties are functionally important in many proteins, including the enzyme adenylate kinase (AK), which undergoes chemically rate-limiting domain motions coupled to substrate binding. Since mesophiles and thermophiles often differ in functionally important motions, we compare coarse-grained simulations of AKmeso and AKthermo as well as several proline and glycine mutational variants designed to interconvert the dynamics. As might be expected, both domain motions and local unfolding motions are reduced in AKthermo relative to AKmeso. In AKthermo, both of these types of motions can be partially shifted toward more flexible AKmeso by heating or by mutating hinge prolines. However, only mutation to highly flexible glycine produces motions like those of AKmeso. Thus, the rate-limiting global transition likely depends on a combination of hinge flexibility and stability within the LID and NMP domains. Finally, this mutagenic framework can inform the rational design of flexibility and allostery in other proteins toward engineering novel biological control systems.

A Grid-enabled web portal for NMR structure refinement with AMBER

MOTIVATION

The typical workflow for NMR structure determination involves collecting thousands of conformational restraints, calculating a bundle of 20-40 conformers in agreement with them and refining the energetics of these conformers. The structure calculation step employs simulated annealing based on molecular dynamics (MD) simulations with very simplified force fields. The value of refining the calculated conformers using restrained MD (rMD) simulations with state-of-art force fields is documented. This refinement however presents various subtleties, from the proper formatting of conformational restraints to the definition of suitable protocols.

RESULTS

We describe a web interface to set up and run calculations with the AMBER package, which we called AMPS-NMR (AMBER-based Portal Server for NMR structures). The interface allows the refinement of NMR structures through rMD. Some predefined protocols are provided for this purpose, which can be personalized; it is also possible to create an entirely new protocol. AMPS-NMR can handle various restraint types. Standard rMD refinement in explicit water of the structures of three different proteins are shown as examples. AMPS-NMR additionally includes a workspace for the user to store different calculations. As an ancillary service, a web interface to AnteChamber is available, enabling the calculation of force field parameters for organic molecules such as ligands in protein-ligand adducts.

AVAILABILITY AND IMPLEMENTATION

AMPS-NMR is embedded within the NMR services of the WeNMR project and is available at http://py-enmr.cerm.unifi.it/access/index/amps-nmr; its use requires registration with a digital certificate.

CONTACT

ivanobertini@cerm.unifi.it

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Weak long-range correlated motions in a surface patch of ubiquitin involved in molecular recognition

Long-range correlated motions in proteins are candidate mechanisms for processes that require information transfer across protein structures, such as allostery and signal transduction. However, the observation of backbone correlations between distant residues has remained elusive, and only local correlations have been revealed using residual dipolar couplings measured by NMR spectroscopy. In this work, we experimentally identified and characterized collective motions spanning four β-strands separated by up to 15 Å in ubiquitin. The observed correlations link molecular recognition sites and result from concerted conformational changes that are in part mediated by the hydrogen-bonding network.

Toward accurate prediction of pKa values for internal protein residues: the importance of conformational relaxation and desolvation energy

Proton uptake or release controls many important biological processes, such as energy transduction, virus replication, and catalysis. Accurate pK(a) prediction informs about proton pathways, thereby revealing detailed acid-base mechanisms. Physics-based methods in the framework of molecular dynamics simulations not only offer pK(a) predictions but also inform about the physical origins of pK(a) shifts and provide details of ionization-induced conformational relaxation and large-scale transitions. One such method is the recently developed continuous constant pH molecular dynamics (CPHMD) method, which has been shown to be an accurate and robust pK(a) prediction tool for naturally occurring titratable residues. To further examine the accuracy and limitations of CPHMD, we blindly predicted the pK(a) values for 87 titratable residues introduced in various hydrophobic regions of staphylococcal nuclease and variants. The predictions gave a root-mean-square deviation of 1.69 pK units from experiment, and there were only two pK(a)'s with errors greater than 3.5 pK units. Analysis of the conformational fluctuation of titrating side-chains in the context of the errors of calculated pK(a) values indicate that explicit treatment of conformational flexibility and the associated dielectric relaxation gives CPHMD a distinct advantage. Analysis of the sources of errors suggests that more accurate pK(a) predictions can be obtained for the most deeply buried residues by improving the accuracy in calculating desolvation energies. Furthermore, it is found that the generalized Born implicit-solvent model underlying the current CPHMD implementation slightly distorts the local conformational environment such that the inclusion of an explicit-solvent representation may offer improvement of accuracy.

Computational exploration of structural hypotheses for an additional sequence in a mammalian mitochondrial protein

BACKGROUND

Proteins involved in mammalian mitochondrial translation, when compared to analogous bacterial proteins, frequently have additional sequence regions whose structural or functional roles are not always clear. For example, an additional short insert sequence in the bovine mitochondrial initiation factor 2 (IF2(mt)) seems sufficient to fulfill the added role of eubacterial initiation factor IF1. Prior to our recent cryo-EM study that showed IF2(mt) to structurally occupy both the IF1 and IF2 binding sites, the spatial separation of these sites, and the short length of the insert sequence, posed ambiguity in whether it could perform the role of IF1 through occupation of the IF1 binding site on the ribosome.

RESULTS

The present study probes how well computational structure prediction methods can a priori address hypothesized roles of such additional sequences by creating quasi-atomic models of IF2(mt) using bacterial IF2 cryo-EM densities (that lack the insert sequences). How such initial IF2(mt) predictions differ from the observed IF2(mt) cryo-EM map and how they can be suitably improved using further sequence analysis and flexible fitting are analyzed.

CONCLUSIONS

By hypothesizing that the insert sequence occupies the IF1 binding site, continuous IF2(mt) models that occupy both the IF2 and IF1 binding sites can be predicted computationally. These models can be improved by flexible fitting into the IF2(mt) cryo-EM map to get reasonable quasi-atomic IF2(mt) models, but the exact orientation of the insert structure may not be reproduced. Specific eukaryotic insert sequence conservation characteristics can be used to predict alternate IF2(mt) models that have minor secondary structure rearrangements but fewer unusually extended linker regions. Computational structure prediction methods can thus be combined with medium-resolution cryo-EM maps to explore structure-function hypotheses for additional sequence regions and to guide further biochemical experiments, especially in mammalian systems where high-resolution structures are difficult to determine.

Free Energy Profiles along Consensus Normal Modes Provide Insight into HIV-1 Protease Flap Opening

Describing biological macromolecular energetics from computer simulations can pose major challenges, and often necessitates enhanced conformational sampling. We describe the calculation of conformational free-energy profiles along carefully chosen collective coordinates: "consensus" normal modes, developed recently as robust alternatives to conventional normal modes. In an application to the HIV-1 protease, we obtain efficient sampling of significant flap opening movements governing inhibitor binding from relatively short simulations, in close correspondence with experimental results.

A conserved asparagine residue in transmembrane segment 1 (TM1) of serotonin transporter dictates chloride-coupled neurotransmitter transport

Na(+)- and Cl(-)-dependent uptake of neurotransmitters via transporters of the SLC6 family, including the human serotonin transporter (SLC6A4), is critical for efficient synaptic transmission. Although residues in the human serotonin transporter involved in direct Cl(-) coordination of human serotonin transport have been identified, the role of Cl(-) in the transport mechanism remains unclear. Through a combination of mutagenesis, chemical modification, substrate and charge flux measurements, and molecular modeling studies, we reveal an unexpected role for the highly conserved transmembrane segment 1 residue Asn-101 in coupling Cl(-) binding to concentrative neurotransmitter uptake.

Insight into the phosphodiesterase mechanism from combined QM/MM free energy simulations

Molecular dynamics simulations employing a combined quantum mechanical and molecular mechanical potential have been carried out to elucidate the reaction mechanism of the hydrolysis of a cyclic nucleotide cAMP substrate by phosphodiesterase 4B (PDE4B). PDE4B is a member of the PDE superfamily of enzymes that play crucial roles in cellular signal transduction. We have determined a two-dimensional potential of mean force (PMF) for the coupled phosphoryl bond cleavage and proton transfer through a general acid catalysis mechanism in PDE4B. The results indicate that the ring-opening process takes place through an S(N)2 reaction mechanism, followed by a proton transfer to stabilize the leaving group. The computed free energy of activation for the PDE4B-catalyzed cAMP hydrolysis is about 13 kcal·mol(-1) and an overall reaction free energy is about -17 kcal·mol(-1), both in accord with experimental results. In comparison with the uncatalyzed reaction in water, the enzyme PDE4B provides a strong stabilization of the transition state, lowering the free energy barrier by 14 kcal·mol(-1). We found that the proton transfer from the general acid residue His234 to the O3' oxyanion of the ribosyl leaving group lags behind the nucleophilic attack, resulting in a shallow minimum on the free energy surface. A key contributing factor to transition state stabilization is the elongation of the distance between the divalent metal ions Zn(2+) and Mg(2+) in the active site as the reaction proceeds from the Michaelis complex to the transition state.

Continuous Constant pH Molecular Dynamics in Explicit Solvent with pH-Based Replica Exchange

A computational tool that offers accurate pKa values and atomically detailed knowledge of protonation-coupled conformational dynamics is valuable for elucidating mechanisms of energy transduction processes in biology, such as enzyme catalysis and electron transfer as well as proton and drug transport. Toward this goal we present a new technique of embedding continuous constant pH molecular dynamics within an explicit-solvent representation. In this technique we make use of the efficiency of the generalized-Born (GB) implicit-solvent model for estimating the free energy of protein solvation while propagating conformational dynamics using the more accurate explicit-solvent model. Also, we employ a pH-based replica exchange scheme to significantly enhance both protonation and conformational state sampling. Benchmark data of five proteins including HP36, NTL9, BBL, HEWL, and SNase yield an average absolute deviation of 0.53 and a root mean squared deviation of 0.74 from experimental data. This level of accuracy is obtained with 1 ns simulations per replica. Detailed analysis reveals that explicit-solvent sampling provides increased accuracy relative to the previous GB-based method by preserving the native structure, providing a more realistic description of conformational flexibility of the hydrophobic cluster, and correctly modeling solvent mediated ion-pair interactions. Thus, we anticipate that the new technique will emerge as a practical tool to capture ionization equilibria while enabling an intimate view of ionization coupled conformational dynamics that is difficult to delineate with experimental techniques alone.

Role of φ29 connector channel loops in late-stage DNA packaging

Double-stranded DNA bacteriophages and their eukaryotic virus counterparts have 12-fold head-tail connector assemblages embedded at a unique capsid vertex. This vertex is the site of assembly of the DNA packaging motor, and the connector has a central channel through which viral DNA passes during genome packaging and subsequent host infection. Crystal structures of connectors from different phages reveal either disordered residues or structured loops that project into the connector channel. Given the proximity to the translocating DNA substrate, these loops have been proposed to play a role in DNA packaging. Previous models have proposed structural motions in either the packaging ATPase or the connector channel loops as the driving force that translocates the DNA into the prohead. Here, we mutate the channel loops of the Bacillus subtilis bacteriophage φ29 connector and show that these loops have no active role in translocation of DNA. Instead, they appear to have an essential function near the end of packaging, acting to retain the packaged DNA in the head in preparation for motor detachment and subsequent tail assembly and virion completion.

The effect of mutations on the alloreactive T cell receptor/peptide-MHC interface structure: a molecular dynamics study

T cells orchestrate adaptive, pathogen-specific immune responses. T cells have a surface receptor (called TCR) whose ligands are complexes (pMHCs) of peptides (derived from pathogens or host proteins) and major histocompatibility complex proteins (MHCs). MHC proteins vary between hosts. During organ transplants, host TCRs interact with peptides present in complex with genetically different MHCs. This usually causes a vigorous immune response: alloreactivity. Studies of alloreactive protein interactions have yielded results that present a puzzle. Some crystallographic studies concluded that the alloreactive TCR/MHC interface is essentially unaffected by changing the TCR peptide-binding region, suggesting that the peptide does not influence the interface. Another biochemical study concluded from mutation data that different peptides can alter the binding interface with the same TCR. To explore the origin of this puzzle, we used molecular dynamics simulations to study the dependence of the TCR/pMHC interface on changes in both the peptide and the TCR. Our simulations show that the footprint of the TCR on the pMHC is insensitive to mutations of the TCR peptide-binding loops, but peptide mutations can make multiple local changes to TCR/pMHC contacts. Therefore, our results demonstrate that the structural and mutation data do not conflict and reveal how subtle, but important, characteristics of the alloreactive TCR/pMHC interface are influenced by the TCR and the peptide.