Kohn-Sham Density Functional Calculations Reveal Proton Wires in the Enolization and Carboxylase Reactions Catalyzed by Rubisco

Ribulose 1,5-bisphosphate (RuBP) carboxylase-oxygenase (Rubisco) plays a fundamental role in the carbon cycle by fixing the atmospheric CO₂ used in photosynthesis. Rubisco is all the more remarkable because it must catalyze some difficult multistep reaction chemistry involving proton transfers within the one active site. In the present study, we have used Kohn-Sham density functional theory at the B3LYP/6-31G* level with basis set superposition error and dispersion corrections (B3LYP-gCP-D3) to examine the possibility that the proton transfers can take place through molecular wires (including active-site water molecules) via the classical Grotthuss proton-shuttle mechanism. The results support an essential role for water molecules found in the crystal structures of Rubisco complexes as facilitators of proton transport in all the rate-limiting (catalytic) reaction steps through a network of short proton wires within the Rubisco active site. We suggest that completion of the initial product turnover (cycle) requires two excess protons produced in the initial carbamylation that is required for Rubisco activation. By use of proton wires, a large number of reaction steps may be accommodated within a single active site without necessitating the input of excessive conformational strain energy arising from the movement of residue side chains into positions where direct protonation of substrates can occur. The involvement of the identified types of proton wires in the kinetic mechanism is capable of providing a unique explanation for various experimental observations, including deuterium isotope effects and the results of site-directed mutagenesis experiments, and may thus provide a realistic solution to the problem of Rubisco's challenging chemistry.

Mechanism of Oxygenase-Pathway Reactions Catalyzed by Rubisco from Large-Scale Kohn-Sham Density Functional Calculations

Ribulose 1,5-bisphosphate carboxylase-oxygenase (Rubisco) is the primary carbon-fixing enzyme in photosynthesis, fixing CO₂ to a 5-carbon sugar, RuBP, in a series of five reactions. However, it also catalyzes an oxygenase reaction by O₂ addition to the same enolized RuBP substrate in an analogous reaction series in the same active site, producing a waste product and loss of photosynthetic efficiency. Starting from RuBP, the reactions are enolization to the enediolate form, addition of CO₂ or O₂ to form the carboxy or peroxo adduct, hydration to form a gemdiolate, scission of the C2-C3 bond of the original RuBP, and stereospecific or nonstereospecific protonation to form two molecules of the 3-carbon PGA product, or one molecule of PGA, one of 2-carbon PG (waste product), and one water molecule. Reducing the loss of efficiency from the oxygenase reaction is an attractive means to increase crop productivity. However, lack of understanding of key aspects of the catalytic mechanisms for both the carboxylase and oxygenase reactions, particularly those involving proton exchanges and roles of water molecules, has stymied efforts at re-engineering Rubisco to reduce losses from the oxygenation reaction. As the stable form of molecular oxygen is the triplet biradical state (³O₂), its reaction with near-universal singlet-state molecules is formally spin forbidden. Although in oxygenase enzymes, ³O₂ activation is usually achieved by one-electron transfers using transition-metal ions or organic cofactors, recently, cofactor-less oxygenases in which the substrate itself is the source of the electron for ³O₂ activation have been identified, but in all such cases an aromatic ring stabilizes the substrate's negative charge. Here we present the first large-scale Kohn-Sham density functional theory study of the reaction mechanism of the Rubisco oxygenase pathway. First, we show that the enediolate substrate complexed to Mg²⁺ and its ligands extends the region for charge delocalization and stabilization of its negative charge to allow formation of a caged biradical enediolate-O₂ complex. Thus, Rubisco is a unique type of oxygenase without precedent in the literature. Second, for the O₂ addition to proceed to the singlet peroxo-adduct intermediate, the system must undergo an intersystem crossing. We found that the presence of protonated LYS334 is required to stabilize this intermediate and that both factors (strongly stabilized anion and protonated LYS334) facilitate a barrier-less activation of ³O₂. This finding supports our recent proposal that deoxygenation, that is, reversal of gas binding, is possible. Third, as neither CO₂ nor O₂ binds to the enzyme, our findings support the proposal from our recent carboxylase study that the observed K_C or K_O (Michaelis-Menten constants) in the steady-state kinetics reflect the respective adducts, carboxy or peroxo. Fourth, after computing hydration pathways with water addition both syn and anti to C3, we found, in contrast to the results of our carboxylation study indicating anti addition, that in the oxygenation reaction only syn-hydration is capable of producing a stable gemdiolate that facilitates the rate-limiting C2-C3 bond scission to final products. Fifth, we propose that an excess proton we previously found was required in the carboxylation reaction for activating the C2-C3 bond scission is utilized in the oxygenation reaction for the required elimination of a water molecule. In summary, despite its oxygenase handicap, Rubisco's success in directing 75% of its substrate through the carboxylation pathway can be considered impressively effective. Although native C3 Rubiscos are in a fix with unwanted activity of ³O₂ hampering its primary carboxylase function, mechanistic differences presented here with findings in our recent carboxylase study for both the gas-addition and subsequent reactions provide some clues as to how creative Rubisco re-engineering may offer a solution to reducing the oxygenase activity.

Revised mechanism of carboxylation of ribulose-1,5-biphosphate by rubisco from large scale quantum chemical calculations

Here, we describe a computational approach for studying enzymes that catalyze complex multi-step reactions and apply it to Ribulose 1,5-bisphosphate carboxylase-oxygenase (Rubisco), the enzyme that fixes atmospheric carbon dioxide within photosynthesis. In the 5-step carboxylase reaction, the substrate Ribulose-1,5-bisphosphate (RuBP) first binds Rubisco and undergoes enolization before binding the second substrate, CO₂ . Hydration of the RuBP.CO₂ complex is followed by CC bond scission and stereospecific protonation. However, details of the roles and protonation states of active-site residues, and sources of protons and water, remain highly speculative. Large-scale computations on active-site models provide a means to better understand this complex chemical mechanism. The computational protocol comprises a combination of hybrid semi-empirical quantum mechanics and molecular mechanics within constrained molecular dynamics simulations, together with constrained gradient minimization calculations using density functional theory. Alternative pathways for hydration of the RuBP.CO₂ complex and associated active-site protonation networks and proton and water sources were investigated. The main findings from analysis of the resulting energetics advocate major revision to existing mechanisms such that: hydration takes place anti to the CO₂ ; both hydration and CC bond scission require early protonation of CO₂ in the RuBP.CO₂ complex; CC bond scission and stereospecific protonation reactions are concerted and, effectively, there is only one stable intermediate, the C3-gemdiolate complex. Our main conclusions for interpreting enzyme kinetic results are that the gemdiolate may represent the elusive Michaelis-Menten-like complex corresponding to the empirical K_m (=K_c ) with turnover to product via bond scission concerted with stereospecific protonation consistent with the observed catalytic rate. © 2018 Wiley Periodicals, Inc.

Directions for Optimization of Photosynthetic Carbon Fixation: RuBisCO's Efficiency May Not Be So Constrained After All

The ubiquitous enzyme Ribulose 1,5-bisphosphate carboxylase-oxygenase (RuBisCO) fixes atmospheric carbon dioxide within the Calvin-Benson cycle that is utilized by most photosynthetic organisms. Despite this central role, RuBisCO's efficiency surprisingly struggles, with both a very slow turnover rate to products and also impaired substrate specificity, features that have long been an enigma as it would be assumed that its efficiency was under strong evolutionary pressure. RuBisCO's substrate specificity is compromised as it catalyzes a side-fixation reaction with atmospheric oxygen; empirical kinetic results show a trend to tradeoff between relative specificity and low catalytic turnover rate. Although the dominant hypothesis has been that the active-site chemistry constrains the enzyme's evolution, a more recent study on RuBisCO stability and adaptability has implicated competing selection pressures. Elucidating these constraints is crucial for directing future research on improving photosynthesis, as the current literature casts doubt on the potential effectiveness of site-directed mutagenesis to improve RuBisCO's efficiency. Here we use regression analysis to quantify the relationships between kinetic parameters obtained from empirical data sets spanning a wide evolutionary range of RuBisCOs. Most significantly we found that the rate constant for dissociation of CO₂ from the enzyme complex was much higher than previous estimates and comparable with the corresponding catalytic rate constant. Observed trends between relative specificity and turnover rate can be expressed as the product of negative and positive correlation factors. This provides an explanation in simple kinetic terms of both the natural variation of relative specificity as well as that obtained by reported site-directed mutagenesis results. We demonstrate that the kinetic behaviour shows a lesser rather than more constrained RuBisCO, consistent with growing empirical evidence of higher variability in relative specificity. In summary our analysis supports an explanation for the origin of the tradeoff between specificity and turnover as due to competition between protein stability and activity, rather than constraints between rate constants imposed by the underlying chemistry. Our analysis suggests that simultaneous improvement in both specificity and turnover rate of RuBisCO is possible.

Optimization of parameters for molecular dynamics simulation using smooth particle-mesh Ewald in GROMACS 4.5

Based on our critique of requirements for performing an efficient molecular dynamics simulation with the particle-mesh Ewald (PME) implementation in GROMACS 4.5, we present a computational tool to enable the discovery of parameters that produce a given accuracy in the PME approximation of the full electrostatics. Calculations on two parallel computers with different processor and communication structures showed that a given accuracy can be attained over a range of parameter space, and that the attributes of the hardware and simulation system control which parameter sets are optimal. This information can be used to find the fastest available PME parameter sets that achieve a given accuracy. We hope that this tool will stimulate future work to assess the impact of the quality of the PME approximation on simulation outcomes, particularly with regard to the trade-off between cost and scientific reliability in biomolecular applications.

Estimating Relative Disulfide Energies: An Accurate Ab Initio Potential Energy Surface

Disulfide torsional energy, a good predictor of disulfide redox potential in proteins, may be estimated by interpolation on a potential energy surface (PES) describing the twisting of diethyl disulfide through its three central dihedral angles. Here we update PES calculations at the M05-2X level of theory with the 6-31G(d) basis set. Although the surface shows no qualitative differences from an earlier MP2(full) PES, energy differences greater than 1 kJ mol–1 were seen for conformations with χ2 between –60° and 30°, or with χ3 below 60° or above 130°. This is particularly significant for highly strained disulfides that are likely to be spontaneously reduced by mechanical means. In benchmarking against the high-level G3X method, M05-2X showed significantly reduced root mean squared deviation compared with MP2(full) (1.0 versus 2.0 kJ mol–1 respectively). Results are incorporated into a web application that calculates relative torsional energies from disulfide dihedral angles (http://www.sbinf.org/applications/pes.html).

Green evolution and dynamic adaptations revealed by genomes of the marine picoeukaryotes Micromonas

Picoeukaryotes are a taxonomically diverse group of organisms less than 2 micrometers in diameter. Photosynthetic marine picoeukaryotes in the genus Micromonas thrive in ecosystems ranging from tropical to polar and could serve as sentinel organisms for biogeochemical fluxes of modern oceans during climate change. These broadly distributed primary producers belong to an anciently diverged sister clade to land plants. Although Micromonas isolates have high 18S ribosomal RNA gene identity, we found that genomes from two isolates shared only 90% of their predicted genes. Their independent evolutionary paths were emphasized by distinct riboswitch arrangements as well as the discovery of intronic repeat elements in one isolate, and in metagenomic data, but not in other genomes. Divergence appears to have been facilitated by selection and acquisition processes that actively shape the repertoire of genes that are mutually exclusive between the two isolates differently than the core genes. Analyses of the Micromonas genomes offer valuable insights into ecological differentiation and the dynamic nature of early plant evolution.

Putative functional role for the invariant aspartate 263 residue of Rhodospirillum rubrum Rubisco

Although aspartate residue D263 of Rhodospirillum rubrum Rubisco is close to the active site and invariant in all reported Rubiscos, its possible functional and structural roles in Rubisco activity have not been investigated. We have mutagenised D263 to several selected amino acids (asparagine, alanine, serine, glutamate, and glutamine) to probe possible roles in facilitating proton movements within the active site and maintaining structural positioning of key active-site groups. The mutants have been characterized by kinetic methods and by differential scanning calorimetry (DSC) to examine the effects of the substitutions on the stability of the folded state. We show that D263 is essential for maintaining effective levels of catalysis with the mutations reducing carboxylation variously by up to 100-fold but having less than 10% effect on the carboxylase/oxygenase specificity of the catalytic reaction. Removing the charge of the residue 263 side chain significantly strengthens binding of the activating (carbamylating) CO(2) molecule. In contrast, a charge on the 263 site has only a small influence on binding of the positively charged Mg(2+) ion, suggesting that the local protein structure provides different shielding of the formal charges on the Mg(2+) ion and the epsilon-lysine group of K191. Interestingly, introduction of an internal cavity (D263S and D263A) and insertion of an extra -CH(2)- group (D263E and D263Q) have opposite effects on catalysis, the former relatively small and the latter much larger, suggesting that the extra side-chain group induces a specific structural distortion that inhibits formation of the transition state. As the DSC results show that the mutations only slightly increase the kinetic stability of the folded state, we conclude that the rate-limiting (activated) step of unfolding involves substantial unfolding of the structure but not in the region of site 263. In summary, interaction of D263 with H287 of a largely electrostatic nature appears critical for maintaining correct positioning of catalytic groups in the active site. The conservation of D263 can thus be accounted for by its contribution to the maintenance of a finely tuned structure in this region abutting the active site.

Identification of the RGG box motif in Shadoo: RNA-binding and signaling roles?

Using comparative genomics and in-silico analyses, we previously identified a new member of the prion-protein (PrP) family, the gene SPRN, encoding the protein Shadoo (Sho), and suggested its functions might overlap with those of PrP. Extended bioinformatics and conceptual biology studies to elucidate Sho’s functions now reveal Sho has a conserved RGG-box motif, a well-known RNA-binding motif characterized in proteins such as FragileX Mental Retardation Protein. We report a systematic comparative analysis of RGG-box containing proteins which highlights the motif’s functional versatility and supports the suggestion that Sho plays a dual role in cell signaling and RNA binding in brain. These findings provide a further link to PrP, which has well-characterized RNA-binding properties.

Redefinition of rubisco carboxylase reaction reveals origin of water for hydration and new roles for active-site residues

Crystallographic, mutagenesis, kinetic, and computational studies on Rubisco over three decades have revealed much about its catalytic mechanism and the role played by several active-site residues. However, key questions remain unanswered. Specific details of the carboxylase and oxygenase mechanisms, required to underpin the rational re-engineering of Rubisco, are still speculative. Here we address critical gaps in knowledge with a definitive comprehensive computational investigation of the mechanism of carboxylase activity at the Rubisco active site. Density functional theory calculations (B3LYP/6-31G(d,p)) were performed on active-site fragment models of a size up to 77 atoms, not previously possible computationally. All amino acid residues suspected to play roles in the acid-base chemistry in the multistep reaction, and interacting directly with the central Mg (2+) atom and the reactive moiety of substrate and intermediates, were included. The results provide a firm basis for us to propose a novel mechanism for the entire sequence of reactions in the carboxylase catalysis and to define precise roles for the active-site residues, singly and in concert. In this mechanism, the carbamylated LYS201 plays a more limited role than previously proposed but is crucial for initiating the reaction by acting as a base in the enolization. We suggest a wider role for HIS294, with involvement in the carboxylation, hydration, and C2-C3 bond-scission steps, consistent with the suggestion of Harpel et al. (1998) but contrary to the consensus view of Cleland et al. (1998). In contrast to the common assumption that the water molecule for the hydration step comes from within the active site, we propose that the Mg-coordinated water is not dissociated at the start of the gas-addition reaction but rather remains coordinated and is used for the hydration of the C3 carbon atom. New roles are also proposed for LYS175, GLU204, and HIS294. The mechanism suggests roles in the gas-addition step for residues in three spatially distinct regions of the active site, HIS294 and LYS334 in the C-terminal domain of the large subunit (LSU), but also hitherto unsuspected roles for a cluster of three residues (ASN123, GLU60, and TYR20) in the N-terminal domain of the partner LSU of the dimer containing the active site. Our new mechanism is supported by existing experimental data, provides new convincing interpretations of previously puzzling data, and allows new insights into mutational strategies for improving Rubisco activity.

Evaluating the stability of disulfide bridges in proteins: a torsional potential energy surface for diethyl disulfide

Disulfide bonds formed by the oxidation of cysteine residues in proteins are the major form of intra- and inter-molecular covalent linkages in the polypeptide chain. To better understand the conformational energetics of this linkage, we have used the MP2(full)/6-31G(d) method to generate a full potential energy surface (PES) for the torsion of the model compound diethyl disulfide (DEDS) around its three critical dihedral angles (χ₂, χ₃, χ₂'). The use of ten degree increments for each of the parameters resulted in a continuous, fine-grained surface. This allowed us to accurately predict the relative stabilities of disulfide bonds in high resolution structures from the Protein Data Bank. The MP2(full) surface showed significant qualitative differences from the PES calculated using the Amber force field. In particular, a different ordering was seen for the relative energies of the local minima. Thus, Amber energies are not reliable for comparison of the relative stabilities of disulfide bonds. Surprisingly, the surface did not show a minimum associated with χ₂ ∼ - 60°, χ₃ ∼ 90, χ₂' ∼ - 60°. This is due to steric interference between Hα atoms. Despite this, significant populations of disulfides were found to adopt this conformation. In most cases this conformation is associated with an unusual secondary structure motif, the cross-strand disulfide. The relative instability of cross-strand disulfides is of great interest, as they have the potential to act as functional switches in redox processes.

Integron-sequestered dihydrofolate reductase: a recently redeployed enzyme

The introduction and wide use of antibacterial drugs has resulted in the emergence of resistant organisms. DfrB dihydrofolate reductase (DHFR) is a bacterial enzyme that is uniquely associated with mobile gene cassettes within integrons, and confers resistance to the drug trimethoprim. This enzyme has intrigued microbiologists since it was discovered more than thirty years ago because of its simple structure, enzymatic inefficiency and its virtual insensitivity to trimethoprim. Here, for the first time, a comprehensive discussion of genetic, evolutionary, structural and functional studies of this enzyme is presented together. This information supports the ideas that DfrB DHFR is a poorly adapted catalyst and has recently been recruited to perform a novel enzymatic activity in response to selective pressure.

The C-type lectin-like domain superfamily

The superfamily of proteins containing C-type lectin-like domains (CTLDs) is a large group of extracellular Metazoan proteins with diverse functions. The CTLD structure has a characteristic double-loop ('loop-in-a-loop') stabilized by two highly conserved disulfide bridges located at the bases of the loops, as well as a set of conserved hydrophobic and polar interactions. The second loop, called the long loop region, is structurally and evolutionarily flexible, and is involved in Ca2+-dependent carbohydrate binding and interaction with other ligands. This loop is completely absent in a subset of CTLDs, which we refer to as compact CTLDs; these include the Link/PTR domain and bacterial CTLDs. CTLD-containing proteins (CTLDcps) were originally classified into seven groups based on their overall domain structure. Analyses of the superfamily representation in several completely sequenced genomes have added 10 new groups to the classification, and shown that it is applicable only to vertebrate CTLDcps; despite the abundance of CTLDcps in the invertebrate genomes studied, the domain architectures of these proteins do not match those of the vertebrate groups. Ca2+-dependent carbohydrate binding is the most common CTLD function in vertebrates, and apparently the ancestral one, as suggested by the many humoral defense CTLDcps characterized in insects and other invertebrates. However, many CTLDs have evolved to specifically recognize protein, lipid and inorganic ligands, including the vertebrate clade-specific snake venoms, and fish antifreeze and bird egg-shell proteins. Recent studies highlight the functional versatility of this protein superfamily and the CTLD scaffold, and suggest further interesting discoveries have yet to be made.

Combining docking and molecular dynamic simulations in drug design

A rational approach is needed to maximize the chances of finding new drugs, and to exploit the opportunities of potential new drug targets emerging from genomic and proteomic initiatives, and from the large libraries of small compounds now readily available through combinatorial chemistry. Despite a shaky early history, computer-aided drug design techniques can now be effective in reducing costs and speeding up drug discovery. This happy outcome results from development of more accurate and reliable algorithms, use of more thoughtfully planned strategies to apply them, and greatly increased computer power to allow studies with the necessary reliability to be performed. Our review focuses on applications and protocols, with the main emphasis on critical analysis of recent studies where docking calculations and molecular dynamics (MD) simulations were combined to dock small molecules into protein receptors. We highlight successes to demonstrate what is possible now, but also point out drawbacks and future directions. The review is structured to lead the reader from the simpler to more compute-intensive methods. Thus, while inexpensive and fast docking algorithms can be used to scan large compound libraries and reduce their size, more accurate but expensive MD simulations can be applied when a few selected ligand candidates remain. MD simulations can be used: during the preparation of the protein receptor before docking, to optimize its structure and account for protein flexibility; for the refinement of docked complexes, to include solvent effects and account for induced fit; to calculate binding free energies, to provide an accurate ranking of the potential ligands; and in the latest developments, during the docking process itself to find the binding site and correctly dock the ligand a priori.

Multiple ligand-binding modes in bacterial R67 dihydrofolate reductase

R67 dihydrofolate reductase (DHFR), a bacterial plasmid-encoded enzyme associated with resistance to the drug trimethoprim, shows neither sequence nor structural homology with the chromosomal DHFR. It presents a highly symmetrical toroidal structure, where four identical monomers contribute to the unique central active-site pore. Two reactants (dihydrofolate, DHF), two cofactors (NADPH) or one of each (R67*DHF*NADPH) can be found simultaneously within the active site, the last one being the reactive ternary complex. As the positioning of the ligands has proven elusive to empirical determination, we addressed the problem from a theoretical perspective. Several potential structures of the ternary complex were generated using the docking programs AutoDock and FlexX. The variability among the final poses, many of which conformed to experimental data, prompted us to perform a comparative scoring analysis and molecular dynamics simulations to assess the stability of the complexes. Analysis of ligand-ligand and ligand-protein interactions along the 4 ns trajectories of eight different structures allowed us to identify important inter-ligand contacts and key protein residues. Our results, combined with published empirical data, clearly suggest that multipe binding modes of the ligands are possible within R67 DHFR. While the pterin ring of DHF and the nicotinamide ring of NADPH assume a stacked endo-conformation at the centre of the pore, probably assisted by V66, Q67 and I68, the tails of the molecules extend towards opposite ends of the cavity, adopting multiple configurations in a solvent rich-environment where hydrogen-bond interactions with K32 and Y69 may play important roles.

The prion protein gene: identifying regulatory signals using marsupial sequence

The function of the prion protein gene (PRNP) and its normal product PrP(C) is elusive. We used comparative genomics as a strategy to understand the normal function of PRNP. As the reliability of comparisons increases with the number of species and increased evolutionary distance, we isolated and sequenced a 66.5 kb BAC containing the PRNP gene from a distantly related mammal, the model Australian marsupial Macropus eugenii (tammar wallaby). Marsupials are separated from eutherians such as human and mouse by roughly 180 million years of independent evolution. We found that tammar PRNP, like human PRNP, has two exons. Prion proteins encoded by the tammar wallaby and a distantly related marsupial, Monodelphis domestica (Brazilian opossum) PRNP contain proximal PrP repeats with a distinct, marsupial-specific composition and a variable number. Comparisons of tammar wallaby PRNP with PRNPs from human, mouse, bovine and ovine allowed us to identify non-coding gene regions conserved across the marsupial-eutherian evolutionary distance, which are candidates for regulatory regions. In the PRNP 3' UTR we found a conserved signal for nuclear-specific polyadenylation and the putative cytoplasmic polyadenylation element (CPE), indicating that post-transcriptional control of PRNP mRNA activity is important. Phylogenetic footprinting revealed conserved potential binding sites for the MZF-1 transcription factor in both upstream promoter and intron/intron 1, and for the MEF2, MyT1, Oct-1 and NFAT transcription factors in the intron(s). The presence of a conserved NFAT-binding site and CPE indicates involvement of PrP(C) in signal transduction and synaptic plasticity.

An efficient Z-score algorithm for assessing sequence alignments

We describe an alternative method for scoring of the pairwise alignment of two biological sequences. Designed to overcome the bias due to the composition of the alignment, it measures the distance (in standard deviations) between the given alignment and the mean value of all other alignments that can be obtained by a permutation of either sequence. We demonstrate that the standard deviation can be calculated efficiently. By concentrating upon the ungapped case, the mean and standard deviation can be calculated exactly and in two steps, the first being O(N) time, where N is the length of the sequence, the second in a fixed number of calculations, i.e., in O(1) time. We argue that this statistic is a more consistent measure than a similarity score based upon a standard scoring matrix. Even in the ungapped case, the statistic proves in many cases to be more accurate than the commonly used (FASTA) (Pearson and Lipman, 1988) gapped Z-score in which the sequence is matched against a random sample of the database. We demonstrate the use of the POZ-score as a secondary filter which screens out several well-known types of false positive, reducing the amount of manual screening to be done by the biologist.

Computational methods for the study of enzymic reaction mechanisms III: A perturbation plus QM/MM approach for calculating relative free energies of protonation

We describe a coupling parameter, that is, perturbation, approach to effectively create and annihilate atoms in the quantum mechanical Hamiltonian within the closed shell restricted Hartree-Fock formalism. This perturbed quantum mechanical atom (PQA) method is combined with molecular mechanics (MM) methods (PQA/MM) within a molecular dynamics simulation, to model the protein environment (MM region) effects that also make a contribution to the overall free energy change. Using the semiempirical PM3 method to model the QM region, the application of this PQA/MM method is illustrated by calculation of the relative protonation free energy of the conserved OD2 (Asp27) and the N5 (dihydrofolate) proton acceptor sites in the active site of Escherichia coli dihydrofolate reductase (DHFR) with the bound nicotinamide adenine dinucleotide phosphate (NADPH) cofactor. For a number of choices for the QM region, the relative protonation free energy was calculated as the sum of contributions from the QM region and the interaction between the QM and MM regions via the thermodynamic integration (TI) method. The results demonstrate the importance of including the whole substrate molecule in the QM region, and the overall protein (MM) environment in determining the relative stabilities of protonation sites in the enzyme active site. The PQA/MM free energies obtained by TI were also compared with those estimated by a less computationally demanding nonperturbative method based on the linear response approximation (LRA). For some choices of QM region, the total free energies calculated using the LRA method were in very close agreement with the PQA/MM values. However, the QM and QM/MM component free energies were found to differ significantly between the two methods.

31P MAS-NMR of human erythrocytes: independence of cell volume from angular velocity

31P magic angle spinning NMR (MAS-NMR) spectra were obtained from suspensions of human red blood cells (RBCs) that contained the cell-volume-sensitive probe molecule, dimethyl methylphosphonate (DMMP). A mathematical representation of the spectral-peak shape, including the separation and width-at-half-height in the 31P NMR spectra, as a function of rotor speed, enabled us to explore the extent to which a change in cell volume would be reflected in the spectra if it occurred. We concluded that a fractional volume change in excess of 3% would have been detected by our experiments. Thus, the experiments indicated that the mean cell volume did not change by this amount even at the highest spinning rate of 7 kHz. The mean cell volume and intracellular 31P line-width were independent of the packing density of the cells and of the initial cell volume. The relationship of these conclusions to other non-NMR studies of pressure effects on cells is noted.

Conformation of prion protein repeat peptides probed by FRET measurements and molecular dynamics simulations

We report the combined use of steady-state fluorescence resonance energy transfer (FRET) experiments and molecular dynamics (MD) simulations to investigate conformational distributions of the prion protein (PrP) repeat system. FRET was used for the first time to probe the distance, as a function of temperature and pH, between a donor Trp residue and an acceptor dansyl group attached to the N-terminus in seven model peptides containing one to three repeats of the second decarepeat of PrP from marsupial possum (PHPGGSNWGQ)nG, and one and two human PrP consensus octarepeats (PHGGGWGQ)nG. In multirepeat peptides, single-Trp mutants were made by replacing other Trp(s) with Phe. As previous work has shown PrP repeats do not adopt a single preferred stable conformation, the FRET values are averages reflecting heterogeneity in the donor-acceptor distances. The T-dependence of the conformational distributions, and derived average dansyl-Trp distances, were obtained directly from MD simulation of the marsupial dansyl-PHPGGSNWGQG peptide. The results show excellent agreement between the FRET and MD T-dependent distances, and demonstrate the remarkable sensitivity and reproducibility of the FRET method in this first-time use for a set of disordered peptides. Based on the results, we propose a model involving cation-pi or pi-pi His-Trp interactions to explain the T- (5-85 degrees C) and pH- (6.0, 7.2) dependencies on distance, with HW i, i + 4 or WH i, i + 4 separations in sequence being more stable than HW i, i + 6 or WH i, i + 6 separations. The model has peptides adopting loosely folded conformations, with dansyl-Trp distances very much less than estimates for fully extended conformations, for example, approximately 16 vs. 33, approximately 21 vs. 69, and approximately 22 vs. 106 A for 1-3 decarepeats, and approximately 14 vs. 25 and approximately 19 vs. 54 A for 1-2 octarepeats, respectively. The study demonstrates the usefulness of combining FRET with MD, a combination reported only once previously. Initial "mapping" of the conformational distribution of flexible peptides by simulation can assist in designing and interpreting experiments using steady-state intensity methods, and indicating how time-resolved or anisotropy methods might be used.

Evolution of vertebrate genes related to prion and Shadoo proteins--clues from comparative genomic analysis

Recent findings of new genes in fish related to the prion protein (PrP) gene PRNP, including our recent report of SPRN coding for Shadoo (Sho) protein found also in mammals, raise issues of their function and evolution. Here we report additional novel fish genes found in public databases, including a duplicated SPRN gene, SPRNB, in Fugu, Tetraodon, carp, and zebrafish encoding the Sho2 protein, and we use comparative genomic analysis to analyze the evolutionary relationships and to infer evolutionary trajectories of the complete data set. Phylogenetic footprinting performed on aligned human, mouse, and Fugu SPRN genes to define candidate regulatory promoter regions, detected 16 conserved motifs, three of which are known transcription factor-binding sites for a receptor and transcription factors specific to or associated with expression in brain. This result and other homology-based (VISTA global genomic alignment; protein sequence alignment and phylogenetics) and context-dependent (genomic context; relative gene order and orientation) criteria indicate fish and mammalian SPRN genes are orthologous and suggest a strongly conserved basic function in brain. Whereas tetrapod PRNPs share context with the analogous stPrP-2-coding gene in fish, their sequences are diverged, suggesting that the tetrapod and fish genes are likely to have significantly different functions. Phylogenetic analysis predicts the SPRN/SPRNB duplication occurred before divergence of fish from tetrapods, whereas that of stPrP-1 and stPrP-2 occurred in fish. Whereas Sho appears to have a conserved function in vertebrate brain, PrP seems to have an adaptive role fine-tuned in a lineage-specific fashion. An evolutionary model consistent with our findings and literature knowledge is proposed that has an ancestral prevertebrate SPRN-like gene leading to all vertebrate PrP-related and Sho-related genes. This provides a new framework for exploring the evolution of this unusual family of proteins and for searching for members in other fish branches and intermediate vertebrate groups.

C-type lectin-like domains in Fugu rubripes

Background

Members of the C-type lectin domain (CTLD) superfamily are metazoan proteins functionally important in glycoprotein metabolism, mechanisms of multicellular integration and immunity. Three genome-level studies on human, C. elegans and D. melanogaster reported previously demonstrated almost complete divergence among invertebrate and mammalian families of CTLD-containing proteins (CTLDcps).

Results

We have performed an analysis of CTLD family composition in Fugu rubripes using the draft genome sequence. The results show that all but two groups of CTLDcps identified in mammals are also found in fish, and that most of the groups have the same members as in mammals. We failed to detect representatives for CTLD groups V (NK cell receptors) and VII (lithostathine), while the DC-SIGN subgroup of group II is overrepresented in Fugu. Several new CTLD-containing genes, highly conserved between Fugu and human, were discovered using the Fugu genome sequence as a reference, including a CSPG family member and an SCP-domain-containing soluble protein. A distinct group of soluble dual-CTLD proteins has been identified, which may be the first reported CTLDcp group shared by invertebrates and vertebrates. We show that CTLDcp-encoding genes are selectively duplicated in Fugu, in a manner that suggests an ancient large-scale duplication event. We have verified 32 gene structures and predicted 63 new ones, and make our annotations available through a distributed annotation system (DAS) server and their sequences as additional files with this paper.

Conclusions

The vertebrate CTLDcp family was essentially formed early in vertebrate evolution and is completely different from the invertebrate families. Comparison of fish and mammalian genomes revealed three groups of CTLDcps and several new members of the known groups, which are highly conserved between fish and mammals, but were not identified in the study using only mammalian genomes. Despite limitations of the draft sequence, the Fugu rubripes genome is a powerful instrument for gene discovery and vertebrate evolutionary analysis. The composition of the CTLDcp superfamily in fish and mammals suggests that large-scale duplication events played an important role in the evolution of vertebrates.

Shadoo, a new protein highly conserved from fish to mammals and with similarity to prion protein

We report evidence from cDNA isolation and expression analysis as well as analyses of genome, expressed sequence tag (EST), cDNA and expression databases for a new gene named SPRN (shadow of prion protein). SPRN comprises two exons, with the open reading frame (ORF) contained within exon 2, and codes for a protein of 130-150 amino acids named Shadoo (Japanese shadow), predicted to be extracellular and GPI-anchored. The SPRN gene was found in fish (zebrafish, Fugu) and mammals (mouse, rat, human). Conservation of order and transcription orientation of two proximal genes between fishes and mammals strongly indicates gene orthology. Sequence comparison shows: a highly conserved N-terminal signal sequence; Arg-rich basic region containing up to six tetrarepeats of consensus XXRG (where X is G, A or S); a hydrophobic region of 20 residues with strong homology to PrP; a less conserved C-terminal domain containing a conserved glycosylation motif; and a C-terminal peptide predicted to be a signal sequence for glycophosphotidylinositol (GPI)-anchor attachment. Fish Shadoos (Sho) show well conserved sequences (identity 54%) over 106 amino acids (zebrafish length), and conservation among the mammalian sequences is very high (identity 81-96%). The fish and mammalian sequences are also well conserved, particularly for zebrafish, to beyond the end of the hydrophobic sequence (identity 41-53%, 78 amino acids, zebrafish length). The overall structure appears closely related to prion proteins (PrPs), although the C-terminal domains of Shos are quite different from those of PrPs, for which conformational changes in mammals are implicated in disease. The structural similarity is particularly interesting given recent reports of three new genes with similarities to PrPs found in Fugu (PrP-like, PrP-461/stPrP-1 and stPrP-2) and other fish, but for which direct evolution to higher vertebrate PrPs is unlikely and for which no other mammalian homologues have been found. Database information indicates expression of SPRN in embryo, brain and retina of mouse and rat, hippocampus of human, and in embryo and retina of zebrafish, and we directly confirmed a strikingly specific expression of the mammalian (human, mouse, rat) transcripts in whole brain. This result together with some common structural features led to the suggestive hypothesis of a possible functional link between mammalian PrP and Sho proteins.

Comparative analysis of structural properties of the C-type-lectin-like domain (CTLD)

The superfamily of proteins containing the C-type-lectin-like domain (CTLD) is a group of abundant extracellular metazoan proteins characterized by evolutionary flexibility and functional versatility. Several CTLDs are also found in parasitic prokaryotes and viruses. The 37 distinct currently available CTLD structures demonstrate significant structural conservation despite low or undetectable sequence similarity. Our aim in this study was to perform an extensive comparative analysis of all available CTLD structures to establish the most conserved structural features of the fold, and to test and extend the early analysis of Drickamer. By implication, these features should be those critical for maintenance of integrity of the fold. By analyzing CTLD structures superimposed by several methods, we have established groups of conserved structural positions involved in fold maintenance but not in ligand binding; these are consistent with the fold's known functional flexibility. In addition to the well-recognized disulfide bridges, groups of conserved residues are involved in hydrophobic interactions stabilizing the core of the fold and the long loop region, and in an alpha2-beta1-beta5 polar interaction. Evaluation of the conclusions of the structure comparison study compared with alignments of all available human, mouse and Caenorhabditis elegans CTLD sequences showed that conservation patterns are preserved throughout the whole CTLD sequence space. Our observations provide an improved understanding of CTLD structure, and will help in identification of new CTLDs and the mechanisms that drive and constrain the coevolution of the structure and function of the fold.

Synthesis of quaternised 2-aminopyrimido[4,5-d]pyrimidin-4(3H)-ones and their biological activity with dihydrofolate reductase

In a program to design and develop mechanism-based compounds active as substrates and inhibitors of dihydrofolate reductase (DHFR), we report the synthesis and physical properties of the 6-methyl- (7), 8-methyl- (8a), and 8-ethyl- (8b) derivatives of the parent 2-aminopyrimido[4,5-d]pyrimidin-4-(3H)-one (6). These compounds are the first members of a class of heterocycles related to 8-alkylpterins (N8-alkyl-2-aminopteridin-4(8H)-ones) (2a-2c), which have been shown to be novel substrates for DHFR. Three methods were developed for the synthesis of target compounds 7, 8a and 8b; however, the optimum yields (1-8%) could not be improved because the products decomposed by ring opening (e.g. to 2,4-diamino-5-methyliminomethylpyrimidin-6(1H)-one (9)) under the reaction conditions. The marked pi-electron deficiency of compounds 7, 8a and 8b is the likely cause for the susceptibility of the quaternised pyrimidine ring in the related cations 10, 15a and 15b, respectively, to add nucleophiles, thus promoting the opening of the pyrimidopyrimidine ring system. 1H-NMR spectroscopic studies of compounds 7, 8a and 8b revealed a fast and reversible covalent hydration of the associated cations across the C7z.sbnd;N8 bond for the N6-methyl derivative 7 and across the N6z.sbnd;C7 bond for the N8-methyl derivative 8a. UV spectroscopic studies of methyl derivatives 7 and 8a as well as the parent heterocycle 6 showed that protonation of the latter occurred at N1, while methylation with iodomethane proceeded at N6 and N8. The basicities of the N-methyl derivatives 7 and 8a (pK(a) ca. 5.5) are similar to those of 8-alkylpterins 2; this is an essential element of the design to promote binding to DHFR in their protonated form. Enzyme kinetics of 7, 8a and 8b with chicken DHFR confirmed our predictions that they are substrates, with apparent K(m) values of 3.8, 0.08, and 0.65 mM, and apparent V(max) values of 0.47, 2.27, and 0.30 nmol L(-1) min(-1) (for enzyme concentration 0.122 micro M), respectively. The parent compound 6 was not a substrate.

Comparison of linear-scaling semiempirical methods and combined quantum mechanical/molecular mechanical methods for enzymic reactions. II. An energy decomposition analysis

QM/MM methods have been developed as a computationally feasible solution to QM simulation of chemical processes, such as enzyme-catalyzed reactions, within a more approximate MM representation of the condensed-phase environment. However, there has been no independent method for checking the quality of this representation, especially for highly nonisotropic protein environments such as those surrounding enzyme active sites. Hence, the validity of QM/MM methods is largely untested. Here we use the possibility of performing all-QM calculations at the semiempirical PM3 level with a linear-scaling method (MOZYME) to assess the performance of a QM/MM method (PM3/AMBER94 force field). Using two model pathways for the hydride-ion transfer reaction of the enzyme dihydrofolate reductase studied previously (Titmuss et al., Chem Phys Lett 2000, 320, 169-176), we have analyzed the reaction energy contributions (QM, QM/MM, and MM) from the QM/MM results and compared them with analogous-region components calculated via an energy partitioning scheme implemented into MOZYME. This analysis further divided the MOZYME components into Coulomb, resonance and exchange energy terms. For the model in which the MM coordinates are kept fixed during the reaction, we find that the MOZYME and QM/MM total energy profiles agree very well, but that there are significant differences in the energy components. Most significantly there is a large change (approximately 16 kcal/mol) in the MOZYME MM component due to polarization of the MM region surrounding the active site, and which arises mostly from MM atoms close to (<10 A) the active-site QM region, which is not modelled explicitly by our QM/MM method. However, for the model where the MM coordinates are allowed to vary during the reaction, we find large differences in the MOZYME and QM/MM total energy profiles, with a discrepancy of 52 kcal/mol between the relative reaction (product-reactant) energies. This is largely due to a difference in the MM energies of 58 kcal/mol, of which we can attribute approximately 40 kcal/mol to geometry effects in the MM region and the remainder, as before, to MM region polarization. Contrary to the fixed-geometry model, there is no correlation of the MM energy changes with distance from the QM region, nor are they contributed by only a few residues. Overall, the results suggest that merely extending the size of the QM region in the QM/MM calculation is not a universal solution to the MOZYME- and QM/MM-method differences. They also suggest that attaching physical significance to MOZYME Coulomb, resonance and exchange components is problematic. Although we conclude that it would be possible to reparameterize the QM/MM force field to reproduce MOZYME energies, a better way to account for both the effects of the protein environment and known deficiencies in semiempirical methods would be to parameterize the force field based on data from DFT or ab initio QM linear-scaling calculations. Such a force field could be used efficiently in MD simulations to calculate free energies.

Copper-induced conformational change in a marsupial prion protein repeat peptide probed using FTIR spectroscopy

We report the first Fourier transform infrared analysis of prion protein (PrP) repeats and the first study of PrP repeats of marsupial origin. Large changes in the secondary structure and an increase in hydrogen bonding within the peptide groups were evident from a red shift of the amide I band by >7 cm(-1) and an approximately five-fold reduction in amide hydrogen-deuterium exchange for peptide interacting with Cu(2+) ions. Changes in the tertiary structure upon copper binding were also evident from the appearance of a new band at 1564 cm(-1), which arises from the ring vibration of histidine. The copper-induced conformational change is pH dependent, and occurs at pH >7.