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

Electrostatic properties of bovine beta-lactoglobulin

Bovine beta-Lactoglobulin (BLG) has been studied for many decades, but only recently structural data have been obtained, making it possible to simulate its molecular properties. In the present study, electrostatic properties of BLG are investigated theoretically using Poisson-Boltzmann calculations and experimentally following pH titration via NMR. Electrostatic properties are determined for several structural models, including an ensemble of NMR structures obtained at low pH. The changes in electrostatic forces upon changes in ionic strength, solvent dielectric constant, and pH are calculated and compared with experiments. pK(a)s are computed for all titratable sites and compared with NMR titration data. The analysis of theoretical and experimental results suggests that (1) there may be more than one binding sites for negatively charged ligands; (2) at low pH the core of the molecule is more compact than observed in the structures obtained via restrained molecular dynamics from NMR data, but loop and terminal regions must be disordered.

Molecular dynamics of the KcsA K(+) channel in a bilayer membrane

Molecular dynamics (MD) simulations of an atomic model of the KcsA K(+) channel embedded in an explicit dipalmitoylphosphatidylcholine (DPPC) phospholipid bilayer solvated by a 150 mM KCl aqueous salt solution are performed and analyzed. The model includes the KcsA K(+) channel, based on the recent crystallographic structure of, Science. 280:69-77), 112 DPPC, K(+) and Cl(-) ions, as well as over 6500 water molecules for a total of more than 40,000 atoms. Three K(+) ions are explicitly included in the pore. Two are positioned in the selectivity filter on the extracellular side and one in the large water-filled cavity. Different starting configurations of the ions and water molecules in the selectivity filter are considered, and MD trajectories are generated for more than 4 ns. The conformation of KcsA is very stable in all of the trajectories, with a global backbone root mean square (RMS) deviation of less than 1.9 A with respect to the crystallographic structure. The RMS atomic fluctuations of the residues surrounding the selectivity filter on the extracellular side of the channel are significantly lower than those on the intracellular side. The motion of the residues with aromatic side chains surrounding the selectivity filter (Trp(67), Trp(68), Tyr(78), and Tyr(82)) is anisotropic with the smallest RMS fluctuations in the direction parallel to the membrane plane. A concerted dynamic transition of the three K(+) ions in the pore is observed, during which the K(+) ion located initially in the cavity moves into the narrow part of the selectivity filter, while the other two K(+) ions move toward the extracellular side. A single water molecule is stabilized between each pair of ions during the transition, suggesting that each K(+) cation translocating through the narrow pore is accompanied by exactly one water molecule, in accord with streaming potential measurements (, Biophys. J. 55:367-371). The displacement of the ions is coupled with the structural fluctuations of Val(76) and Gly(77), in the selectivity filter, as well as the side chains of Glu(71), Asp(80), and Arg(89), near the extracellular side. Thus the mechanical response of the channel structure at distances as large as 10-20 A from the ions in the selectivity filter appears to play an important role in the concerted transition.

Substrate recognition by the Lyn protein-tyrosine kinase. NMR structure of the immunoreceptor tyrosine-based activation motif signaling region of the B cell antigen receptor

The immunoreceptor tyrosine-based activation motif (ITAM) plays a central role in transmembrane signal transduction in hematopoietic cells by mediating responses leading to proliferation and differentiation. An initial signaling event following activation of the B cell antigen receptor is phosphorylation of the CD79a (Ig-alpha) ITAM by Lyn, a Src family protein-tyrosine kinase. To elucidate the structural basis for recognition between the ITAM substrate and activated Lyn kinase, the structure of an ITAM-derived peptide bound to Lyn was determined using exchange-transferred nuclear Overhauser NMR spectroscopy. The bound substrate structure has an irregular helix-like character. Docking based on the NMR data into the active site of the closely related Lck kinase strongly favors ITAM binding in an orientation similar to binding of cyclic AMP-dependent protein kinase rather than that of insulin receptor tyrosine kinase. The model of the complex provides a rationale for conserved ITAM residues, substrate specificity, and suggests that substrate binds only the active conformation of the Src family tyrosine kinase, unlike the ATP cofactor, which can bind the inactive form.

Pathways of ligand clearance in acetylcholinesterase by multiple copy sampling

The clearance of seven different ligands from the deeply buried active-site of Torpedo californica acetylcholinesterase is investigated by combining multiple copy sampling molecular dynamics simulations, with the analysis of protein-ligand interactions, protein motion and the electrostatic potential sampled by the ligand copies along their journey outwards. The considered ligands are the cations ammonium, methylammonium, and tetramethylammonium, the hydrophobic methane and neopentane, and the anionic product acetate and its neutral form, acetic acid. We find that the pathways explored by the different ligands vary with ligand size and chemical properties. Very small ligands, such as ammonium and methane, exit through several routes. One involves the main exit through the mouth of the enzyme gorge, another is through the so-called back door near Trp84, and a third uses a side door at a direction of approximately 45 degrees to the main exit. The larger polar ligands, methylammonium and acetic acid, leave through the main exit, but the bulkiest, tetramethylammonium and neopentane, as well as the smaller acetate ion, remain trapped in the enzyme gorge during the time of the simulations. The pattern of protein-ligand contacts during the diffusion process is highly non-random and differs for different ligands. A majority is made with aromatic side-chains, but classical H-bonds are also formed. In the case of acetate, but not acetic acid, the anionic and neutral form, respectively, of one of the reaction products, specific electrostatic interactions with protein groups, seem to slow ligand motion and interfere with protein flexibility; protonation of the acetate ion is therefore suggested to facilitate clearance. The Poisson-Boltzmann formalism is used to compute the electrostatic potential of the thermally fluctuating acetylcholinesterase protein at positions actually visited by the diffusing ligand copies. Ligands of different charge and size are shown to sample somewhat different electrostatic potentials during their migration, because they explore different microscopic routes. The potential along the clearance route of a cation such as methylammonium displays two clear minima at the active and peripheral anionic site. We find moreover that the electrostatic energy barrier that the cation needs to overcome when moving between these two sites is small in both directions, being of the order of the ligand kinetic energy. The peripheral site thus appears to play a role in trapping inbound cationic ligands as well as in cation clearance, and hence in product release.

Experiments and simulations of ion-enhanced interfacial chemistry on aqueous NaCl aerosols

A combination of experimental, molecular dynamics, and kinetics modeling studies is applied to a system of concentrated aqueous sodium chloride particles suspended in air at room temperature with ozone, irradiated at 254 nanometers to generate hydroxyl radicals. Measurements of the observed gaseous molecular chlorine product are explainable only if reactions at the air-water interface are dominant. Molecular dynamics simulations show the availability of substantial amounts of chloride ions for reaction at the interface, and quantum chemical calculations predict that in the gas phase chloride ions will strongly attract hydroxl radicals. Model extrapolation to the marine boundary layer yields daytime chlorine atom concentrations that are in good agreement with estimates based on field measurements of the decay of selected organics over the Southern Ocean and the North Atlantic. Thus, ion-enhanced interactions with gases at aqueous interfaces may play a more generalized and important role in the chemistry of concentrated inorganic salt solutions than was previously recognized.

Diffusion of water molecules in crystalline beta-cyclodextrin hydrates

To understand the rapid diffusion mechanism of water molecules in the crystal lattice of hydrated beta-cyclodextrin (beta-CD), molecular dynamics (MD) simulations of crystalline beta-CD were performed at five different relative humidities corresponding to hydration states ranging from beta-CD-9.4H2O to beta-CD-12.3H2O, and in aqueous solution. The trajectories for the crystalline beta-CD hydrates had lengths of 4 ns each, whereas the simulation in aqueous solution extended to 2 ns. Transport of water molecules in the crystal was characterized in terms of a spatially varying diffusion constant and the main direction of diffusion, which were evaluated using newly developed algorithms. The main diffusion pathway winds through the cavities of adjacent doughnut shaped beta-CD molecules and is slightly slanted with respect to the crystallographic b-axis. Water molecules outside the beta-CD cavities have access to the main diffusion pathway. The diffusion constant for transport of water molecules along the main pathway calculated from the MD simulation data adopts 1/30 of the value in bulk water at room temperature. This is in agreement with estimates that can be made from experimental data on the adjustment of a beta-CD crystal to changes in relative humidity.

Protein-induced membrane disorder: a molecular dynamics study of melittin in a dipalmitoylphosphatidylcholine bilayer

A molecular dynamics simulation of melittin in a hydrated dipalmitoylphosphatidylcholine (DPPC) bilayer was performed. The 19, 000-atom system included a 72-DPPC phospholipid bilayer, a 26-amino acid peptide, and more than 3000 water molecules. The N-terminus of the peptide was protonated and embedded in the membrane in a transbilayer orientation perpendicular to the surface. The simulation results show that the peptide affects the lower (intracellular) layer of the bilayer more strongly than the upper (extracellular) layer. The simulation results can be interpreted as indicating an increased level of disorder and structural deformation for lower-layer phospholipids in the immediate vicinity of the peptide. This conclusion is supported by the calculated deuterium order parameters, the observed deformation at the intracellular interface, and an increase in fractional free volume. The upper layer was less affected by the embedded peptide, except for an acquired tilt relative to the bilayer normal. The effect of melittin on the surrounding membrane is localized to its immediate vicinity, and its asymmetry with respect to the two layers may result from the fact that it is not fully transmembranal. Melittin's hydrophilic C-terminus anchors it at the extracellular interface, leaving the N-terminus "loose" in the lower layer of the membrane. In general, the simulation supports a role for local deformation and water penetration in melittin-induced lysis. As for the peptide, like other membrane-embedded polypeptides, melittin adopts a significant 25 degree tilt relative to the membrane normal. This tilt is correlated with a comparable tilt of the lipids in the upper membrane layer. The peptide itself retains an overall helical structure throughout the simulation (with the exception of the three N-terminal residues), adopting a 30 degree intrahelical bend angle.

Electrostatic stress in catalysis: structure and mechanism of the enzyme orotidine monophosphate decarboxylase

Orotidine 5'-monophosphate decarboxylase catalyzes the conversion of orotidine 5'-monophosphate to uridine 5'-monophosphate, the last step in biosynthesis of pyrimidine nucleotides. As part of a Structural Genomics Initiative, the crystal structures of the ligand-free and the6-azauridine 5'-monophosphate-complexed forms have been determined at 1.8 and 1.5 A, respectively. The protein assumes a TIM-barrel fold with one side of the barrel closed off and the other side binding the inhibitor. A unique array of alternating charges (Lys-Asp-Lys-Asp) in the active site prompted us to apply quantum mechanical and molecular dynamics calculations to analyze the relative contributions of ground state destabilization and transition state stabilization to catalysis. The remarkable catalytic power of orotidine 5'-monophosphate decarboxylase is almost exclusively achieved via destabilization of the reactive part of the substrate, which is compensated for by strong binding of the phosphate and ribose groups. The computational results are consistent with a catalytic mechanism that is characterized by Jencks's Circe effect.

Correlation between changes in nuclear magnetic resonance order parameters and conformational entropy: molecular dynamics simulations of native and denatured staphylococcal nuclease

Recent work has suggested that changes in NMR order parameters may quantitatively reflect changes in the conformational entropy of a protein ensemble. The extent of the mathematical relationship between local entropy changes as seen by NMR order parameters and the full protein entropy change is a complex issue. As a step towards a fuller understanding of this problem, molecular dynamics calculations of both native and denatured staphylococcal nuclease were performed. The N-H bond vector motion, in both explicit and implicit solvent, was analyzed to estimate local and global entropy changes. The calculated N-H bond vector order parameters from simulation agreed on average with experimental values for both native and denatured structures. However, the inverted-U profile of order parameters versus residue number observed experimentally for denatured nuclease was only partially reproduced by simulation of compact denatured structures. Comparisons made across the full set of simulations revealed a correlation between the N-H order parameter-based conformational entropy change and the total quasiharmonic-based conformational entropy change between the native and denatured structures. The calculations showed that about 25% of the total entropy change was reflected by changes in simulated S2 values. This result suggests that NMR-derived order parameters may be used to provide a reasonable estimate of the total conformational entropy change on protein folding.

Molecular dynamics study of the nature and origin of retinal's twisted structure in bacteriorhodopsin

The planarity of the polyene chain of the retinal chromophore in bacteriorhodopsin is studied using molecular dynamics simulation techniques and applying different force-field parameters and starting crystal structures. The largest deviations from a planar structure are observed for the C(13)==C(14) and C(15)==N(16) double bonds in the retinal Schiff base structure. The other dihedral angles along the polyene chain of the chromophore, although having lower torsional barriers in some cases, do not significantly deviate from the planar structure. The results of the simulations of different mutants of the pigment show that, among the studied amino acids of the binding pocket, the side chain of Trp-86 has the largest impact on the planarity of retinal, and the mutation of this amino acid to alanine leads to chromophore planarity. Deletion of the methyl C(20), removal of a water molecule hydrogen-bonded to H(15), or mutation of other amino acids to alanine did not show any significant influence on the distortion of the chromophore. The results from the present study suggest the importance of the bulky residue of Trp-86 in the isomerization process, in both ground and excited states of the chromophore, and in fine-tuning of the pK(a) of the retinal protonated Schiff base in bacteriorhodopsin. The dark adaptation of the pigment and the last step of the bacteriorhodopsin photocycle imply low barriers against the rotation of the double bonds in the Schiff base region. The twisted double bonds found in the present study are consistent with the proposed mechanism of these ground state isomerization events.

Vibrational population relaxation of carbon monoxide in the heme pocket of photolyzed carbonmonoxy myoglobin: comparison of time-resolved mid-IR absorbance experiments and molecular dynamics simulations

The vibrational energy relaxation of carbon monoxide in the heme pocket of sperm whale myoglobin was studied by using molecular dynamics simulation and normal mode analysis methods. Molecular dynamics trajectories of solvated myoglobin were run at 300 K for both the delta- and epsilon-tautomers of the distal His-64. Vibrational population relaxation times of 335 +/- 115 ps for the delta-tautomer and 640 +/- 185 ps for the epsilon-tautomer were estimated by using the Landau-Teller model. Normal mode analysis was used to identify those protein residues that act as the primary "doorway" modes in the vibrational relaxation of the oscillator. Although the CO relaxation rates in both the epsilon- and delta-tautomers are similar in magnitude, the simulations predict that the vibrational relaxation of the CO is faster in the delta-tautomer with the distal His playing an important role in the energy relaxation mechanism. Time-resolved mid-IR absorbance measurements were performed on photolyzed carbonmonoxy hemoglobin (Hb(13)CO). From these measurements, a T(1) time of 600 +/- 150 ps was determined. The simulation and experimental estimates are compared and discussed.

Structural correlates for enhanced stability in the E2 DNA-binding domain from bovine papillomavirus

Papillomaviral E2 proteins participate in viral DNA replication and transcriptional regulation. We have solved the solution structure of the DNA-binding domain of the E2 protein from bovine papillomavirus (BPV-1). The structure calculation used 2222 distance and 158 dihedral angle restraints for the homodimer (202 residues in total), which were derived from homonuclear and heteronuclear multidimensional nuclear magnetic resonance (NMR) spectroscopic data. The root-mean-square deviation for structured regions of the monomer when superimposed to the average is 0.73 +/- 0.10 A for backbone atoms and 1.42 +/- 0.16 A for heavy atoms. The 101 residue construct used in this study (residues 310-410) is about 4.5 kcal/mol more stable than a minimal domain comprising the C-terminal 85 amino acid residues (residues 326-410). The structure of the core domain contained within BPV-1 E2 is similar to the corresponding regions of other papilloma viral E2 proteins. Here, however, the extra N-terminal 16 residues form a flap that covers a cavity at the dimer interface and play a role in DNA binding. Interactions between residues in the N-terminal extension and the core domain correlate with the greater stability of the longer form of the protein relative to the minimal domain.

Structural determinants of trypsin affinity and specificity for cationic inhibitors

The binding free energies of four inhibitors to bovine beta-trypsin are calculated. The inhibitors use either ornithine, lysine, or arginine to bind to the S1 specificity site. The electrostatic contribution to binding free energy is calculated by solving the finite difference Poisson-Boltzmann equation, the contribution of nonpolar interactions is calculated using a free energy-surface area relationship and the loss of conformational entropy is estimated both for trypsin and ligand side chains. Binding free energy values are of a reasonable magnitude and the relative affinity of the four inhibitors for trypsin is correctly predicted. Electrostatic interactions are found to oppose binding in all cases. However, in the case of ornithine- and lysine-based inhibitors, the salt bridge formed between their charged group and the partially buried carboxylate of Asp189 is found to stabilize the complex. Our analysis reveals how the molecular architecture of the trypsin binding site results in highly specific recognition of substrates and inhibitors. Specifically, partially burying Asp189 in the inhibitor-free enzyme decreases the penalty for desolvation of this group upon complexation. Water molecules trapped in the binding interface further stabilize the buried ion pair, resulting in a favorable electrostatic contribution of the ion pair formed with ornithine and lysine side chains. Moreover, all side chains that form the trypsin specificity site are partially buried, and hence, relatively immobile in the inhibitor-free state, thus reducing the entropic cost of complexation. The implications of the results for the general problem of recognition and binding are considered. A novel finding in this regard is that like charged molecules can have electrostatic contributions to binding that are more favorable than oppositely charged molecules due to enhanced interactions with the solvent in the highly charged complex that is formed.

Modelling the catalytic reaction in human aldose reductase

Aldose reductase (ALR2) has received considerable attention due to its possible link to long-term diabetic complications. Although crystal structures and kinetic data reveal important aspects of the reaction mechanism, details of the catalytic step are still unclear. In this paper a computer simulation study is presented that utilizes the hybrid quantum mechanical and molecular mechanical (QM-MM) potential to elucidate the nature of the hydride and proton transfer steps in the reduction of D-glyceraldehyde by ALR2. Several reaction pathways were investigated in two models with either Tyr48 or protonated His110+ acting as the potential proton donor in the active site. Calculations show that the substrate binds to ALR2 through hydrogen bonds in an orientation that facilitates the stereospecific catalytic step in both models. It is established that in the case that His110 is present in the protonated form in the native complex, it is the energetically favored proton donor compared with Tyr48 in the active pocket with neutral His110. The reaction mechanisms in the different models are discussed based on structural and energetic considerations.

Active site dynamics of the HhaI methyltransferase: insights from computer simulation

A molecular dynamics study was performed on the DNA methyltransferase M. Hha I in a ternary complex with DNA and AdoMet in solution. Methylation involves addition of the Cys81 sulfhydryl anion to the 6-position of Cyt18, followed by a nucleophilic attack of the resultant carbanion at C5 on the AdoMet methyl group. It was found in this simulation that the distances between the sulfhydryl group (SG) of Cys81 to the C6 of Cyt18 (SG-C6) and methyl carbon (CH3) of AdoMet to the C5 of cytosine (CH3-C5) are dependent on the dihedral angle chi (O4'-C1'-N1-C2) of the nucleotide. When the chi angle of Cyt18 is low (< -80 degrees), the SG-C6 and CH3-C5 distances are large. A high chi angle (> -80 degrees) for the target cytosine residue reduces the distances for both SG-C6 and CH3-C5, and the angles formed between the cytosine ring and AdoMet correspond well to values for the transition state structures formed during methylation of cytosine from ab initio calculations. Two possible proton sources for protonation of N3 of the cytosine residue upon formation of the covalent intermediate were found in the simulation. The protonated amine group of AdoMet could provide a proton via a water bridge, or Arg163 could also be the source of the proton for N3 via a water bridge. The simulation provides insights into how the H5 of cytosine could go from the active site into solvent. Conserved residues Asn304 and Gln82 stabilize a water network within the active site of M. Hha I which provides a route for H5 to diffuse into bulk solvent. An initially distant water molecule was able to diffuse into the active site of the enzyme and replace a position of a crystallographic water molecule in close proximity to the C5 of cytosine. The movement of this water molecule showed that a channel exists between Gln82 and the AdoMet in M. Hha I which allows both water and protons to easily gain access to the active site of the enzyme.

The influence of helix morphology on co-operative polyamide backbone conformational flexibility in peptide nucleic acid complexes

A systematic analysis of peptide nucleic acid (PNA) complexes deposited in the Protein Data Bank has been carried out using a set of contiguous atom torsion angle definitions. The analysis is complemented by molecular mechanics adiabatic potential energy calculations on hybrid PNA-nucleic acid model systems. Hitherto unobserved correlations in the values of the (alpha and epsilon) dihedral angles flanking the backbone secondary amide bond are found. This dihedral coupling forms the basis of a PNA backbone conformation classification scheme. Six conformations are thus characterised in experimental structures. Helix morphology is found to exert a significant influence on backbone conformation and flexibility: Watson-Crick PNA strands in complexes with DNA and RNA, that possess A-like base-pair stacking, adopt backbone conformations distinct from those in PNA.DNA-PNA triplex and PNA-PNA duplex P-helix forms. Solvation effects on Watson-Crick PNA backbone conformation in heterotriplexes are discussed and the possible involvement of inter-conformational transitions and dihedral angle uncoupling in asymmetric heteroduplex base-pair breathing is suggested.

Free energy calculations and molecular dynamics simulations of wild-type and variants of the DNA-EcoRI complex

Molecular dynamics simulations and free energy calculations of the wild-type EcoRI-DNA complex and several variants have been performed in aqueous solvent. In general, he theoretical estimations of the free energy differences (DeltaDeltaA) qualitatively agree well with the corresponding experimental data. The modifications which were experimentally found unfavorable compared to the wild-type complex were also found to be so in theoretical estimates. The mutant where the amino group of the base Ade(6) was replaced by a hydrogen atom eliminating one H-bond between the DNA and the protein, was experimentally found to be more stable than the wild-type complex. It was speculated that the modification also caused a structural relaxation in the DNA making DeltaDeltaA favorable. Our theoretical estimate yields a positive DeltaDeltaA in this case, but the difference is small, and no significant local structural relaxation was observed. The major H-bonds between the DNA and the protein in the wild-type complex are found to be maintained in the different mutants although the specific and non-specific interaction energies between the interacting the DNA bases and the protein residues are different in different mutants. The interaction pattern of the other nearby nucleotides are significantly influenced by each modification. Thus, the alteration of the non-specific interactions may also play an indirect role in determining the specificity of the complex. The interaction of the Gua(4) of the DNA with the protein is found to be most sensitive to any alteration in the recognition site. Because Gua(4) is the nucleotide closest to the scissile bond, this extra sensitivity seems to play an important role in altering the functional activity of the complex.