Vibrational subsystem analysis: A method for probing free energies and correlations in the harmonic limit

A new vibrational subsystem analysis (VSA) method is presented for coupling global motion to a local subsystem while including the inertial effects of the environment. The premise of the VSA method is a partitioning of a system into a smaller region of interest and a usually larger part referred to as environment. This method allows the investigation of local-global coupling, a more accurate estimation of vibrational free energy contribution for parts of a large system, and the elimination of the "tip effect" in elastic network model calculations. Additionally, the VSA method can be used as a probe of specific degrees of freedom that may contribute to free energy differences. The VSA approach can be employed in many ways, but it will likely be most useful for estimating activation free energies in QM/MM reaction path calculations. Four examples are presented to demonstrate the utility of this method.

Factors governing helix formation in peptides confined to carbon nanotubes

The effect of confinement on the stability and dynamics of peptides and proteins is relevant in the context of a number of problems in biology and biotechnology. We have examined the stability of different helix-forming sequences upon confinement to a carbon nanotube using Langevin dynamics simulations of a coarse-grained representation of the polypeptide chain. We show that the interplay of several factors that include sequence, solvent conditions, strength (lambda) of nanotube-peptide interactions, and the nanotube diameter (D) determines confinement-induced stability of helicies. In agreement with predictions based on polymer theory, the helical state is entropically stabilized for all sequences when the interaction between the peptide and the nanotube is weakly hydrophobic and D is small. However, there is a strong sequence dependence as the strength of the lambda increases. For an amphiphilic sequence, the helical stability increases with lambda, whereas for polyalanine the diagram of states is a complex function of lambda and D. In addition, decreasing the size of the "hydrophobic patch" lining the nanotube, which mimics the chemical heterogeneity of the ribosome tunnel, increases the helical stability of the polyalanine sequence. Our results provide a framework for interpreting a number of experiments involving the structure formation of peptides in the ribosome tunnel as well as transport of biopolymers through nanotubes.

Backbone relaxation coupled to the ionization of internal groups in proteins: a self-guided Langevin dynamics study

Pathways of structural relaxation triggered by ionization of internal groups in staphylococcal nuclease (SNase) were studied through multiple self-guided Langevin dynamics (SGLD) simulations. Circular dichroism, steady-state Trp fluorescence, and nuclear magnetic resonance spectroscopy have shown previously that variants of SNase with internal Glu, Asp, and Lys at positions 66 or 92, and Arg at position 66, exhibit local reorganization or global unfolding when the internal ionizable group is charged. Except for Arg-66, these internal ionizable groups have unusual pKa values and are neutral at physiological pH. The structural trends observed in the simulations are in general agreement with experimental observations. The I92D variant, which unfolds globally upon ionization of Asp-92, in simulations often exhibits extensive hydration of the protein core, and sometimes also significant perturbations of the beta-barrel. In the crystal structure of the V66R variant, the beta1 strand from the beta-barrel is domain-swapped; in the simulations, the beta1 strand is sometimes partially released. The V66K variant, which in solutions shows reorganization of six residues at the C-terminus of helix alpha1 and perturbations in the beta-barrel structure, exhibits fraying of three residues of helix alpha1 in one simulation, and perturbations and partial unfolding of three beta-strands in a few other simulations. In sharp contrast, very small structural changes were observed in simulations of the wild-type protein. The simulations indicate that charging of internal groups frequently triggers penetration of water into the protein interior. The pKa values of Asp-92 and Arg-66 calculated with continuum methods on SGLD-relaxed structures reached the normal values in most simulations. Detailed analysis of accuracy and performance of SGLD demonstrates that SGLD outperforms LD in sampling of alternative protein conformations without loss of the accuracy and level of detail characteristic of regular LD.

Using the isotropic periodic sum method to calculate long-range interactions of heterogeneous systems

Isotropic periodic sum (IPS) is a method for the calculation of long-range interactions in molecular simulation based on the homogeneity of simulation systems. Three IPS models, 3D IPS, 2D IPS, and 1D IPS have been developed for three common types of homogeneous systems. Based on the fact that 3D IPS can well describe the long-range interactions of a heterogeneous system if a local region larger than the homogeneity scale is used, this work presents a method based on 3D IPS to calculate long-range interactions for all kinds of simulation systems, including homogeneous, heterogeneous, and finite systems. Unlike the original 3D IPS method that uses a local region defined by the cutoff distance, this method uses a local region larger than that defined by the cutoff distance to reach the homogeneity scale. To efficiently calculate interactions within such a large local region, this method split long-range interactions into two parts, a cutoff part and a long-range part. The cutoff part is calculated by summing over atom pairs within a cutoff range (about 10 A), and the long-range part is calculated using the discrete fast Fourier transform (DFFT) technique. This method is applied to electrostatic and van der Waals interactions for both periodic and non-periodic systems. Example simulations demonstrate that this method can accurately and efficiently calculate long-range interactions for molecular simulation.

A new model defines the minimal set of polymorphism in HLA-DQ and -DR that determines susceptibility and resistance to autoimmune diabetes

Background

The mechanism underlying autoimmune diabetes has been difficult to define. There is a strong genetic contribution and numerous studies associate the major histocompatibility complex, especially the class II region, with predisposition or resistance. However, how these molecules are implicated remains obscure.

Presentation of the hypothesis

We have supplemented structural analysis with computational biophysical and sequence analyses and propose an heuristic for distinguishing between human leukocyte antigen molecules that predispose to insulin dependent diabetes mellitus and those that are protective. Polar residues at both β37 and β9 suffice to distinguish accurately between class II alleles that predispose to type 1 diabetes and those that do not. The electrostatic potential within the peptide binding pocket exerts a strong influence on diabetogenic epitopes with basic residues. Diabetes susceptibility alleles are predicted to bind autoantigens strongly with tight affinity, prolonged association and altered cytokine expression profile. Protective alleles bind moderately, and neutral alleles poorly or not at all. Non-Asp β57 is a modifier that supplements disease risk but only in the presence of the polymorphic, polar pair at β9 and β37. The nature of β37 determines resistance on one hand, and susceptibility or dominant protection on the other.

Conclusion

The proposed ideas are illustrated with structural, functional and population studies from the literature. The hypothesis, in turn, rationalizes their results. A plausible mechanism of immune mediated diabetes based on binding affinity and peptide kinetics is discussed. The number of the polymorphic markers present correlates with onset of disease and severity. The molecular elucidation of disease susceptibility and resistance paves the way for risk prediction, treatment and prevention of disease based on analogue peptides.

Reviewers

This article was reviewed by Eugene V. Koonin, Michael Lenardo, Hossam Ashour, and Bhagirath Singh. For the full reviews, please go to the Reviewers' comments section.

Open science grid study of the coupling between conformation and water content in the interior of a protein

Computational grids are a promising resource for modeling complex biochemical processes such as protein folding, penetration of gases or water into proteins, or protein structural rearrangements coupled to ligand binding. We have enabled the molecular dynamics program CHARMM to run on the Open Science Grid. The implementation is general, flexible, easily modifiable for use with other molecular dynamics programs and other grids and automated in terms of job submission, monitoring, and resubmission. The usefulness of grid computing was demonstrated through the study of hydration of the Glu-66 side chain in the interior of protein staphylococcal nuclease. Multiple simulations started with and without two internal water molecules shown crystallographically to be associated with the side chain of Glu-66 yielded two distinct populations of rotameric states of Glu-66 that differed by as much as 20%. This illustrates how internal water molecules can bias protein conformations. Furthermore, there appeared to be a temporal correlation between dehydration of the side chain and conformational transitions of Glu-66. This example demonstrated how difficult it is to get convergence even in the relatively simple case of a side chain oscillating between two conformations. With grid computing, we also benchmarked the self-guided Langevin dynamics method against the Langevin dynamics method traditionally used for temperature control in molecular dynamics simulations and showed that the two methods yield comparable results.

Multiscale methods for macromolecular simulations

In this article we review the key modeling tools available for simulating biomolecular systems. We consider recent developments and representative applications of mixed quantum mechanics/molecular mechanics (QM/MM), elastic network models (ENMs), coarse-grained molecular dynamics, and grid-based tools for calculating interactions between essentially rigid protein assemblies. We consider how the different length scales can be coupled, both in a sequential fashion (e.g. a coarse-grained or grid model using parameterization from MD simulations), and via concurrent approaches, where the calculations are performed together and together control the progression of the simulation. We suggest how the concurrent coupling approach familiar in the context of QM/MM calculations can be generalized, and describe how this has been done in the CHARMM macromolecular simulation package.

A computational characterization of boron-oxygen multiple bonding in HN=CH-CH=CH-NH-BO

Structures, relative energies, and bonding characteristics for various conformers of 3-imino-N-(oxoboryl)prop-1-en-1-amine, HN=CH-CH=CH-NH-BO, and the corresponding borocycle (-HN=CH-CH=CH-NH-B-)O are discussed using results from second-order Møller-Plesset (MP2) perturbation theory with the Dunning-Woon correlation-consistent cc-pVDZ, aug-cc-pVDZ, and cc-pVTZ basis sets. These MP2 results are compared to those from computationally efficient density functional theory (DFT) calculations using the LDA, PBE, TPSS, BLYP, B3LYP, BVP86, OLYP, O3LYP, and PBE1PBE functionals in conjunction with the economical Pople-type 6-311++G(d,p) basis set to evaluate the suitability of these DFT/6-311++G(d,p) levels for use with larger boron-containing systems. The effects of an aqueous environment were incorporated into the calculations using COSMO methodology. The calculated boron-oxygen bond lengths, orbital compositions, and bond orders in all the (acyclic) HN=CH-CH=CH-NH-BO conformers were consistent with the presence of a boron-oxygen triple bond, similar to that found in H-BO and H2N-BO. The (-HN=CH-CH=CH-NH-B-)O borocycle is predicted to be planar (C2v symmetry), and it is approximately 30 kcal/mol lower in energy than any of the (acyclic) HN=CH-CH=CH-NH-BO conformers; the boron-oxygen bond in this borocycle has significant double bond character, a bonding scheme for which there has been only one experimental structure reported in the literature (Vidovic, D. ; et al. J. Am. Chem. Soc. 2005, 127, 4566- 4569).

CHARMMing: a new, flexible web portal for CHARMM

A new web portal for the CHARMM macromolecular modeling package, CHARMMing (CHARMM interface and graphics, http://www.charmming.org), is presented. This tool provides a user-friendly interface for the preparation, submission, monitoring, and visualization of molecular simulations (i.e., energy minimization, solvation, and dynamics). The infrastructure used to implement the web application is described. Two additional programs have been developed and integrated with CHARMMing: GENRTF, which is employed to define structural features not supported by the standard CHARMM force field, and a job broker, which is used to provide a portable method for using grid and cluster computing with CHARMMing. The use of the program is described with three proteins: 1YJP , 1O1O , and 1UFY . Source code is provided allowing CHARMMing to be downloaded, installed, and used by supercomputing centers and research groups that have a CHARMM license. Although no software can replace a scientist's own judgment and experience, CHARMMing eases the introduction of newcomers to the molecular modeling discipline by providing a graphical method for running simulations.

Molecular dynamics simulations of perylene and tetracene librations: comparison with femtosecond upconversion data

In a prior manuscript by Xu et al. [Xu, J.; Shen, X.; Knutson, J. R. J. Phys. Chem. A 2003, 107, 8383], time-resolved fluorescence emission anisotropy measurements were performed on perylene and tetracene in hexadecane using an upconversion technique with approximately 100 fs resolution. The anisotropy transients contained previously unseen decay terms of approximately 300 fs. In perylene, their amplitude corresponded to the " r(o) defect" that has gathered interest over decades. We ascribed this term to a predominantly in-plane libration. In this manuscript, we present molecular dynamics simulations for the motions of perylene and tetracene using the CHARMm molecular dynamics program (version c29b2). Both rotational correlation functions contain subpicosecond decay terms that resemble experimental anisotropy decays. It was suggested that the r(o) defect might arise from excited-state distortions of perylene, so we conducted quantum mechanical calculations to show that such distortion does not significantly displace the oscillators. We compare the case of perylene, with a strongly allowed singlet emission transition, to that of the weakly allowed tetracene transition. In perylene, motion alone can explain subpicosecond anisotropy decay, while tetracene decay also contains vibrational coupling terms, as previously reported by Sarkar et al. [Sarkar, N.; Takeuchi, S.; Tahara, T. J. Phys. Chem. A 1999, 103, 4808].

Pathways and populations: stereoelectronic insights into the exocyclic torsion of 5-(hydroxymethyl)tetrahydropyran

High level ab initio computations in vacuum and with the IEFPCM implicit solvent model are carried out on 5-(hydroxymethyl)tetrahydropyran to investigate the effects of water on the exocyclic torsional surface. Rotamer populations evaluated from the omega(C-C-C-O), theta(C-C-C-O) solvent surface agree almost quantitatively with experimental values for the closely related methyl 4-deoxy-alpha-D-xylohexopyranoside. Potentials of mean force obtained from the two surfaces show substantial solvent stabilization of the TG (omega = 180 +/- 60 degrees) rotamer and the barriers at omega= 120 and 240 degrees but solvent destabilization at the cis barrier (omega = 0 degrees). Natural bond orbital analyses indicate that energetics of these effects are largely explained by overstabilization of the vacuum GT (omega= 60 +/- 60 degrees) and GG (omega = 300 +/- 60 degrees) rotamers. Solvent stabilization of theta conformations provides entropic stabilization.

Molecular modeling of nucleic acid structure

This unit is the first in a series of four units covering the analysis of nucleic acid structure by molecular modeling. This unit provides an overview of computer simulation of nucleic acids. Topics include the static structure model, computational graphics and energy models, generation of an initial model, and characterization of the overall three-dimensional structure.

Molecular modeling of nucleic acid structure: setup and analysis

The last in a set of units by these authors, this unit addresses some important remaining questions about molecular modeling of nucleic acids. It describes how to choose an appropriate molecular mechanics force field; how to set up and equilibrate the system for accurate simulation of a nucleic acid in an explicit solvent by molecular dynamics or Monte Carlo simulation; and how to analyze molecular dynamics trajectories.

Extended polypeptide linkers establish the spatial architecture of a pyruvate dehydrogenase multienzyme complex

Icosahedral pyruvate dehydrogenase (PDH) enzyme complexes are molecular machines consisting of a central E2 core decorated by a shell of peripheral enzymes (E1 and E3) found localized at a distance of approximately 75-90 A from the core. Using a combination of biochemical, biophysical, and cryo-electron microscopic techniques, we show here that the gap between the E2 core and the shell of peripheral enzymes is maintained by the flexible but extended conformation adopted by 60 linker polypeptides that radiate outwards from the inner E2 core, irrespective of the E1 or E3 occupancy. The constancy of the gap is thus not due to protein-protein interactions in the outer protein shell. The extended nature of the E2 inner-linker regions thereby creates the restricted annular space in which the lipoyl domains of E2 that carry catalytic intermediates shuttle between E1, E2, and E3 active sites, while their conformational flexibility facilitates productive encounters.

Langevin network model of myosin

Langevin mode theory and the coarse-grained elastic network model (ENM) for proteins are combined to yield the Langevin network model (LNM). Hydrodynamic radii of 6 A were assigned to each alpha-carbon on the basis of matching experimental translational and rotational diffusion constants of lysozyme, myoglobin, and hemoglobin with those calculated using a rigid body bead model with hydrodynamic interactions described by the Rotne-Prager tensor. LNM analysis of myosin II indicates that all ENM-like modes are overdamped at water viscosities. The low-frequency LNM modes in the pre-power stroke structure (PDB code: 1VOM) are substantially less mixed than the corresponding modes of the post-power stroke structure (1Q5G). Results from a four-bead model of the myosin "lever arm" indicate that coupling between modes increases as the array departs from linearity and are consistent with the results for 1VOM and 1Q5G. The decay times for all overdamped Langevin modes are shorter than the calculated rotational tumbling times found for lysozyme and myosin.

Collective and noncollective models of NMR relaxation in lipid vesicles and multilayers

NMR (13)C spin lattice relaxation (1/T(1)) rates of dipalmitoylphosphatidylcholine (DPPC) bilayers obtained from molecular dynamics simulations of 72 and 288 lipids are compared with each other, with experimental values from large liposomes obtained by magic angle spinning, and with previously published experimental data from small vesicles. The experimental results for multilayers and vesicles at the same frequencies differ only slightly. The simulation results indicate that T(1) relaxation in the 15.1 to 201.2 MHz carbon frequency range and up to 100 A length scale is dominated by fast isomerizations and slower lipid wobble (D perpendicular approximately 2.5 x 10(8) s(-1)). Rotational diffusion about the lipid long axis (described by D(parallel)) does not make a substantial contribution to the T(1). Modifications to the acyl chain torsional potential energy function used for the simulations substantially improve agreement with experiment.

Dimers of boroglycine and methylamine boronic acid: a computational comparison of the relative importance of dative versus hydrogen bonding

Boronic acids are widely used in materials science, pharmacology, and the synthesis of biologically active compounds. In this Article, geometrical structures and relative energies of dimers of boroglycine, H2N-CH2-B(OH)2, and its constitutional isomer H3C-NH-B(OH)2, were computed using second-order Møller-Plesset perturbation theory and density functional theory; Dunning-Woon correlation-consistent cc-pVDZ, aug-cc-pVDZ, cc-pVTZ, and aug-cc-pVTZ basis sets were employed for the MP2 calculations, and the Pople 6-311++G(d,p) basis set was employed for a majority of the DFT calculations. Effects of an aqueous environment were incorporated into the results using PCM and COSMO-RS methodology. The lowest-energy conformer of the H2N-CH2-B(OH)2 dimer was a six-membered ring structure (chair conformation; Ci symmetry) with two intermolecular B:N dative-bonds; it was 14.0 kcal/mol lower in energy at the MP2/aug-cc-pVDZ computational level than a conformer with the classic eight-centered ring structure (Ci symmetry) in which the boroglycine monomers are linked by a pair of H-O...H bonds. Compared to the results of MP2 calculations with correlation-consistent basis sets, DFT calculations using the PBE1PBE and TPSS functionals with the 6-311++G(d,p) basis set were significantly better at predicting relative conformational energies of the H2N-CH2-B(OH)2 and H3C-NH-B(OH)2 dimers than corresponding calculations using the BLYP, B3LYP, OLYP, and O3LYP functionals, particularly with respect to dative-bonded structures.

Protein conformational transitions explored by mixed elastic network models

We develop a mixed elastic network model (MENM) to study large-scale conformational transitions of proteins between two (or more) known structures. Elastic network potentials for the beginning and end states of a transition are combined, in effect, by adding their respective partition functions. The resulting effective MENM energy function smoothly interpolates between the original surfaces, and retains the beginning and end structures as local minima. Saddle points, transition paths, potentials of mean force, and partition functions can be found efficiently by largely analytic methods. To characterize the protein motions during a conformational transition, we follow "transition paths" on the MENM surface that connect the beginning and end structures and are invariant to parameterizations of the model and the mathematical form of the mixing scheme. As illustrations of the general formalism, we study large-scale conformation changes of the motor proteins KIF1A kinesin and myosin II. We generate possible transition paths for these two proteins that reveal details of their conformational motions. The MENM formalism is computationally efficient and generally applicable even for large protein systems that undergo highly collective structural changes.

Protein recognition motifs: design of peptidomimetics of helix surfaces

Helices represent one of the most common recognition motifs in proteins. The design of nonpeptidic scaffolds, such as the 3,2',2''-tris-substituted terphenyl, that can imitate the side-chain orientation along one face of an alpha-helix potentially provides an effective means to modulate helix-recognition functions. Here, based on theoretical arguments, we described novel alpha-helix mimetics which are more effective than the terphenyl at constraining the aryl-aryl torsion angles to those associated with structures suitable for mimicking the alpha-helical twist for side-chain orientation and for superimposing the side chains of residues i, i + 3 or i + 4, i + 7 when compared with the alpha-beta side-chain vectors of the regular alpha-helix with an improved root mean square deviation (RMSD) of approximately 0.5 A. In addition, this study suggests that rotamer distributions around the C(alpha)--C(beta) bonds of these helix mimetics are similar to those of alpha-helices, except that these rotamer distributions show an approximately 60 degrees shift compared to those of alpha-helices when the mimetic axis is superimposed upon the helix axis. This change in rotamer orientation complicates mimicry of the helix surface.

Ab initio modeling of glycosyl torsions and anomeric effects in a model carbohydrate: 2-ethoxy tetrahydropyran

A range of ab initio calculations were carried out on the axial and equatorial anomers of the model carbohydrate 2-ethoxy tetrahydropyran to evaluate the level of theory required to accurately evaluate the glycosyl dihedral angle and the anomeric ratio. Vacuum CCSD(T)/CBS extrapolations at the global minimum yield DeltaE = E(equatorial) - E(axial) = 1.42 kcal/mol. When corrected for solvent (by the IEFPCM model), zero-point vibrations and entropy, DeltaG(298) = 0.49 kcal/mol, in excellent agreement with the experimental value of 0.47 +/- 0.3 kcal/mol. A new additivity scheme, the layered composite method (LCM), yields DeltaE to within 0.1 kcal/mol of the CCSD(T)/CBS result at a fraction of the computer requirements. Anomeric ratios and one-dimensional torsional surfaces generated by LCM and the even more efficient MP2/cc-pVTZ level of theory are in excellent agreement, indicating that the latter is suitable for force-field parameterization of carbohydrates. Hartree-Fock and density functional theory differ from CCSD(T)/CBS for DeltaE by approximately 1 kcal/mol; they show similar deviations in torsional surfaces evaluated from LCM. A comparison of vacuum and solvent-corrected one- and two-dimensional torsional surfaces indicates the equatorial form of 2-ethoxy tetrahydropyran is more sensitive to solvent than the axial.