Search

PDZ domains are found in many signaling proteins. One of their functions is to provide scaffolds for forming membrane-associated protein complexes by binding to the carboxyl termini of their partners. PDZ domains are thought also to play a signal transduction role by propagating the information that binding has occurred to remote sites. In this study, a molecular dynamics (MD) simulation-based approach, referred to as an interaction correlation analysis, is applied to the PDZ2 domain to identify the possible signal transduction pathways. A residue correlation matrix is constructed from the interaction energy correlations between all residue pairs obtained from the MD simulations. Two continuous interaction pathways, starting at the ligand binding pocket, are identified by a hierarchical clustering analysis of the residue correlation matrix. One pathway is mainly localized at the N-terminal side of helix alpha1 and the adjacent C-terminus of loop beta1-beta2. The other pathway is perpendicular to the central beta-sheet and extends toward the side of PDZ2 domain opposite to the ligand binding pocket. The results complement previous studies based on multiple sequence analysis, NMR, and MD simulations. Importantly, they reveal the energetic origin of the long-range coupling. The PDZ2 results, as well as the earlier rhodopsin analysis, show that the interaction correlation analysis is a robust approach for determining pathways of intramolecular signal transduction., ((c) 2008 Wiley-Liss, Inc.)

A calculation of the binding free energy for the dimerization of insulin has been performed using the molecular mechanics-generalized Born surface area approach. The calculated absolute binding free energy is -11.9 kcal/mol, in approximate agreement with the experimental value of -7.2 kcal/mol. The results show that the dimerization is mainly due to nonpolar interactions. The role of the hydrogen bonds between the 2 monomers appears to give the direction of the interactions. A per-atom decomposition of the binding free energy has been performed to identify the residues contributing most to the self association free energy. Residues B24-B26 are found to make the largest favorable contributions to the dimerization. Other residues situated at the interface between the 2 monomers were found to make favorable but smaller contributions to the dimerization: Tyr B16, Val B12, and Pro B28, and to an even lesser extent, Gly B23. The energy decomposition on a per-residue basis is in agreement with experimental alanine scanning data. The results obtained from a single trajectory (i.e., the dimer trajectory is also used for the monomer analysis) and 2 trajectories (i.e., separate trajectories are used for the monomer and dimer) are similar., ((c) 2005 Wiley-Liss, Inc.)

Possible insulin binding sites for D-glucose have been investigated theoretically by docking and molecular dynamics (MD) simulations. Two different docking programs for small molecules were used; Multiple Copy Simultaneous Search (MCSS) and Solvation Energy for Exhaustive Docking (SEED) programs. The configurations resulting from the MCSS search were evaluated with a scoring function developed to estimate the binding free energy. SEED calculations were performed using various values for the dielectric constant of the solute. It is found that scores emphasizing non-polar interactions gave a preferential binding site in agreement with that inferred from recent fluorescence and NMR NOESY experiments. The calculated binding affinity of -1.4 to -3.5 kcal/mol is within the measured range of -2.0 +/- 0.5 kcal/mol. The validity of the binding site is suggested by the dynamical stability of the bound glucose when examined with MD simulations with explicit solvent. Alternative binding sites were found in the simulations and their relative stabilities were estimated. The motions of the bound glucose during molecular dynamics simulations are correlated with the motions of the insulin side chains that are in contact with it and with larger scale insulin motions. These results raise the question of whether glucose binding to insulin could play a role in its activity. The results establish the complementarity of molecular dynamics simulations and normal mode analyses with the search for binding sites proposed with small molecule docking programs., (Copyright 2004 Wiley-Liss, Inc.)

The effect of the mutation Arg 96 to His on the stability of bacteriophage T4 lysozyme has been previously studied by calorimetric experiments, X-ray crystallography, and free energy simulation techniques. The experimental and calculated values for the difference between the free energy of denaturation of the mutant and the wild type are in reasonable agreement. However, the two approaches led to different explanations for the loss in stability. To analyze the differences, a series of refinements based on the crystallographic data were performed, a number of aspects of the simulations were reexamined, and continuum electrostatic calculations were done to complement the latter. The results of those comparisons provide a better understanding of the origin of the free energy difference in this mutant. Furthermore, they show the importance of the combined use of simulations and crystallography for interpreting the effects of mutations on the energetics of the system., (Copyright 2004 Wiley-Liss, Inc.)

Although most side-chain torsion angles correspond to low-energy rotameric positions, deviations occur with significant frequency. One striking example arises in Trp residues, which have an important role in stabilizing protein structures because of their size and mixed hydrophobic/hydrophilic character. Ten percent of Trp side-chains have unexplained conformations with chi(2) near 0 degrees instead of the expected 90 degrees. The current work is a structural and energetic analysis of these conformations. It is shown that many Trp residues with these orientations are stabilized by three-center carbon-donor hydrogen bonds of the form C-H...X...H-C, where X is a polar hydrogen-bond acceptor in the environment of the side-chain. The bridging hydrogen bonds occur both within the Trp side-chain and between the side-chain and the local protein backbone. Free energy maps of an isolated Trp residue in an explicit water environment show a minimum corresponding to the off-rotamer peak observed in the crystallographic data. Bridging carbon-donor hydrogen bonds are also shown to stabilize on-rotamer Trp conformations, and similar bridging hydrogen bonds also stabilize some off-rotamer Asp conformations. The present results suggest a previously unrecognized role for three-center carbon-donor hydrogen bonds in protein structures and support the view that the off-rotamer Trp side-chain orientations are real rather than artifacts of crystallographic refinements. Certain of the off-rotamer Trp conformations appear to have a functional role., (Copyright 2004 Wiley-Liss, Inc.)

Delineation of the structural properties of transition states is key to deriving models for protein folding. Here we describe the structures of the transition states of the bacterial immunity proteins Im7 and Im9 obtained by all-atom molecular dynamics simulations with phi value restraints derived from protein engineering experiments. This pair of proteins is of special interest because, at pH 7 and 10 degrees C, Im7 folds via an intermediate while Im9 folds with a two-state transition. The structures of the transition states for Im7 and Im9, together with their radii of gyration and distances from the native state, are similar. The typical distance between any two members of the transition state ensemble of both proteins is large, with that of Im9 nearly twice that of Im7. Thus, a broad range of structures make up the transition state ensembles of these proteins. The ensembles satisfy the set of rather low phi values and yet are consistent with high beta(T) values (> 0.85 for both proteins). For both Im7 and Im9 the inter-helical angles are highly variable in the transition state ensembles, although the native contacts between helices I and IV are well conserved. By measuring the distribution of the accessible surface area for each residue we show that the hydrophobic residues that are buried in the native state remain buried in the transition state, corresponding to a hydrophobic collapse to a relatively ordered globule. The data provide new insights into the structural properties of the transition states of these proteins at an atomic level of detail and show that molecular dynamics simulations with phi value restraints can significantly enhance the knowledge of the transition state ensembles (TSE) provided by the experimental phi values alone., (Copyright 2003 Wiley-Liss, Inc.)

Gō-type models, which include only native contact interactions in the energy function, are being used increasingly to describe the protein folding reaction. To investigate the validity of such models, we determine the role of native and non-native interactions along folding and unfolding pathways. For this purpose, we use a molecular mechanics energy function with an implicit solvation model (an effective energy function or potential of mean force) that can be expressed in a pairwise decomposable form. We find that for the native state and a wide range of other configurations, the contact energy is an accurate description, in part due to the cancellation of non-zero contributions from more distant residues. However, significant errors in the energy are introduced for non-native structures if the energy is calculated from the native contacts alone. Non-native contacts tend to make a significant contribution, particularly for molten globules and collapsed states along the unfolding pathways. The implication of these results for the use of Gō-type models in studies of protein folding are discussed., (Copyright 2002 Wiley-Liss, Inc.)

The R-factor and R-free are commonly used to measure the quality of protein models obtained in X-ray crystallography. Well-refined protein structures usually have R-factors in the range of 20-25%, whereas intrinsic errors in the experimental data are usually around 5%. We use molecular dynamics simulations to perform a self-consistent analysis by which we determine the major factors contributing to large values of protein R-factors. The analysis shows that significant R-factor values can arise from the use of isotropic B-factors to model anisotropic protein motions and from coordinate errors. Even in the absence of coordinate errors, the use of isotropic B-factors can cause the R-factors to be around 10%; for coordinate errors smaller than 0.2 A, the two errors types make similar contributions. The inaccuracy of the energy function used and multistate protein dynamics are unlikely to make significant contributions to the large R-factors., (Copyright 2002 Wiley-Liss, Inc.)

The structures of two species of potato carboxypeptidase inhibitor with nonnative disulfide bonds were determined by molecular dynamics simulations in explicit solvent using disulfide bond constraints that have been shown to work for the native species. Ten structures were determined; five for scrambled A (disulfide bonds between Cys8-Cys27, Cys12-Cys18, and Cys24-Cys34) and five for the scrambled C (disulfide bonds Cys8-Cys24, Cys12-Cys18, and Cys27-Cys34). The two scrambled species were both more solvent exposed than the native structure; the scrambled C species was more solvent exposed and less compact than the scrambled A species. Analysis of the loop regions indicates that certain loops in scrambled C are more nativelike than in scrambled A. These factors, combined with the fact that scrambled C has one native disulfide bond, may contribute to the observed faster conversion to the native structure from scrambled C than from scrambled A. Results from the PROCHECK program using the standard parameter database and a database specially constructed for small, disulfide-rich proteins indicate that the 10 scrambled structures have correct stereochemistry. Further, the results show that a characteristic feature of small, disulfide-rich proteins is that they score poorly using the standard PROCHECK parameter database. Proteins 2000;40:482-493., (Copyright 2000 Wiley-Liss, Inc.)

The dynamic ligand design (DLD) algorithm, an automated method for the creation of novel ligands, links up small functional groups that have been placed in energetically favorable positions in the binding site of a target molecule. The positions and orientations of the small functional groups can be determined using the multi-copy simultaneous search approach (MCSS) or experimental data. In this work the original DLD methodology is extended by using a modified version of the pseudo-potential energy function. A novel simulated annealing protocol is presented for optimizing the pseudo-potential energy of ligands in the binding site; the protocol is expected to be applicable to other optimization problems. The utility of the method is illustrated by designing an inhibitor for endothiapepsin. The binding affinity of the inhibitor is assessed using a thermodynamic cycle that decomposes the binding free energy into a sum of translational, rotational, configurational, hydrophobic, and electrostatic contributions. The calculations suggest that the designed molecule will bind endothiapepsin with high affinity., (Copyright 2000 Wiley-Liss, Inc.)

The Multiple Copy Simultaneous Search (MCSS) methodology for finding energetically favorable positions and orientations of small functional groups in a binding site is extended to include flexibility of the target. This makes possible the finding of novel minima not present in a fixed structure and so extends the diversity of inhibitors that can be constructed starting with the MCSS procedure. Quenched molecular dynamics is used to generate energetically favorable positions and orientations of the functional groups in the field of a flexible protein. The method is applied to the viral protein HIV-1 protease with methanol and methyl ammonium as a test case. If the protein is quenched with many copies of functional groups randomly distributed in the binding site, the resulting minima have ligand-protein interaction energies that are, on average, less favorable than those obtained with standard MCSS. This is a consequence of the renormalized potential function employed in the Locally Enhanced Sampling (LES) approximation. However, local optimizations of existing MCSS minima with a flexible protein results in lower energy minima in regions of the protein that are of particular interest. Their use in constructing a consensus protein model for ligand design is discussed.

A primary and a secondary neural network are applied to secondary structure and structural class prediction for a database of 681 non-homologous protein chains. A new method of decoding the outputs of the secondary structure prediction network is used to produce an estimate of the probability of finding each type of secondary structure at every position in the sequence. In addition to providing a reliable estimate of the accuracy of the predictions, this method gives a more accurate Q3 (74.6%) than the cutoff method which is commonly used. Use of these predictions in jury methods improves the Q3 to 74.8%, the best available at present. On a database of 126 proteins commonly used for comparison of prediction methods, the jury predictions are 76.6% accurate. An estimate of the overall Q3 for a given sequence is made by averaging the estimated accuracy of the prediction over all residues in the sequence. As an example, the analysis is applied to the target beta-cryptogein, which was a difficult target for ab initio predictions in the CASP2 study; it shows that the prediction made with the present method (62% of residues correct) is close to the expected accuracy (66%) for this protein. The larger database and use of a new network training protocol also improve structural class prediction accuracy to 86%, relative to 80% obtained previously. Secondary structure content is predicted with accuracy comparable to that obtained with spectroscopic methods, such as vibrational or electronic circular dichroism and Fourier transform infrared spectroscopy.

A Gaussian solvent-exclusion model for the solvation free energy is developed. It is based on theoretical considerations and parametrized with experimental data. When combined with the CHARMM 19 polar hydrogen energy function, it provides an effective energy function (EEF1) for proteins in solution. The solvation model assumes that the solvation free energy of a protein molecule is a sum of group contributions, which are determined from values for small model compounds. For charged groups, the self-energy contribution is accounted for primarily by the exclusion model. Ionic side-chains are neutralized, and a distance-dependent dielectric constant is used to approximate the charge-charge interactions in solution. The resulting EEF1 is subjected to a number of tests. Molecular dynamics simulations at room temperature of several proteins in their native conformation are performed, and stable trajectories are obtained. The deviations from the experimental structures are similar to those observed in explicit water simulations. The calculated enthalpy of unfolding of a polyalanine helix is found to be in good agreement with experimental data. Results reported elsewhere show that EEF1 clearly distinguishes correctly from incorrectly folded proteins, both in static energy evaluations and in molecular dynamics simulations and that unfolding pathways obtained by high-temperature molecular dynamics simulations agree with those obtained by explicit water simulations. Thus, this energy function appears to provide a realistic first approximation to the effective energy hypersurface of proteins.

The folding ability of a heteropolymer model for proteins subject to Monte Carlo dynamics on a simple cubic lattice is shown to be strongly correlated with the stability of the native state. We consider a number of estimates of the stability that can be determined without simulation, including the energy gap between the native state and the structurally dissimilar part of the spectrum (Z score) and, for sequences with fully compact native states, the gap in energy between the native and first excited fully compact states. These estimates are found to be more robust predictors of folding ability than a parameter sigma that requires simulation for its evaluation: sigma = 1 - Tf/Ttheta, where Tf is the temperature at which the fluctuation of an order parameter is at its maximum and Ttheta is the temperature at which the specific heat is at its maximum. We show that the interpretation of Ttheta as the collapse transition temperature is not correct in general and that the correlation between sigma and the folding ability arises from the fact that sigma is related to the energy gap (Z score).

We investigate the folding of a 125-bead heteropolymer model for proteins subject to Monte Carlo dynamics on a simple cubic lattice. Detailed study of a few sequences revealed a folding mechanism consisting of a rapid collapse followed by a slow search for a stable core that served as the transition state for folding to a near-native intermediate. Rearrangement from the intermediate to the native state slowed folding further because it required breaking native-like local structure between surface monomers so that those residues could condense onto the core. We demonstrate here the generality of this mechanism by a statistical analysis of a 200 sequence database using a method that employs a genetic algorithm to pick the sequence attributes that are most important for folding and an artificial neural network to derive the corresponding functional dependence of folding ability on the chosen sequence attributes [quantitative structure-property relationships (QSPRs)]. QSPRs that use three sequence attributes yielded substantially more accurate predictions than those that use only one. The results suggest that efficient search for the core is dependent on both the native state's overall stability and its amount of kinetically accessible, cooperative structure, whereas rearrangement from the intermediate is facilitated by destabilization of contacts between surface monomers. Implications for folding and design are discussed.

Several methods for including the conformational flexibility of proteins in the calculation of titration curves are compared. The methods use the linearized Poisson-Boltzmann equation to calculate the electrostatic free energies of solvation and are applied to bovine pancreatic trypsin inhibitor (BPTI) and hen egg-white lysozyme (HEWL). An ensemble of conformations is generated by a molecular dynamics simulation of the proteins with explicit solvent. The average titration curve of the ensemble is calculated in three different ways: an average structure is used for the pKa calculation; the electrostatic interaction free energies are averaged and used for the pKa calculation; and the titration curve for each structure is calculated and the curves are averaged. The three averaging methods give very similar results and improve the pKa values to approximately the same degree. This suggests, in contrast to implications from other work, that the observed improvement of pKa values in the present studies is due not to averaging over an ensemble of structures, but rather to the generation of a single properly averaged structure for the pKa calculation.

A combinatorial ligand design approach based on the multiple copy simultaneous search (MCSS) method and a simple scheme for joining MCSS functional group sites was applied to the binding pocket of P3/Sabin poliovirus and rhinovirus 14. The MCSS method determines where specific functional (chemical) groups have local potential energy minima in the binding site. Before the virus application, test calculations were run to determine the optimal set of input parameters to be used in evaluating the MCSS results. The MCSS minima are analysed and selected minima are connected with (CH2)n linkers to form candidate ligands, whose structures are optimized in the binding site. Estimates of the binding strength were made for the ligands and compared with those for known drugs. The results indicate that the proposed ligands should bind to P3/Sabin poliovirus at least as well as the best of the existing drugs, and that they should also bind to P1/Mahoney poliovirus and rhinovirus 14. A detailed comparison of the poliovirus and rhinovirus binding pockets and an analysis of drug binding specificity is presented.

An automated method for the dynamic ligand design (DLD) for a binding site of known structure is described. The method can be used for the creation of de novo ligands and for the modification of existing ligands. The binding site is saturated with atoms (sp3 carbon atoms in the present implementation) that form molecules under the influence of a potential function that joins atoms to each other with the correct stereochemistry. The resulting molecules are linked to precomputed functional group minimum energy positions in the binding site. The generalized potential function allows atoms to sample a continuous parameter space that includes the Cartesian coordinates and their occupancy and type, e.g., the method allows change of an sp3 carbon into an sp2 carbon or oxygen. A parameter space formulated in this way can then be sampled and optimized by a variety of methods. In this work, molecules are generated by use of a Monte Carlo simulated annealing algorithm. The DLD method is illustrated by its application to the binding site of FK506 binding protein (FKBP), an immunophilin. De novo ligands are designed and modification of the immunosuppressant drug FK506 are suggested. The results demonstrate that the dynamic ligand design approach can automatically construct ligands which complement both the shape and charge distribution of the binding site.

We evaluate 3D models of human nucleoside diphosphate kinase, mouse cellular retinoic acid binding protein I, and human eosinophil neurotoxin that were calculated by MODELLER, a program for comparative protein modeling by satisfaction of spatial restraints. The models have good stereochemistry and are at least as similar to the crystallographic structures as the closest template structures. The largest errors occur in the regions that were not aligned correctly or where the template structures are not similar to the correct structure. These regions correspond predominantly to exposed loops, insertions of any length, and non-conserved side chains. When a template structure with more than 40% sequence identity to the target protein is available, the model is likely to have about 90% of the mainchain atoms modeled with an rms deviation from the X-ray structure of approximately 1 A, in large part because the templates are likely to be that similar to the X-ray structure of the target. This rms deviation is comparable to the overall differences between refined NMR and X-ray crystallography structures of the same protein.

Force field parameters that use a combination of Lennard-Jones and electrostatic interactions are developed for divalent zinc and tested in solution and protein simulations. It is shown that the parameter set gives free energies of solution in good agreement with experiment. Molecular dynamics simulations of carboxypeptidase A and carbonic anhydrase are performed with these zinc parameters and the CHARMM 22 beta all-atom parameter set. The structural results are as accurate as those obtained in published simulations that use specifically bonded models for the zinc ion and the AMBER force field. The inclusion of longer-range electrostatic interactions by use of the Extended Electrostatics model is found to improve the equilibrium conformation of the active site It is concluded that the present parameter set, which permits different coordination geometries and ligand exchange for the zinc ion, can be employed effectively for both solution and protein simulations of zinc-containing systems.

A theoretical analysis is made of the decomposition into contributions from individual interactions of the free energy calculated by thermodynamic integration. It is demonstrated that such a decomposition, often referred to as "component analysis," is meaningful, even though it is a function of the integration path. Moreover, it is shown that the path dependence can be used to determine the relation of the contribution of a given interaction to the state of the system. To illustrate these conclusions, a simple transformation (Cl- to Br- in aqueous solution) is analyzed by use of the Reference Interaction Site Model-Hypernetted Chain Closure integral equation approach; it avoids the calculational difficulties of macromolecular simulation while retaining their conceptual complexity. The difference in the solvation free energy between chloride and bromide is calculated, and the contributions of the Lennard-Jones and electrostatic terms in the potential function are analyzed by the use of suitably chosen integration paths. The model is also used to examine the path dependence of individual contributions to the double free energy differences (delta delta G or delta delta A) that are often employed in free energy simulations of biological systems. The alchemical path, as contrasted with the experimental path, is shown to be appropriate for interpreting the effects of mutations on ligand binding and protein stability. The formulation is used to obtain a better understanding of the success of the Poisson-Boltzmann continuum approach for determining the solvation properties of polar and ionic systems.

A program (HOOK) is described for generating potential ligands that satisfy the chemical and steric requirements of the binding region of a macromolecule. Functional group sites with defined positions and orientations are derived from known ligand structures or the multicopy simulation search (MCSS) method (Miranker, A., Karplus, M. Proteins 11:29-34, 1991). HOOK places molecular "skeletons" from a database into the protein binding region by making bonds between sites ("hooks") on the skeleton and functional groups. The nonpolar interactions with the binding region of candidate molecules are assessed by use of a simplified van der Waals potential. The method is illustrated by constructing ligands for the sialic acid binding site of the hemagglutinin from the influenza A virus and the active site of chloramphenicol acetyltransferase. Aspects of the HOOK program that lead to a highly efficient search of 10(5) or more skeletons for binding to 10(2) or more functional group minima are outlined.

The motions of water molecules, the acyl moiety, the catalytic triad, and the oxyanion binding site of acyl-chymotrypsin were studied by means of a stochastic boundary molecular dynamics simulation. A water molecule that could provide the nucleophilic OH- for the deacylation stage of the catalysis was found to be trapped between the imidazole ring of His-57 and the carbonyl carbon of the acyl group. It makes a hydrogen bond with the N epsilon 2 of His-57 and is held in place through a network of hydrogen-bonded water molecules in the active site. The water molecule was found as close as 2.8 A to the carbonyl carbon. This appears to be due to the constraints imposed by nonbonded interaction in the active site. Configurations were found in which one hydrogen of the trapped water shared a bifurcated hydrogen bond with His-57-N epsilon 2 and Ser-195-O gamma, with the water oxygen very close to the carbonyl carbon. The existence of such a water molecule suggests that large movement of the His-57 imidazole ring between positions suitable for providing general-base catalyzed assistance and for providing general-acid catalyzed assistance may not be required during the reaction. The simulation indicates that the side chains of residues involved in catalysis (i.e., His-57, Ser-195, and Asp-102) are significantly less flexible than other side chains in the protein. The 40% reduction in rms fluctuations is consistent with a comparable reduction calculated from the temperature factors obtained in the X-ray crystallographic data of gamma-chymotrypsin. The greater rigidity of active site residues seems to result from interconnected hydrogen bonding networks among the residues and between the residues and the solvent water in the active site.

The vibrational entropy of native BPTI, with three disulfide bonds, was determined by use of normal mode calculations and compared with that of folded variants having either one less disulfide bond or lacking a peptide bond at the trypsin-reactive site. Favorable contributions to the free energy of 2.5-5.1 kcal/mol at 300 K were calculated for the reduction of disulfide bonds in the folded state, whereas no favorable contribution was found for the hydrolysis of the peptide bond cleaved by trypsin. This is on the order of the effect of disulfides in the unfolded state. The implications of these results for the stabilization of a folded protein by the introduction of cross-links are discussed.

A simplified model of a polypeptide chain is used to study the dynamics of the beta-sheet-coil transition. Each amino acid residue is treated as a single quasiparticle in an effective potential that approximates the potential of mean force in solution. The model is used to study the equilibrium and dynamic aspects of the sheet-coil transition. Systems studied include ones with both strands free to move (two-strand sheet), and ones with either strand fixed in position (multistrand sheet). The equilibrium properties examined include sheet-coil equilibrium constants and their dependence on chain position. Dynamic properties are investigated by a stochastic simulation of the Brownian motion of the chain in its solvent surroundings. Time histories of the dihedral angles and residue-residue cross-strand distances are used to study the behavior of the sheet structure. Auto- and cross-correlation functions are calculated from the time histories with relaxation times of tens to hundreds of picoseconds. Sheet-coil rate constants of tens of ns-1 were found for the fixed strand cases.

Conformational disorder in crystal structures of ribonuclease-A and crambin is studied by including two independent structures in least-squares optimizations against X-ray data. The optimizations are carried out by X-ray restrained molecular dynamics (simulated annealing refinement) and by conventional least-squares optimization. Starting from two identical structures, the optimizations against X-ray data lead to significant deviations between the two, with rms backbone displacements of 0.45 A for refinement of ribonuclease at 1.53 A resolution, and 0.31 A for crambin at 0.945 A. More than 15 independent X-ray restrained molecular dynamics runs have been carried out for ribonuclease, and the displacements between the resulting structures are highly reproducible for most atoms. These include residues with two or more conformations with significant dihedral angle differences and alternative hydrogen bonding, as well as groups of residues that undergo displacements that are suggestive of rigid-body librations. The crystallographic R-values obtained are approximately 13%, as compared to 15.3% for a comparable refinement with a single structure. Least-squares optimization without an intervening restrained molecular dynamics stage is sufficient to reproduce most of the observed displacements. Similar results are obtained for crambin, where the higher resolution of the X-ray data allows for refinement of unconstrained individual anisotropic temperature factors. These are shown to be correlated with the displacements in the two-structure refinements.

A method is described for identifying collective motions in proteins from molecular dynamics trajectories or normal mode simulations. The method makes use of the covariances of atomic positional fluctuations. It is illustrated by an analysis of the bovine pancreatic trypsin inhibitor. Comparison of the covariance and cross-correlation matrices shows that the relative motions have many similar features in the different simulations. Many regions of the protein, especially regions of secondary structure, move in a correlated manner. Anharmonic effects, which are included in the molecular dynamics simulations but not in the normal analysis, are of some importance in determining the larger scale collective motions, but not the more local fluctuations. Comparisons of molecular dynamics simulations in the present and absence of solvent indicate that the environment is of significance for the long-range motions.

A new method is proposed for determining energetically favorable positions and orientations for functional groups on the surface of proteins with known three-dimensional structure. From 1,000 to 5,000 copies of a functional group are randomly placed in the site and subjected to simultaneous energy minimization and/or quenched molecular dynamics. The resulting functionality maps of a protein receptor site, which can take account of its flexibility, can be used for the analysis of protein ligand interactions and rational drug design. Application of the method to the sialic acid binding site of the influenza coat protein, hemagglutinin, yields functional group minima that correspond with those of the ligand in a cocrystal structure.

Positional probability density functions (pdf) for the atomic fluctuations are determined from a molecular dynamics simulation for hen egg-white lysozyme. Most atoms are found to have motions that are highly anisotropic but only slightly anharmonic. The largest deviations from harmonic motion are in the direction of the largest rms fluctuations in the local principal axis frame. Backbone atoms tend to be more nearly harmonic than sidechain atoms. The atoms with the largest anharmonicities tend to have pdfs with multiple peaks, each of which is close to harmonic. Several model pdfs are evaluated on the basis of how well they fit probability densities from the dynamics simulations when parameterized in terms of the moments of the distribution. Gram-Charlier and Edgeworth perturbation expansions, which have been successful in describing the motions of small molecules in crystals, are shown to be inadequate for the distributions found in the dynamics of proteins. Multipeaked distribution functions are found to be more appropriate.

A method for the prediction of hydrogen positions in proteins is presented. The method is based on the knowledge of the heavy atom positions obtained, for instance, from X-ray crystallography. It employs an energy minimization limited to the environment of the hydrogen atoms bound to a common heavy atom or to a single water molecule. The method is not restricted to proteins and can be applied without modification to nonpolar hydrogens and to nucleic acids. The method has been applied to the neutron diffraction structures of trypsin, ribonuclease A, and bovine pancreatic trypsin inhibitor. A comparison of the constructed and the observed hydrogen positions shows few deviations except in situations in which several energetically similar conformations are possible. Analysis of the potential energy of rotation of Lys amino and Ser, Thr, Tyr hydroxyl groups reveals that the conformations of lowest intrinsic torsion energies are statistically favored in both the crystal and the constructed structures.

The uncertainties in the refined parameters for a 1.5-A X-ray structure of carbon-monoxy (FeII) myoglobin are estimated by combining energy minimization with least-squares refinement against the X-ray data. The energy minimizations, done without reference to the X-ray data, provide perturbed structures which are used to restart conventional X-ray refinement. The resulting refined structures have the same, or better, R-factor and stereochemical parameters as the original X-ray structure, but deviate from it by 0.13 A rms for the backbone atoms and 0.31 A rms for the sidechain atoms. Atoms interacting with a disordered sidechain, Arg 45 CD3, are observed to have larger positional uncertainties. The uncertainty in the B-factors, within the isotropic harmonic motion approximation, is estimated to be 15%. The resulting X-ray structures are more consistent with the energy parameters used in simulations.

The relation between protein secondary structure and internal motions was examined by using molecular dynamics to calculate positional fluctuations of individual helix, beta-sheet, and loop structural elements in free and substrate-bound hen egg-white lysozyme. The time development of the fluctuations revealed a general correspondence between structure and dynamics; the fluctuations of the helices and beta-sheets converged within the 101 psec period of the simulation and were lower than average in magnitude, while the fluctuations of the loop regions were not converged and were mostly larger than average in magnitude. Notable exceptions to this pattern occurred in the substrate-bound simulation. A loop region (residues 101-107) of the active site cleft had significantly reduced motion due to interactions with the substrate. Moreover, part of a loop and a 3(10) helix (residues of 67-88) not in contact with the substrate showed a marked increase in fluctuations. That these differences in dynamics of free and substrate-bound lysozyme did not result simply from sampling errors was established by an analysis of the variations in the fluctuations of the two halves of the 101 psec simulation of free lysozyme. Concerted transitions of four to five mainchain phi and psi angles between dihedral wells were shown to be responsible for large coordinate shifts in the loops. These transitions displaced six or fewer residues and took place either abruptly, in 1 psec or less, or with a diffusive character over 5-10 psec. Displacements of rigid secondary structures involved longer timescale motions in bound lysozyme; a 0.5 A rms change in the position of a helix occurred over the 55 psec simulation period. This helix reorientation within the protein appears to be a response to substrate binding. There was little correlation between the solvent accessible surface area and the dynamics of the different structural elements.