NWChem: Past, present, and future.

Specialized computational chemistry packages have permanently reshaped the landscape of chemical and materials science by providing tools to support and guide experimental efforts and for the prediction of atomistic and electronic properties. In this regard, electronic structure packages have played a special role by using first-principle-driven methodologies to model complex chemical and materials processes. Over the past few decades, the rapid development of computing technologies and the tremendous increase in computational power have offered a unique chance to study complex transformations using sophisticated and predictive many-body techniques that describe correlated behavior of electrons in molecular and condensed phase systems at different levels of theory. In enabling these simulations, novel parallel algorithms have been able to take advantage of computational resources to address the polynomial scaling of electronic structure methods. In this paper, we briefly review the NWChem computational chemistry suite, including its history, design principles, parallel tools, current capabilities, outreach, and outlook.

Expression regulation by a methyl-CpG binding domain in an E. coli based, cell-free TX-TL system.

Cytosine methylation plays an important role in the epigenetic regulation of eukaryotic gene expression. The methyl-CpG binding domain (MBD) is common to a family of eukaryotic transcriptional regulators. How MBD, a stretch of about 80 amino acids, recognizes CpGs in a methylation dependent manner, and as a function of sequence, is only partly understood. Here we show, using an Escherichia coli cell-free expression system, that MBD from the human transcriptional regulator MeCP2 performs as a specific, methylation-dependent repressor in conjunction with the BDNF (brain-derived neurotrophic factor) promoter sequence. Mutation of either base flanking the central CpG pair changes the expression level of the target gene. However, the relative degree of repression as a function of MBD concentration remains unaltered. Molecular dynamics simulations that address the DNA B fiber ratio and the handedness reveal cooperative transitions in the promoter DNA upon MBD binding that correlate well with our experimental observations. We suggest that not only steric hindrance, but also conformational changes of the BDNF promoter as a result of MBD binding are required for MBD to act as a specific inhibitory element. Our work demonstrates that the prokaryotic transcription machinery can reproduce features of epigenetic mammalian transcriptional regulatory elements.

IMP2/p62 induces genomic instability and an aggressive hepatocellular carcinoma phenotype.

Hepatocellular carcinoma (HCC) represents the third leading cause of cancer-related deaths and commonly develops in inflammatory environments. The IGF2 mRNA-binding protein IMP2-2/IGF2BP2-2/p62 was originally identified as an autoantigen in HCC. Aim of this study was to investigate a potential pathophysiological role of p62 in hepatocarcinogenesis. Human HCC tissue showed overexpression of IMP2, which strongly correlated with the fetal markers AFP and DLK1/Pref-1/FA-1 and was particularly elevated in tumors with stem-like features and hypervascularization. Molecular classification of IMP2-overexpressing tumors revealed an aggressive phenotype. Livers of mice overexpressing the IMP2 splice variant p62 highly expressed the stem cell marker DLK1 and secreted DLK1 into the blood. p62 was oncogenic: diethylnitrosamine (DEN)-treated p62 transgenic mice exhibited a higher tumor incidence and multiplicity than wild types. Tumors of transgenics showed a more aggressive and stem-like phenotype and displayed more oncogenic chromosomal aberrations determined with aCGH analysis. DEN-treated p62 transgenic mice exhibited distinct signs of inflammation, such as inflammatory cytokine expression and oxidative stress markers, that is, thiobarbituric acid-reactive substance (TBARS) levels. Reactive oxygen species (ROS) production was elevated in HepG2 cells, which either overexpressed p62 or were treated with DLK1. p62 induced this ROS production by a DLK1-dependent induction and activation of the small Rho-GTPase RAC1, activating NADPH oxidase and being overexpressed in human HCC. Our data indicate that p62/IMP2 promotes hepatocarcinogenesis by an amplification of inflammation.

Functional classification of membrane transporters and channels based on filtered TM/non-TM amino acid composition.

Membrane transporters catalyze the transport of small solute molecules across biological barriers such as lipid bilayer membranes. As the experimental annotation of which proteins transport which substrates is incomplete it is highly desirable to develop computational methods that can assist in the classification and substrate annotation of putative membrane transport proteins. Here, we determined the similarity of membrane transporter sequences annotated in the Transport Classification Database (Saier et al., Nucleic Acids Res 2006, 34, D181-D186) and Arabidopsis thaliana membrane transporters annotated in the database Aramemnon (Schwacke et al., Plant Physiol 2003, 131, 16-26). The similarity measure was based on the amino acid composition either considering the full sequences or separately in the transmembrane (TM) and external parts of the sequences. We considered four different substrate sets and three different subfamilies and tried to classify the given proteins into these classes. Family or substrate prediction based on the simple amino acid frequency had an average accuracy of 76%. The differentiation between TM and non-TM regions led to an improved accuracy of 80% on average.

Interfacing Brownian dynamics simulations.

Starting from the flux of particles in a Brownian dynamics simulation we derive boundary conditions, which allow us (i) to couple a Brownian dynamics calculation to a reservoir of particles of a given density, i.e., setting up constant density boundary conditions, and (ii) to build an interface between Brownian dynamics and a diffusional treatment of adjacent simulation volumes. With these algorithms it is sometimes possible to dramatically reduce the system size--and therefore the necessary resources--of multiparticle Brownian dynamics calculations. In this paper we give one-dimensional examples which illustrate potential applications and savings.

Brownian dynamics simulations of simplified cytochrome c molecules in the presence of a charged surface.

Simulations were performed for up to 150 simplified spherical horse heart cytochrome c molecules in the presence of a charged surface, which serves as an approximate model for a lipid membrane. Screened electrostatic and short-ranged attractive as well as repulsive van der Waals forces for interparticle and particle-membrane interactions are utilized in the simulations. At a distance from the membrane, where particle-membrane interactions are negligible, the simulation is coupled to a noninteraction continuum analogous to a heat bath [Geyer et al., J. Chem. Phys. 120, 4573 (2004)]. From the particles' density profiles perpendicular to the planar surface binding isotherms are derived and compared to experimental results [Heimburg et al. (1999)]. Using a negatively charged structureless membrane surface a saturation effect was found for relatively large particle concentrations. Since biological membranes often contain membrane proteins, we also studied the influence of additional charges on our model membrane mimicking bacterial reaction centers. We find that the onset of the saturation occurs for much lower concentrations and is sensitive to the detailed implementation. Therefore we suggest that local distortion of membrane planarity (undulation), or lipid demixing, or the presence of charged integral membrane proteins create preferential binding sites on the membrane. Only then do we observe saturation at physiological concentrations.

Re-face stereospecificity of methylenetetrahydromethanopterin and methylenetetrahydrofolate dehydrogenases is predetermined by intrinsic properties of the substrate.

Four different dehydrogenases are known that catalyse the reversible dehydrogenation of N5,N10-methylenetetrahydromethanopterin (methylene-H4MPT) or N5,N10-methylenetetrahydrofolate (methylene-H4F) to the respective N5,N10-methenyl compounds. Sequence comparison indicates that the four enzymes are phylogenetically unrelated. They all catalyse the Re-face-stereospecific removal of the pro-R hydrogen atom of the coenzyme's methylene group. The Re-face stereospecificity is in contrast to the finding that in solution the pro-S hydrogen atom of methylene-H4MPT and of methylene-H4F is more reactive to heterolytic cleavage. For a better understanding we determined the conformations of methylene-H4MPT in solution and when enzyme-bound by using NMR spectroscopy and semiempirical quantum mechanical calculations. For the conformation free in solution we find an envelope conformation for the imidazolidine ring, with the flap at N10. The methylene pro-S C-H bond is anticlinal and the methylene pro-R C-H bond is synclinal to the lone electron pair of N10. Semiempirical quantum mechanical calculations of heats of formation of methylene-H4MPT and methylene-H4F indicate that changing this conformation into an activated one in which the pro-S C-H bond is antiperiplanar, resulting in the preformation of the leaving hydride, would require a deltadeltaH(f) of +53 kJ mol-1 for methylene-H4MPT and of +51 kJ mol-1 for methylene-H4F. This is almost twice the energy required to force the imidazolidine ring in the enzyme-bound conformation of methylene-H4MPT (+29 kJ mol-1) or of methylene-H4F (+35 kJ mol-1) into an activated conformation in which the pro-R hydrogen atom is antiperiplanar to the lone electron pair of N10. The much lower energy for pro-R hydrogen activation thus probably predetermines the Re-face stereospecificity of the four dehydrogenases. Results are also presented explaining why the chemical reduction of methenyl-H4MPT+ and methenyl-H4F+ with NaBD4 proceeds Si-face-specific, in contrast to the enzyme-catalysed reaction.

Multicopy molecular dynamics simulations suggest how to reconcile crystallographic and product formation data for camphor enantiomers bound to cytochrome P-450cam.

Multiple ligand binding modes are possible in many enzyme active sites; their presence in cytochrome P450cam (P450cam) is evident from crystallographic studies of the binding of thiocamphor and phenylimidazoles. Here, we use multicopy molecular dynamics simulations to compare the binding modes of (1R)- and (1S)-camphor in the active site of P450cam. Simulations with (1R)-camphor, the natural substrate, serve to calibrate our protocol: 19 out of 20 copies of (1R)-camphor converged to coordinates very close to those observed for (1R)-camphor in its crystallographic complex with P450cam during the simulations. Simulations with the (1S)-camphor enantiomer showed greater mobility of the substrate, consistent with spectroscopic data, and resulted in 3 major binding modes. One of these is similar to the major conformation (of the two conformations assigned) in a recently determined crystal structure, but this conformation is not correctly oriented for regiospecific hydroxylation at C-5. The simulations, however, provide evidence for reorientation of (1S)-camphor upon formation of the reactive Fe-O intermediate to an orientation suitable for hydroxylation. The simulations thus permit rationalisation of the apparent inconsistency between the crystal structure and the reaction products.

Phosphoryl transfer by a concerted reaction mechanism in UMP/CMP-kinase.

The reaction mechanism of phosphoryl transfer catalyzed by UMP/CMP-kinase from Dictyostelium discoideum was investigated by semiempirical AM1 molecular orbital computations of an active site model system derived from crystal structures that contain a transition state analog or a bisubstrate inhibitor. The computational results suggest that the nucleoside monophosphate must be protonated for the forward reaction while it is unprotonated in the presence of aluminium fluoride, a popular transition state analog for phosphoryl transfer reactions. Furthermore, a compactification of the active site model system during the reaction and for the corresponding complex containing AlF3 was observed. For the active site residues that are part of the LID domain, conformational flexibility during the reaction proved to be crucial. On the basis of the calculations, a concerted phosphoryl transfer mechanism is suggested that involves the synchronous shift of a proton from the monophosphate to the transferred PO3-group. The proposed mechanism is thus analogous to the phosphoryl transfer mechanism in cAMP-dependent protein kinase that phosphorylates the hydroxyl groups of serine residues.

Dynamical properties of fasciculin-2.

Fasciculin-2 (FAS2) is a potent protein inhibitor of the hydrolytic enzyme acetylcholinesterase. A 2-ns isobaric-isothermal ensemble molecular dynamics simulation of this toxin was performed to examine the dynamic structural properties which may play a role in this inhibition. Conformational fluctuations of the FAS2 protein were examined by a variety of techniques to identify flexible residues and determine their characteristic motion. The tips of the toxin "finger" loops and the turn connecting loops I and II were found to fluctuate, while the rest of the protein remained fairly rigid throughout the simulation. Finally, the structural fluctuations were compared to NMR data of fluctuations on a similar timescale in a related three-finger toxin. The molecular dynamics results were in good qualitative agreement with the experimental measurements. Proteins 1999;36:447-453.

Molecular dynamics of mouse acetylcholinesterase complexed with huperzine A.

Two molecular dynamics simulations were performed for a modeled complex of mouse acetylcholinesterase liganded with huperzine A (HupA). Analysis of these simulations shows that HupA shifts in the active site toward Tyr 337 and Phe 338, and that several residues in the active site area reach out to make hydrogen bonds with the inhibitor. Rapid fluctuations of the gorge width are observed, ranging from widths that allow substrate access to the active site, to pinched structures that do not allow access of molecules as small as water. Additional openings or channels to the active site are found. One opening is formed in the side wall of the active site gorge by residues Val 73, Asp 74, Thr 83, Glu 84, and Asn 87. Another opening is formed at the base of the gorge by residues Trp 86, Val 132, Glu 202, Gly 448, and Ile 451. Both of these openings have been observed separately in the Torpedo californica form of the enzyme. These channels could allow transport of waters and ions to and from the bulk solution.

Shedding light on the dark and weakly fluorescent states of green fluorescent proteins.

Recent experiments on various similar green fluorescent protein (GFP) mutants at the single-molecule level and in solution provide evidence of previously unknown short- and long-lived "dark" states and of related excited-state decay channels. Here, we present quantum chemical calculations on cis-trans photoisomerization paths of neutral, anionic, and zwitterionic GFP chromophores in their ground and first singlet excited states that explain the observed behaviors from a common perspective. The results suggest that favorable radiationless decay channels can exist for the different protonation states along these isomerizations, which apparently proceed via conical intersections. These channels are suggested to rationalize the observed dramatic reduction of fluorescence in solution. The observed single-molecule fast blinking is attributed to conversions between the fluorescent anionic and the dark zwitterionic forms whereas slow switching is attributed to conversions between the anionic and the neutral forms. The predicted nonadiabatic crossings are seen to rationalize the origins of a variety of experimental observations on a common basis and may have broad implications for photobiophysical mechanisms in GFP.

Determinants of ligand binding to cAMP-dependent protein kinase.

Protein kinases are essential for the regulation of cellular growth and metabolism. Since their dysfunction leads to debilitating diseases, they represent key targets for pharmaceutical research. The rational design of kinase inhibitors requires an understanding of the determinants of ligand binding to these proteins. In the present study, a theoretical model based on continuum electrostatics and a surface-area-dependent nonpolar term is used to calculate binding affinities of balanol derivatives, H-series inhibitors, and ATP analogues toward the catalytic subunit of cAMP-dependent protein kinase (cAPK or protein kinase A). The calculations reproduce most of the experimental trends and provide insight into the driving forces responsible for binding. Nonpolar interactions are found to govern protein-ligand affinity. Hydrogen bonds represent a negligible contribution, because hydrogen bond formation in the complex requires the desolvation of the interacting partners. However, the binding affinity is decreased if hydrogen-bonding groups of the ligand remain unsatisfied in the complex. The disposition of hydrogen-bonding groups in the ligand is therefore crucial for binding specificity. These observations should be valuable guides in the design of potent and specific kinase inhibitors.

Influence of key residues on the reaction mechanism of the cAMP-dependent protein kinase.

The reaction mechanism of the catalytic phosphoryl transfer of cAMP-dependent protein kinase (cAPK) was investigated by semi-empirical AM1 molecular orbital computations of an active site model system derived from the crystal structure of the catalytic subunit of the enzyme. The activation barrier is calculated as 20.7 kcal mol(-1) and the reaction itself to be exothermic by 12.2 kcal mol(-1). The active site residue Asp166, which was often proposed to act as a catalytic base, does not accept a proton in any of the reaction steps. Instead, the hydroxyl hydrogen of serine is shifted to the simultaneously transferred phosphate group of ATP. Although the calculated transition state geometry indicates an associative phosphoryl transfer, no concentration of negative charge is found. To study the influence of protein mutations on the reaction mechanism, we compared two-dimensional energy hypersurfaces of the protein kinase wild-type model and a corresponding mutant in which Asp166 was replaced by alanine. Surprisingly, they show similar energy profiles despite the experimentally known decrease of catalytic activity for corresponding mutants. Furthermore, a model structure was examined, where the charged NH3 group of Lys168 was replaced by a neutral methyl group. The energetic hypersurface of this hypothetical mutant shows two possible pathways for phosphoryl transfer, which both require significantly higher activation energies than the other systems investigated, while the energetic stabilization of the reaction product is similar in all systems. As the position of the amino acid side chains and the substrate peptide is virtually unchanged in all model systems, our results suggest that the exchange of Asp166 by other amino acid is less important to the phosphoryl transfer itself, but crucial to maintain the configuration of the active site in vivo. The positively charged side chain of Lys168, however, is necessary to stabilize the intermediate reaction states, particularly the side chain of the substrate peptide.

Hydration energy landscape of the active site cavity in cytochrome P450cam.

Hydration of protein cavities influences protein stability, dynamics, and function. Protein active sites usually contain water molecules that, upon ligand binding, are either displaced into bulk solvent or retained to mediate protein-ligand interactions. The contribution of water molecules to ligand binding must be accounted for to compute accurate values of binding affinities. This requires estimation of the extent of hydration of the binding site. However, it is often difficult to identify the water molecules involved in the binding process when ligands bind on the surface of a protein. Cytochrome P450cam is, therefore, an ideal model system because its substrate binds in a buried active site, displacing partially disordered solvent, and the protein is well characterized experimentally. We calculated the free energy differences for having five to eight water molecules in the active site cavity of the unliganded enzyme from molecular dynamics simulations by thermodynamic integration employing a three-stage perturbation scheme. The computed free energy differences between the hydration states are small (within 12 kJ mol-1) but distinct. Consistent with the crystallographic determination and studies employing hydrostatic pressure, we calculated that, although ten water molecules could in principle occupy the volume of the active site, occupation by five to six water molecules is thermodynamically most favorable.

Kinase conformations: a computational study of the effect of ligand binding.

Protein function is often controlled by ligand-induced conformational transitions. Yet, in spite of the increasing number of three-dimensional crystal structures of proteins in different conformations, not much is known about the driving forces of these transitions. As an initial step toward exploring the conformational and energetic landscape of protein kinases by computational methods, intramolecular energies and hydration free energies were calculated for different conformations of the catalytic domain of cAMP-dependent protein kinase (cAPK) with a continuum (Poisson) model for the electrostatics. Three protein kinase crystal structures for ternary complexes of cAPK with the peptide inhibitor PKI(5-24) and ATP or AMP-PNP were modeled into idealized intermediate and open conformations. Concordant with experimental observation, we find that the binding of PKI(5-24) is more effective in stabilizing the closed and intermediate forms of cAPK than ATP. PKI(5-24) seems to drive the final closure of the active site cleft from intermediate to closed state because ATP does not distinguish between these two states. Binding of PKI(5-24) and ATP is energetically additive.

Structural changes in cytochrome P-450cam effected by the binding of the enantiomers (1R)-camphor and (1S)-camphor.

A comparative study of the enantiomeric substrate [(1R)-camphor- and (1S)-camphor)-bound cytochrome P-450cam concerns the spin-state equilibrium, substrate dissociation, the thermal unfolding of the protein structure, and the subconformer equilibria observed in the infrared spectra of the carbon monoxide (CO) complex of cytochrome P-450cam. The behavior of the different conformational equilibria in dependence on temperature, pressure, pH-value, cosolvent, and cation binding led us to suggest that (1S)-camphor is more loosely and less optimally bound in the heme pocket, which facilitates the access of solvent molecules into the heme-iron environment. The spin reaction volume difference measured using the high pressure technique is smaller by 16 +/- 9 cm3/mol for (1S)-camphor-bound P-450cam compared to the (1R)-camphor-bound P-450cam, which might indicate a higher water content in the protein and in the heme environment in the (1S)-camphor complex. The half-transition temperature of the thermal unfolding of 53.8 degrees C for the (1S)-camphor-bound oxidized cytochrome P-450cam is one degree lower than the value for the (1R)-camphor-bound protein (54.8 degrees C). In the reduced, CO-bound form of cytochrome P-450cam at 290 K the (1S)-camphor complex reveals another CO stretch vibration population distribution with slightly higher frequencies [1940.2 cm-1 (major band) and 1946.3 cm-1 (minor band)] compared to the (1R)-camphor complex [1939.7 cm-1 (major band) and 1930 cm-1 (minor band)]. A loosening of the contact between the iron-bound CO ligand and amino acids of the I-helix, probably induced by compensating effects of the increased water content, is suggested. Assuming the carbon monoxide complex as a model for the dioxygen complex, the more loosened binding of (1S)-camphor, therefore the increased water accessibility, and the weaker contact of the iron ligand to the I-helix might explain the higher amount of uncoupling of the cytochrome P-450 reaction cycle compared to that when (1R)-camphor is used as substrate.

Improved binding of cytochrome P450cam substrate analogues designed to fill extra space in the substrate binding pocket.

Cytochrome P450cam catalyzes the 5-exo-hydroxylation of camphor. Camphor analogues were designed to fill an empty region of the substrate binding pocket with the expectation that they would bind more tightly than camphor itself due to increased van der Waals interactions with the protein and the displacement of any solvent occupying this site. A series of compounds (endo-borneol methyl ether, endo-borneol propyl ether, endo-borneol allyl ether and endo-borneol dimethyl allyl ether) were synthesized with substituents at the camphor carbonyl oxygen. The spin conversion and thermodynamic properties of this series of compounds were measured for wild type and Y96F mutant cytochrome P450cam and were interpreted in the context of molecular dynamics simulations of the camphor analogues in the P450 binding site and in solution. Compounds with a 3-carbon chain substituent were predicted to match the size of the unoccupied region most optimally and thus bind best. Consistent with this prediction, the borneol allyl ether binds to cytochrome P450cam with highest affinity with a Kd = 0.6 +/- 0.1 microM (compared to a Kd = 1.7 +/- 0.2 microM for camphor under the same experimental conditions). Binding of the camphor analogues to the Y96F mutant is much enhanced over the binding of camphor, indicating that hydrogen bonding plays a less important role in binding of these analogues. Binding enthalpies calculated from the simulations, taking all solvent contributions into account, agree very well with experimental binding enthalpies. Binding affinity is not however correlated with the calculated binding enthalpy because the binding of the substrate analogues is characterized by enthalpy-entropy compensation. The new compounds are useful probes for further studies of the mechanism of cytochrome P450cam due to their high binding affinities and high spin properties.

Thermodynamics of water mediating protein-ligand interactions in cytochrome P450cam: a molecular dynamics study.

Ordered water molecules are observed by crystallography and nuclear magnetic resonance to mediate protein-ligand interactions. Here, we examine the energetics of hydrating cavities formed at protein-ligand interfaces using molecular dynamics simulations. The free energies of hydrating two cavities in the active site of two liganded complexes of cytochrome P450cam were calculated by multiconfigurational thermodynamic integration. The complex of cytochrome P450cam with 2-phenyl-imidazole contains a crystallographically well defined water molecule mediating hydrogen bonds between the protein and the inhibitor. We calculate that this water molecule is stabilized by a binding free energy of -11.6 +/- kJ/mol. The complex of cytochrome P450cam with its natural substrate, camphor, contains a cavity that is empty in the crystal structure although a water molecule in it could make a hydrogen bond to camphor. Here, solvation of this cavity is calculated to be unfavorable by +15.8 +/- 5.0 kJ/mol. The molecular dynamics simulations can thus distinguish a hydrated interfacial cavity from an empty one. They also provide support for the notion that protein-ligand complexes can accommodate empty interfacial cavities and that such cavities are likely to be unhydrated unless more than one hydrogen bond can be made to a water molecule in the cavity.