Defining Coarse-Grained Representations of Large Biomolecules and Biomolecular Complexes from Elastic Network Models

Coarse-grained (CG) models of large biomolecular complexes enable simulations of these systems over long timescales that are not accessible for atomistic molecular dynamics (MD) simulations. A systematic methodology, called essential dynamics coarse-graining (ED-CG), has been developed for defining coarse-grained sites in a large biomolecule. The method variationally determines the CG sites so that key dynamic domains in the protein are preserved in the CG representation. The original ED-CG method relies on a principal component analysis (PCA) of a MD trajectory. However, for many large proteins and multi-protein complexes such an analysis may not converge or even be possible. This work develops a new ED-CG scheme using an elastic network model (ENM) of the protein structure. In this procedure, the low-frequency normal modes obtained by ENM are used to define dynamic domains and to define the CG representation accordingly. The method is then applied to several proteins, such as the HIV-1 CA protein dimer, ATP-bound G-actin, and the Arp2/3 complex. Numerical results show that ED-CG with ENM (ENM-ED-CG) is much faster than ED-CG with PCA because no MD is necessary. The ENM-ED-CG models also capture functional essential dynamics of the proteins almost as well as those using full MD with PCA. Therefore, the ENM-ED-CG method may be better suited to coarse-grain a very large biomolecule or biomolecular complex that is too computationally expensive to be simulated by conventional MD, or when a high resolution atomic structure is not even available.     

Conformational sampling by NMR solution structures calculated with the program DIANA evaluated by comparison with long-time molecular dynamics calculations in explicit water

The NMR solution structure of bovine pancreatic trypsin inhibitor (BPTI) obtained by distance geometry calculations with the program DIANA is compared with groups of conformers generated by molecular dynamics (MD) simulations in explicit water at ambient temperature and pressure. The MD simulations started from a single conformer and were free or restrained either by the experimental NOE distance restraints or by time-averaged restraints; the groups of conformers were collected either in 10 ps intervals during 200 ps periods of simulation, or in 50 ps intervals during a 1 ns period of simulation. Overall, these comparisons show that the standard protein structure determination protocol with the program DIANA provides a picture of the protein structure that is in agreement with MD simulations using “realistic” potential functions over a nanosecond timescale. For well-constrained molecular regions there is a trend in the free MD simulation of duration 1 ns that the sampling of the conformation space is slightly increased relative to the DIANA calculations. In contrast, for surface-exposed side-chains that are less extensively constrained by the NMR data, the DIANA conformers tend to sample larger regions of conformational space than conformers selected from any of the MD trajectories. Additional insights into the behavior of surface side-chains come from comparison of the MD runs of 200 ps or 1 ns duration. In this time range the sampling of conformation space by the protein surface depends strongly on the length of the simulation, which indicates that significant side-chain transitions occur on the nanosecond timescale and that much longer simulations will be needed to obtain statistically significant data on side-chain dynamics.     

Brownian and essential dynamics studies of the HIV-1 integrase catalytic domain

The three-dimensional structure of the active site region of the enzyme HIV-1 integrase is not unambiguously known. This region includes a flexible peptide loop that cannot be well resolved in crystallographic determinations. Here we present two different computational approaches with different levels of resolution and on different time-scales to understand this flexibility and to analyze the dynamics of this part of the protein. We have used molecular dynamics simulations with an atomic model to simulate the region in a realistic and reliable way for 1 ns. It is found that parts of the loop wind up after 300 ps to extend an existing helix. This indicates that the helix is longer than in the earlier crystal structures that were used as basis for this study. Very recent crystal data confirms this finding, underlining the predictive value of accurate MD simulations. Essential dynamics analysis of the MD trajectory yields an anharmonic motion of this loop. We have supplemented the MD data with a much lower resolution Brownian dynamics simulation of 600 ns length. It provides ideas about the slow-motion dynamics of the loop. It is found that the loop explores a conformational space much larger than in the MD trajectory, leading to a "gating"-like motion with respect to the active site.     

Brownian and Essential Dynamics Studies of the HIV-1 Integrase Catalytic Domain

Abstract

The three-dimensional structure of the active site region of the enzyme HIV-1 integrase is not unambiguously known. This region includes a flexible peptide loop that cannot be well resolved in crystallographic determinations. Here we present two different computional approaches with different levels of resolution and on different time-scales to understand this flexibility and to analyze the dynamics of this part of the protein. We have used molecular dynamics simulations with an atomic model to simulate the region in a realistic and reliable way for 1 ns. It is found that parts of the loop wind up after 300 ps to extend an existing helix. This indicates that the helix is longer than in the earlier crystal structures that were used as basis for this study. Very recent crystal data confirms this finding, underlining the predictive value of accurate MD simulations. Essential dynamics analysis of the MD trajectory yields an anharmonic motion of this loop. We have supplemented the MD data with a much lower resolution Brownian dynamics simulation of 600 ns length. It provides ideas about the slow-motion dynamics of the loop. It is found that the loop explores a conformational space much larger than in the MD trajectory, leading to a “gating”-like motion with respect to the active site.     

Molecular structures of human argininosuccinate synthetase pseudogenes. Evolutionary and mechanistic implications

In the human genome there is one expressed gene for argininosuccinate synthetase and 14 pseudogenes. A cDNA coding for human argininosuccinate synthetase was used to screen a human genomic library. Twenty-five unique genomic clones were isolated and extensively characterized. At least seven clones represented processed argininosuccinate synthetase pseudogenes that lost the introns in the expressed gene. Restriction mapping demonstrated that these processed pseudogenes were located in distinct regions of the human genome. Complete nucleotide sequences of two processed pseudogenes, psi AS-1 and psi AS-3, and a partial sequence of psi AS-7 were determined. Both psi AS-1 and psi AS-3 had an adenine-rich region at their 3' end and were flanked by distinct imperfect direct repeats. A comparison of these pseudogene sequences to that of the cDNA demonstrated that psi AS-1 and psi AS-3 were 93% homologous to the cDNA, whereas psi AS-7 was 89% homologous to the cDNA. Therefore, it is estimated that psi AS-1 and psi AS-3 were created 10-11 million years ago, whereas psi AS-7 arose approximately 21 million years ago. We have estimated the evolutionary rate for the expressed argininosuccinate synthetase gene based on the sequences of psi AS-1 and psi AS-3. These data indicate that the expressed argininosuccinate synthetase gene is evolving at a rate similar to that of the beta-globin gene and much faster than the alpha-tubulin gene. Furthermore, a comparison of the sequences of psi AS-1 and psi AS-3 suggests the possibility that these pseudogenes arose from a common intermediate.     

A structure-based simulation approach for electron paramagnetic resonance spectra using molecular and stochastic dynamics simulations

Electron paramagnetic resonance (EPR) spectroscopy using site-directed spin-labeling is an appropriate technique to analyze the structure and dynamics of flexible protein regions as well as protein-protein interactions under native conditions. The analysis of a set of protein mutants with consecutive spin-label positions leads to the identification of secondary and tertiary structure elements. In the first place, continuous-wave EPR spectra reflect the motional freedom of the spin-label specifically linked to a desired site within the protein. EPR spectra calculations based on molecular dynamics (MD) and stochastic dynamics simulations facilitate verification or refinement of predicted computer-aided models of local protein conformations. The presented spectra simulation algorithm implies a specialized in vacuo MD simulation at 600 K with additional restrictions to sample the entire accessible space of the bound spin-label without large temporal effort. It is shown that the distribution of spin-label orientations obtained from such MD simulations at 600 K agrees well with the extrapolated motion behavior during a long timescale MD at 300 K with explicit water. The following potential-dependent stochastic dynamics simulation combines the MD data about the site-specific orientation probabilities of the spin-label with a realistic rotational diffusion coefficient yielding a set of trajectories, each more than 700 ns long, essential to calculate the EPR spectrum. Analyses of a structural model of the loop between helices E and F of bacteriorhodopsin are illustrated to demonstrate the applicability and potentials of the reported simulation approach. Furthermore, effects on the motional freedom of bound spin-labels induced by solubilization of bacteriorhodopsin with Triton X-100 are examined.     

Molecular dynamics of synthetic leucine-serine ion channels in a phospholipid membrane

Molecular dynamics calculations were carried out on models of two synthetic leucine-serine ion channels: a tetrameric bundle with sequence (LSLLLSL)(3)NH(2) and a hexameric bundle with sequence (LSSLLSL)(3)NH(2). Each protein bundle is inserted in a palmitoyloleoylphosphatidylcholine bilayer membrane and solvated by simple point charge water molecules inside the pore and at both mouths. Both systems appear to be stable in the absence of an electric field during the 4 ns of molecular dynamics simulation. The water motion in the narrow pore of the four-helix bundle is highly restricted and may provide suitable conditions for proton transfer via a water wire mechanism. In the wider hexameric pore, the water diffuses much more slowly than in bulk but is still mobile. This, along with the dimensions of the pore, supports the observation that this peptide is selective for monovalent cations. Reasonable agreement of predicted conductances with experimentally determined values lends support to the validity of the simulations.     

Fast assignment of protein structures to sequences using the intermediate sequence library PDB-ISL

MOTIVATION: For large-scale structural assignment to sequences, as in computational structural genomics, a fast yet sensitive sequence search procedure is essential. A new approach using intermediate sequences was tested as a shortcut to iterative multiple sequence search methods such as PSI-BLAST. RESULTS: A library containing potential intermediate sequences for proteins of known structure (PDB-ISL) was constructed. The sequences in the library were collected from a large sequence database using the sequences of the domains of proteins of known structure as the query sequences and the program PSI-BLAST. Sequences of proteins of unknown structure can be matched to distantly related proteins of known structure by using pairwise sequence comparison methods to find homologues in PDB-ISL. Searches of PDB-ISL were calibrated, and the number of correct matches found at a given error rate was the same as that found by PSI-BLAST. The advantage of this library is that it uses pairwise sequence comparison methods, such as FASTA or BLAST2, and can, therefore, be searched easily and, in many cases, much more quickly than an iterative multiple sequence comparison method. The procedure is roughly 20 times faster than PSI-BLAST for small genomes and several hundred times for large genomes. AVAILABILITY: Sequences can be submitted to the PDB-ISL servers at http://stash.mrc-lmb.cam.ac.uk/PDB_ISL/ or http://cyrah.ebi.ac.uk:1111/Serv/PDB_ISL/ and can be downloaded from ftp://ftp.ebi.ac.uk/pub/contrib/jong/PDB_+ ++ISL/ CONTACT: sat@mrc-lmb.cam.ac.uk and jong@ebi.ac.uk     

Evolution of Germline-Limited Sequences in Two Populations of the Ciliate Chilodonella uncinata

Ciliates are microbial eukaryotes that separate their nuclear functions into a germline micronucleus and a somatic macronucleus. During development of the macronucleus the genome undergoes a series of reorganization events that includes the precise excision of intervening DNA. Here, we determine the architecture of four loci in the micronuclear and macronuclear genomes of the ciliate Chilodonella uncinata and compare the levels of variation in micronuclear-limited sequences to macronuclear destined sequences at two of these loci. We find that within a population, germline-limited sequences are evolving at the same rate as other putatively neutral sites, but between populations germline-limited sequences are accumulating mutations at a much faster rate than other sites. We also find evidence of macronuclear recombination and incomplete elimination of intervening DNA, which result in increased diversity in the macronuclear genome. Our results support the assertion that the unusual genomic features of ciliates can result in rapid and unpredicted patterns of diversification.     

Relocation or duplication of the helix A sequence of T4 lysozyme causes only modest changes in structure but can increase or decrease the rate of folding

In T4 lysozyme, helix A is located at the amino terminus of the sequence but is associated with the C-terminal domain in the folded structure. To investigate the implications of this arrangement for the folding of the protein, we first created a circularly permuted variant with a new amino terminus at residue 12. In effect, this moves the sequence corresponding to helix A from the N- to the C-terminus of the molecule. The protein crystallized nonisomorphously with the wild type but has a very similar structure, showing that the unit consisting of helix A and the C-terminal domain can be reconstituted from a contiguous polypeptide chain. The protein is less stable than the wild type but folds slightly faster. We then produced a second variant in which the helix A sequence was appended at the C-terminus (as in the first variant), but was also restored at the N-terminus (as in the wild type). This variant has two helix A sequences, one at the N-terminus and the other at the C-terminus, each of which can compete for the same site in the folded protein. The crystal structure shows that it is the N-terminal sequence that folds in a manner similar to that of the wild type, whereas the copy at the C-terminus is forced to loop out. The stability of this protein is much closer to that of the wild type, but its rate of folding is significantly slower. The reduction in rate is attributed to the presence of the two identical sequence segments which compete for a single, mutually exclusive, site.     

Improving the accuracy of protein stability predictions with multistate design using a variety of backbone ensembles

Multistate computational protein design (MSD) with backbone ensembles approximating conformational flexibility can predict higher quality sequences than single‐state design with a single fixed backbone. However, it is currently unclear what characteristics of backbone ensembles are required for the accurate prediction of protein sequence stability. In this study, we aimed to improve the accuracy of protein stability predictions made with MSD by using a variety of backbone ensembles to recapitulate the experimentally measured stability of 85 Streptococcal protein G domain β1 sequences. Ensembles tested here include an NMR ensemble as well as those generated by molecular dynamics (MD) simulations, by Backrub motions, and by PertMin, a new method that we developed involving the perturbation of atomic coordinates followed by energy minimization. MSD with the PertMin ensembles resulted in the most accurate predictions by providing the highest number of stable sequences in the top 25, and by correctly binning sequences as stable or unstable with the highest success rate (≈90%) and the lowest number of false positives. The performance of PertMin ensembles is due to the fact that their members closely resemble the input crystal structure and have low potential energy. Conversely, the NMR ensemble as well as those generated by MD simulations at 500 or 1000 K reduced prediction accuracy due to their low structural similarity to the crystal structure. The ensembles tested herein thus represent on‐ or off‐target models of the native protein fold and could be used in future studies to design for desired properties other than stability. Proteins 2014; 82:771–784. © 2013 Wiley Periodicals, Inc.     

Temperature and urea induced denaturation of the TRP‐cage mini protein TC5b: A simulation study consistent with experimental observations

The effects of temperature and urea denaturation (6M urea) on the dominant structures of the 20‐residue Trp‐cage mini‐protein TC5b are investigated by molecular dynamics simulations of the protein at different temperatures in aqueous and in 6M urea solution using explicit solvent degrees of freedom and the GROMOS force‐field parameter set 45A3. In aqueous solution at 278 K, TC5b is stable throughout the 20 ns of MD simulation and the trajectory structures largely agree with the NMR‐NOE atom–atom distance data available. Raising the temperature to 360 K and to 400 K, the protein denatures within 22 ns and 3 ns, showing that the denaturation temperature is well below 360 K using the GROMOS force field. This is 40–90 K lower than the denaturation temperatures observed in simulations using other much used protein force fields. As the experimental denaturation temperature is about 315 K, the GROMOS force field appears not to overstabilize TC5b, as other force fields and the use of continuum solvation models seem to do. This feature may directly stem from the GROMOS force‐field parameter calibration protocol, which primarily involves reproduction of condensed phase thermodynamic quantities such as energies, densities, and solvation free energies of small compounds representative for protein fragments. By adding 6M urea to the solution, the onset of denaturation is observed in the simulation, but is too slow to observe a particular side‐chain side‐chain contact (Trp6‐Ile4) that was experimentally observed to be characteristic for the denatured state. Interestingly, using temperature denaturation, the process is accelerated and the experimental data are reproduced.     

Variations on probabilistic suffix trees: statistical modeling and prediction of protein families

MOTIVATION: We present a method for modeling protein families by means of probabilistic suffix trees (PSTs). The method is based on identifying significant patterns in a set of related protein sequences. The patterns can be of arbitrary length, and the input sequences do not need to be aligned, nor is delineation of domain boundaries required. The method is automatic, and can be applied, without assuming any preliminary biological information, with surprising success. Basic biological considerations such as amino acid background probabilities, and amino acids substitution probabilities can be incorporated to improve performance. RESULTS: The PST can serve as a predictive tool for protein sequence classification, and for detecting conserved patterns (possibly functionally or structurally important) within protein sequences. The method was tested on the Pfam database of protein families with more than satisfactory performance. Exhaustive evaluations show that the PST model detects much more related sequences than pairwise methods such as Gapped-BLAST, and is almost as sensitive as a hidden Markov model that is trained from a multiple alignment of the input sequences, while being much faster.     

Electronic excitation transfer in the complex of lumazine protein with bacterial bioluminescence intermediates

Fluorescence dynamics measurements have been made on the bioluminescence reaction intermediates using Photobacterium leiognathi, Vibrio fischeri, and Vibrio harveyi luciferases, both alone and in mixtures with Photobacterium phosphoreum lumazine protein. Each luciferase produces a "fluorescent transient" intermediate on reaction with the bioluminescence substrates, FMNH2, tetradecanal, and O2, and all have a fluorescence quantum yield about 0.3, with a predominant lifetime around 10 ns. The P. leiognathi luciferase fluorescent transient has a rotational correlation time of 79 ns at 2 degrees C, as expected for the rotational diffusion of a 77-kDa macromolecule. In the presence of lumazine protein however a faster correlation time of about 3 ns predominates. This rapid channel of anisotropy loss is attributed to energy transfer from the flavin intermediate bound on the luciferase to the lumazine ligand, reflects the presence of protein-protein complexation, and is greatest in the case of P. leiognathi, but not at all for V. fischeri. This fact is consistent with the strong influence of lumazine protein on the bioluminescence reaction of P. leiognathi, and not at all with V. fischeri. The rate of energy transfer is of order 10(9) s-1, much greater than the 10(8) s-1 fluorescence rate of the donor. Thus the bioluminescence excitation of lumazine protein could occur by a similar photophysical mechanism of interprotein energy transfer from a chemically excited fluorescent transient donor to the lumazine acceptor.     

An efficient algorithm for the identification of structured motifs in DNA promoter sequences

We propose a new algorithm for identifying cis-regulatory modules in genomic sequences. The proposed algorithm, named RISO, uses a new data structure, called box-link, to store the information about conserved regions that occur in a well-ordered and regularly spaced manner in the data set sequences. This type of conserved regions, called structured motifs, is extremely relevant in the research of gene regulatory mechanisms since it can effectively represent promoter models. The complexity analysis shows a time and space gain over the best known exact algorithms that is exponential in the spacings between binding sites. A full implementation of the algorithm was developed and made available online. Experimental results show that the algorithm is much faster than existing ones, sometimes by more than four orders of magnitude. The application of the method to biological data sets shows its ability to extract relevant consensi.     

Isolation of an active gene encoding human hnRNP protein A1. Evidence for alternative splicing

Heterogeneous nuclear ribonucleoprotein (hnRNP) core protein A1 is a major component of mammalian hnRNP 40 S particles. We describe the structure of an active A1 gene and report on the partial characterization of the A1 gene family. About 30 A1-specific sequences are present per haploid human genome: 15 such sequences were isolated from a human genomic DNA library. Many corresponded to pseudogenes of the processed type but by applying a selection for actively transcribed regions we isolated an active A1 gene. The gene spans a region of 4.6 x 10(3) base-pairs and it is split into ten exons that encode the 320 amino acid residues of the protein. The amino acid sequence derived from the exon sequences is identical with that deduced from cDNA and reported for the protein. One intron exactly separates the two structural domains that constitute the protein. Each of the two RNA-binding domains in protein A1 is encoded by one exon. Experimental evidence indicates that the A1 gene can encode for more than one protein by alternative splicing. The gene is preceded by a strong promoter that contains at least two CCAAT boxes and two possible Sp1 binding sites, but it lacks a TATA box.     

Mechanism of Interaction between the General Anesthetic Halothane and a Model Ion Channel Protein, II: Fluorescence and Vibrational Spectroscopy Using a Cyanophenylalanine Probe

We demonstrate that cyano-phenylalanine (Phe_CN) can be utilized to probe the binding of the inhalational anesthetic halothane to an anesthetic-binding, model ion channel protein hbAP-Phe_CN. The Trp to Phe_CN mutation alters neither the α-helical conformation nor the 4-helix bundle structure. The halothane binding properties of this Phe_CN mutant hbAP-Phe_CN, based on fluorescence quenching, are consistent with those of the prototype, hbAP1. The dependence of fluorescence lifetime as a function of halothane concentration implies that the diffusion of halothane in the nonpolar core of the protein bundle is one-dimensional. As a consequence, at low halothane concentrations, the quenching of the fluorescence is dynamic, whereas at high concentrations the quenching becomes static. The 4-helix bundle structure present in aqueous detergent solution and at the air-water interface, is preserved in multilayer films of hbAP-Phe_CN, enabling vibrational spectroscopy of both the protein and its nitrile label (-CN). The nitrile groups' stretching vibration band shifts to higher frequency in the presence of halothane, and this blue-shift is largely reversible. Due to the complexity of this amphiphilic 4-helix bundle model membrane protein, where four Phe_CN probes are present adjacent to the designed cavity forming the binding site within each bundle, all contributing to the infrared absorption, molecular dynamics (MD) simulation is required to interpret the infrared results. The MD simulations indicate that the blue-shift of -CN stretching vibration induced by halothane arises from an indirect effect, namely an induced change in the electrostatic protein environment averaged over the four probe oscillators, rather than a direct interaction with the oscillators. hbAP-Phe_CN therefore provides a successful template for extending these investigations of the interactions of halothane with the model membrane protein via vibrational spectroscopy, using cyano-alanine residues to form the anesthetic binding cavity.