BANMOKI: a searchable database of homology-based 3D models and their electrostatic properties of five bacterial nucleoside monophosphate kinase families

The nucleoside monophosphate kinases (NMPK) are important enzymes that control the ratio of mono- and di-phosphate nucleosides and participate in gene regulation and signal transduction in the cell. However, despite their importance only several 3D structures were experimentally determined in contrast to the wealth of sequences available for each of the NMPK families. To fill this gap we present a Web-based database containing structural models for all proteins of the five bacterial nucleoside monophosphate kinase (bNMPK) families. The models were computed by means of homology-based approach using a few experimentally determined bNMPK structures. The database also contains pK(a) values and their components calculated for the homology-based 3D models, which is a unique feature of the database. The BActerial Nucleoside MOnophosphate KInases (BANMOKI) database is freely accessible (http://www.ces.clemson.edu/compbio/banmoki) and offers an easy user-friendly interface for browsing, searching and downloading content of the database. The users can investigate, using the searching tools of the database, the properties of the bNMP kinases in respect to sequence composition, electrostatic interactions and structural differences.     

Acid-base equilibria in rhodopsin: dependence of the protonation state of glu134 on its environment

Glutamic acid E134 in rhodopsin is part of a highly conserved triad, D(E)RY, located near the cytoplasmic lipid/water interface in transmembrane helix 3 of G protein-coupled receptors (GPCRs). A large body of experimental evidence suggests that the protonation of E134 plays a role in the mechanism of activation of rhodopsin and other GPCRs as well. For E134 to change its protonation state, its pK(a) value must shift from values below physiological pH to higher values. Because of the proximity of the triad to the lipid/water interface, it was hypothesized that a change in solvent around E134 from water to lipid could induce such a shift in pK(a). To test this hypothesis, the pK(a) values of the titratable amino acid residues in rhodopsin have been calculated and the change in solvent around E134 was modeled by shifting the position of the lipid/water interface. The approach used to carry out the pK(a) calculations takes into account the partial immersion of transmembrane proteins in lipid. Qualitative experimental evidence is available for several residues regarding their likely protonation state in rhodopsin at or near physiological pH. Comparison of the calculated pK(a) values with these experimental findings shows good agreement between the two. Notably, glutamic acids E122 and E181 were found to be protonated. The pK(a) values were then calculated for a range of lipid/water interface positions. Although the surrounding solvent of several titratable residues changed from water to lipid in this range, leading to pK(a) shifts in most cases, only for E134 would the shift lead to a change in protonation state at physiological pH. Thus, our results show that the protonation state of E134 is particularly sensitive to its environment. This sensitivity together with the location of E134 near the actual position of the lipid/water interface could be a strategic element in the mechanism of activation of rhodopsin.     

Optimizing the hydrogen-bond network in Poisson-Boltzmann equation-based pK(a) calculations

pK(a) calculation methods that are based on finite difference solutions to the Poisson-Boltzmann equation (FDPB) require that energy calculations be performed for a large number of different protonation states of the protein. Normally, the differences between these protonation states are modeled by changing the charges on a few atoms, sometimes the differences are modeled by adding or removing hydrogens, and in a few cases the positions of these hydrogens are optimized locally. We present an FDPB-based pK(a) calculation method in which the hydrogen-bond network is globally optimized for every single protonation state used. This global optimization gives a significant improvement in the accuracy of calculated pK(a) values, especially for buried residues. It is also shown that large errors in calculated pK(a) values are often due to structural artifacts induced by crystal packing. Optimization of the force fields and parameters used in pK(a) calculations should therefore be performed with X-ray structures that are corrected for crystal artifacts.     

Calculation of pK(a) in proteins with the microenvironment modulated-screened coulomb potential

The MM-SCP has been applied to predict pK(a) values of titratable residues in wild type and mutants of staphylococcal nuclease (SNase). The calculations were based on crystal structures made available by the Garcia-Moreno Laboratory. In the mutants, mostly deeply buried hydrophobic residues were replaced with ionizable residues, and thus their pK(a) values could be measured and calculated using various methods. The data set used here consisted of a set of WT SNase for which His pK(a) at several ionic strengths had been measured, a set of mutants for which measured pK(a) were available and a set of 11 mutants for which the measured pK(a) were not known at the time of calculation. For this latter set, blind predictions were submitted to the protein pK(a) cooperative, 2009 workshop at Telluride, where the results of the blind predictions were discussed (the RMSD of the submitted set was 1.10 pH units). The calculations on the structures with known pK(a) indicated that in addition to weaknesses of the method, structural issues were observed that led to larger errors (>1) in pK(a) predictions. For example, different crystallographic conditions or steric clashes can lead to differences in the local environment around the titratable residue, which can produce large differences in the calculated pK(a) . To gain further insight into the reliability of the MM-SCP, pK(a) of an extended set of 54 proteins belonging to several structural classes were carried out. Here some initial results from this study are reported to help place the SNase results in the appropriate context.     

Predicting extreme pKa shifts in staphylococcal nuclease mutants with constant pH molecular dynamics

Accurate computational methods of determining protein and nucleic acid pK(a) values are vital to understanding pH-dependent processes in biological systems. In this article, we use the recently developed method constant pH molecular dynamics (CPHMD) to explore the calculation of highly perturbed pK(a) values in variants of staphylococcal nuclease (SNase). Simulations were performed using the replica exchange (REX) protocol for improved conformational sampling with eight temperature windows, and yielded converged proton populations in a total sampling time of 4 ns. Our REX-CPHMD simulations resulted in calculated pK(a) values with an average unsigned error (AUE) of 0.75 pK units for the acidic residues in Δ + PHS, a hyperstable variant of SNase. For highly pK(a)-perturbed SNase mutants with known crystal structures, our calculations yielded an AUE of 1.5 pK units and for those mutants based on modeled structures an AUE of 1.4 pK units was found. Although a systematic underestimate of pK shifts was observed in most of the cases for the highly perturbed pK mutants, correlations between conformational rearrangement and plasticity associated with the mutation and error in pK(a) prediction was not evident in the data. This study further extends the scope of electrostatic environments explored using the REX-CPHMD methodology and suggests that it is a reliable tool for rapidly characterizing ionizable amino acids within proteins even when modeled structures are employed.     

pKa's of ionizable groups in proteins: atomic detail from a continuum electrostatic model

A macroscopic electrostatic model is used to calculate the pKa values of the titratable groups in lysozyme. The model makes use of detailed structural information and treats solvation self-energies and interactions arising from permanent partial charges and titratable charges. Both the tetragonal and triclinic crystal structures are analyzed. Half of the experimentally observed pKa shifts (11 out of 21) are well reproduced by calculations for both structures; this includes the unusually high pKa of Glu 35 in the active site. For more than half the pKa's (13 out of 21), there is a large difference (1-3.3 pK units) between the results from the two structures. Many of these correspond to the titrating groups for which the calculations are in error. Since for an ionic strength of 0.1 M the Debye screening between titratable groups leads to a very high effective dielectric constant (the average value for all pairs of titrating groups is approximately 900), near-neighbor interactions dominate the pKa perturbations. Thus, the pKa values are very sensitive to the details of the local protein conformation, and it is likely that side-chain mobility has an important role in determining the observed pKa shifts.     

pK(a) calculations for class C beta-lactamases: the role of Tyr-150

The Poisson-Boltzmann method was used to compute the pK(a) values of titratable residues in a set of class C beta-lactamases. In these calculations, the pK(a) of the phenolic group of residue Tyr150 is the only one to stand out with an abnormally low value of 8.3, more than one pK(a) unit lower than the measured reference value for tyrosine in solution. Other important residues of the catalytic pocket, such as the conserved Lys67, Lys315, His314, and Glu272 (hydrogen-bonded to the ammonium group of Lys315), display normal protonation states at neutral pH. pK(a) values were also computed in catalytically impaired beta-lactamase mutants. Comparisons between the relative k(cat) values and the Tyr150 pK(a) value in these mutants revealed a striking correlation. In active enzymes, this pK(a) value is always lower than the solution reference value while it is close to normal in inactive enzymes. These results thus support the hypothesis that the phenolate form of Tyr150 is responsible for the activation of the nucleophilic serine. The possible roles of Lys67 and Lys315 during catalysis are also discussed.     

Comparison of Composer and ORCHESTRAR

Although the number of known protein structures is increasing, the number of protein sequences without determined structures is still much larger. Three-dimensional (3D) protein structure information helps in the understanding of functional mechanisms, but solving structures by X-ray crystallography or NMR is often a lengthy and difficult process. A relatively fast way of determining a protein's 3D structure is to construct a computer model using homologous sequence and structure information. Much work has gone into algorithms that comprise the ORCHESTRAR homology modeling program in the SYBYL software package. This novel homology modeling tool combines algorithms for modeling conserved cores, variable regions, and side chains. The paradigm of using existing knowledge from multiple templates and the underlying protein environment knowledgebase is used in all of these algorithms, and will become even more powerful as the number of experimentally derived protein structures increases. To determine how ORCHESTRAR compares to Composer (a broadly used, but an older tool), homology models of 18 proteins were constructed using each program so that a detailed comparison of each step in the modeling process could be carried out. Proteins modeled include kinases, dihydrofolate reductase, HIV protease, and factor Xa. In almost all cases ORCHESTRAR produces models with lower root-mean-squared deviation (RMSD) values when compared with structures determined by X-ray crystallography or NMR. Moreover, ORCHESTRAR produced a homology model for three target sequences where Composer failed to produce any. Data for RMSD comparisons between structurally conserved cores, structurally variable regions, side-chain conformations are presented, as well as analyses of active site and protein-protein interface configurations.     

Homology modelling and molecular dynamics simulations: comparative studies of human aquaporin-1

The structures of the mammalian water transport protein Aqp1 and of its bacterial homologue GlpF enables us to test whether homology models can be used to explore relationships between structure, dynamics and function in mammalian transport proteins. Molecular dynamics simulations (totalling almost 40 ns) were performed starting from: the X-ray structure of Aqp1; a homology model of Aqp1 based on the GlpF structure; and intermediate resolution structures of Aqp1 derived from electron microscopy. Comparisons of protein RMSDs vs. time suggest that the homology models are of comparable conformational stability to the X-ray structure, whereas the intermediate resolution structures exhibit significant conformation drift. For simulations based on the X-ray structure and on homology models, the flexibility profile vs. residue number correlates well with the crystallographic B-values for each residue. In the simulations based on intermediate resolution structures, mobility of the highly conserved NPA loops is substantially higher than in the simulations based on the X-ray structure or the homology models. Pore radius profiles remained relatively constant in the X-ray and homology model simulations but showed substantial fluctuations (reflecting the higher NPA loop mobility) in the intermediate resolution simulations. The orientation of the dipoles of water molecules within the pore is of key importance in maintaining low proton permeability through Aqp1. This property seems to be quite robust to the starting model used in the simulation. These simulations suggest that homology models based on bacterial homologues may be used to derive functionally relevant information on the structural dynamics of mammalian transport proteins.     

Energetics of a possible proton exit pathway for water oxidation in photosystem II

The crystal structure of photosystem II (PSII) at 3.0-A resolution suggests that titratable residues on the lumenal side of D1/D2 and PsbO form a polar channel, which might serve as a proton exit pathway associated with water oxidation on the Mn-cluster. With full account of protein environment, we calculated the pK(a) of these residues by solving the linearized Poisson-Boltzmann equation. Along the prospective proton channel, the calculated pK(a) of titratable residues (namely via D1-Asp61, D1-Glu65, D2-Glu312, D2-Lys317 D1-Asp59, D1-Arg64, PsbO-Arg152, and PsbO-Asp224) monotonically increase from the Mn-cluster to the lumenal bulk side. We suggest that these residues form the exit pathway guiding protons, which are released at the Mn-cluster as a product of water oxidation, in an exergonic process out of PSII. Upon the S2 to S3 transition, CP43-Arg357 showed a dramatic deprotonation of ca. one H(+), suggesting that this residue is coupled to the redox states of the Mn-cluster and the tyrosine Y(Z). The calculated pK(a) values of 4.2-4.4 for D2-Glu312 and those of approximately 8-10.9 for D1-Asp59 and D1-Arg64 are indicative of the experimentally determined pK(a) values for inhibition of S-state transitions. Upon removal of the atomic coordinates of PsbO, the pK(a) of these residues are dramatically affected, indicating a significant role of PsbO in tuning the pK(a) of those residues in the proton exit pathway.     

Electrostatic pKa computations in proteins: role of internal cavities

The solvent accessible surface area (SASA) algorithm is conventionally used to characterize protein surfaces in electrostatic energy computations of proteins. Unfortunately, it often fails to find narrow cavities inside a protein. As a consequence pK(a) computations based on this algorithm perform badly. In this study a new cavity-algorithm is introduced, which solves this problem and provides improved pK(a) values. The procedure is applied to 20 pK(a) values of titratable groups introduced as point mutations in SNase variants, where crystal structures are available. The computations of these pK(a)s are particular challenging, since they are placed in a rather hydrophobic environment. For nine mutants, where the titratable residue is in contact with a large cavity, the RMSD(pKa) between computed and measured pK(a) values is 2.04, which is a considerable improvement as compared to the original results obtained with Karlsberg(+) (http://agknapp.chemie.fu-berlin.de/karlsberg/) that yielded an RMSD(pKa) of 8.8. However, for 11 titratable residues the agreement with experiments remains poor (RMSD(pKa) = 6.01). Considering 15 pK(a)s of SNase, which are in a more conventional less hydrophobic protein environment, the RMSD(pKa) is 2.1 using the SASA-algorithm and 1.7 using the new cavity-algorithm. The agreement is reasonable but less good than what one would expect from the general performance of Karlsberg(+) indicating that SNase belongs to the more difficult proteins with respect to pK(a) computations. We discuss the possible reasons for the remaining discrepancies between computed and measured pK(a)s.     

The ionization state of D37 in E. coli porin OmpF and the nature of conductance fluctuations in D37 mutants

Understanding the flow of ions through E. coli porin outer membrane protein F (OmpF) requires knowledge of the charge state of all titratable residues located along the permeation pathway. Earlier theoretical studies proved successful in the calculation of the pK values of most residues. The (apparent) pK of Asp37 (D37), on the other hand, appeared rather sensitive to the (unknown) protein dielectric used. We addressed the protonation state of D37 experimentally by replacing D37 with a (neutral) valine. This D37V mutant expressed reduced cation selectivity, in agreement with the view that D37 in wild-type (WT) OmpF is fully ionized, i.e., deprotonated. The introduction of a (positively charged) arginine at position 37 evoked current fluctuations. Similar behavior was observed in the D37K mutant and the cysteine mutants D37C-MTSEA and D37C-MTSET. Nontitratable [2-(trimethylammonium)ethyl]-methanethiosulfonate (MTSET) carries a permanent and pH-independent charge of 1e, implying that the fluctuations of the D37C-MTSET mutant do not represent (de)protonation reactions of MTSET. We therefore conclude that these fluctuations reflect transitions between conformational substates evoked by structural instabilities due to the positive charge at that particular position in the pore lumen. Based on the similarities between D37C-MTSET fluctuations and those seen in the other mutants, notably D37K, the underlying mechanism of these fluctuations may be (essentially) the same in all four mutants studied.     

Homology modeling of the human 5-HT1A, 5-HT2A, D1, and D2 receptors: model refinement with molecular dynamics simulations and docking evaluation

5-HT(1A) serotonin and D1 dopamine receptor agonists have been postulated to be able to improve negative and cognitive impairment symptoms of schizophrenia, while partial agonists and antagonists of the D2 and 5-HT(2A) receptors have been reported to be effective in reducing positive symptoms. There is therefore a need for well-defined homology models for the design of more selective antipsychotic agents, since no three-dimensional (3D) crystal structures of these receptors are currently available. In this study, homology models were built based on the high-resolution crystal structure of the β(2)-adrenergic receptor (2RH1) and further refined via molecular dynamics simulations in a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) lipid bilayer system with the GROMOS96 53A6 united atom force field. Docking evaluations with representative agonists and antagonists using AutoDock 4.2 revealed binding modes in agreement with experimentally determined site-directed mutagenesis data and significant correlations between the computed and experimental pK (i) values. The models are also able to distinguish between antipsychotic agents with different selectivities and binding affinities for the four receptors, as well as to differentiate active compounds from decoys. Hence, these human 5-HT(1A), 5-HT(2A), D1 and D2 receptor homology models are capable of predicting the activities of novel ligands, and can be used as 3D templates for antipsychotic drug design and discovery.     

Measuring the successes and deficiencies of constant pH molecular dynamics: a blind prediction study

A constant pH molecular dynamics method has been used in the blind prediction of pK(a) values of titratable residues in wild type and mutated structures of the Staphylococcal nuclease (SNase) protein. The predicted values have been subsequently compared to experimental values provided by the laboratory of García-Moreno. CpHMD performs well in predicting the pK(a) of solvent-exposed residues. For residues in the protein interior, the CpHMD method encounters some difficulties in reaching convergence and predicting the pK(a) values for residues having strong interactions with neighboring residues. These results show the need to accurately and sufficiently sample conformational space in order to obtain pK(a) values consistent with experimental results.     

The K-pathway revisited: a computational study on cytochrome c oxidase

Cytochrome c oxidase contains two established proton-conducting structures, the D- and K-pathways. The role of the K-pathway appears to be to conduct the first two protons to be used in water formation, which are taken up on reduction of the oxidized enzyme. Previous computational work has suggested that Lys(I)-319 is neutral over a large pH range and in various redox states. We have constructed oxidase models in different redox states using quantum-chemically derived charge parameters for the redox metal centers. The protonation behaviour of titratable sites in the two-subunit enzyme was defined by continuum electrostatics. The calculations reported here show substantial protonation of Lys(I)-319 at neutral pH once the stable X-ray crystallographic water molecule found immediately next to it is treated explicitly. The immediate structure of the Lys(I)-319 environment is independent of redox state, but the pK(a) value of this residue changes with the redox state of the binuclear heme a3/Cu(B) site whenever that change is electrically uncompensated. Lys(I)-319 is also found to interact electrostatically with the conserved residue Glu(II)-62 in subunit II. These results are discussed in relation to the role of the K-pathway in oxidase function.     

Molecular dynamics effects on protein electrostatics

Electrostatic calculations have been carried out on a number of structural conformers of tuna cytochrome c. Conformers were generated using molecular dynamics simulations with a range of solvent simulating, macroscopic dielectric formalisms, and one solvent model that explicitly included solvent water molecules. Structures generated using the lowest dielectric models were relatively tight, with side chains collapsed on the surface, while those from the higher dielectric models had more internal and external fluidity, with surface side chains exploring a fuller range of conformational space. The average structure generated with the explicitly solvated model corresponded most closely with the crystal structure. Individual pK values, overall titration curves, and electrostatic potential surfaces were calculated for average structures and structures along each simulation. Differences between structural conformers within each simulation give rise to substantial changes in calculated local electrostatic interactions, resulting in pK value fluctuations for individual sites in the protein that vary by 0.3-2.0 pK units from the calculated time average. These variations are due to the thermal side chain reorientations that produce fluctuations in charge site separations. Properties like overall titration curves and pH dependent stability are not as sensitive to side chain fluctuations within a simulation, but there are substantial effects between simulations due to marked differences in average side chain behavior. These findings underscore the importance of proper dielectric formalism in molecular dynamics simulations when used to generate alternate solution structures from a crystal structure, and suggest that conformers significantly removed from the average structure have altered electrostatic properties that may prove important in episodic protein properties such as catalysis.     

pKa predictions with a coupled finite difference Poisson-Boltzmann and Debye-Hückel method

Modeling charge interactions is important for understanding many aspects of biological structure and function, and continuum methods such as Finite Difference Poisson-Boltzmann (FDPB) are commonly employed. Calculations of pH-dependence have identified separate populations; surface groups that can be modeled with a simple Debye-Hückel (DH) model, and buried groups, with stronger resultant interactions that are dependent on detailed conformation. This observation led to the development of a combined FDPB and DH method for pK(a) prediction (termed FD/DH). This study reports application of this method to ionizable groups, including engineered buried charges, in staphylococcal nuclease. The data had been made available to interested research groups before publication of mutant structures and/or pK(a) values. Overall, FD/DH calculations perform as intended with low ΔpK(a) values for surface groups (RMSD between predicted and experimental pK(a) values of 0.74), and much larger ΔpK(a) values for the engineered internal groups, with RMSD = 1.64, where mutant structures were known and RMSD = 1.80, where they were modeled. The weaker resultant interactions of the surface groups are determined mostly by charge-charge interactions. For the buried groups, R(2) for correlation between predicted and measured ΔpK(a) values is 0.74, despite the high RMSDs. Charge-charge interactions are much less important, with the R(2) value for buried group ΔpK(a) values increasing to 0.80 when the term describing charge desolvation alone is used. Engineered charge burial delivers a relatively uniform, nonspecific effect, in terms of pK(a) . How the protein environment adapts in atomic detail to deliver this resultant effect is still an open question.     

NMR structure of the conserved hypothetical protein TM0487 from Thermotoga maritima: implications for 216 homologous DUF59 proteins

The NMR structure of the conserved hypothetical protein TM0487 from Thermotoga maritima represents an alpha/beta-topology formed by the regular secondary structures alpha1-beta1-beta2-alpha2-beta3-beta4-alpha3- beta5-3(10)-alpha4, with a small anti-parallel beta-sheet of beta-strands 1 and 2, and a mixed parallel/anti-parallel beta-sheet of beta-strands 3-5. Similar folds have previously been observed in other proteins, with amino acid sequence identity as low as 3% and a variety of different functions. There are also 216 sequence homologs of TM0487, which all have the signature sequence of domains of unknown function 59 (DUF59), for which no three-dimensional structures have as yet been reported. The TM0487 structure thus presents a platform for homology modeling of this large group of DUF59 proteins. Conserved among most of the DUF59s are 13 hydrophobic residues, which are clustered in the core of TM0487. A putative active site of TM0487 consisting of residues D20, E22, L23, T51, T52, and C55 is conserved in 98 of the 216 DUF59 sequences. Asp20 is buried within the proposed active site without any compensating positive charge, which suggests that its pK(a) value may be perturbed. Furthermore, the DUF59 family includes ORFs that are part of a conserved chromosomal group of proteins predicted to be involved in Fe-S cluster metabolism.     

The pH dependence of red cell membrane transport of titratable anions studied by NMR spectroscopy

The effects of varying extracellular pH on the rates of uptake of titratable anions by human erythrocytes under conditions of constant intracellular pH have been determined for a series of highly related anions, the phosphate "analogs." These compounds are simply substituted phosphorus oxyacids, differing in the number and acidity of titratable protons: phosphate (HPO4(2-), pKa 6.8); phosphite (HPO3(2-), pKa 6.4); hypophosphite (H2PO2-); methylphosphonate ((CH3)PO3(2-), pKa 7.4); dimethylphosphinate ((CH3)2PO2-); fluorophosphate [PO3F2-, pKa 4.7); and thiophosphate (HSPO3(2-), pKa 5.5). Suspensions of intact, Cl(-)-loaded erythrocytes (intracellular pH, 7.2) were incubated at 37 degrees C in isotonic buffers (pH 4-8) containing 60 mM phosphate analog for specified time intervals, whereupon influx was halted by the addition of 1 mM 4-acetamido-4'-isothiocyanatostilbene-2,2'-disulfonic acid (SITS), an inhibitor of anion exchange. The intracellular anion concentrations were determined from 31P or 19F nuclear magnetic resonance spectra from the erythrocyte suspensions. The influx rates for the titratable phosphate analogs exhibited bimodal pH dependence, reaching maximal levels at pH values that increased with increasing anion pK. This pH-dependent behavior is consistent with a transport channel that contains a titratable regulatory site which interacts with the translocated anion. Based upon the Henderson-Hasselbalch equation, the probability that a titratable anion will have an electric charge of equal magnitude to that of the titratable carrier is highest at a pH value exactly midway between the pK of the regulatory site and that of the anion. The pH maxima observed for the phosphate analogs indicate a pK for this site of 5.5 at 37 degrees C. Intracellular pH changes associated with influx indicated that transport of the "fast" anion phosphite is largely in monoionized form. Intracellular pH changes associated with transport of slow anions were predominantly determined by partial ionic equilibrium effects and did not indicate the ionization state of the transported anion.