All-atom empirical potential for molecular modeling and dynamics studies of proteins

Continuum and atomistic modeling of ion partitioning into a peptide nanotube

Continuum and atomistic descriptions of the partitioning of ions into a self-assembled (D,L)-octapeptide nanotube, cyclo[-(L-Ala-D-Ala)(4)-], are presented. Perturbation free energy calculations, including Ewald electrostatics, are used to estimate the electrostatic component of the excess free energy of charging Li(+), Na(+), Rb(+), and Cl(minus sign) ions inside the nanotube. The radial density and orientational distribution of water around the ion is calculated for the ion at two different positions inside the tube; it is seen that the calculated distributions are sensitive to the location of the ions. Two different continuum electrostatic models are formulated to describe the ion solvation inside the nanotube. When enhanced orientational structuring of water dipoles is evidenced, explicitly including the first solvation shell as part of the low dielectric nanotube environment provides good agreement with molecular dynamics simulations. When water orientational structuring is as in the reference bulk solvent, we find that treating the first shell water explicitly or as a high dielectric continuum leads to similar results. These results are discussed, and their importance for continuum electrostatic modeling of ion channels are highlighted.

Excitons in a photosynthetic light-harvesting system: a combined molecular dynamics, quantum chemistry, and polaron model study

The dynamics of pigment-pigment and pigment-protein interactions in light-harvesting complexes is studied with an approach that combines molecular dynamics simulations with quantum chemistry calculations and a polaron model analysis. The molecular dynamics simulation of light-harvesting (LH) complexes was performed on an 87 055 atom system comprised of a LH-II complex of Rhodospirillum molischianum embedded in a lipid bilayer and surrounded with appropriate water layers. For each of the 16 B850 bacteriochlorophylls (BChls), we performed 400 ab initio quantum chemistry calculations on geometries that emerged from the molecular dynamical simulations, determining the fluctuations of pigment excitation energies as a function of time. From the results of these calculations we construct a time-dependent Hamiltonian of the B850 exciton system from which we determine within linear response theory the absorption spectrum. Finally, a polaron model is introduced to describe both the excitonic and coupled phonon degrees of freedom by quantum mechanics. The exciton-phonon coupling that enters into the polaron model, and the corresponding phonon spectral function, are derived from the molecular dynamics and quantum chemistry simulations. The model predicts that excitons in the B850 BChl ring are delocalized over five pigments at room temperature. Also, the polaron model permits the calculation of the absorption and circular dichroism spectra of the B850 excitons from the sole knowledge of the autocorrelation function of the excitation energies of individual BChls, which is readily available from the combined molecular dynamics and quantum chemistry simulations. The obtained results are found to be in good agreement with the experimentally measured absorption and circular dichroism spectra.

Molecular mechanisms of calcium and magnesium binding to parvalbumin

Molecular dynamics simulations have been used to investigate the relationship between the coordinating residues of the EF-hand calcium binding loop of parvalbumin and the overall plasticity and flexibility of the protein. The first simulation modeled the transition from Ca(2+) to Mg(2+) coordination by varying the van der Waals parameters for the bound metal ions. The glutamate at position 12 could be accurately and reversibly seen to be a source of selective bidentate ligation of Ca(2+) in the simulations. A second simulation correlated well with the experimental observation that an E101D substitution at EF loop position 12 results in a dramatically less tightly bound monodentate Ca(2+) coordination by aspartate. A final set of simulations investigated Ca(2+) binding in the E101D mutant loop in the presence of applied external forces designed to impose bidentate coordination. The results of these simulations illustrate that the aspartate is capable of attaining a suitable orientation for bidentate coordination, thus implying that it is the inherent rigidity of the loop that prevents bidentate coordination in the parvalbumin E101D mutant.

A structural model for force regulated integrin binding to fibronectin's RGD-synergy site

The synergy site on fibronectin's FN-III(9) module, located approximately 32 A away from the RGD-loop on FN-III(10), greatly enhances integrin alpha(5)beta(1) mediated cell binding. Since fibronectin is exposed to mechanical forces acting on the extracellular matrix in vivo, we used steered molecular dynamics to study how mechanical stretching of FN-III(9-10) affects the relative distance between these two synergistic sites. Our simulations predict the existence of an intermediate state prior to unfolding. In this state, the synergy-RGD distance is increased from 32 A to approximately 55 A, while the conformations of both sites remain unperturbed. This distance is too large for both sites to co-bind the same receptor, as indicated by experiments that confirm that increasing the length of the linker chain between FN-III(9) and FN-III(10) reduces alpha(5)beta(1) binding. Our simulations thus suggest that increased alpha(5)beta(1)-binding attributed to the synergy site, along with the associated downstream cell-signaling events, can be turned off mechanically by stretching FN-III(9-10) into this intermediate state. The potential physiological implications are discussed.

Temperature dependence of protein dynamics: computer simulation analysis of neutron scattering properties

The temperature dependence of the internal dynamics of an isolated protein, bovine pancreatic trypsin inhibitor, is examined using normal mode analysis and molecular dynamics (MD) simulation. It is found that the protein exhibits marked anharmonic dynamics at temperatures of approximately 100-120 K, as evidenced by departure of the MD-derived average mean square displacement from that of the harmonic model. This activation of anharmonic dynamics is at lower temperatures than previously detected in proteins and is found in the absence of solvent molecules. The simulation data are also used to investigate neutron scattering properties. The effects are determined of instrumental energy resolution and of approximations commonly used to extract mean square displacement data from elastic scattering experiments. Both the presence of a distribution of mean square displacements in the protein and the use of the Gaussian approximation to the dynamic structure factor lead to quantified underestimation of the mean square displacement obtained.

New space warping method for the simulation of large-scale macromolecular conformational changes

A space warping method, facilitating the modeling of large-scale conformational changes in mesoscopic systems, is presented. The method uses a set of "global (or collective) coordinates" that capture overall behavior, in conjunction with the set of atomic coordinates. Application of the space warping method to energy minimization is discussed. Several simulations where the method is used to determine the energy minimizing structures of simple central force systems are analyzed. Comparing the results and behavior of the space warping method to simulations involving atomic coordinates only, it is found that the space warping method scales better with system size and also finds lower minima when the potential energy surface has multiple minima. It is shown that the transformation of [Ala16]+ in vacuo from linear to globular is captured efficiently using the space warping method.

Imaging the electrostatic potential of transmembrane channels: atomic probe microscopy of OmpF porin

The atomic force microscope (AFM) was used to image native OmpF porin and to detect the electrostatic potential generated by the protein. To this end the OmpF porin trimers from Escherichia coli was reproducibly imaged at a lateral resolution of approximately 0.5 nm and a vertical resolution of approximately 0.1 nm at variable electrolyte concentrations of the buffer solution. At low electrolyte concentrations the charged AFM probe not only contoured structural details of the membrane protein surface but also interacted with local electrostatic potentials. Differences measured between topographs recorded at variable ionic strength allowed mapping of the electrostatic potential of OmpF porin. The potential map acquired by AFM showed qualitative agreement with continuum electrostatic calculations based on the atomic OmpF porin embedded in a lipid bilayer at the same electrolyte concentrations. Numerical simulations of the experimental conditions showed the measurements to be reproduced quantitatively when the AFM probe was included in the calculations. This method opens a novel avenue to determine the electrostatic potential of native protein surfaces at a lateral resolution better than 1 nm and a vertical resolution of approximately 0.1 nm.

Role of protein flexibility in enzymatic catalysis: quantum mechanical-molecular mechanical study of the deacylation reaction in class A beta-lactamases

We present a theoretical study of a mechanism for the hydrolysis of the acyl-enzyme complex formed by a class A beta-lactamase (TEM1) and an antibiotic (penicillanate), as a part of the process of antibiotic's inactivation by this type of enzymes. In the presented mechanism the carboxylate group of a particular residue (Glu166) activates a water molecule, accepting one of its protons, and afterward transfers this proton directly to the acylated serine residue (Ser70). In our study we employed a quantum mechanics (AM1)-molecular mechanics partition scheme (QM/MM) where all the atoms of the system were allowed to relax. For this purpose we used the GRACE procedure in which part of the system is used to define the Hessian matrix while the rest is relaxed at each step of the stationary structures search. By use of this computational scheme, the hydrolysis of the acyl-enzyme is described as a three-step process: The first step corresponds to the proton transfer from the hydrolytic water molecule to the carboxylate group of Glu166 and the subsequent formation of a tetrahedral adduct as a consequence of the attack of this activated water molecule to the carbonyl carbon atom of the beta-lactam. In the second step, the acyl-enzyme bond is broken, obtaining a negatively charged Ser70. In the last step this residue is protonated by means of a direct proton transfer from Glu166. The large mobility of Glu166, a residue that is placed in a Ohms-loop, is essential to facilitate this mechanism. The geometry of the acyl-enzyme complex shows a large distance between Glu166 and Ser70 and thus, if protein coordinates were kept frozen during the reaction path, it would be difficult to get a direct proton transfer between these two residues. This computational study shows how a flexible treatment suggests the feasibility of a mechanism that could have been discounted on the basis of crystallographic positions.

Docking multiple conformations of a flexible ligand into a protein binding site using NMR restraints

A method is described for docking a large, flexible ligand using intra-ligand conformational restraints from exchange-transferred NOE (etNOE) data. Numerous conformations of the ligand are generated in isolation, and a subset of representative conformations is selected. A crude model of the protein-ligand complex is used as a template for overlaying the selected ligand structures, and each complex is conformationally relaxed by molecular mechanics to optimize the interaction. Finally, the complexes were assessed for structural quality. Alternative approaches are described for the three steps of the method: generation of the initial docking template; selection of a subset of ligand conformations; and conformational sampling of the complex. The template is generated either by manual docking using interactive graphics or by a computational grid-based search of the binding site. A subset of conformations from the total number of peptides calculated in isolation is selected based on either low energy and satisfaction of the etNOE restraints, or a cluster analysis of the full set. To optimize the interactions in the complex, either a restrained Monte Carlo-energy minimization (MCM) protocol or a restrained simulated annealing (SA) protocol were used. This work produced 53 initial complexes of which 8 were assessed in detail. With the etNOE conformational restraints, all of the approaches provide reasonable models. The grid-based approach to generate an initial docking template allows a large volume to be sampled, and as a result, two distinct binding modes were identified for a fifteen-residue peptide binding to an enzyme active site.

Structural basis of RasGRP binding to high-affinity PKC ligands

The Ras guanyl releasing protein RasGRP belongs to the CDC25 class of guanyl nucleotide exchange factors that regulate Ras-related GTPases. These GTPases serve as switches for the propagation and divergence of signaling pathways. One interesting feature of RasGRP is the presence of a C-terminal C1 domain, which has high homology to the PKC C1 domain and binds to diacylglycerol (DAG) and phorbol esters. RasGRP thus represents a novel, non-kinase phorbol ester receptor. In this paper, we investigate the binding of indolactam(V) (ILV), 7-(n-octyl)-ILV, 8-(1-decynyl)benzolactam(V) (benzolactam), and 7-methoxy-8-(1-decynyl)benzolactam(V) (methoxylated benzolactam) to RasGRP through both experimental binding assays and molecular modeling studies. The binding affinities of these lactams to RasGRP are within the nanomolar range. Homology modeling was used to model the structure of the RasGRP C1 domain (C1-RasGRP), which was subsequently used to model the structures of C1-RasGRP in complex with these ligands and phorbol 13-acetate using a computational docking method. The structural model of C1-RasGRP exhibits a folding pattern that is nearly identical to that of C1b-PKCdelta and is comprised of three antiparallel-strand beta-sheets capped against a C-terminal alpha-helix. Two loops A and B comprising residues 8-12 and 21-27 form a binding pocket that has some positive charge character. The ligands phorbol 13-acetate, benzolactam, and ILV are recognized by C1-RasGRP through a number of hydrogen bonds with loops A and B. In the models of C1-RasGRP in complex with phorbol 13-acetate, benzolactam, and ILV, common hydrogen bonds are formed with two residues Thr12 and Leu21, whereas other hydrogen bond interactions are unique for each ligand. Furthermore, our modeling results suggest that the shallower insertion of ligands into the binding pocket of C1-RasGRP compared to C1b-PKCdelta may be due to the presence of Phe rather than Leu at position 20 in C1-RasGRP. Taken together, our experimental and modeling studies provide us with a better understanding of the structural basis of the binding of PKC ligands to the novel phorbol ester receptor RasGRP.

Tetraethylammonium binding to the outer mouth of the KcsA potassium channel: implications for ion permeation

Extracellular tetraethylammonium (TEA+) inhibits the current carried out by K+ ions in potassium channels. Structural models of wild-type (WT) and Y82C KcsA K+ channel/TEA+ complexes are here built using docking procedures, electrostatics calculations and molecular dynamics simulations. The calculations are based on the structure determined by Doyle et al. (11) Our calculations suggest that in WT, the TEA+ cation turns binds at the outer mouth of the selectivity filter, stabilized by electrostatic and hydrophobic interactions with the four Tyr82 side chains. Replacement of Tyr82 with Cys causes a decrease of the affinity of the cation for the channel, consistently with the available site-directed mutagenesis data (16). An MD simulation in which K+ replaces TEA+ provides evidence that TEA+ binding site can accommodate a potassium ion, in agreement with the high-resolution structure recently reported by Zhou et al. (20)

Dynamic and structural analysis of isotropically distributed molecular ensembles

An efficient new method is presented for the characterization of motional correlations derived from a set of protein structures without requiring the separation of overall and internal motion. In this method, termed isotropically distributed ensemble (IDE) analysis, each structure is represented by an ensemble of isotropically distributed replicas corresponding to the situation found in an isotropic protein solution. This leads to a covariance matrix of the cartesian atomic positions with elements proportional to the ensemble average of scalar products of the position vectors with respect to the center of mass. Diagonalization of the covariance matrix yields eigenmodes and amplitudes that describe concerted motions of atoms, including overall rotational and intramolecular dynamics. It is demonstrated that this covariance matrix naturally distinguishes between "rigid" and "mobile" parts without necessitating a priori selection of a reference structure and an atom set for the orientational alignment process. The method was applied to the analysis of a 5-ns molecular dynamics trajectory of native ubiquitin and a 40-ns trajectory of a partially folded state of ubiquitin. The results were compared with essential dynamics analysis. By taking advantage of the spherical symmetry of the IDE covariance matrix, more than a 10-fold speed up is achieved for the computation of eigenmodes and mode amplitudes. IDE analysis is particularly suitable for studying the correlated dynamics of flexible and large molecules.

The ionization state and the conformation of Glu-71 in the KcsA K(+) channel

The side chain of Glu-71 of the KcsA K(+) channel, an important residue in the vicinity of the selectivity filter, was not resolved in the crystallographic structure of Doyle et al. (Doyle, D. A., J. M. Cabral, R. A. Pfuetzner, A. Kuo, J. M. Gulbis, S. L. Cohen, B. T. Chait, and R. MacKinnon. 1998. Science. 280:69-77). Its atomic coordinates are undetermined and its ionization state is unknown. For meaningful theoretical and computational studies of the KcsA K(+) channel, it is essential to address questions about the conformation and the ionization state of this residue in detail. In previous MD simulations in which the side chain of Glu-71 is protonated and forming a strong hydrogen bond with Asp-80 it was observed that the channel did not deviate significantly from the crystallographic structure (Bernèche, S., and B. Roux. 2000. Biophys. J. 78:2900-2917). In contrast, we show here that the structure of the selectivity filter of the KcsA channel is significantly disrupted when these side chains are fully ionized on each of the four monomers. To further resolve questions about the ionization state of Glu-71 we calculated the pK(a) value of this residue using molecular dynamics free energy simulations (MD/FES) with a fully flexible system including explicit solvent and membrane and finite-difference Poisson-Boltzmann (PB) continuum electrostatics. It is found that the pK(a) of Glu-71 is shifted by approximately +10 pK(a) units. These results strongly suggest that Glu-71 is protonated under normal conditions.

Side-chain modeling with an optimized scoring function

Modeling side-chain conformations on a fixed protein backbone has a wide application in structure prediction and molecular design. Each effort in this field requires decisions about a rotamer set, scoring function, and search strategy. We have developed a new and simple scoring function, which operates on side-chain rotamers and consists of the following energy terms: contact surface, volume overlap, backbone dependency, electrostatic interactions, and desolvation energy. The weights of these energy terms were optimized to achieve the minimal average root mean square (rms) deviation between the lowest energy rotamer and real side-chain conformation on a training set of high-resolution protein structures. In the course of optimization, for every residue, its side chain was replaced by varying rotamers, whereas conformations for all other residues were kept as they appeared in the crystal structure. We obtained prediction accuracy of 90.4% for chi(1), 78.3% for chi(1 + 2), and 1.18 A overall rms deviation. Furthermore, the derived scoring function combined with a Monte Carlo search algorithm was used to place all side chains onto a protein backbone simultaneously. The average prediction accuracy was 87.9% for chi(1), 73.2% for chi(1 + 2), and 1.34 A rms deviation for 30 protein structures. Our approach was compared with available side-chain construction methods and showed improvement over the best among them: 4.4% for chi(1), 4.7% for chi(1 + 2), and 0.21 A for rms deviation. We hypothesize that the scoring function instead of the search strategy is the main obstacle in side-chain modeling. Additionally, we show that a more detailed rotamer library is expected to increase chi(1 + 2) prediction accuracy but may have little effect on chi(1) prediction accuracy.

Protonation and stability of the globular domain of influenza virus hemagglutinin

A partial dissociation of the HA1 subunits of influenza virus hemagglutinin (HA) is considered to be the initial step of conformational changes of the HA ectodomain leading to a membrane fusion active conformation (L. Godley, J. Pfeifer, D. Steinhauer, B. Ely, G. Shaw, R. Kaufman, E. Suchanek, C. Pabo, J.J. Skehel, D.C. Wiley, and S. Wharton, 1992, Cell 68:635-645; G.W. Kemble, D.L.Bodian, J. Rose, I.A. Wilson, and J.M. White, 1992, J. Virol. 66:4940-4950). Here, we explore a mechanism that provides an understanding of the physical and chemical basis for such dissociation and relies on two essential observations. First, based on the x-ray structure of HA from X31 (I.A. Wilson, J.J. Skehel, and D.C. Wiley, 1981, Nature 289:366-373), and by employing techniques of molecular modeling, we show that the protonation of the HA1 subunits is enhanced at the conditions known to trigger conformational changes of the HA ectodomain. Second, we found that the dependence of the calculated relative degree of protonation of the HA1 domain on temperature and pH is similar to that observed experimentally for the conformational change of HA assessed by proteinase K sensitivity. We suggest that at the pH-temperature conditions typical for the conformational change of HA and membrane fusion, dissociation of the HA1 subunits is caused by the enhanced protonation of the HA1 subunits leading to an increase in the positive net charge of these subunits and, in turn, to a weakened attraction between them.

Combined ab initio/empirical approach for optimization of Lennard-Jones parameters for polar-neutral compounds

The study of small functionalized organic molecules in aqueous solution is a useful step toward gaining a basic understanding of the behavior of biomolecular systems in their native aqueous environment. Interest in studying amines and fluorine-substituted compounds has risen from their intrinsic physicochemical properties and their prevalence in biological and pharmaceutical compounds. In the present study, a previously developed approach which optimizes Lennard-Jones (LJ) parameters via the use of rare gas atoms combined with the reproduction of experimental condensed phase properties was extended to polar-neutral compounds. Compounds studied included four amines (ammonia, methylamine, dimethylamine, and trimethylamine) and three fluoroethanes (1-fluoroethane, 1,1-difluoroethane, and 1,1,1-trifluoroethane). The resulting force field yielded heats of vaporization and molecular volumes in excellent agreement with the experiment, with average differences less than 1%. The current amine CHARMM parameters successfully reproduced experimental aqueous solvation data where methylamine is more hydrophilic than ammonia, with hydrophobicity increasing with additional methylation on the nitrogen. For both the amines and fluoroethanes the parabolic relationship of the extent of methylation or fluorination, respectively, to the heats of vaporization were reproduced by the new parameters. The present results are also discussed with respect to the impact of parameterization approach to molecular details obtained from computer simulations and to the unique biological properties of fluorine in pharmaceutical compounds.

X-ray structure of a ClC chloride channel at 3.0 A reveals the molecular basis of anion selectivity

The ClC chloride channels catalyse the selective flow of Cl- ions across cell membranes, thereby regulating electrical excitation in skeletal muscle and the flow of salt and water across epithelial barriers. Genetic defects in ClC Cl- channels underlie several familial muscle and kidney diseases. Here we present the X-ray structures of two prokaryotic ClC Cl- channels from Salmonella enterica serovar typhimurium and Escherichia coli at 3.0 and 3.5 A, respectively. Both structures reveal two identical pores, each pore being formed by a separate subunit contained within a homodimeric membrane protein. Individual subunits are composed of two roughly repeated halves that span the membrane with opposite orientations. This antiparallel architecture defines a selectivity filter in which a Cl- ion is stabilized by electrostatic interactions with alpha-helix dipoles and by chemical coordination with nitrogen atoms and hydroxyl groups. These findings provide a structural basis for further understanding the function of ClC Cl- channels, and establish the physical and chemical basis of their anion selectivity.

Combined QM/MM study of the opsin shift in bacteriorhodopsin

Combined quantum mechanical and molecular mechanical (QM/MM) calculations and molecular dynamics simulations of bacteriorhodopsin (bR) in the membrane matrix have been carried out to determine the factors that make significant contributions to the opsin shift. We found that both solvation and interactions with the protein significantly shifts the absorption maximum of the retinal protonated Schiff base, but the effects are much more pronounced in polar solvents such as methanol, acetonitrile, and water than in the protein environment. The differential solvatochromic shifts of PSB in methanol and in bR leads to a bathochromic shift of about 1800 cm(-1). Because the combined QM/MM configuration interaction calculation is essentially a point charge model, this contribution is attributed to the extended point-charge model of Honig and Nakanishi. The incorporation of retinal in bR is accompanied by a change in retinal conformation from the 6-s-cis form in solution to the 6-s-trans configuration in bR. The extension of the pi-conjugated system further increases the red-shift by 2400 cm(-1). The remaining factors are due to the change in dispersion interactions. Using an estimate of about 1000 cm(-1) in the dispersion contribution by Houjou et al., we obtained a theoretical opsin shift of 5200 cm(-1) in bR, which is in excellent agreement with the experimental value of 5100 cm(-1). Structural analysis of the PSB binding site revealed the specific interactions that make contributions to the observed opsin shift. The combined QM/MM method used in the present study provides an opportunity to accurately model the photoisomerization and proton transfer reactions in bR.

Relation between sequence and structure of HIV-1 protease inhibitor complexes: a model system for the analysis of protein flexibility

The flexibility of different regions of HIV-1 protease was examined by using a database consisting of 73 X-ray structures that differ in terms of sequence, ligands or both. The root-mean-square differences of the backbone for the set of structures were shown to have the same variation with residue number as those obtained from molecular dynamics simulations, normal mode analyses and X-ray B-factors. This supports the idea that observed structural changes provide a measure of the inherent flexibility of the protein, although specific interactions between the protease and the ligand play a secondary role. The results suggest that the potential energy surface of the HIV-1 protease is characterized by many local minima with small energetic differences, some of which are sampled by the different X-ray structures of the HIV-1 protease complexes. Interdomain correlated motions were calculated from the structural fluctuations and the results were also in agreement with molecular dynamics simulations and normal mode analyses. Implications of the results for the drug-resistance engendered by mutations are discussed briefly.

Calculation of protein ionization equilibria with conformational sampling: pK(a) of a model leucine zipper, GCN4 and barnase

The use of conformational ensembles provided by nuclear magnetic resonance (NMR) experiments or generated by molecular dynamics (MD) simulations has been regarded as a useful approach to account for protein motions in the context of pK(a) calculations, yet the idea has been tested occasionally. This is the first report of systematic comparison of pK(a) estimates computed from long multiple MD simulations and NMR ensembles. As model systems, a synthetic leucine zipper, the naturally occurring coiled coil GCN4, and barnase were used. A variety of conformational averaging and titration curve-averaging techniques, or combination thereof, was adopted and/or modified to investigate the effect of extensive global conformational sampling on the accuracy of pK(a) calculations. Clustering of coordinates is proposed as an approach to reduce the vast diversity of MD ensembles to a few structures representative of the average electrostatic properties of the system in solution. Remarkable improvement of the accuracy of pK(a) predictions was achieved by the use of multiple MD simulations. By using multiple trajectories the absolute error in pK(a) predictions for the model leucine zipper was reduced to as low as approximately 0.25 pK(a) units. The validity, advantages, and limitations of explicit conformational sampling by MD, compared with the use of an average structure and a high internal protein dielectric value as means to improve the accuracy of pK(a) calculations, are discussed.

Exploratory studies of ab initio protein structure prediction: multiple copy simulated annealing, AMBER energy functions, and a generalized born/solvent accessibility solvation model

A theoretical and computational approach to ab initio structure prediction for polypeptides in water is described and applied to selected amino acid sequences for testing and preliminary validation. The method builds systematically on the extensive efforts applied to parameterization of molecular dynamics (MD) force fields, employs an empirically well-validated continuum dielectric model for solvation, and an eminently parallelizable approach to conformational search. The effective free energy of polypeptide chains is estimated from AMBER united atom potential functions, with internal degrees of freedom for both backbone and amino acid side chains explicitly treated. The hydration free energy of each structure is determined using the Generalized Born/Solvent Accessibility (GBSA) method, modified and reparameterized to include atom types consistent with the AMBER force field. The conformational search procedure employs a multiple copy, Monte Carlo simulated annealing (MCSA) protocol in full torsion angle space, applied iteratively on sets of structures of progressively lower free energy until a prediction of a structure with lowest effective free energy is obtained. Calibration tests for the effective energy function and search algorithm are performed on the alanine dipeptide, selected protein crystal structures, and united atom decoys on barnase, crambin, and six examples from the Rosetta set. Specific demonstration cases of the method are provided for the 8-mer sequence of Ala residues, a 12-residue peptide with longer side chains QLLKKLLQQLKQ, a de novo designed 16 residue peptide of sequence (AAQAA)3Y, a 15-residue sequence with a beta sheet motif, GEWTWDATKTFTVTE, and a 36 residue small protein, Villin headpiece. The Ala 8-mer readily formed an alpha-helix. An alpha-helix structure was predicted for the 16-mer, consistent with observed results from IR and CD spectroscopy and with the pattern in psi/straight phi angles of known protein structures. The predicted structure for the 12-mer, composed of a mix of helix and less regular elements of secondary structure, lies 2.65 A RMS from the observed crystal structure. Structure prediction for the 8-mer beta-motif resulted in form 4.50 A RMS from the crystal geometry. For Villin, the predicted native form is very close to the crystal structure, RMS values of 3.5 A (including sidechains), and 1.01 A (main chain only). The methodology permits a detailed analysis of the molecular forces which dominate various segments of the predicted folding trajectory. Analysis of the results in terms of internal torsional, electrostatic and van der Waals and the electrostatic and non-electrostatic contributions to hydration, including the hydrophobic effect, is presented.

Use of MCSS to design small targeted libraries: application to picornavirus ligands

Computational methods were used to design structure-based combinatorial libraries of antipicornaviral capsid-binding ligands. The multiple copy simultaneous search (MCSS) program was employed to calculate functionality maps for many diverse functional groups for both the poliovirus and rhinovirus capsid structures in the region of the known drug binding pocket. Based on the results of the MCSS calculations, small combinatorial libraries consisting of 10s or 100s of three-monomer compounds were designed and synthesized. Ligand binding was demonstrated by a noncell-based mass spectrometric assay, a functional immuno-precipitation assay, and crystallographic analysis of the complexes of the virus with two of the candidate ligands. The P1/Mahoney poliovirus strain was used in the experimental studies. A comparison showed that the MCSS calculations had correctly identified the observed binding site for all three monomer units in one ligand and for two out of three in the other ligand. The correct central monomer position in the second ligand was reproduced in calculations in which the several key residues lining the pocket were allowed to move. This study validates the computational methodology. It also illustrates that subtle changes in protein structure can lead to differences in docking results and points to the importance of including target flexibility, as well as ligand flexibility, in the design process.

Conformational preferences and vibrational frequency distributions of short peptides in relation to multidimensional infrared spectroscopy

Molecular dynamics simulations of the structural distributions and the associated amide-I vibrational modes are carried out for dialanine peptide in water and carbon tetrachloride. The various manifestations in nonlinear-infrared spectroscopic experiments of the distributions of conformations of solvated dialanine are examined. The two-dimensional infrared (2D-IR) spectrum of dialanine exhibits the coupling between the amide oscillators and the correlations of the frequency fluctuations. An internally hydrogen-bonded conformation exists in CCl(4) but not in H(2)O where two externally hydrogen-bonded forms are preferred. Simulations of solvated dialanine show how the 2D-IR spectra expose the underlying structural distributions and dynamics that are not deducible from linear-infrared spectra. In H(2)O the 2D-IR shows cross-peaks from large coupling in the alpha-helical conformer and an elongated higher frequency diagonal peak, reflecting the broader distribution of structures for the more flexible acetyl end. In CCl(4), the computed cross-peak portion of the 2D-IR shows evidence of two amide-I transitions in the high-frequency region which are not apparent from the diagonal peak profile. The vibrational frequency inhomogeneity of the amide-I band arises from fluctuations of the instantaneous normal modes of these conformers rather than the shifts induced by hydrogen bonding. The simulation shows that there are correlations between fluctuations of the acetyl and amino end frequencies in H(2)O that arise from mechanical coupling and not from hydrogen bonding at the two ends of the molecule. The angular relationships between the two amide units which also show up in 2D-IR were computed, and spectral manifestations of them are discussed. The simulations also permit a calculation of the rate of energy transfer from one side of the molecule to the other. From these calculations, 2D-IR spectroscopy in conjunction with simulations is seen to be a promising tool for determining dynamics of structure changes in dipeptides.

Characterization of the Saccharomyces cerevisiae homolog of the melatonin rhythm enzyme arylalkylamine N-acetyltransferase (EC 2.3.1.87)

Arylalkylamine N-acetyltransferase (AANAT, serotonin N-acetyltransferase, EC ) plays a unique transduction role in vertebrate physiology by converting information about day and night into a hormonal signal: melatonin. Only vertebrate members of the AANAT family have been functionally characterized. Here a putative AANAT from Saccharomyces cerevisiae (scAANAT) was studied to determine whether it possessed the catalytic activity of the vertebrate enzyme. scAANAT is 47% similar to ovine AANAT, but lacks the regulatory N- and C-terminal flanking regions conserved in all vertebrate AANATs. It was found to have enzyme activity generally typical for AANAT family members, although the substrate preference pattern was somewhat broader, the specific activity was lower, and the pH optimum was higher. Deletion of scAANAT reduced arylalkylamine acetylation by S. cerevisiae extracts, indicating that scAANAT contributes significantly to this process. The scAANAT sequence conformed to the three-dimensional structure of ovine AANAT catalytic core; however, an important structural element (loop 1) was found to be shorter and to lack a proline involved in substrate binding. These differences could explain the lower specific activity of scAANAT, because of the importance of loop 1 in catalysis. Data base analysis revealed the presence of putative AANATs in other fungi but not in the nearly complete genomes of Drosophila melanogaster or Caenorhabditis elegans. These studies indicate that the catalytic and kinetic characteristics of fungal and vertebrate enzymes can be considered to be generally similar, although some differences exist that appear to be linked to changes in one structural element. Perhaps the most striking difference is that fungal AANATs lack the regulatory domains of the vertebrate enzyme, which appear to be essential for the regulatory role the enzyme plays in photochemical transduction.

Refined crystallographic structure of Pseudomonas aeruginosa exotoxin A and its implications for the molecular mechanism of toxicity

Exotoxin A of Pseudomonas aeruginosa asserts its cellular toxicity through ADP-ribosylation of translation elongation factor 2, predicated on binding to specific cell surface receptors and intracellular trafficking via a complex pathway that ultimately results in translocation of an enzymatic activity into the cytoplasm. In early work, the crystallographic structure of exotoxin A was determined to 3.0 A resolution, revealing a tertiary fold having three distinct structural domains; subsequent work has shown that the domains are individually responsible for the receptor binding (domain I), transmembrane targeting (domain II), and ADP-ribosyl transferase (domain III) activities, respectively. Here, we report the structures of wild-type and W281A mutant toxin proteins at pH 8.0, refined with data to 1.62 A and 1.45 A resolution, respectively. The refined models clarify several ionic interactions within structural domains I and II that may modulate an obligatory conformational change that is induced by low pH. Proteolytic cleavage by furin is also obligatory for toxicity; the W281A mutant protein is substantially more susceptible to cleavage than the wild-type toxin. The tertiary structures of the furin cleavage sites of the wild-type and W281 mutant toxins are similar; however, the mutant toxin has significantly higher B-factors around the cleavage site, suggesting that the greater susceptibility to furin cleavage is due to increased local disorder/flexibility at the site, rather than to differences in static tertiary structure. Comparison of the refined structures of full-length toxin, which lacks ADP-ribosyl transferase activity, to that of the enzymatic domain alone reveals a salt bridge between Arg467 of the catalytic domain and Glu348 of domain II that restrains the substrate binding cleft in a conformation that precludes NAD+ binding. The refined structures of exotoxin A provide precise models for the design and interpretation of further studies of the mechanism of intoxication.

Mutations in PTPN11, encoding the protein tyrosine phosphatase SHP-2, cause Noonan syndrome

Noonan syndrome (MIM 163950) is an autosomal dominant disorder characterized by dysmorphic facial features, proportionate short stature and heart disease (most commonly pulmonic stenosis and hypertrophic cardiomyopathy). Webbed neck, chest deformity, cryptorchidism, mental retardation and bleeding diatheses also are frequently associated with this disease. This syndrome is relatively common, with an estimated incidence of 1 in 1,000-2,500 live births. It has been mapped to a 5-cM region (NS1) [corrected] on chromosome 12q24.1, and genetic heterogeneity has also been documented. Here we show that missense mutations in PTPN11 (MIM 176876)-a gene encoding the nonreceptor protein tyrosine phosphatase SHP-2, which contains two Src homology 2 (SH2) domains-cause Noonan syndrome and account for more than 50% of the cases that we examined. All PTPN11 missense mutations cluster in interacting portions of the amino N-SH2 domain and the phosphotyrosine phosphatase domains, which are involved in switching the protein between its inactive and active conformations. An energetics-based structural analysis of two N-SH2 mutants indicates that in these mutants there may be a significant shift of the equilibrium favoring the active conformation. This implies that they are gain-of-function changes and that the pathogenesis of Noonan syndrome arises from excessive SHP-2 activity.

Helix-coil transition of PrP106-126: molecular dynamic study

A set of 34 molecular dynamic (MD) simulations totaling 305 ns of simulation time of the prion protein-derived peptide PrP106-126 was performed with both explicit and implicit solvent models. The objective of these simulations is to investigate the relative stability of the alpha-helical conformation of the peptide and the mechanism for conversion from the helix to a random-coil structure. At neutral pH, the wild-type peptide was found to lose its initial helical structure very fast, within a few nanoseconds (ns) from the beginning of the simulations. The helix breaks up in the middle and then unwinds to the termini. The spontaneous transition into the random coil structure is governed by the hydrophobic interaction between His(111) and Val(122). The A117V mutation, which is linked to GSS disease, was found to destabilize the helix conformation of the peptide significantly, leading to a complete loss of helicity approximately 1 ns faster than in the wild-type. Furthermore, the A117V mutant exhibits a different mechanism for helix-coil conversion, wherein the helix begins to break up at the C-terminus and then gradually to unwind towards the N-terminus. In most simulations, the mutation was found to speed up the conversion through an additional hydrophobic interaction between Met(112) and the mutated residue Val(117), an interaction that did not exist in the wild-type peptide. Finally, the beta-sheet conformation of the wild-type peptide was found to be less stable at acidic pH due to a destabilization of the His(111)-Val(122), since at acidic pH this histidine is protonated and is unlikely to participate in hydrophobic interaction.

Molecular mechanics and dynamics of biomolecules using a solvent continuum model

An easy implementation of molecular mechanics and molecular dynamics simulation using a continuum solvent model is presented that is particularly suitable for biomolecular simulations. The computation of solvation forces is made using the linear Poisson-Boltzmann equation (polar contribution) and the solvent-accessible surface area approach (nonpolar contribution). The feasibility of the methodology is demonstrated on a small protein and a small DNA hairpin. Although the parameters employed in this model must be refined to gain reliability, the performance of the method, with a standard choice of parameters, is comparable with results obtained by explicit water simulations. Copyright 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1830-1842, 2001

Effective atom volumes for implicit solvent models: comparison between Voronoi volumes and minimum fluctuation volumes

An essential element of implicit solvent models, such as the generalized Born method, is a knowledge of the volume associated with the individual atoms of the solute. Two approaches for determining atomic volumes for the generalized Born model are described; one is based on Voronoi polyhedra and the other, on minimizing the fluctuations in the overall volume of the solute. Volumes to be used with various parameter sets for protein and nucleic acids in the CHARMM force field are determined from a large set of known structures. The volumes resulting from the two different approaches are compared with respect to various parameters, including the size and solvent accessibility of the structures from which they are determined. The question of whether to include hydrogens in the atomic representation of the solute volume is examined. Copyright 2001 John Wiley & Sons, Inc. J Comput Chem 22: 1857-1879, 2001

Theoretical study of the conformation of the lipoamide arm in a mutant H protein

The lipoamide arm of the H protein plays a pivotal role in the catalytic cycle of the glycine decarboxylase complex (GDC) by being successively methylamine loaded (Hmet), reduced (Hred), and oxidized (Hox). In a previous study, we calculated free-energy surfaces as a function of the lipoamide arm position of the three forms of the wild-type protein and found close agreement with the available experimental data. Our simulations, together with crystallographic and NMR data, showed that the methylamine-loaded arm is locked in a cavity by interaction with Ser12, Glu14, and Asp67. In this work, we investigate the behavior of the methylamine-loaded form of a mutant H protein (HEA) where Glu14 has been replaced by Ala. We find that the arm can still be held in the cavity but that the energy barrier to release of the arm is halved from approximately 40 kcal mol(-1) for Hmet to approximately 12 kcal mol(-1) for HEA. To compensate for the loss of Glu14, the methylamine group shifts toward Ser66 in the mutant form. These results provide a structural basis for the equilibrium between the loaded and the unloaded forms of the arm observed by Gueguen et al. (Gueguen et al., J Biol Chem 1999;274:26344-26352) in HEA.

Contryphan-Vn: a novel peptide from the venom of the Mediterranean snail Conus ventricosus

The isolation, purification, and biochemical characterization of the novel peptide Contryphan-Vn, extracted from the venom of the Mediterranean marine snail Conus ventricosus, is reported. Contryphan-Vn is the first Conus peptide described from a vermivorous species and the first purified from the venom of the single Mediterranean Conus species. The amino acid sequence of Contryphan-Vn is As with other contryphans, Contryphan-Vn contains a d-tryptophan residue, is amidated at the C-terminus, and maintains the five-residue intercystine loop size. However, Contryphan-Vn differs from the known contryphans by the insertion of the Asp residue at position 2, by the lack of hydroxylation of Pro(4), and, remarkably, by the presence of the basic residue Lys(6) within the intercystine loop. Although the biological function(s) of contryphans is still unknown, these characteristics suggest distinct molecular target(s) and/or function(s) for Contryphan-Vn.

Dielectric relaxation in an enzyme active site: molecular dynamics simulations interpreted with a macroscopic continuum model

Dielectric relaxation plays an important role in many chemical processes in proteins, including acid-base titration, ligand binding, and charge transfer reactions. Its complexity makes experimental characterization difficult, and so, theoretical approaches are valuable. The comparison of molecular dynamics free energy simulations with simpler models such as a dielectric continuum model is especially useful for obtaining qualitative insights. We have analyzed a charge insertion process that models deprotonation or mutation of an important side chain in the active site of the enzyme aspartyl-tRNA synthetase. Complexes with the substrate aspartate and the analogue asparagine were studied. The resulting dielectric relaxation was found to involve both ligand and side chain rearrangements in the active site and to account for a large part of the overall charging free energy. With the continuum model, charge insertion is performed along a two-step pathway: insertion into a static environment, followed by relaxation of the environment. These correspond to different physical processes and require different protein dielectric constants. A low value of approximately 1 is needed for the static step, consistent with the parametrization of the molecular mechanics charge set used. A value of 3-6 (depending on the exact insertion site and the nature of the ligand) is needed to describe the dielectric relaxation step. This moderate value indicates that, for this system, the local protein polarizability in the active site is within at most a factor of 2 of that expected at nonspecific positions in a protein interior.

Electrostatic contributions to protein-protein interactions: fast energetic filters for docking and their physical basis

The methods of continuum electrostatics are used to calculate the binding free energies of a set of protein-protein complexes including experimentally determined structures as well as other orientations generated by a fast docking algorithm. In the native structures, charged groups that are deeply buried were often found to favor complex formation (relative to isosteric nonpolar groups), whereas in nonnative complexes generated by a geometric docking algorithm, they were equally likely to be stabilizing as destabilizing. These observations were used to design a new filter for screening docked conformations that was applied, in conjunction with a number of geometric filters that assess shape complementarity, to 15 antibody-antigen complexes and 14 enzyme-inhibitor complexes. For the bound docking problem, which is the major focus of this paper, native and near-native solutions were ranked first or second in all but two enzyme-inhibitor complexes. Less success was encountered for antibody-antigen complexes, but in all cases studied, the more complete free energy evaluation was able to identify native and near-native structures. A filter based on the enrichment of tyrosines and tryptophans in antibody binding sites was applied to the antibody-antigen complexes and resulted in a native and near-native solution being ranked first and second in all cases. A clear improvement over previously reported results was obtained for the unbound antibody-antigen examples as well. The algorithm and various filters used in this work are quite efficient and are able to reduce the number of plausible docking orientations to a size small enough so that a final more complete free energy evaluation on the reduced set becomes computationally feasible.

The mechanism of glycerol conduction in aquaglyceroporins

BACKGROUND

The E. coli glycerol facilitator, GlpF, selectively conducts glycerol and water, excluding ions and charged solutes. The detailed mechanism of the glycerol conduction and its relationship to the characteristic secondary structure of aquaporins and to the NPA motifs in the center of the channel are unknown.

RESULTS

Molecular dynamics simulations of GlpF reveal spontaneous glycerol and water conduction driven, on a nanosecond timescale, by thermal fluctuations. The bidirectional conduction, guided and facilitated by the secondary structure, is characterized by breakage and formation of hydrogen bonds for which water and glycerol compete. The conduction involves only very minor changes in the protein structure, and cooperativity between the GlpF monomers is not evident. The two conserved NPA motifs are strictly linked together by several stable hydrogen bonds and their asparagine side chains form hydrogen bonds with the substrates passing the channel in single file.

CONCLUSIONS

A complete conduction of glycerol through the GlpF was deduced from molecular dynamics simulations, and key residues facilitating the conduction were identified. The nonhelical parts of the two half-membrane-spanning segments expose carbonyl groups towards the channel interior, establishing a curve-linear pathway. The conformational stability of the NPA motifs is important in the conduction and critical for selectivity. Water and glycerol compete in a random manner for hydrogen bonding sites in the protein, and their translocations in single file are correlated. The suggested conduction mechanism should apply to the whole family.