CHARMM: the biomolecular simulation program

The 1.4 A crystal structure of the ArsD arsenic metallochaperone provides insights into its interaction with the ArsA ATPase

Arsenic is a carcinogen that tops the Superfund list of hazardous chemicals. Bacterial resistance to arsenic is facilitated by ArsD, which delivers As(III) to the ArsA ATPase, the catalytic subunit of the ArsAB pump. Here we report the structure of the arsenic metallochaperone ArsD at 1.4 A and a model for its binding of metalloid. There are two ArsD molecules in the asymmetric unit. The overall structure of the ArsD monomer has a thioredoxin fold, with a core of four beta-strands flanked by four alpha-helices. Based on data from structural homologues, ArsD was modeled with and without bound As(III). ArsD binds one arsenic per monomer coordinated with the three sulfur atoms of Cys12, Cys13, and Cys18. Using this structural model, an algorithm was used to dock ArsD and ArsA. The resulting docking model provides testable predictions of the contact points of the two proteins and forms the basis for future experiments.

Using fit functions in computational dielectric spectroscopy

This work deals with the development of an appropriate set of fit functions for describing dielectric spectra based on simulated raw data. All these fit functions are of exponential character with properly chosen cofunctions. The type of the cofunctions is different for translation, rotation and their coupling. As an alternative to multiexponential fits we also discuss Kohlrausch-Williams-Watts functions. Since the corresponding Fourier-Laplace series for these stretched exponentials has severe convergence problems, we represent their Fourier-Laplace spectrum as a Havriliak-Negami expression with properly chosen parameters. A general relation between the parameter of the Kohlrausch-Williams-Watts and the Havriliak-Negami parameters is given. The set of fit functions is applied to the concrete simulation of the hydrated ionic liquid 1-ethyl-3-methyl-imidazolium triflate with H(2)O. The systematic variation of the water mole fraction permits to study the gradual transition from a neutral molecular liquid to molecular ionic liquids.

Comparison of solid-state dipolar couplings and solution relaxation data provides insight into protein backbone dynamics

Analyses of solution (15)N relaxation data and solid-state (1)H(N)-(15)N dipolar couplings from a small globular protein, alpha-spectrin SH3 domain, produce a surprisingly similar pattern of order parameters. This result suggests that there is little or no ns-mus dynamics throughout most of the sequence and, in particular, in the structured portion of the backbone. At the same time, evidence of ns-mus motions is found in the flexible loops and termini. These findings, corroborated by the MD simulations of alpha-spectrin SH3 in a hydrated crystalline environment and in solution, are consistent with the picture of protein dynamics that has recently emerged from the solution studies employing residual dipolar couplings.

Role of electrostatics in modulating hydrophobic interactions and barriers to hydrophobic assembly

Hydrophobic effects continue to be an active area of research due to implications for a wide range of physicochemical phenomena. Molecular dynamics simulations have been used extensively in the study of such effects using various water potential models, with few studies addressing the differences between models. In particular, studies considering the explicit treatment of water polarizability are underrepresented in the literature. We present results from molecular dynamics simulations that systematically compare the dependence of large-scale hydrophobic effects on the water model. We consider three common nonpolarizable models (SPC/E, TIP3P, and TIP4P) and two common polarizable models (TIP4P-FQ and SWM4-NDP). Results highlight the similarities and differences of the different water models in the vicinity of two large hydrophobic plates. In particular, profiles of average density, density fluctuations, orientation, and hydrogen bonding show only minor differences among the water models studied. However, the potential of mean force for the hydrophobe dimerization is significantly reduced in the polarizable water systems. TIP4P-FQ shows the deepest minimum of approximately -54(+/-3) kcal/mol compared to -40(+/-3), -40(+/-2), -42(+/-3), and -45(+/-5) kcal/mol for TIP4P, TIP3P, SPC/E, and SWM4-NDP (all relative to the dissociated state). We discuss the relationship between hydrophobic association and the strength of water-water interactions in the liquid phase. Results suggest that models treating polarizability (both implicitly and explicitly) influence a stronger driving force toward hydrophobic assembly. Implications of these results, as well as prospectives on future work, are discussed.

Molecular switch controlling the binding of anionic bile acid conjugates to human apical sodium-dependent bile acid transporter

The human apical sodium-dependent bile acid transporter (hASBT) may serve as a prodrug target for oral drug absorption. Synthetic, biological, NMR, and computational approaches identified the structure-activity relationships of mono- and dianionic bile acid conjugates for hASBT binding. Experimental data combined with a conformationally sampled pharmacophore/QSAR modeling approach (CSP-SAR) predicted that dianionic substituents with intramolecular hydrogen bonding between hydroxyls on the cholane skeleton and the acid group on the conjugate's aromatic ring increased conjugate hydrophobicity and improved binding affinity. Notably, the model predicted the presence of a conformational molecular switch, where shifting the carboxylate substituent on an aromatic ring by a single position controlled binding affinity. Model validation was performed by effectively shifting the spatial location of the carboxylate by inserting a methylene adjacent to the aromatic ring, resulting in the predicted alteration in binding affinity. This work illustrates conformation as a determinant of ligand physiochemical properties and ligand binding affinity to a biological transporter.

Population distribution of flexible molecules from maximum entropy analysis using different priors as background information: application to the Φ, Ψ-conformational space of the α-(1-->2)-linked mannose disaccharide present in N- and O-linked glycoproteins

The conformational space available to the flexible molecule α-D-Manp-(1-->2)-α-D-Manp-OMe, a model for the α-(1-->2)-linked mannose disaccharide in N- or O-linked glycoproteins, is determined using experimental data and molecular simulation combined with a maximum entropy approach that leads to a converged population distribution utilizing different input information. A database survey of the Protein Data Bank where structures having the constituent disaccharide were retrieved resulted in an ensemble with >200 structures. Subsequent filtering removed erroneous structures and gave the database (DB) ensemble having three classes of mannose-containing compounds, viz., N- and O-linked structures, and ligands to proteins. A molecular dynamics (MD) simulation of the disaccharide revealed a two-state equilibrium with a major and a minor conformational state, i.e., the MD ensemble. These two different conformation ensembles of the disaccharide were compared to measured experimental spectroscopic data for the molecule in water solution. However, neither of the two populations were compatible with experimental data from optical rotation, NMR (1)H,(1)H cross-relaxation rates as well as homo- and heteronuclear (3)J couplings. The conformational distributions were subsequently used as background information to generate priors that were used in a maximum entropy analysis. The resulting posteriors, i.e., the population distributions after the application of the maximum entropy analysis, still showed notable deviations that were not anticipated based on the prior information. Therefore, reparameterization of homo- and heteronuclear Karplus relationships for the glycosidic torsion angles Φ and Ψ were carried out in which the importance of electronegative substituents on the coupling pathway was deemed essential resulting in four derived equations, two (3)J(COCC) and two (3)J(COCH) being different for the Φ and Ψ torsions, respectively. These Karplus relationships are denoted JCX/SU09. Reapplication of the maximum entropy analysis gave excellent agreement between the MD- and DB-posteriors. The information entropies show that the current reparametrization of the Karplus relationships constitutes a significant improvement. The Φ(H) torsion angle of the disaccharide is governed by the exo-anomeric effect and for the dominating conformation Φ(H) = -40 degrees and Ψ(H) = 33 degrees. The minor conformational state has a negative Ψ(H) torsion angle; the relative populations of the major and the minor states are approximately 3 : 1. It is anticipated that application of the methodology will be useful to flexible molecules ranging from small organic molecules to large biomolecules.

Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types

A significant modification to the additive all-atom CHARMM lipid force field (FF) is developed and applied to phospholipid bilayers with both choline and ethanolamine containing head groups and with both saturated and unsaturated aliphatic chains. Motivated by the current CHARMM lipid FF (C27 and C27r) systematically yielding values of the surface area per lipid that are smaller than experimental estimates and gel-like structures of bilayers well above the gel transition temperature, selected torsional, Lennard-Jones and partial atomic charge parameters were modified by targeting both quantum mechanical (QM) and experimental data. QM calculations ranging from high-level ab initio calculations on small molecules to semiempirical QM studies on a 1,2-dipalmitoyl-sn-phosphatidylcholine (DPPC) bilayer in combination with experimental thermodynamic data were used as target data for parameter optimization. These changes were tested with simulations of pure bilayers at high hydration of the following six lipids: DPPC, 1,2-dimyristoyl-sn-phosphatidylcholine (DMPC), 1,2-dilauroyl-sn-phosphatidylcholine (DLPC), 1-palmitoyl-2-oleoyl-sn-phosphatidylcholine (POPC), 1,2-dioleoyl-sn-phosphatidylcholine (DOPC), and 1-palmitoyl-2-oleoyl-sn-phosphatidylethanolamine (POPE); simulations of a low hydration DOPC bilayer were also performed. Agreement with experimental surface area is on average within 2%, and the density profiles agree well with neutron and X-ray diffraction experiments. NMR deuterium order parameters (S(CD)) are well predicted with the new FF, including proper splitting of the S(CD) for the aliphatic carbon adjacent to the carbonyl for DPPC, POPE, and POPC bilayers. The area compressibility modulus and frequency dependence of (13)C NMR relaxation rates of DPPC and the water distribution of low hydration DOPC bilayers also agree well with experiment. Accordingly, the presented lipid FF, referred to as C36, allows for molecular dynamics simulations to be run in the tensionless ensemble (NPT), and is anticipated to be of utility for simulations of pure lipid systems as well as heterogeneous systems including membrane proteins.

Update of the CHARMM all-atom additive force field for lipids: validation on six lipid types

A significant modification to the additive all-atom CHARMM lipid force field (FF) is developed and applied to phospholipid bilayers with both choline and ethanolamine containing head groups and with both saturated and unsaturated aliphatic chains. Motivated by the current CHARMM lipid FF (C27 and C27r) systematically yielding values of the surface area per lipid that are smaller than experimental estimates and gel-like structures of bilayers well above the gel transition temperature, selected torsional, Lennard-Jones and partial atomic charge parameters were modified by targeting both quantum mechanical (QM) and experimental data. QM calculations ranging from high-level ab initio calculations on small molecules to semiempirical QM studies on a 1,2-dipalmitoyl-sn-phosphatidylcholine (DPPC) bilayer in combination with experimental thermodynamic data were used as target data for parameter optimization. These changes were tested with simulations of pure bilayers at high hydration of the following six lipids: DPPC, 1,2-dimyristoyl-sn-phosphatidylcholine (DMPC), 1,2-dilauroyl-sn-phosphatidylcholine (DLPC), 1-palmitoyl-2-oleoyl-sn-phosphatidylcholine (POPC), 1,2-dioleoyl-sn-phosphatidylcholine (DOPC), and 1-palmitoyl-2-oleoyl-sn-phosphatidylethanolamine (POPE); simulations of a low hydration DOPC bilayer were also performed. Agreement with experimental surface area is on average within 2%, and the density profiles agree well with neutron and X-ray diffraction experiments. NMR deuterium order parameters (S(CD)) are well predicted with the new FF, including proper splitting of the S(CD) for the aliphatic carbon adjacent to the carbonyl for DPPC, POPE, and POPC bilayers. The area compressibility modulus and frequency dependence of (13)C NMR relaxation rates of DPPC and the water distribution of low hydration DOPC bilayers also agree well with experiment. Accordingly, the presented lipid FF, referred to as C36, allows for molecular dynamics simulations to be run in the tensionless ensemble (NPT), and is anticipated to be of utility for simulations of pure lipid systems as well as heterogeneous systems including membrane proteins.

How does bone sialoprotein promote the nucleation of hydroxyapatite? A molecular dynamics study using model peptides of different conformations

Bone sialoprotein (BSP) is a highly phosphorylated, acidic, noncollagenous protein in bone matrix. Although BSP has been proposed to be a nucleator of hydroxyapatite (Ca(5)(PO(4))(3)OH), the major mineral component of bone, no detailed mechanism for the nucleation process has been elucidated at the atomic level to date. In the present work, using a peptide model, we apply molecular dynamics (MD) simulations to study the conformational effect of a proposed nucleating motif of BSP (a phosphorylated, acidic, 10 amino-acid residue sequence) on controlling the distributions of Ca(2+) and inorganic phosphate (Pi) ions in solution, and specifically, we explore whether a nucleating template for orientated hydroxyapatite could be formed in different peptide conformations. Both the alpha-helical conformation and the random coil structure have been studied, and inorganic solutions without the peptide are simulated as reference. Ca(2+) distributions around the peptide surface and interactions between Ca(2+) and Pi in the presence of the peptide are examined in detail. From the MD simulations, although in some cases for the alpha-helical conformation, we observe that a Ca(2+) equilateral triangle forms around the surface of peptide, which matches the distribution of Ca(2+) ions on the (001) face of the hydroxyapatite crystal, we do not consistently find a stable nucleating template formation in general for either the helical conformation or the random coil structure. Therefore, independent of conformations, the BSP nucleating motif is more likely to help nucleate an amorphous calcium phosphate cluster, which ultimately converts to crystalline hydroxyapatite.

Complete phenotypic recovery of an Alzheimer's disease model by a quinone-tryptophan hybrid aggregation inhibitor

The rational design of amyloid oligomer inhibitors is yet an unmet drug development need. Previous studies have identified the role of tryptophan in amyloid recognition, association and inhibition. Furthermore, tryptophan was ranked as the residue with highest amyloidogenic propensity. Other studies have demonstrated that quinones, specifically anthraquinones, can serve as aggregation inhibitors probably due to the dipole interaction of the quinonic ring with aromatic recognition sites within the amyloidogenic proteins. Here, using in vitro, in vivo and in silico tools we describe the synthesis and functional characterization of a rationally designed inhibitor of the Alzheimer's disease-associated β-amyloid. This compound, 1,4-naphthoquinon-2-yl-L-tryptophan (NQTrp), combines the recognition capacities of both quinone and tryptophan moieties and completely inhibited Aβ oligomerization and fibrillization, as well as the cytotoxic effect of Aβ oligomers towards cultured neuronal cell line. Furthermore, when fed to transgenic Alzheimer's disease Drosophila model it prolonged their life span and completely abolished their defective locomotion. Analysis of the brains of these flies showed a significant reduction in oligomeric species of Aβ while immuno-staining of the 3^rd instar larval brains showed a significant reduction in Aβ accumulation. Computational studies, as well as NMR and CD spectroscopy provide mechanistic insight into the activity of the compound which is most likely mediated by clamping of the aromatic recognition interface in the central segment of Aβ. Our results demonstrate that interfering with the aromatic core of amyloidogenic peptides is a promising approach for inhibiting various pathogenic species associated with amyloidogenic diseases. The compound NQTrp can serve as a lead for developing a new class of disease modifying drugs for Alzheimer's disease.

PARP1 ADP-ribosylates lysine residues of the core histone tails

The chromatin-associated enzyme PARP1 has previously been suggested to ADP-ribosylate histones, but the specific ADP-ribose acceptor sites have remained enigmatic. Here, we show that PARP1 covalently ADP-ribosylates the amino-terminal histone tails of all core histones. Using biochemical tools and novel electron transfer dissociation mass spectrometric protocols, we identify for the first time K13 of H2A, K30 of H2B, K27 and K37 of H3, as well as K16 of H4 as ADP-ribose acceptor sites. Multiple explicit water molecular dynamics simulations of the H4 tail peptide into the catalytic cleft of PARP1 indicate that two stable intermolecular salt bridges hold the peptide in an orientation that allows K16 ADP-ribosylation. Consistent with a functional cross-talk between ADP-ribosylation and other histone tail modifications, acetylation of H4K16 inhibits ADP-ribosylation by PARP1. Taken together, our computational and experimental results provide strong evidence that PARP1 modifies important regulatory lysines of the core histone tails.

Molecular description of flexibility in an antibody combining site

Mature antibodies (Abs) that are exquisitely specific for virtually any foreign molecule may be produced by affinity maturation of naïve (or germline) Abs. However, the finite number of germline Abs available suggests that, in contrast to mature Abs, germline Abs must be broadly polyspecific so that they are able to recognize a wide range of ligands. Thus, affinity maturation must play a role in mediating Ab specificity. One biophysical property that distinguishes polyspecificity from specificity is protein flexibility; a flexible combining site is able to adopt different conformations that recognize different foreign molecules (or antigens), while a rigid combining site is locked into a conformation that is specific for a given antigen. Recent studies (Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 8821-8826) have examined, at the atomic level, the structural properties that mediate changes in flexibility at four stages of affinity maturation in the 4-4-20 Ab. These studies employed molecular dynamics simulations to reveal a network of residue interactions that mediate the flexibility changes accompanying maturation. The flexibility of the Ab combining sites in these molecular systems was originally measured using three-pulse photon echo spectroscopy (3PEPS). The present investigation extends this work by providing a concrete link between structural properties of the Ab molecules and features of the spectroscopic measurements used to characterize their flexibility. Results obtained from the simulations are in good qualitative agreement with the experimental measurements and indicate that the spectroscopic signal is sensitive to protein dynamics distributed throughout the entire combining site. Thus, the simulations provide a molecular-level interpretation of the changes induced by affinity maturation of the Ab. The results suggest that 3PEPS spectroscopy in combination with molecular dynamics simulations can provide a detailed description of protein dynamics and, in this case, how it is evolved for biological function.

Molecular dynamics simulation of bovine pancreatic ribonuclease A-CpA and transition state-like complexes

The mechanisms of enzymes are intimately connected with their overall structure and dynamics in solution. Experimentally, it is considerably challenging to provide detailed atomic level information about the conformational events that occur at different stages along the chemical reaction path. Here, theoretical tools may offer new potential insights that complement those obtained from experiments that may not yield an unambiguous mechanistic interpretation. In this study, we apply molecular dynamics simulations of bovine pancreatic ribonuclease A, an archetype ribonuclease, to study the conformational dynamics, structural relaxation, and differential solvation that occur at discrete stages of the transesterification and cleavage reaction. Simulations were performed with explicit solvation with rigorous electrostatics and utilize recently developed molecular mechanical force field parameters for transphosphorylation and hydrolysis transition state analogues. Herein, we present results for the enzyme complexed with the dinucleotide substrate cytidilyl-3',5'-adenosine (CpA) in the reactant, and transphosphorylation and hydrolysis transition states. A detailed analysis of active site structures and hydrogen-bond patterns is presented and compared. The integrity of the overall backbone structure is preserved in the simulations and supports a mechanism whereby His12 stabilizes accumulating negative charge at the transition states through hydrogen-bond donation to the nonbridge oxygens. Lys41 is shown to be highly versatile along the reaction coordinate and can aid in the stabilization of the dianionic transition state, while being poised to act as a general acid catalyst in the hydrolysis step.

On the orientation of the catalytic dyad in aspartic proteases

The recent re-refinement of the X-ray structure of apo plasmepsin II from Plasmodium falciparum suggests that the two carboxylate groups in the catalytic dyad are noncoplanar, (Robbins et al., Acta Crystallogr D Biol Crystallogr 2009;65: 294-296) in remarkable contrast with the vast majority of structures of aspartic proteases. Here, evidence for the noncoplanarity of the catalytic aspartates is provided by analysis of multiple explicit water molecular dynamics (MD) simulations of plasmepsin II, human beta-secretase, and HIV-protease. In the MD runs of plasmepsin II, the angle between the planes of the two carboxylates of the catalytic dyad is almost always in the range 60 degrees -120 degrees , in agreement with the perpendicular orientation in the re-refined X-ray structure. The noncoplanar arrangement is prevalent also in the beta-secretase simulations, as well as in the runs with the inhibitor-bound proteases. Quantum-mechanics calculations provide further evidence that before catalysis the noncoplanar arrangement is favored energetically in eukaryotic aspartic proteases. Remarkably, the coplanar orientation of the catalytic dyad is observed in MD simulations of HIV-protease at 100 K but not at 300 K, which indicates that the noncoplanar arrangement is favored by conformational entropy. This finding suggests that the coplanar orientation in the crystal structures of apo aspartic proteases is promoted by the very low temperature used for data collection (usually around 100 K).

Evolution: a guide to perturb protein function and networks

Protein interactions give rise to networks that control cell fate in health and disease; selective means to probe these interactions are therefore of wide interest. We discuss here Evolutionary Tracing (ET), a comparative method to identify protein functional sites and to guide experiments that selectively block, recode, or mimic their amino acid determinants. These studies suggest, in principle, a scalable approach to perturb individual links in protein networks.

The sodium/galactose symporter crystal structure is a dynamic, not so occluded state

The recent elucidation of the sodium/galactose symporter structure from the Vibrio parahaemolyticus bacterium, vSGLT, has revealed a similarity in the core architecture with transporters from different gene families, including the leucine transporter (LeuT). Even though several transporters sharing this core have been structurally determined over the past few years, vSGLT is the only one crystallized in the substrate-bound inward-facing conformation so far. In this study, we report the first insight into the dynamics and coordination of the galactose (Gal) and proposed Na+ ion in vSGLT using a series of molecular dynamics simulations with a total time of about 0.1 micros. Our study reveals new residues, not observed in the crystal structure, which closely interact with the Na(+) ion or the substrate for extended times, and shows that in the crystallized conformation, a Na+ ion placed at the site equivalent to Na2 in LeuT can escape into the intracellular (IC) space in the absence of external forces. We have identified the highly conserved Asp189 as a likely binding residue on the pathway of Na(+) into the IC cavity. The release of Gal, on the other hand, requires the rotation of the side chain of the inner hydrophobic gate, Tyr263, without a significant change in vSGLT backbone conformation. Our simulations further show that the crystal structure represents but one accessible binding pose of Gal and Na+ among an ensemble of microstates, and that the Gal undergoes versatile alternate interactions at the binding pocket.

Coarse-grained model of nucleic acid bases

Atomistic simulations of nucleic acids are prohibitively expensive and, consequently, reduced models of these compounds are of great interest in the field. In this work, we propose a physics-based coarse-grained model of nucleic-acid bases in which each base is represented by several (3-5) interaction centers. van der Waals interactions are modeled by Lennard-Jones spheres with a 12-6 potential energy function. The charge distribution is modeled by a set of electric dipole moments located at the centers of the Lennard-Jones spheres. The method for computing the Lennard-Jones parameters, electric dipole moments (their magnitude and orientation) and positions of the interaction centers is described. Several models with different numbers of interaction centers were tested. The model with three-center cytosine, four-center guanine, four-center thymine, and five-center adenine satisfactorily reproduces the canonical Watson-Crick hydrogen bonding and stacking interaction energies of the all-atom AMBER model. The computation time with the coarse-grained model is reduced seven times compared with that of the all-atom model.

Fine grained sampling of residue characteristics using molecular dynamics simulation

In a fine-grained computational analysis of protein structure, we investigated the relationships between a residue's backbone conformations and its side-chain packing as well as conformations. To produce continuous distributions in high resolution, we ran molecular dynamics simulations over a set of protein folds (dynameome). In effect, the dynameome dataset samples not only the states well represented in the PDB but also the known states that are not well represented in the structural database. In our analysis, we characterized the mutual influence among the backbone phi,psi angles with the first side-chain torsion angles (chi(1)) and the volumes occupied by the side-chains. The dependencies of these relationships on side-chain environment and amino acids are further explored. We found that residue volumes exhibit dependency on backbone 2 degrees structure conformation: side-chains pack more densely in extended beta-sheet than in alpha-helical structures. As expected, residue volumes on the protein surface were larger than those in the interior. The first side-chain torsion angles are found to be dependent on the backbone conformations in agreement with previous studies, but the dynameome dataset provides higher resolution of rotamer preferences based on the backbone conformation. All three gauche(-), gauche(+), and trans rotamers show different patterns of phi,psi dependency, and variations in chi(1) value are skewed from their canonical values to relieve the steric strains. By demonstrating the utility of dynameomic modeling on the native state ensemble, this study reveals details of the interplay among backbone conformations, residue volumes and side-chain conformations.

Current trends in antimicrobial agent research: chemo- and bioinformatics approaches

Databases and chemo- and bioinformatics tools that contain genomic, proteomic and functional information have become indispensable for antimicrobial drug research. The combination of chemoinformatics tools, bioinformatics tools and relational databases provides means of analyzing, linking and comparing online search results. The development of computational tools feeds on a diversity of disciplines, including mathematics, statistics, computer science, information technology and molecular biology. The computational approach to antimicrobial agent discovery and design encompasses genomics, molecular simulation and dynamics, molecular docking, structural and/or functional class prediction, and quantitative structure-activity relationships. This article reviews progress in the development of computational methods, tools and databases used for organizing and extracting biological meaning from antimicrobial research.

The impact of viral RNA on the association rates of capsid protein assembly: bacteriophage MS2 as a case study

A large number of single-stranded RNA viruses, which form a major class of all viruses, co-assemble their protein container and their genomic material. The multiple roles of the viral genome in this process are presently only partly understood. Recent experimental results indicate that RNA, in addition to its function as a repository for genetic information, could play important functional roles during the assembly of the viral protein containers. An investigation of the impact of genomic RNA on the association of the protein subunits may therefore provide further insights into the mechanism of virus assembly. We study here the impact of viral RNA on the association rates of the capsid proteins during virus assembly. As a case study, we consider the viral capsid of bacteriophage MS2, which is formed from 60 asymmetric (AB) and 30 symmetric (CC) protein dimers. Using Brownian dynamics simulations, we investigate the effect of the binding of an RNA stem-loop (the translational repressor) on the association rates of the capsid protein dimers. Our analysis shows that translational repressor binding results in self-association of AB dimers being inhibited, whilst association of AB with CC dimers is greatly enhanced. This provides an explanation for experimental results in which an alternating assembly pattern of AB and CC dimer addition to the growing assembly intermediate has been observed to be the dominant mode of assembly. The presence of the RNA hence dramatically decreases the number of dominant assembly pathways and thereby reduces the complexity of the self-assembly process of these viruses.

Protein dynamics investigated by inherent structure analysis

Molecular dynamics (MD) simulations provide essential information about the thermodynamics and dynamics of proteins. To construct the free-energy surface from equilibrium trajectories, it is necessary to group the individual snapshots in a meaningful way. The inherent structures (IS) are shown to provide an appropriate discretization of the trajectory and to avoid problems that can arise in clustering algorithms that have been employed previously. The IS-based approach is illustrated with a 30-ns room temperature "native" state MD simulation of a 10-residue peptide in a beta-hairpin conformation. The transitions between the IS are used to construct a configuration space network from which a one-dimensional free-energy profile is obtained with the mincut method. The results demonstrate that the IS approach is useful and that even for this simple system, there exists a nontrivial organization of the native state into several valleys separated by barriers as high as 3 kcal/mol. Further, by introducing a coarse-grained network, it is demonstrated that there are multiple pathways connecting the valleys. This scenario is hidden when the snapshots of the trajectory are used directly with rmsd clustering to compute the free-energy profile. Application of the IS approach to the native state of the PDZ2 signaling domain indicates its utility for the study of biologically relevant systems.

Interaction of the α2A domain of integrin with small collagen fragments

We here present a detailed study of the ligand-receptor interactions between single and triple-helical strands of collagen and the α2A domain of integrin (α2A), providing valuable new insights into the mechanisms and dynamics of collagen-integrin binding at a sub-molecular level. The occurrence of single and triple-helical strands of the collagen fragments was scrutinized with atom force microscopy (AFM) techniques. Strong interactions of the triple-stranded fragments comparable to those of collagen can only be detected for the 42mer triple-helical collagen-like peptide under study (which contains 42 amino acid residues per strand) by solid phase assays as well as by surface plasmon resonance (SPR) measurements. However, changes in NMR signals during titration and characteristic saturation transfer difference (STD) NMR signals are also detectable when α2A is added to a solution of the 21mer single-stranded collagen fragment. Molecular dynamics (MD) simulations employing different sets of force field parameters were applied to study the interaction between triple-helical or single-stranded collagen fragments with α2A. It is remarkable that even single-stranded collagen fragments can form various complexes with α2A showing significant differences in the complex stability with identical ligands. The results of MD simulations are in agreement with the signal alterations in our NMR experiments, which are indicative of the formation of weak complexes between single-stranded collagen and α2A in solution. These results provide useful information concerning possible interactions of α2A with small collagen fragments that are of relevance to the design of novel therapeutic A-domain inhibitors.

Binding of imidazole to the heme of cytochrome c1 and inhibition of the bc1 complex from Rhodobacter sphaeroides: I. Equilibrium and modeling studies

We have used imidazole (Im) and N-methylimidazole (MeIm) as probes of the heme-binding cavity of membrane-bound cytochrome (cyt) c(1) in detergent-solubilized bc(1) complex from Rhodobacter sphaeroides. Imidazole binding to cyt c(1) substantially lowers the midpoint potential of the heme and fully inhibits bc(1) complex activity. Temperature dependences showed that binding of Im (K(d) approximately 330 microM, 25 degrees C, pH 8) is enthalpically driven (DeltaH(0) = -56 kJ/mol, DeltaS(0) = -121 J/mol/K), whereas binding of MeIm is 30 times weaker (K(d) approximately 9.3 mM) and is entropically driven (DeltaH(0) = 47 kJ/mol, DeltaS(0)(o) = 197 J/mol/K). The large enthalpic and entropic contributions suggest significant structural and solvation changes in cyt c(1) triggered by ligand binding. Comparison of these results with those obtained previously for soluble cyts c and c(2) suggested that Im binding to cyt c(1) is assisted by formation of hydrogen bonds within the heme cleft. This was strongly supported by molecular dynamics simulations of Im adducts of cyts c, c(2), and c(1), which showed hydrogen bonds formed between the N(delta)H of Im and the cyt c(1) protein, or with a water molecule sequestered with the ligand in the heme cleft.

Slow folding of cross-linked alpha-helical peptides due to steric hindrance

The folding process of a 16-residue alpha-helical peptide with an azobenzene cross-linker (covalently bound to residues Cys3 and Cys14) is investigated by 50 molecular dynamics simulations of 4 micros each. The folding kinetics at 281 K show a stretched exponential behavior but become simpler and much faster when a distance restraint is used to emulate a nonbulky cross-linker. The free-energy basin of the helical state is divided into two subbasins by a barrier that separates helical conformations with opposite orientations of the Arg10 side chain with respect to the azobenzene cross-linker. In contrast, such barrier is not present in the helical basin of the peptide with the nonbulky cross-linker, which folds with speed similar to the unrestrained peptide. These results indicate that the cross-linker slows down folding because of steric hindrance rather than its restraining effect on the two ends of the helical segment.

Function inferences from a molecular structural model of bacterial ParE toxin

Toxin-antitoxin (TA) systems contribute to plasmid stability by a mechanism that relies on the differential stabilities of the toxin and antitoxin proteins and leads to the killing of daughter bacteria that did not receive a plasmid copy at the cell division. ParE is the toxic component of a TA system that constitutes along with RelE an important class of bacterial toxin called RelE/ParE superfamily. For ParE toxin, no crystallographic structure is available so far and rare in vitro studies demonstrated that the target of toxin activity is E. coli DNA gyrase. Here, a 3D Model for E. coli ParE toxin by molecular homology modeling was built using MODELLER, a program for comparative modeling. The Model was energy minimized by CHARMM and validated using PROCHECK and VERIFY3D programs. Resulting Ramachandran plot analysis it was found that the portion residues failing into the most favored and allowed regions was 96.8%. Structural similarity search employing DALI server showed as the best matches RelE and YoeB families. The Model also showed similarities with other microbial ribonucleases but in a small score. A possible homologous deep cleft active site was identified in the Model using CASTp program. Additional studies to investigate the nuclease activity in members of ParE family as well as to confirm the inhibitory replication activity are needed. The predicted Model allows initial inferences about the unexplored 3D structure of the ParE toxin and may be further used in rational design of molecules for structure-function studies.

A repulsive electrostatic mechanism for protein export through the type III secretion apparatus

Many Gram-negative bacteria initiate infections by injecting effector proteins into host cells through the type III secretion apparatus, which is comprised of a basal body, a needle, and a tip. The needle channel is formed by the assembly of a single needle protein. To explore the export mechanisms of MxiH needle protein through the needle of Shigella flexneri, an essential step during needle assembly, we have performed steered molecular dynamics simulations in implicit solvent. The trajectories reveal a screwlike rotation motion during the export of nativelike helix-turn-helix conformations. Interestingly, the channel interior with excessive electronegative potential creates an energy barrier for MxiH to enter the channel, whereas the same may facilitate the ejection of the effectors into host cells. Structurally known basal regions and ATPase underneath the basal region also have electronegative interiors. Effector proteins also have considerable electronegative potential patches on their surfaces. From these observations, we propose a repulsive electrostatic mechanism for protein translocation through the type III secretion apparatus. Based on this mechanism, the ATPase activity and/or proton motive force could be used to energize the protein translocation through these nanomachines. A similar mechanism may be applicable to macromolecular channels in other secretion systems or viruses through which proteins or nucleic acids are transported.

Versatile substrate protein recognition mechanism of the eukaryotic chaperonin CCT

Group II chaperonins, found in eukaryotic and archaeal organisms, recognize substrate proteins through diverse mechanisms that involve either hydrophobic- or electrostatic-dominated interactions. This action is distinct from the universal substrate recognition mechanism of group I chaperonins, which bind a wide spectrum of non-native proteins primarily through hydrophobic interactions. We use computational approaches to pinpoint the substrate protein binding sites of the gamma-subunit of the eukaryotic chaperonin CCT and to identify its interactions with the stringent substrate beta-tubulin. Protein-protein docking methods reveal intrinsic binding sites of CCT comprising a helical (HL) region, homologous to the GroEL-binding site, and the helical protrusion (HP) region. We performed molecular dynamics simulations of the solvated CCTgamma apical domain, beta-tubulin peptide-CCTgamma complexes, and isolated beta-tubulin peptides. We find that tubulin binds to CCTgamma through an extensive interface that spans both the HL region and the HP region. HL interactions involve both hydrophobic and electrostatic contacts, while binding to the HP region is stabilized almost exclusively by a salt bridge network. On the basis of additional simulations of a beta-tubulin-CCTgamma complex that involves a reduced interface, centered onto the HP region, we conclude that this salt bridge network is the minimal stabilizing interaction required. Strong conservation of the charged amino acids that participate in the salt bridge network, Arg306 and Glu271, indicates a general mechanism across the nonidentical CCT subunits and group II chaperonins.

Engineering of substrate selectivity for tissue factor.factor VIIa complex signaling through protease-activated receptor 2

The complex of factor VIIa (FVIIa) with tissue factor (TF) triggers coagulation by recognizing its macromolecular substrate factors IX (FIX) and X (FX) predominantly through extended exosite interactions. In addition, TF mediates unique cell-signaling properties in cancer, angiogenesis, and inflammation that involve proteolytic cleavage of protease-activated receptor 2 (PAR2). PAR2 is cleaved by FVIIa in the binary TF.FVIIa complex and by FXa in the ternary TF.FVIIa.FXa complex, but physiological roles of these signaling complexes are incompletely understood. In a screen of FVIIa protease domain mutants, three variants (Q40A, Q143N, and T151S) activated macromolecular coagulation substrates and supported signaling of the ternary TF.FVIIa-Xa complex normally but were severely impaired in binary TF.FVIIa.PAR2 signaling. The residues identified were located in the model-predicted S2' pocket of FVIIa, and complementary PAR2 P2' Leu-38 replacements demonstrated that the P2' side chain was indeed crucial for PAR2 cleavage by TF.FVIIa. In addition, PAR2 was activated more efficiently by FVIIa T99Y, consistent with further contributions from the S2 subsite. The P2 residue preference of FVIIa and FXa predicted additional PAR2 mutants that were efficiently activated by TF.FVIIa but resistant to cleavage by the alternative PAR2 activator FXa. Thus, contrary to the paradigm of exosite-assisted cleavage of PAR1 by thrombin, the cofactor-associated protease FVIIa recognizes PAR2 predominantly by catalytic cleft interactions. Furthermore, the delineated molecular details of this substrate interaction enabled protein engineering of protease-selective PAR2 receptors that will aid further studies to dissect the roles of TF signaling complexes in vivo.

Computational design and biological testing of highly cytotoxic colchicine ring A modifications

Microtubules are the primary target for many anti-cancer drugs, the majority of which bind specifically to beta-tubulin. The existence of several beta-tubulin isotypes, coupled with their varied expression in normal and cancerous cells provides a platform upon which to construct selective chemotherapeutic agents. We have examined five prevalent human beta-tubulin isotypes and identified the colchicine-binding site as the most promising for drug design based on specificity. Using this binding site as a template, we have designed several colchicine derivatives and computationally probed them for affinity to the beta-tubulin isotypes. These compounds were synthesized and subjected to cytotoxicity assays to determine their effectiveness against several cancerous cell lines. We observed a correlation between computational-binding predictions and experimentally determined IC(50) values, demonstrating the utility of computational screening in the design of more effective colchicine derivatives. The most promising derivative exhibited an IC(50) approximately threefold lower than values previously reported for either colchicine or paclitaxel, demonstrating the utility of computational design and assessment of binding to tubulin.

Mechanical unfolding of an ankyrin repeat protein

Ankryin repeat proteins comprise tandem arrays of a 33-residue, predominantly alpha-helical motif that stacks roughly linearly to produce elongated and superhelical structures. They function as scaffolds mediating a diverse range of protein-protein interactions, and some have been proposed to play a role in mechanical signal transduction processes in the cell. Here we use atomic force microscopy and molecular-dynamics simulations to investigate the natural 7-ankyrin repeat protein gankyrin. We find that gankyrin unfolds under force via multiple distinct pathways. The reactions do not proceed in a cooperative manner, nor do they always involve fully stepwise unfolding of one repeat at a time. The peeling away of half an ankyrin repeat, or one or more ankyrin repeats, occurs at low forces; however, intermediate species are formed that are resistant to high forces, and the simulations indicate that in some instances they are stabilized by nonnative interactions. The unfolding of individual ankyrin repeats generates a refolding force, a feature that may be more easily detected in these proteins than in globular proteins because the refolding of a repeat involves a short contraction distance and incurs a low entropic cost. We discuss the origins of the differences between the force- and chemical-induced unfolding pathways of ankyrin repeat proteins, as well as the differences between the mechanics of natural occurring ankyrin repeat proteins and those of designed consensus ankyin repeat and globular proteins.

Human TLRs 10 and 1 share common mechanisms of innate immune sensing but not signaling

TLRs are central receptors of the innate immune system that drive host inflammation and adaptive immune responses in response to invading microbes. Among human TLRs, TLR10 is the only family member without a defined agonist or function. Phylogenetic analysis reveals that TLR10 is most related to TLR1 and TLR6, both of which mediate immune responses to a variety of microbial and fungal components in cooperation with TLR2. The generation and analysis of chimeric receptors containing the extracellular recognition domain of TLR10 and the intracellular signaling domain of TLR1, revealed that TLR10 senses triacylated lipopeptides and a wide variety of other microbial-derived agonists shared by TLR1, but not TLR6. TLR10 requires TLR2 for innate immune recognition, and these receptors colocalize in the phagosome and physically interact in an agonist-dependent fashion. Computational modeling and mutational analysis of TLR10 showed preservation of the essential TLR2 dimer interface and lipopeptide-binding channel found in TLR1. Coimmunoprecipitation experiments indicate that, similar to TLR2/1, TLR2/10 complexes recruit the proximal adaptor MyD88 to the activated receptor complex. However, TLR10, alone or in cooperation with TLR2, fails to activate typical TLR-induced signaling, including NF-kappaB-, IL-8-, or IFN-beta-driven reporters. We conclude that human TLR10 cooperates with TLR2 in the sensing of microbes and fungi but possesses a signaling function distinct from that of other TLR2 subfamily members.

A second generation distributed point polarizable water model

A distributed point polarizable model (DPP2) for water, with explicit terms for charge penetration, induction, and charge transfer, is introduced. The DPP2 model accurately describes the interaction energies in small and large water clusters and also gives an average internal energy per molecule and radial distribution functions of liquid water in good agreement with experiment. A key to the success of the model is its accurate description of the individual terms in the n-body expansion of the interaction energies.

Probing the folding of mini-protein Beta3s by two-dimensional infrared spectroscopy; simulation study

We propose to use infrared coherent two-dimensional correlation spectroscopy (2DCS) to characterize the folding mechanism of the mini-protein Beta3s. In this study Beta3s was folded by molecular dynamics (MD) simulation and intermediate conformational ensembles were identified. The one and two-dimensional correlation spectrum was calculated for the intermediate and native states of the mini-protein. A direct structure-spectra relationship was determined by analysis of conformational properties and specific residue contributions. We identified the structural origin of diagonal and off-diagonal peaks in the 2DCS spectra for the native and intermediate conformational ensembles in the folding mechanism. This work supports the implementation of computational techniques in conjunction with experimental 2DCS to study the folding mechanism of proteins. In addition to exploring the folding mechanism the work presented here can be applied in combination with experiment to refine and validate current molecular dynamics force fields.

PACS Codes: 87.15.Cc, 87.15.hm, 87.15.hp

Ion selectivity of alpha-hemolysin with a beta-cyclodextrin adapter. I. Single ion potential of mean force and diffusion coefficient

The alpha-hemolysin (alphaHL) is a self-assembling exotoxin that binds to the membrane of a susceptible host cell and causes its death. Experimental studies show that electrically neutral beta-cyclodextrin (betaCD) can insert into the alphaHL channel and significantly increase its anion selectivity. To understand how betaCD can affect ion selectivity, molecular dynamics simulations and potential of mean force (PMF) calculations are carried out for different alphaHL channels with and without the betaCD adapter. A multiscale approach based on the generalized solvent boundary potential is used to reduce the size of the simulated system. The PMF profiles reveal that betaCD has no anion selectivity by itself but can increase the Cl(-) selectivity of the alphaHL channel when lodged into the pore lumen. Analysis shows that betaCD causes a partial desolvation of ions and affects the orientation of nearby charged residues. The ion selectivity appears to result from increased electrostatic interaction between the ion and the channel due to a reduction in dielectric shielding by the solvent. These observations suggest a reasonable explanation of the ion selectivity and provide important information for further ion channel modification.

Pi release from myosin: a simulation analysis of possible pathways

The release of phosphate (Pi) is an important element in actomyosin function and has been shown to be accelerated by the binding of myosin to actin. To provide information about the structural elements important for Pi release, possible escape pathways from various isolated myosin II structures have been determined by molecular dynamics simulations designed for studying such slow processes. The residues forming the pathways were identified and their role was evaluated by mutant simulations. Pi release is slow in the pre-powerstroke structure, an important element in preventing the powerstroke prior to actin binding, and is much more rapid for Pi modeled into the post-rigor and rigor-like structures. The previously proposed backdoor route is dominant in the pre-powerstroke and post-rigor states, whereas a different path is most important in the rigor-like state. This finding suggests a mechanism for the actin-activated acceleration of Pi release.

Ion selectivity of alpha-hemolysin with beta-cyclodextrin adapter. II. Multi-ion effects studied with grand canonical Monte Carlo/Brownian dynamics simulations

In a previous study of ion selectivity of alpha-hemolysin (alphaHL) in complex with beta-cyclodextrin (betaCD) adapter, we calculated the potential of mean force (PMF) and characterized the self-diffusion coefficients of isolated K(+) and Cl(-) ions using molecular dynamics simulations (Y. Luo et al., "Ion Selectivity of alpha-Hemolysin with beta-Cyclodextrin Adapter: I. Single Ion Potential of Mean Force and Diffusion Coefficient"). In the present effort, these results pertaining to single isolated ions in the wide aqueous pore are extended to take into account multi-ion effects. The grand canonical Monte Carlo/Brownian dynamics (GCMC/BD) algorithm is used to simulate ion currents through the wild-type alphaHL ion channel, as well as two engineered alphaHL mutants, with and without the cyclic oligosaccaride betaCD lodged in the lumen of the pore. The GCMC/BD current-voltage curves agree well with experimental results and show that betaCD increases the anion selectivity of alphaHL. Comparisons between multi-ion PMFs from GCMC/BD simulations and single-ion PMFs demonstrate that multi-ion effects and pore shape are crucial for explaining this behavior. It is concluded that the narrow betaCD adapter increases the anion selectivity of alphaHL because it reduces the pore radius locally, which decreases the ionic screening and the dielectric shielding of the strong electrostatic field induced by a nearby ring of positively charged alphaHL side chains.

Current status of the AMOEBA polarizable force field

Molecular force fields have been approaching a generational transition over the past several years, moving away from well-established and well-tuned, but intrinsically limited, fixed point charge models toward more intricate and expensive polarizable models that should allow more accurate description of molecular properties. The recently introduced AMOEBA force field is a leading publicly available example of this next generation of theoretical model, but to date, it has only received relatively limited validation, which we address here. We show that the AMOEBA force field is in fact a significant improvement over fixed charge models for small molecule structural and thermodynamic observables in particular, although further fine-tuning is necessary to describe solvation free energies of drug-like small molecules, dynamical properties away from ambient conditions, and possible improvements in aromatic interactions. State of the art electronic structure calculations reveal generally very good agreement with AMOEBA for demanding problems such as relative conformational energies of the alanine tetrapeptide and isomers of water sulfate complexes. AMOEBA is shown to be especially successful on protein-ligand binding and computational X-ray crystallography where polarization and accurate electrostatics are critical.

Molecular dynamics simulations show that conformational selection governs the binding preferences of imatinib for several tyrosine kinases

Tyrosine kinases transmit cellular signals through a complex mechanism, involving their phosphorylation and switching between inactive and active conformations. The cancer drug imatinib binds tightly to several homologous kinases, including Abl, but weakly to others, including Src. Imatinib specifically targets the inactive, so-called "DFG-out" conformation of Abl, which differs from the preferred, "DFG-in" conformation of Src in the orientation of a conserved Asp-Phe-Gly (DFG) activation loop. However, recent x-ray structures showed that Src can also adopt the DFG-out conformation and uses it to bind imatinib. The Src/Abl-binding free energy difference can thus be decomposed into two contributions. Contribution i measures the different protein-imatinib interactions when either kinase is in its DFG-out conformation. Contribution ii depends on the ability of imatinib to select or induce this conformation, i.e. on the relative stabilities of the DFG-out and DFG-in conformations of each kinase. Neither contribution has been measured experimentally. We use molecular dynamics simulations to show that contribution i is very small, 0.2 +/- 0.6 kcal/mol; imatinib interactions are very similar in the two kinases, including long range electrostatic interactions with the imatinib positive charge. Contribution ii, deduced using the experimental binding free energy difference, is much larger, 4.4 +/- 0.9 kcal/mol. Thus, conformational selection, easy in Abl, difficult in Src, underpins imatinib specificity. Contribution ii has a simple interpretation; it closely approximates the stability difference between the DFG-out and DFG-in conformations of apo-Src. Additional calculations show that conformational selection also governs the relative binding of imatinib to the kinases c-Kit and Lck. These results should help clarify the current framework for engineering kinase signaling.

CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields

The widely used CHARMM additive all-atom force field includes parameters for proteins, nucleic acids, lipids, and carbohydrates. In the present article, an extension of the CHARMM force field to drug-like molecules is presented. The resulting CHARMM General Force Field (CGenFF) covers a wide range of chemical groups present in biomolecules and drug-like molecules, including a large number of heterocyclic scaffolds. The parametrization philosophy behind the force field focuses on quality at the expense of transferability, with the implementation concentrating on an extensible force field. Statistics related to the quality of the parametrization with a focus on experimental validation are presented. Additionally, the parametrization procedure, described fully in the present article in the context of the model systems, pyrrolidine, and 3-phenoxymethylpyrrolidine will allow users to readily extend the force field to chemical groups that are not explicitly covered in the force field as well as add functional groups to and link together molecules already available in the force field. CGenFF thus makes it possible to perform "all-CHARMM" simulations on drug-target interactions thereby extending the utility of CHARMM force fields to medicinally relevant systems.