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

Energetics of ion conduction through the K+ channel

K+ channels are transmembrane proteins that are essential for the transmission of nerve impulses. The ability of these proteins to conduct K+ ions at levels near the limit of diffusion is traditionally described in terms of concerted mechanisms in which ion-channel attraction and ion-ion repulsion have compensating effects, as several ions are moving simultaneously in single file through the narrow pore. The efficiency of such a mechanism, however, relies on a delicate energy balance-the strong ion-channel attraction must be perfectly counterbalanced by the electrostatic ion-ion repulsion. To elucidate the mechanism of ion conduction at the atomic level, we performed molecular dynamics free energy simulations on the basis of the X-ray structure of the KcsA K+ channel. Here we find that ion conduction involves transitions between two main states, with two and three K+ ions occupying the selectivity filter, respectively; this process is reminiscent of the 'knock-on' mechanism proposed by Hodgkin and Keynes in 1955. The largest free energy barrier is on the order of 2-3 kcal mol-1, implying that the process of ion conduction is limited by diffusion. Ion-ion repulsion, although essential for rapid conduction, is shown to act only at very short distances. The calculations show also that the rapidly conducting pore is selective.

Photochemically induced electron transfer (PET) catalyzed radical cyclization: a practical method for inducing structural changes in peptides by formation of cyclic amino acid derivatives

A new radical cyclization reaction of unsaturated amino acid derivatives is presented. The reaction is induced by photoelectron transfer (PET) catalysis and proceeds, in comparison to commonly applied methods, under mild, nonoxidizing, and nontoxic conditions in neutral medium. This type of radical cyclization reaction can be used in peptide chemistry for inducing structural changes in peptides.

Uracil-DNA glycosylase acts by substrate autocatalysis

In humans, uracil appears in DNA at the rate of several hundred bases per cell each day as a result of misincorporation of deoxyuridine (dU) or deamination of cytosine. Four enzymes that catalyse the hydrolysis of the glycosylic bond of dU in DNA to yield an apyridiminic site as the first step in base excision repair have been identified in the human genome. The most efficient and well characterized of these uracil-DNA glycosylases is UDG (also known as UNG and present in almost all known organisms), which excises U from single- or double-stranded DNA and is associated with DNA replication forks. We used a hybrid quantum-mechanical/molecular-mechanical (QM/MM) approach to determine the mechanism of catalysis by UDG. In contrast to the concerted associative mechanism proposed initially, we show here that the reaction proceeds in a stepwise dissociative manner. Cleavage of the glycosylic bond yields an intermediate comprising an oxocarbenium cation and a uracilate anion. Subsequent attack by a water molecule and transfer of a proton to D145 result in the products. Surprisingly, the primary contribution to lowering the activation energy comes from the substrate, rather than from the enzyme. This 'autocatalysis' derives from the burial and positioning of four phosphate groups that stabilize the rate-determining transition state. The importance of these phosphates explains the residual activity observed for mutants that lack key residues. A corresponding catalytic mechanism could apply to the DNA glycosylases TDG and SMUG1, which belong to the same structural superfamily as UDG.

The energetics of off-rotamer protein side-chain conformations

Non-rotameric ("off-rotamer") conformations are commonly observed for the side-chains of protein crystal structures. This study examines whether such conformations are real or artifactual by comparing the energetics of on and off-rotamer side-chain conformations calculated with the CHARMM energy function. Energy-based predictions of side-chain orientation are carried out by rigid-geometry mapping in the presence of the fixed protein environment for 1709 non-polar side-chains in 24 proteins for which high-resolution (2.0 A or better) structures are available. For on-rotamer conformations, 97.6 % are correctly predicted; i.e. they correspond to the absolute minima of their local side-chain energy maps (generally to within 10 degrees or less). By contrast, for the observed off-rotamer side-chain conformations, 63.8 % are predicted correctly. This difference is statistically significant (P<0.001) and suggests that while most of the observed off-rotamer conformations are real, many of the erroneously predicted ones are likely to be artifacts of the X-ray refinements. Probabilities for off-rotamer conformations of the non-polar side-chains are calculated to be 5.0-6.1 % by adaptive umbrella-sampled molecular dynamics trajectories of individual amino acid residues in vacuum and in the presence of an average protein or aqueous dielectric environment. These results correspond closely to the 5.7 % off-rotamer fraction predicted by the rigid-geometry mapping studies. Since these values are about one-half of the 10.2 % off-rotamer fraction observed in the X-ray structures, they support the conclusion that many of the latter are artifacts. In both the rigid-geometry mapping and the molecular dynamics studies, the discrepancies between the predicted and observed fractions of off-rotamer conformations are largest for leucine residues (approximately 6 % versus 16.6 %). The simulations for the isolated amino acid residues indicate that the real off-rotamer frequency of 5-6 % is consistent with the internal side-chain and local side-chain-backbone energetics and does not originate from shifts due to the protein. The present results suggest that energy-based rotation maps can be used to find side-chain positional artifacts that appear in crystal structures based on refinements in the 2 A resolution range.

Simulated refolding of stretched titin immunoglobulin domains

Steered molecular dynamics (SMD) is used to investigate forced unfolding and spontaneous refolding of immunoglobulin I27, a domain of the muscle protein titin. Previous SMD simulations revealed the events leading to stretch-induced unfolding of I27, the rupture of hydrogen bonds bridging beta-strands A and B, and those bridging beta-strands A' and G, the latter rupture occurring at an extension of approximately 15 A and preceding the complete unfolding. Simulations are now used to study the refolding of partially unfolded I27 domains. The results reveal that stretched domains with ruptured interstrand hydrogen bonds shrink along the extension direction. Two types of refolding patterns are recognized: for separated beta-strands A' and G, in most simulations five of the six hydrogen bonds between A' and G stably reformed in 2 ns, whereas for separated beta-strands A and B hydrogen bonds seldom reformed in eight 2-ns simulations. The mechanical stability of the partially refolded intermediates has been tested by re-stretching.

Structure and dynamics of calcium-activated calmodulin in solution

Two 4-ns molecular dynamics simulations of calcium loaded calmodulin in solution have been performed, using both standard nonbonded cutoffs and Ewald summation to treat electrostatic interactions. Our simulation results are generally consistent with solution experimental studies of calmodulin structure and dynamics, including NMR, cross-linking, fluorescence and x-ray scattering. The most interesting result of the molecular dynamics simulations is the detection of large-scale structural fluctuations of calmodulin in solution. The globular N- and C-terminal domains tend to move approximately like rigid bodies, with fluctuations of interdomain distances within a 7 A range and of interdomain angles by up to 60 deg. Essential dynamics analysis indicates that the three dominant types of motion involve bending of the central helix in two perpendicular planes and a twist in which the domains rotate in opposite directions around the central helix. In the more realistic Ewald trajectory the protein backbone remains mostly within a 2-3 A root-mean-square distance from the crystal structure, the secondary structure within the domains is conserved and middle part of the central helix becomes disordered. The central helix itself exhibits limited fluctuations, with its bend angle exploring the 0-50 degrees range and the end-to-end distance falling in 39-43 A. The results of the two simulations were similar in many respects. However, the cutoff trajectory exhibited a larger deviation from the crystal, loss of several helical hydrogen bonds in the N-terminal domain and lack of structural disorder in the central helix.

Ligand diffusion in the catalase from Proteus mirabilis: a molecular dynamics study

The role of the channels and cavities present in the catalase from Proteus mirabilis (PMC) was investigated using molecular dynamics (MD) simulations. The reactant and products of the reaction, H(2)O(2) -->1/2 O(2) + H(2)O, catalyzed by the enzyme were allowed to diffuse to and from the active site. Dynamic fluctuations in the structure are found necessary for the opening of the major channel, identified in the X-ray model, which allows access to the active site. This channel is the only pathway to the active site observed during the dynamics, and both the products and reactant use it. H(2)O and O(2) are also detected in a cavity defined by the heme and Ser196, which could play an important role during the reaction. Free energy profiles of the ligands diffusing through the major channel indicate that the barriers to ligand diffusion are less than 20 kJ mol(-1) for each of the species. It is not clear from our study that minor channels play a role for access to the protein active site or to the protein surface.

Functional group placement in protein binding sites: a comparison of GRID and MCSS

One approach to combinatorial ligand design begins by determining optimal locations (i.e., local potential energy minima) for functional groups in the binding site of a target macromolecule. MCSS and GRID are two methods, based on significantly different algorithms, which are used for this purpose. A comparison of the two methods for the same functional groups is reported. Calculations were performed for nonpolar and polar functional groups in the internal hydrophobic pocket of the poliovirus capsid protein, and on the binding surface of the src SH3 domain. The two approaches are shown to agree qualitatively; i.e., the global characteristics of the functional group maps generated by MCSS and GRID are similar. However, there are significant differences in the relative interaction energies of the two sets of minima, a consequence of the different functional form used to evaluate polar interactions (electrostatics and hydrogen bonding) in the two methods. The single sphere representation used by GRID affords only positional information, supplemented by the identification of hydrogen bonding interactions. By contrast, the multi-atom representation of most MCSS groups yields in both positional and orientational information. The two methods are most similar for small functional groups, while for larger functional groups MCSS yields results consistent with GRID but superior in detail. These results are in accord with the somewhat different purposes for which the two methods were developed. GRID has been used mainly to introduce functionalities at specific positions in lead compounds, in which case the orientation is predetermined by the structure of the latter. The orientational information provided by MCSS is important for its use in the de novo design of large, multi-functional ligands, as well as for improving lead compounds.

pK(a) Calculations suggest storage of an excess proton in a hydrogen-bonded water network in bacteriorhodopsin

Calculations of protonation states and pK(a) values for the ionizable groups in the resting state of bacteriorhodopsin have been carried out using the recently available 1.55 A resolution X-ray crystallographic structure. The calculations are in reasonable agreement with the available experimental data for groups on or near the ion transport chain (the retinal Schiff base; Asp85, 96, 115, 212, and Arg82). In contrast to earlier studies using lower-resolution structural data, this agreement is achieved without manipulations of the crystallographically determined heavy-atom positions or ad hoc adjustments of the intrinsic pK(a) of the Schiff base. Thus, the theoretical methods used provide increased reliability as the input structural data are improved. Only minor effects on the agreement with experiment are found with respect to methodological variations, such as single versus multi-conformational treatment of hydrogen atom placements, or retaining the crystallographically determined internal water molecules versus treating them as high-dielectric cavities. The long-standing question of the identity of the group that releases a proton to the extracellular side of the membrane during the L-to-M transition of the photocycle is addressed by including as pH-titratable sites not only Glu204 and Glu194, residues near the extracellular side that have been proposed as the release group, but also an H(5)O(2)(+) molecule in a nearby cavity. The latter represents the recently proposed storage of the release proton in an hydrogen-bonded water network. In all calculations where this possibility is included, the proton is stored in the H(5)O(2)(+) rather than on either of the glutamic acids, thus establishing the plausibility on theoretical grounds of the storage of the release proton in bacteriorhodopsin in a hydrogen-bonded water network. The methods used here may also be applicable to other proteins that may store a proton in this way, such as the photosynthetic reaction center and cytochrome c oxidase.

Radially softening diffusive motions in a globular protein

Molecular dynamics simulation, quasielastic neutron scattering and analytical theory are combined to characterize diffusive motions in a hydrated protein, C-phycocyanin. The simulation-derived scattering function is in approximate agreement with experiment and is decomposed to determine the essential contributions. It is found that the geometry of the atomic motions can be modeled as diffusion in spheres with a distribution of radii. The time dependence of the dynamics follows stretched exponential behavior, reflecting a distribution of relaxation times. The average side chain and backbone dynamics are quantified and compared. The dynamical parameters are shown to present a smooth variation with distance from the core of the protein. Moving outward from the center of the protein there is a progressive increase of the mean sphere size, accompanied by a narrowing and shifting to shorter times of the relaxation time distribution. This smooth, "radially softening" dynamics may have important consequences for protein function. It also raises the possibility that the dynamical or "glass" transition with temperature observed experimentally in proteins might be depth dependent, involving, as the temperature decreases, progressive freezing out of the anharmonic dynamics with increasing distance from the center of the protein.

Role of salt bridges in homeodomains investigated by structural analyses and molecular dynamics simulations

Homeodomains are a class of helix-turn-helix DNA-binding protein motifs that play an important role in the control of cellular development in eukaryotes. They fold in a three alpha-helix structural module, where the third helix is the recognition helix that fits into the major groove of DNA. Structural analysis of the members of the homeodomain family led to the identification of interactions likely to stabilize the protein domains. Linking the helices pairwise, three salt bridges were found to be well preserved within the family. Also well conserved were two cation-pi interactions between aromatic and positively charged side chains. To analyze the structural role of the salt bridges, molecular dynamics simulations (MD) were carried out on the wild-type homeodomain from the Drosophila paired protein (1fjl) and on three mutants, which lack one or two salt bridges and mimic natural mutations in other homeodomains. Analysis of the trajectories revealed only small structural rearrangements of the three helices in all MD simulations, thereby suggesting that the salt bridges have no essential stabilizing role at room temperature, but rather might be important for improving thermostability. The latter hypothesis is supported by a good correlation between the melting midpoint temperatures of several homeodomains and the number of salt bridges and cation-pi interactions that connect secondary structures.

Quantum mechanics simulation of protein dynamics on long timescale

Protein structure and dynamics are the keys to a wide range of problems in biology. In principle, both can be fully understood by using quantum mechanics as the ultimate tool to unveil the molecular interactions involved. Indeed, quantum mechanics of atoms and molecules have come to play a central role in chemistry and physics. In practice, however, direct application of quantum mechanics to protein systems has been prohibited by the large molecular size of proteins. As a consequence, there is no general quantum mechanical treatment that not only exceeds the accuracy of state-of-the-art empirical models for proteins but also maintains the efficiency needed for extensive sampling in the conformational space, a requirement mandated by the complexity of protein systems. Here we show that, given recent developments in methods, a general quantum mechanical-based treatment can be constructed. We report a molecular dynamics simulation of a protein, crambin, in solution for 350 ps in which we combine a semiempirical quantum-mechanical description of the entire protein with a description of the surrounding solvent, and solvent-protein interactions based on a molecular mechanics force field. Comparison with a recent very high-resolution crystal structure of crambin (Jelsch et al., Proc Natl Acad Sci USA 2000;102:2246-2251) shows that geometrical detail is better reproduced in this simulation than when several alternate molecular mechanics force fields are used to describe the entire system of protein and solvent, even though the structure is no less flexible. Individual atomic charges deviate in both directions from "canonical" values, and some charge transfer is found between the N and C-termini. The capability of simulating protein dynamics on and beyond the few hundred ps timescale with a demonstrably accurate quantum mechanical model will bring new opportunities to extend our understanding of a range of basic processes in biology such as molecular recognition and enzyme catalysis.

Modeling of the metallo-beta-lactamase from B. fragilis: structural and dynamic effects of inhibitor binding

The structure and dynamics of an inhibitor-bound complex of the metallo-beta-lactamase from Bacteroides fragilis are studied by using molecular dynamics. A search of the conformational space was performed to obtain three distinct models of the complex, which were then subjected to solvated molecular dynamics. A solvated molecular dynamics study of the apo protein was performed to serve as a baseline for comparison with the bound simulations. We find loop conformation changes due to binding as well as a decrease in flexibility of the protein as a whole and especially in the major loop of the beta-lactamase. We report the structural and dynamical features of the inhibitor-bound and apo models, as well as experimentally measurable quantities, which should be capable of distinguishing the two binding modes we have determined.

Molecular dynamics study of aquaporin-1 water channel in a lipid bilayer

The aquaporin-1 water channel was modeled in a palmitoyl-oleoyl-phosphatidyl-choline lipid bilayer, by means of molecular dynamics simulations. Interaction of the protein with the membrane and inter-monomer interactions were analyzed. Structural features of the channel important for its biological function, including the Asn-Pro-Ala (NPA) motifs, and the diffusion of water molecules into the channels, were investigated. Simulations revealed the formation of single file water inside the channels for certain relative positions of the NPA motifs.

Solution structure of Co(III)-bleomycin-OOH bound to a phosphoglycolate lesion containing oligonucleotide: implications for bleomycin-induced double-strand DNA cleavage

Bleomycin (BLM) is an antitumor antibiotic that is used clinically. Its major cause of cytotoxicity is thought to be related to BLM's ability to cause double-strand (ds) DNA cleavage. A single molecule of BLM appears to cleave both strands of DNA in the presence of its required cofactors Fe(2+) and oxygen without dissociating from the helix. A mechanism for this process has been proposed based on a model structure of the hydroperoxide of Co(III)-BLM (CoBLM) bound sequence-specifically to an intact duplex containing a GTAC site, a hot spot for ds cleavage [Vanderwall, D. E., Lui, S. M., Wu, W., Turner, C. J., Kozarich, J. W., and Stubbe, J. (1997) Chem. Biol. 4, 373-387]. In this paper, we present a structural model for the second cleavage event. Two-dimensional NMR spectroscopy and molecular modeling were carried out to study CoBLM bound to d(CCAAAGXACTGGG).d(CCCAGTACTTTGG), where X represents a 3'-phosphoglycolate lesion next to a 5'-phosphate. Assignments of 729 NOEs, including 51 between the drug and the DNA and 126 within the BLM molecule, have been made. These NOEs in addition to 96 dihedral angle constraints have been used to obtain a well-defined structural model for this complex. The model reveals that the bithiazole tail is partially intercalated between the T19 and the A20 of the duplex and that the metal binding domain is poised for abstraction of the T19 H4' in the minor groove. The modeling further reveals that the predominant conformation of the bithiazole protons is trans. Two cis conformations of these protons are also observed, and ROESY experiments provide evidence for interconversion of all of these forms. The relationship of these observations to the model for ds cleavage is presented.

Dynamics of hydration in hen egg white lysozyme

We investigate the hydration dynamics of a small globular protein, hen egg-white lysozyme. Extensive simulations (two trajectories of 9 ns each) were carried out to identify the time-scales and mechanism of water attachment to this protein. The location of the surface and integral water molecules in lysozyme was also investigated. Three peculiar temporal scales of the hydration dynamics can be discerned: two among these, with sub-nanosecond mean residence time, tau(w), are characteristic of surface hydration water; the slower time-scale (tau(w) approximately 2/3 ns) is associated with buried water molecules in hydrophilic pores and in superficial clefts. The computed tau(w) values in the two independent runs fall in a similar range and are consistent with each other, thus adding extra weight to our result. The tau(w) of surface water obtained from the two independent trajectories is 20 and 24 ps. In both simulations only three water molecules are bound to lysozyme for the entire length of the trajectories, in agreement with nuclear magnetic relaxation dispersion estimates. Locations other than those identified in the protein crystal are found to be possible for these long-residing water molecules. The dynamics of the hydration water molecules observed in our simulations implies that each water molecule visits a multitude of residues during the lifetime of its bound with the protein. The number of residues seen by a single water molecule increases with the time-scale of its residence time and, on average, is equal to one only for the water molecules with shorter residence time. Thus, tau(w) values obtained from inelastic neutron scattering and based on jump-diffusion models are likely not to account for the contribution of water molecules with longer residence time.

Biomolecular simulations: recent developments in force fields, simulations of enzyme catalysis, protein-ligand, protein-protein, and protein-nucleic acid noncovalent interactions

Computer modeling has been developed and widely applied in studying molecules of biological interest. The force field is the cornerstone of computer simulations, and many force fields have been developed and successfully applied in these simulations. Two interesting areas are (a) studying enzyme catalytic mechanisms using a combination of quantum mechanics and molecular mechanics, and (b) studying macromolecular dynamics and interactions using molecular dynamics (MD) and free energy (FE) calculation methods. Enzyme catalysis involves forming and breaking of covalent bonds and requires the use of quantum mechanics. Noncovalent interactions appear ubiquitously in biology, but here we confine ourselves to review only noncovalent interactions between protein and protein, protein and ligand, and protein and nucleic acids.

Extending the accuracy limits of prediction for side-chain conformations

Current techniques for the prediction of side-chain conformations on a fixed backbone have an accuracy limit of about 1.0-1.5 A rmsd for core residues. We have carried out a detailed and systematic analysis of the factors that influence the prediction of side-chain conformation and, on this basis, have succeeded in extending the limits of side-chain prediction for core residues to about 0.7 A rmsd from native, and 94 % and 89 % of chi(1) and chi(1+2 ) dihedral angles correctly predicted to within 20 degrees of native, respectively. These results are obtained using a force-field that accounts for only van der Waals interactions and torsional potentials. Prediction accuracy is strongly dependent on the rotamer library used. That is, a complete and detailed rotamer library is essential. The greatest accuracy was obtained with an extensive rotamer library, containing over 7560 members, in which bond lengths and bond angles were taken from the database rather than simply assuming idealized values. Perhaps the most surprising finding is that the combinatorial problem normally associated with the prediction of the side-chain conformation does not appear to be important. This conclusion is based on the fact that the prediction of the conformation of a single side-chain with all others fixed in their native conformations is only slightly more accurate than the simultaneous prediction of all side-chain dihedral angles.

Probing the origins of increased activity of the E22Q "Dutch" mutant Alzheimer's beta-amyloid peptide

The amyloid peptide congener A beta(10--35)-NH(2) is simulated in an aqueous environment in both the wild type (WT) and E22Q "Dutch" mutant forms. The origin of the noted increase in deposition activity resulting from the Dutch mutation is investigated. Multiple nanosecond time scale molecular dynamics trajectories were performed and analyzed using a variety of measures of the peptide's average structure, hydration, conformational fluctuations, and dynamics. The results of the study support the conclusions that 1) the E22Q mutant and WT peptide are both stable in "collapsed coil" conformations consistent with the WT structure of, J. Struct. Biol. 130:130--141); 2) the E22Q peptide is more flexible in solution, supporting early claims that its equilibrium structural fluctuations are larger than those of the WT peptide; and 3) the local E22Q mutation leads to a change in the first solvation layer in the region of the peptide's "hydrophobic patch," resulting in a more hydrophobic solvation of the mutant peptide. The simulation results support the view that the noted increase in activity due to the Dutch mutation results from an enhancement of the desolvation process that is an essential step in the aggregation of the peptide.

Extracellular blockade of K(+) channels by TEA: results from molecular dynamics simulations of the KcsA channel

TEA is a classical blocker of K(+) channels. From mutagenesis studies, it has been shown that external blockade by TEA is strongly dependent upon the presence of aromatic residue at Shaker position 449 which is located near the extracellular entrance to the pore (Heginbotham, L., and R. MacKinnon. 1992. Neuron. 8:483-491). The data suggest that TEA interacts simultaneously with the aromatic residues of the four monomers. The determination of the 3-D structure of the KcsA channel using X-ray crystallography (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) has raised some issues that remain currently unresolved concerning the interpretation of these observations. In particular, the center of the Tyr82 side chains in KcsA (corresponding to position 449 in Shaker) forms a square of 11.8-A side, a distance which is too large to allow simultaneous interactions of a TEA molecule with the four aromatic side chains. In this paper, the external blockade by TEA is explored by molecular dynamics simulations of an atomic model of KcsA in an explicit phospholipid bilayer with aqueous salt solution. It is observed, in qualitative accord with the experimental results, that TEA is stable when bound to the external side of the wild-type KcsA channel (with Tyr82), but is unstable when bound to a mutant channel in which the tyrosine residue has been substituted by a threonine. The free energy profile of TEA relative to the pore is calculated using umbrella sampling simulations to characterize quantitatively the extracellular blockade. It is found, in remarkable agreement with the experiment, that the TEA is more stably bound by 2.3 kcal/mol to the channel with four tyrosine residues. In the case of the wild-type KcsA channel, TEA (which has the shape of a flattened oblate spheroid) acts as an ideal plug blocking the pore. In contrast, it is considerably more off-centered and tilted in the case of the mutant channel. The enhanced stability conferred by the tyrosine residues does not arise from Pi-cation interactions, but appears to be due to differences in the hydration structure of the TEA. Finally, it is shown that the experimentally observed voltage dependence of TEA block, which is traditionally interpreted in terms of the physical position of the TEA along the axis of the pore, must arise indirectly via coupling with the ions in the pore.

Interaction of hypericin with serum albumins: surface-enhanced Raman spectroscopy, resonance Raman spectroscopy and molecular modeling study

Surface-enhanced Raman spectroscopy, resonance Raman spectroscopy and molecular modeling were employed to study the interaction of hypericin (Hyp) with human (HSA), rat (RSA) and bovine (BSA) serum albumins. The identification of the binding site of Hyp in serum albumins as well as the structural model for Hyp/HSA complex are presented. The interactions mainly reflect: (1) a change of the strength of H bonding at the N1-H site of Trp; (2) a change of the Trp side-chain conformation; (3) a change of the hydrophobicity of the Trp environment; and (4) a formation of an H-bond between the carbonyl group of Hyp and a proton donor in HSA and RSA which leads to a protonated-like carbonyl in Hyp. Our results indicate that Hyp is rigidly bound in IIA subdomain of HSA close to Trp214 (distance 5.12 A between the centers of masses). In the model presented the carbonyl group of Hyp is hydrogen bonded to Asn458. Two other candidates for hydrogen bonds have been identified between the bay-region hydroxyl group of Hyp and the carbonyl group of the Trp214 peptidic link and between the peri-region hydroxyl group of Hyp and the Asn458 carbonyl group. It is shown that the structures of the Hyp/HSA and Hyp/RSA complexes are similar to, and in some aspects different from, those found for the Hyp/BSA complex. The role of aminoacid sequence in the IIA subdomains of HSA, RSA and BSA is discussed to explain the observed differences.

Sequence determination of reduction potentials by cysteinyl hydrogen bonds and peptide pipoles in [4Fe-4S] ferredoxins

A sequence determinant of reduction potentials is reported for bacterial [4Fe-4S]-type ferredoxins. The residue that is four residues C-terminal to the fourth ligand of either cluster is generally an alanine or a cysteine. In five experimental ferredoxin structures, the cysteine has the same structural orientation relative to the nearest cluster, which is stabilized by the SH...S bond. Although such bonds are generally considered weak, indications that Fe-S redox site sulfurs are better hydrogen-bond acceptors than most sulfurs include the numerous amide NH...S bonds noted by Adman and our quantum mechanical calculations. Furthermore, electrostatic potential calculations of 11 experimental ferredoxin structures indicate that the extra cysteine decreases the reduction potential relative to an alanine by approximately 60 mV, in agreement with experimental mutational studies. Moreover, the decrease in potential is due to a shift in the polar backbone stabilized by the SH...S bond rather than to the slightly polar cysteinyl side chain. Thus, these cysteines can "tune" the reduction potential, which could optimize electron flow in an electron transport chain. More generally, hydrogen bonds involving sulfur can be important in protein structure/function, and mutations causing polar backbone shifts can alter electrostatics and thus affect redox properties or even enzymatic activity of a protein.

Atomic charges in molecular mechanical force fields: a theoretical insight

Based on quantum theory of the Coulombic interactions between a molecule and its surrounding molecules, a theoretical derivation is presented to obtain an atomic charge model. The charge model shows that the appropriate atomic charge, used for example in molecular mechanical force fields, is simply the average value of the electrostatic potential (ESP) derived charge for an isolated molecule and that obtained for the molecule in the interaction system. Computational scheme to calculate the ESP derived atomic charges of a molecule in the interaction system is presented. The method is applied to two common liquids with associated hydrogen bonds: water and methanol. The obtained theoretical atomic charges are similar to those found in the common interaction potential models, such as: SPC, TIP3P, OPLS, etc. However, for methanol they differ considerably from those obtained using the restrained ESP method.

Free energy decomposition of protein-protein interactions

A free energy decomposition scheme has been developed and tested on antibody-antigen and protease-inhibitor binding for which accurate experimental structures were available for both free and bound proteins. Using the x-ray coordinates of the free and bound proteins, the absolute binding free energy was computed assuming additivity of three well-defined, physical processes: desolvation of the x-ray structures, isomerization of the x-ray conformation to a nearby local minimum in the gas-phase, and subsequent noncovalent complex formation in the gas phase. This free energy scheme, together with the Generalized Born model for computing the electrostatic solvation free energy, yielded binding free energies in remarkable agreement with experimental data. Two assumptions commonly used in theoretical treatments; viz., the rigid-binding approximation (which assumes no conformational change upon complexation) and the neglect of vdW interactions, were found to yield large errors in the binding free energy. Protein-protein vdW and electrostatic interactions between complementary surfaces over a relatively large area (1400--1700 A(2)) were found to drive antibody-antigen and protease-inhibitor binding.

Reorientational eigenmode dynamics: a combined MD/NMR relaxation analysis method for flexible parts in globular proteins

An approach is presented for the interpretation of heteronuclear NMR spin relaxation data in mobile protein parts in terms of reorientational eigenmode dynamics. The method is based on the covariance matrix of the spatial functions of the nuclear spin interactions that cause relaxation expressed as spherical harmonics of rank 2. The approach was applied to characterize the dynamics of a loop region of ubiquitin. The covariance matrix was determined from a conformational ensemble generated by a 5 ns molecular dynamics simulation. It was found that the time correlation functions of the dominant eigenmodes decay in good approximation with a single correlation time. From the reorientational eigenmodes, their eigenvalues, and correlation times, NMR relaxation data were calculated in accordance with Bloch-Wangsness-Redfield relaxation theory and directly compared with experimental (15)N relaxation parameters. Using a fitting procedure, agreement between calculated and experimental data was improved significantly by adjusting eigenvalues and correlation times of the dominant modes. The presented procedure provides detailed information on correlated reorientational dynamics of flexible parts in globular proteins. The covariance matrix was linked to the covariance matrix of backbone dihedral angle fluctuations, allowing one to study the motional behavior of these degrees of freedom on nano- and subnanosecond time scales.

Substrate conformational transitions in the active site of chorismate mutase: their role in the catalytic mechanism

Chorismate mutase acts at the first branch-point of aromatic amino acid biosynthesis and catalyzes the conversion of chorismate to prephenate. The results of molecular dynamics simulations of the substrate in solution and in the active site of chorismate mutase are reported. Two nonreactive conformers of chorismate are found to be more stable than the reactive pseudodiaxial chair conformer in solution. It is shown by QM/MM molecular dynamics simulations, which take into account the motions of the enzyme, that when these inactive conformers are bound to the active site, they are rapidly converted to the reactive chair conformer. This result suggests that one contribution of the enzyme is to bind the more prevalent nonreactive conformers and transform them into the active form in a step before the chemical reaction. The motion of the reactive chair conformer in the active site calculated by using the QM/MM potential generates transient structures that are closer to the transition state than is the stable CHAIR conformer.

Structure, dynamics and electrostatics of the active site of glutaredoxin 3 from Escherichia coli: comparison with functionally related proteins

The chemistry of active-site cysteine residues is central to the activity of thiol-disulfide oxidoreductases of the thioredoxin superfamily. In these reactions, a nucleophilic thiolate is required, but the associated pK(a) values differ vastly in the superfamily, from less than 4 in DsbA to greater than 7 in Trx. The factors that stabilize this thiolate are, however, not clearly established. The glutaredoxins (Grxs), which are members of this superfamily, contain a Cys-Pro-Tyr-Cys motif in their active site. In reduced Grxs, the pK(a) of the N-terminal active-site nucleophilic cysteine residue is lowered significantly, and the stabilization of the corresponding thiolate is expected to influence the redox potential of these enzymes. Here, we use a combination of long molecular dynamics (MD) simulations, pK(a) calculations, and experimental investigations to derive the structure and dynamics of the reduced active site from Escherichia coli Grx3, and investigate the factors that stabilize the thiolate. Several different MD simulations converged toward a consensus conformation for the active-site cysteine residues (Cys11 and Cys14), after a number of local conformational changes. Key features of the model were tested experimentally by measurement of NMR scalar coupling constants, and determination of pK(a) values of selected residues. The pK(a) values of the Grx3 active-site residues were calculated during the MD simulations, and support the underlying structural model. The structure of Grx3, in combination with the pK(a) calculations, indicate that the pK(a) of the N-terminal active-site cysteine residue in Grx3 is intermediate between that of its counterpart in DsbA and Trx. The pK(a) values in best agreement with experiment are obtained with a low (<4) protein dielectric constant. The calculated pK(a) values fluctuate significantly in response to protein dynamics, which underscores the importance of the details of the underlying structures when calculating pK(a) values. The thiolate of Cys11 is stabilized primarily by direct hydrogen bonding with the amide protons of Tyr13 and Cys14 and the thiol proton of Cys14, rather than by long-range interactions from charged groups or from a helix macrodipole. From the comparison of reduced Grx3 with other members of the thioredoxin superfamily, a unifying theme for the structural basis of thiol pK(a) differences in this superfamily begins to emerge.

A barley polyamine oxidase isoform with distinct structural features and subcellular localization

Two cDNAs encoding polyamine oxidase (PAO) isoforms (BPAO1 and BPAO2) and the corresponding gene copies were isolated from barley cultivar Aura. Gene organization is not conserved between these two nonallelic coding sequences. Both precursor proteins include a cleavable N-terminal leader of 25 amino acids. N-terminal sequencing of PAO purified from barley seedlings reveals a unique amino-acid sequence corresponding to the BPAO2 N-terminus as predicted from the corresponding cDNA. BPAO2 has been purified, characterized and compared to maize PAO (MPAO), the best characterized member of this enzyme class. The two proteins show different pH optima for catalytic activity, Km and Vmax values with spermidine and spermine as substrates. Molecular modelling of BPAO2 reveals the same global fold as in MPAO. However, substitution of the active site residue Phe403 by a tyrosine, provides a rationale for the different catalytic properties of the two enzymes. In barley leaves PAO-specific activity is higher in isolated mesophyll protoplasts than in the extracellular fluids, whereas in maize the reverse is true. The C-terminus of BPAO2 shows homology with the endoplasmic reticulum retention signal that might be responsible for the subcellular localization observed. We conclude that BPAO2 is a symplastic PAO in barley mesophyll cells. Production of BPAO2 mRNA and the corresponding protein is induced by light, and has a different pattern of accumulation in leaves and coleoptiles.

Lipid-mediated interactions between intrinsic membrane proteins: dependence on protein size and lipid composition

The present study is an application of an approach recently developed by the authors for describing the structure of the hydrocarbon chains of lipid-bilayer membranes (LBMs) around embedded protein inclusions ( Biophys. J. 79:2867-2879). The approach is based on statistical mechanical integral equation theories developed for the study of dense liquids. First, the configurations extracted from molecular dynamics simulations of pure LBMs are used to extract the lateral density-density response function. Different pure LBMs composed of different lipid molecules were considered: dioleoyl phosphatidylcholine (DOPC), palmitoyl-oleoyl phosphatidylcholine (POPC), dipalmitoyl phosphatidylcholine (DPPC), and dimyristoyl phosphatidylcholine (DMPC). The results for the lateral density-density response function was then used as input in the integral equation theory. Numerical calculations were performed for protein inclusions of three different sizes. For the sake of simplicity, protein inclusions are represented as hard smooth cylinders excluding the lipid hydrocarbon core from a small cylinder of 2.5 A radius, corresponding roughly to one aliphatic chain, a medium cylinder of 5 A radius, corresponding to one alpha-helix, and a larger cylinder of 9 A radius, representing a small protein such as the gramicidin channel. The lipid-mediated interaction between protein inclusions was calculated using a closed-form expression for the configuration-dependent free energy. This interaction was found to be repulsive at intermediate range and attractive at short range for two small cylinders in POPC, DPPC, and DMPC bilayers, whereas it oscillates between attractive and repulsive values in DOPC bilayers. For medium size cylinders, it is again repulsive at intermediate range and attractive at short range, but for every model LBM considered here. In the case of a large cylinder, the lipid-mediated interaction was shown to be repulsive for both short and long ranges for the DOPC, POPC, and DPPC bilayers, whereas it is again repulsive and attractive for DMPC bilayers. The results indicate that the packing of the hydrocarbon chains around protein inclusions in LBMs gives rise to a generic (i.e., nonspecific) lipid-mediated interaction which favors the association of two alpha-helices and depends on the lipid composition of the membrane.

Conformational dynamics of an alanine dipeptide analog: an ab initio molecular dynamics study

An ab initio molecular dynamics (MD) simulation technique employing the Born-Oppenheimer approach in the framework of a Gaussian implementation of Kohn-Sham density functional theory is used to study the gas-phase conformational dynamics of an alanine dipeptide analog. It is found that conformational transformation between C5 and C7(eq) occurs on the picosecond time scale. Classical MD simulations using most of the popular force fields do not yield a transition even after nanoseconds. An analysis is given of the difference, for this small gas-phase system, between ab initio MD and traditional MD simulation using force fields.

Influence of the heme pocket conformation on the structure and vibrations of the Fe-CO bond in myoglobin: a QM/MM density functional study

The influence of the distal pocket conformation on the structure and vibrations of the heme-CO bond in carbonmonoxy myoglobin (MbCO) is investigated by means of hybrid QM/MM calculations based on density functional theory combined with a classical force field. It is shown that the heme-CO structure (QM treated) is quite rigid and not influenced by the distal pocket conformation (MM treated). This excludes any relation between FeCO distortions and the different CO absorptions observed in the infrared spectra of MbCO (A states). In contrast, both the CO stretch frequency and the strength of the CO...His64 interaction are very dependent on the orientation and tautomerization state of His64. Our calculations indicate that the CO...N(epsilon) type of approach does not contribute to the A states, whereas the CO...H-N(epsilon) interaction is the origin of the A(1) and A(3) states, the His64 residue being protonated at N(epsilon). The strength of the CO...His64 interaction is quantified, in comparison with the analogous O(2)...His64 interaction and with the observed changes in the CO stretch frequency. Additional aspects of the CO...His64 interaction and its biological implications are discussed.

Artificial cytochrome b: computer modeling and evaluation of redox potentials

We generated atomic coordinates of an artificial protein that was recently synthesized to model the central part of the native cytochrome b (Cb) subunit consisting of a four-helix bundle with two hemes. Since no X-ray structure is available, the structural elements of the artificial Cb were assembled from scratch using all known chemical and structural information available and avoiding strain as much as possible. Molecular dynamics (MD) simulations applied to this model protein exhibited root-mean-square deviations as small as those obtained from MD simulations starting with the crystal structure of the native Cb subunit. This demonstrates that the modeled structure of the artificial Cb is relatively rigid and strain-free. The model structure of the artificial Cb was used to determine the redox potentials of the two hemes by calculating the electrostatic energies from the solution of the linearized Poisson-Boltzmann equation (LPBE). The calculated redox potentials agree within 20 meV with the experimentally measured values. The dependence of the redox potentials of the hemes on the protein environment was analyzed. Accordingly, the total shift in the redox potentials is mainly due to the low dielectric medium of the protein, the protein backbone charges, and the salt bridges formed between the arginines and the propionic acid groups of the hemes. The difference in the shift of the redox potentials is due to the interactions with the hydrophilic side chains and the salt bridges formed with the propionic acids of the hemes. For comparison and to test the computational procedure, the redox potentials of the two hemes in the native Cb from the cytochrome bc(1) (Cbc(1)) complex were also calculated. Also in this case the computed redox potentials agree well with experiments.

Model-free approach to the dynamic interpretation of residual dipolar couplings in globular proteins

The effects of internal motions on residual dipolar NMR couplings of proteins partially aligned in a liquid-crystalline environment are analyzed using a 10 ns molecular dynamics (MD) computer simulation of ubiquitin. For a set of alignment tensors with different orientations and rhombicities, MD-averaged dipolar couplings are determined and subsequently interpreted for different scenarios in terms of effective alignment tensors, average orientations of dipolar vectors, and intramolecular reorientational vector distributions. Analytical relationships are derived that reflect similarities and differences between motional scaling of dipolar couplings and scaling of dipolar relaxation data (NMR order parameters). Application of the self-consistent procedure presented here to dipolar coupling measurements of biomolecules aligned in different liquid-crystalline media should allow one to extract in a "model-free" way average orientations of dipolar vectors and specific aspects of their motions.

Proposed lipocalin fold for apolipoprotein M based on bioinformatics and site-directed mutagenesis

Apolipoprotein M (apoM) is a novel apolipoprotein that is predominantly present in high-density lipoprotein. Sensitive sequence searches, threading and comparative model building experiments revealed apoM to be structurally related to the lipocalin protein family. In a 3D model, characterized by an eight-stranded anti-parallel beta-barrel, a segment including Asn135 could adopt a closed or open conformation. Using site-directed mutagenesis, we demonstrated Asn135 in wild-type apoM to be glycosylated, suggesting that the segment is solvent exposed. ApoM displays two strong acidic patches of potential functional importance, one around the N-terminus and the other next to the opening of the beta-barrel.

An optimized protein kinase C activating diacylglycerol combining high binding affinity (Ki) with reduced lipophilicity (log P)

A small, focused combinatorial library encompassing all possible permutations of acyl branched alkyl chains-small and large, saturated and unsaturated-was generated from the active diacylglycerol enantiomer (S-DAG) to help identify the analogue with the highest binding affinity (lowest Ki) for protein kinase C (PK-C) combined with the minimum lipophilicity (log P). The selected ligand (3B) activated PK-C more effectively than sn-1,2-dioctanoylglycerol (diC8) despite being 1.4 log units more hydrophilic. Compound 3B indeed represents the most potent, hydrophilic DAG ligand to date. With the help of a green fluorescent protein (GFP)-tagged PK-Calpha, 3B was able to translocate the full length protein to the membrane with an optimal dose of 100 microM in CHO-K1 cells, while diC8 failed to achieve translocation even at doses 3-fold higher. Molecular modeling of 3B into an empty C1b domain of PK-Cdelta clearly showed the existence of a preferred binding orientation. In addition, molecular dynamic simulations suggest that binding discrimination could result from a favorable van der Waals (VDW) interaction between the large, branched sn-1 acyl group of 3B and the aromatic rings of Trp252 (PK-Cdelta) or Tyr252 (PK-Calpha). The DAG analogue of 3B in which the acyl groups are reversed (2C) showed a decrease in binding affinity reflecting the capacity of PK-C to effectively discriminate between alternative orientations of the acyl chains.

Predicting sequences and structures of MHC-binding peptides: a computational combinatorial approach

Peptides bound to MHC molecules on the surface of cells convey critical information about the cellular milieu to immune system T cells. Predicting which peptides can bind an MHC molecule, and understanding their modes of binding, are important in order to design better diagnostic and therapeutic agents for infectious and autoimmune diseases. Due to the difficulty of obtaining sufficient experimental binding data for each human MHC molecule, computational modeling of MHC peptide-binding properties is necessary. This paper describes a computational combinatorial design approach to the prediction of peptides that bind an MHC molecule of known X-ray crystallographic or NMR-determined structure. The procedure uses chemical fragments as models for amino acid residues and produces a set of sequences for peptides predicted to bind in the MHC peptide-binding groove. The probabilities for specific amino acids occurring at each position of the peptide are calculated based on these sequences, and these probabilities show a good agreement with amino acid distributions derived from a MHC-binding peptide database. The method also enables prediction of the three-dimensional structure of MHC-peptide complexes. Docking, linking, and optimization procedures were performed with the XPLOR program [1].