A Nonorthogonal Configuration Interaction Approach to Singlet Fission in Perylenediimide Compounds

Perylenediimide molecules constitute a family of chromophores that undergo singlet fission, a process in which an excited singlet state converts into lower energy triplets on two neighboring molecules, potentially increasing the efficiency of organic solar cells. Here, the nonorthogonal configuration interaction method is applied to study the effect of the different crystal packing of various perylenediimide derivatives on the relative energies of the singlet and triplet states, the intermolecular electronic couplings, and the relative rates for singlet fission. The analysis of the wave functions and electronic couplings reveals that charge transfer states play an important role in the singlet fission mechanism. Dimer conformations where the PDI molecules are at large displacements along the long axis and short on the short axis are posed as the most favorable for singlet fission. The role of the substituent at the imide group has been inspected concluding that, although it has no effect in the energies, for some conformations it significantly influences the electronic couplings, and therefore, replacing this substituent with hydrogen may introduce artifacts in the computational modeling of the PDI molecules.

GronOR: Scalable and Accelerated Nonorthogonal Configuration Interaction for Molecular Fragment Wave Functions

GronOR is a program package for nonorthogonal configuration interaction calculations. Electronic wave functions are constructed in terms of antisymmetrized products of multiconfiguration molecular fragment wave functions. The computational complexity of the nonorthogonal methodologies implemented in GronOR applied to large molecular assemblies requires a design that takes full advantage of massively parallel supercomputer architectures and accelerator technologies. This work describes the implementation strategy and resulting performance characteristics. In addition to parallelization and acceleration, the software development strategy includes aspects of fault resiliency and heterogeneous computing. The program was designed for large-scale supercomputers but also runs effectively on small clusters and workstations for small molecular systems. GronOR is available as open source to the scientific community.

On the role of dynamic electron correlation in non-orthogonal configuration interaction with fragments

Two different approaches have been implemented to include the effect of dynamic electron correlation in the Non-Orthogonal Configuration Interaction for Fragments (NOCI-F) method. The first is based on shifting the diagonal matrix elements of the NOCI matrix, while the second incorporates the dynamic correlation explicitly in the fragment wave functions used to construct the many-electron basis functions of the NOCI. The two approaches are illustrated for the calculation of the electronic coupling relevant in singlet fission and the coupling of spin moments in organic radicals. Comparison of the calculated diabatic couplings, the NOCI energies and wave functions shows that dynamic electron correlation is not only efficiently but also effectively incorporated by the shifting approach and can largely affect the coupling between electronic states. Also, it brings the NOCI coupling of the spin moments in close agreement with benchmark calculations.

From NWChem to NWChemEx: Evolving with the Computational Chemistry Landscape

Since the advent of the first computers, chemists have been at the forefront of using computers to understand and solve complex chemical problems. As the hardware and software have evolved, so have the theoretical and computational chemistry methods and algorithms. Parallel computers clearly changed the common computing paradigm in the late 1970s and 80s, and the field has again seen a paradigm shift with the advent of graphical processing units. This review explores the challenges and some of the solutions in transforming software from the terascale to the petascale and now to the upcoming exascale computers. While discussing the field in general, NWChem and its redesign, NWChemEx, will be highlighted as one of the early codesign projects to take advantage of massively parallel computers and emerging software standards to enable large scientific challenges to be tackled.

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.

GronOR: Massively parallel and GPU-accelerated non-orthogonal configuration interaction for large molecular systems

GronOR is a program package for non-orthogonal configuration interaction calculations for an electronic wave function built in terms of anti-symmetrized products of multi-configuration molecular fragment wave functions. The two-electron integrals that have to be processed may be expressed in terms of atomic orbitals or in terms of an orbital basis determined from the molecular orbitals of the fragments. The code has been specifically designed for execution on distributed memory massively parallel and Graphics Processing Unit (GPU)-accelerated computer architectures, using an MPI+OpenACC/OpenMP programming approach. The task-based execution model used in the implementation allows for linear scaling with the number of nodes on the largest pre-exascale architectures available, provides hardware fault resiliency, and enables effective execution on systems with distinct central processing unit-only and GPU-accelerated partitions. The code interfaces with existing multi-configuration electronic structure codes that provide optimized molecular fragment orbitals, configuration interaction coefficients, and the required integrals. Algorithm and implementation details, parallel and accelerated performance benchmarks, and an analysis of the sensitivity of the accuracy of results and computational performance to thresholds used in the calculations are presented.

Retention of conformational entropy upon calmodulin binding to target peptides is driven by transient salt bridges

Calmodulin (CaM) is a highly flexible calcium-binding protein that mediates signal transduction through an ability to differentially bind to highly variable binding sequences in target proteins. To identify how binding affects CaM motions, and its relationship to conformational entropy and target peptide sequence, we have employed fully atomistic, explicit solvent molecular dynamics simulations of unbound CaM and CaM bound to five different target peptides. The calculated CaM conformational binding entropies correlate with experimentally derived conformational entropies with a correlation coefficient R(2) of 0.95. Selected side-chain interactions with target peptides restrain interhelical loop motions, acting to tune the conformational entropy of the bound complex via widely distributed CaM motions. In the complex with the most conformational entropy retention (CaM in complex with the neuronal nitric oxide synthase binding sequence), Lys-148 at the C-terminus of CaM forms transient salt bridges alternating between Glu side chains in the N-domain, the central linker, and the binding target. Additional analyses of CaM structures, fluctuations, and CaM-target interactions illuminate the interplay between electrostatic, side chain, and backbone properties in the ability of CaM to recognize and discriminate against targets by tuning its conformational entropy, and suggest a need to consider conformational dynamics in optimizing binding affinities.

Force-Field Development and Molecular Dynamics of [NiFe] Hydrogenase

Classical molecular force-field parameters describing the structure and motion of metal clusters in [NiFe] hydrogenase enzymes can be used to compare the dynamics and thermodynamics of [NiFe] under different oxidation, protonation, and ligation circumstances. Using density functional theory (DFT) calculations of small model clusters representative of the active site and the proximal, medial, and distal Fe/S metal centers and their attached protein side chains, we have calculated classical force-field parameters for [NiFe] in reduced and oxidized states, including internal coordinates, force constants, and atom-centered charges. Derived force constants revealed that cysteinate ligands bound to the metal ions are more flexible in the Ni-B active site, which has a bridging hydroxide ligand, than in the Ni-C active site, which has a bridging hydride. Ten nanosecond all-atom, explicit-solvent MD simulations of [NiFe] hydrogenase in oxidized and reduced catalytic states established the stability of the derived force-field parameters in terms of Cα and metal cluster fluctuations. Average active site structures from the protein MD simulations are consistent with [NiFe] structures from the Protein Data Bank, suggesting that the derived force-field parameters are transferrable to other hydrogenases beyond the structure used for testing. A comparison of experimental H2-production rates demonstrated a relationship between cysteinate side chain rotation and activity, justifying the use of a fully dynamic model of [NiFe] metal cluster motion.

Characterization of the outer membrane protein OprF of Pseudomonas aeruginosa in a lipopolysaccharide membrane by computer simulation

The N-terminal domain of outer membrane protein OprF of Pseudomonas aeruginosa forms a membrane spanning eight-stranded antiparallel beta-barrel domain that folds into a membrane channel with low conductance. The structure of this protein has been modeled after the crystal structure of the homologous protein OmpA of Escherichia coli. A number of molecular dynamics simulations have been carried out for the homology modeled structure of OprF in an explicit molecular model for the rough lipopolysaccharide (LPS) outer membrane of P. aeruginosa. The structural stability of the outer membrane model as a result of the strong electrostatic interactions compared with simple lipid bilayers is restricting both the conformational flexibility and the lateral diffusion of the porin in the membrane. Constricting side-chain interactions within the pore are similar to those found in reported simulations of the protein in a solvated lipid bilayer membrane. Because of the strong interactions between the loop regions of OprF and functional groups in the saccharide core of the LPS, the entrance to the channel from the extracellular space is widened compared with the lipid bilayer simulations in which the loops are extruding in the solvent. The specific electrostatic signature of the LPS membrane, which results in a net intrinsic dipole across the membrane, is found to be altered by the presence of OprF, resulting in a small electrically positive patch at the position of the channel.

Electronic Coupling between Heme Electron-Transfer Centers and Its Decay with Distance Depends Strongly on Relative Orientation

A method for calculating the electron-transfer matrix element V(RP) using density functional theory Kohn-Sham orbitals is presented and applied to heme dimers of varying relative orientation. The electronic coupling decays with increased iron separation according to V(RP) = V(0)(RP)exp(-beta r/2) with a distance dependence parameter beta approximately 2 A(-1) for hemes with parallel porphyrins and either 1.1 or 4.0 A(-1) when the porphyrin planes are perpendicular, depending on the alignment of the iron d(pi) orbital. These findings are used to interpret the observed orientation of the hemes in tetraheme redox proteins such as Flavocytochrome c(3) fumarate reductase (Ifc(3), PDB code 1QJD) of Shewanella frigidimarina, another flavocytochrome from the same bacterium (Fcc(3), 1E39) and a small tetraheme cytochrome of Shewanella oneidensis strain MR1 (1M1P). Our results show that shifting and rotating the hemes controls the adiabaticity of the three electron hopping steps.

Molecular docking of balanol to dynamics snapshots of protein kinase A

Even if the structure of a receptor has been determined experimentally, it may not be a conformation to which a ligand would bind when induced fit effects are significant. Molecular docking using such a receptor structure may thus fail to recognize a ligand to which the receptor can bind with reasonable affinity. Here, we examine one way to alleviate this problem by using an ensemble of receptor conformations generated from a molecular dynamics simulation for molecular docking. Two molecular dynamics simulations were conducted to generate snapshots for protein kinase A: one with the ligand bound, the other without. The ligand, balanol, was then docked to conformations of the receptors presented by these trajectories. The Lamarckian genetic algorithm in Autodock [Goodsell et al. J Mol Recognit 1996;9(1):1-5; Morris et al. J Comput Chem 1998;19(14):1639-1662] was used in the docking. Three ligand models were used: rigid, flexible, and flexible with torsional potentials. When the snapshots were taken from the molecular dynamics simulation of the protein-ligand complex, the correct docking structure could be recovered easily by the docking algorithm in all cases. This was an easier case for challenging the docking algorithm because, by using the structure of the protein in a protein-ligand complex, one essentially assumed that the protein already had a pocket to which the ligand can fit well. However, when the snapshots were taken from the ligand-free protein simulation, which is more useful for a practical application when the structure of the protein-ligand complex is not known, several clusters of structures were found. Of the 10 docking runs for each snapshot, at least one structure was close to the correctly docked structure when the flexible-ligand models were used. We found that a useful way to identify the correctly docked structure was to locate the structure that appeared most frequently as the lowest energy structure in the docking experiments to different snapshots.

Characterization of the active site of DNA polymerase beta by molecular dynamics and quantum chemical calculation

It is well established that the fully formed polymerase active site of the DNA repair enzyme, polymerase beta (pol beta), including two bound Mg2+ cations and the nucleoside triphosphate (dNTP) substrate, exists at only one point in the catalytic cycle just prior to the chemical nucleotidyl transfer step. The structure of the active conformation has been the subject of much interest as it relates to the mechanism of the chemical step and also to the question of fidelity assurance. Although crystal structures of ternary pol beta-(primer-template) DNA-dNTP complexes have provided the main structural features of the active site, they are necessarily incomplete due to intentional alterations (e.g., removal of the 3'OH groups from primer and substrate) needed to obtain a structure from midcycle. Working from the crystal structure closest to the fully formed active site [Protein Data Bank (PDB) code: 1bpy], two molecular dynamics (MD) simulations of the solvated ternary complex were performed: one with the missing 3'OHs restored, via modeling, to the primer and substrate, and the other without restoration of the 3'OHs. The results of the simulations, together with ab initio optimizations on simplified active-site models, indicate that the missing primer 3'OH in the crystal structure is responsible for a significant perturbation in the coordination sphere of the catalytic cation and allow us to suggest several corrections and additions to the active-site structure as observed by crystallography. In addition, the calculations help to resolve questions raised regarding the protonation states of coordinating ligands.

8-Oxoguanine enhances bending of DNA that favors binding to glycosylases

Molecular dynamics (MD) simulations were carried out on the DNA oligonucleotide GGGAACAACTAG:CTAGTTGTTCCC in its native form and with guanine in the central G(19):C(6) base pair replaced by 8-oxoguanine (8oxoG). A box of explicit water molecules was used for solvation, and Na(+) counterions were added to neutralize the system. The direction and magnitude of global bending were assessed by a technique used previously to analyze simulations of DNA containing a thymine dimer. The presence of 8oxoG did not greatly affect the magnitude of DNA bending; however, bending into the major groove was significantly more probable when 8oxoG replaced G(19). Crystal structures of glycosylases bound to damaged-DNA substrates consistently show a sharp bend into the major groove at the damage site. We conclude that changes in bending dynamics that assist the formation of this kink are a part of the mechanism by which glycosylases of the base excision repair pathway recognize the presence of 8oxoG in DNA.

Characterization of electronic structure and properties of a Bis(histidine) heme model complex

Ferric and ferrous hemes, such as those present in electron transfer proteins, often have low-lying spin states that are very close in energy. To explore the relationship between spin state, geometry, and cytochrome electron transfer, we investigate, using density functional theory, the relative energies, electronic structure, and optimized geometries for a high- and low-spin ferric and ferrous heme model complex. Our model consists of an iron-porphyrin axially ligated by two imidazoles, which model the interaction of a heme with histidine residues. Using the B3LYP hybrid functional, we found that, in the ferric model heme complex, the doublet is lower in energy than the sextet by 8.4 kcal/mol and the singlet ferrous heme is 6.7 kcal/mol more stable than the quintet. The difference between the high-spin ferric and ferrous model heme energies yields an adiabatic electron affinity (AEA) of 5.24 eV, and the low-spin AEA is 5.17 eV. Both values are large enough to ensure electron trapping, and electronic structure analysis indicates that the iron d(pi) orbital is involved in the electron transfer between hemes. Mössbauer parameters calculated to verify the B3LYP electronic structure correlate very well with experimental values. Isotropic hyperfine coupling constants for the ligand nitrogen atoms were also evaluated. The optimized geometries of the ferric and ferrous hemes are consistent with structures from X-ray crystallography and reveal that the iron-imidazole distances are significantly longer in the high-spin hemes, which suggests that the protein environment, modeled here by the imidazoles, plays an important role in regulating the spin state. Iron-imidazole dissociation energies, force constants, and harmonic frequencies were calculated for the ferric and ferrous low-spin and high-spin hemes. In both the ferric and the ferrous cases, a single imidazole ligand is more easily dissociated from the high-spin hemes.

Molecular basis for microbial adhesion to geochemical surfaces: computer simulation of Pseudomonas aeruginosa adhesion to goethite

The adhesion of Pseudomonas aeruginosa to the goethite mineral is investigated using classical molecular simulation. A fragment model for goethite has been integrated into a fully atomistic membrane model. Properties for the resulting system are evaluated for a 1.5-ns simulation in the isothermal-isobaric ensemble. The response of the membrane to the presence of the mineral is investigated. Radial distribution functions are used to present an average picture of the hydrogen bonding. Orientational vectors, assigned to the saccharide groups, reveal the extent of the mineral's perturbations on the membrane. Significant structural changes were observed for the outermost saccharide groups, several of which rotate to form hydrogen bonds with the mineral surface. The structure of the inner core, and the corresponding integrity of the membrane, is maintained. The mineral surface dehydrates slightly in the presence of the membrane as saccharide hydroxyl groups compete with water molecules for hydrogen-bonding sites on its surface.

Molecular structure of the outer bacterial membrane of Pseudomonas aeruginosa via classical simulation

A detailed structural analysis has been performed of the outer bacterial membrane of Pseudomonas aeruginosa using a parameterized classical simulation model (R. D. Lins and T. P. Straatsma, Biophysical Journal, 2001, Vol. 81, pp. 1037-1046) with modest modifications. The structural analysis of the membrane is presented and newly discovered characteristics of the membrane are discussed. Simulations indicate that the relative contribution of different ligands to calcium ion coordination varies across the membrane, while maintaining a constant average coordination number of 6.1. Water penetrates the surface of the membrane to a depth of about 30 A. The hydration of ions and phosphate groups is shown to depend on location within the membrane. A measure of saccharide residue orientation is defined and average orientations are presented. Saccharide residues possess varying degrees of motion with a trend of greater mobility at the membrane surface. However, their motion is limited and even in the membrane outer core region the average structure appears fairly rigid over a period of 1 ns.

Internal dynamics and ionization states of the macrophage migration inhibitory factor: comparison between wild-type and mutant forms

The macrophage migration inhibitory factor (MIF) is a cytokine that shares a common structural architecture and catalytic strategy with three isomerases: 4-oxalocrotonate tautomerase, 5-carboxymethyl-2-hydroxymuconate isomerase, and D-dopachrome tautomerase. A highly conserved N-terminal proline acts as a base-acid during the proton transfer reaction catalyzed by these enzymes. Such unusual catalytic strategy appears to be possible only due to the N-terminal proline pK(a) shifted to 5.0-6.0 units. Mutations of this residue result in a significant decrease of the catalytic activity of MIF. Two hypotheses have been proposed to explain the catalytic inefficiency of MIF: the lower basicity of primary amines with regard to secondary ones and the increased flexibility resulting from the replacement of a proline by residues like glycine. To investigate that, we have performed molecular dynamics simulations of MIF wild-type and its mutant P1G, as well as calculated the protonation properties of several mutant forms. It was found that the N-terminal glycine does not show larger fluctuations compared to proline, but the former residue is more exposed to the solvent throughout the simulations. The apparent pK(a) of these residues displays very little change (as expected from the structural rigidity of MIF) and is not significantly affected by the surrounding ionizable residues. Instead, the hydrophobic character of the active site seems to be the main factor in determining the pKa of the N-terminal residue and the catalytic efficiency of MIF.

Revisiting the structural flexibility of the complex p21(ras)-GTP: the catalytic conformation of the molecular switch II

The hydrolysis of GTP in p21(ras) triggers conformational changes that regulate the ras/ERK signaling pathway. An important active site residue is Gln61, which has been found to be mutated in 30% of human tumors. The dynamics of the active site conformation is studied by using molecular dynamics simulation of two independent structures of the GTP-bound uncomplexed enzyme. Two distinct conformations of the enzyme are observed, in which the side-chain residue Gln61 is in different orientations. Essential dynamics analysis is used to describe the essential motions in the transition between the two conformations. Results are compared with earlier simulations of p21(ras) and its complex with GTPase activating protein p21-GAP.

Computer simulation of the rough lipopolysaccharide membrane of Pseudomonas aeruginosa

Lipopolysaccharides (LPSs) form the major constituent of the outer membrane of Gram-negative bacteria, and are believed to play a key role in processes that govern microbial metal binding, microbial adsorption to mineral surfaces, and microbe-mediated oxidation/reduction reactions at the bacterial exterior surface. A computational modeling capability is being developed for the study of geochemical reactions at the outer bacterial envelope of Gram-negative bacteria. A molecular model for the rough LPS of Pseudomonas aeruginosa has been designed based on experimentally determined structural information. An electrostatic model was developed based on Hartree-Fock SCF calculations of the complete LPS molecule to obtain partial atomic charges. The exterior of the bacterial membrane was assembled by replication of a single LPS molecule and a single phospholipid molecule. Molecular dynamics simulations of the rough LPS membrane of P. aeruginosa were carried out and trajectories were analyzed for the energetic and structural factors that determine the role of LPS in processes at the cell surface.

The arginine finger of RasGAP helps Gln-61 align the nucleophilic water in GAP-stimulated hydrolysis of GTP

The Ras family of GTPases is a collection of molecular switches that link receptors on the plasma membrane to signaling pathways that regulate cell proliferation and differentiation. The accessory GTPase-activating proteins (GAPs) negatively regulate the cell signaling by increasing the slow intrinsic GTP to GDP hydrolysis rate of Ras. Mutants of Ras are found in 25-30% of human tumors. The most dramatic property of these mutants is their insensitivity to the negative regulatory action of GAPs. All known oncogenic mutants of Ras map to a small subset of amino acids. Gln-61 is particularly important because virtually all mutations of this residue eliminate sensitivity to GAPs. Despite its obvious importance for carcinogenesis, the role of Gln-61 in the GAP-stimulated GTPase activity of Ras has remained a mystery. Our molecular dynamics simulations of the p21ras-p120GAP-GTP complex suggest that the local structure around the catalytic region can be different from that revealed by the x-ray crystal structure. We find that the carbonyl oxygen on the backbone of the arginine finger supplied in trans by p120GAP (Arg-789) interacts with a water molecule in the active site that is forming a bridge between the NH(2) group of the Gln-61 and the gamma-phosphate of GTP. Thus, Arg-789 may play a dual role in generating the nucleophile as well as stabilizing the transition state for PO bond cleavage.

Similarities in the HIV-1 and ASV integrase active sites upon metal cofactor binding

The HIV-1 integrase, which is essential for viral replication, catalyzes the insertion of viral DNA into the host chromosome thereby recruiting host cell machinery into making viral proteins. It represents the third main HIV enzyme target for inhibitor design, the first two being the reverse transcriptase and the protease. We report here a fully hydrated 2 ns molecular dynamics simulation performed using parallel NWChem3.2.1 with the AMBER95 force field. The HIV-1 integrase catalytic domain previously determined by crystallography (1B9D) and modeling including two Mg(2+) ions placed into the active site based on an alignment against an ASV integrase structure containing two divalent metals (1VSH), was used as the starting structure. The simulation reveals a high degree of flexibility in the region of residues 140-149 even in the presence of a second divalent metal ion and a dramatic conformational change of the side chain of E152 when the second metal ion is present. This study shows similarities in the behavior of the catalytic residues in the HIV-1 and ASV integrases upon metal binding. The present simulation also provides support to the hypothesis that the second metal ion is likely to be carried into the HIV-1 integrase active site by the substrate, a strand of DNA.

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.

Mouse acetylcholinesterase unliganded and in complex with huperzine A: a comparison of molecular dynamics simulations

A 1 ns molecular dynamics simulation of unliganded mouse acetylcholinesterase (AChE) is compared to a previous simulation of mouse AChE complexed with huperzine A (HupA). Several common features are observed. In both simulations, the active site gorge fluctuates in size during the 1 ns trajectory and is completely pinched off several times. Many of the residues in the gorge that formed hydrogen bonds with HupA in the simulation of the complex now form hydrogen bonds with other protein residues and water molecules in the gorge. The opening of a "backdoor" entrance to the active site that was found in the simulation of the complex is also observed in the unliganded simulation. Differences between the two simulations include overall lower structural rms deviations for residues in the gorge in the unliganded simulation, a smaller diameter of the gorge in the absence of HupA, and the disappearance of a side channel that was frequently present in the liganded simulation. The differences between the two simulations can be attributed, in part, to the interaction of AChE with HupA.

Molecular dynamics studies on the HIV-1 integrase catalytic domain

The HIV-1 integrase, which is essential for viral replication, catalyzes the insertion of viral DNA into the host chromosome, thereby recruiting host cell machinery into making viral proteins. It represents the third main HIV enzyme target for inhibitor design, the first two being the reverse transcriptase and the protease. Two 1-ns molecular dynamics simulations have been carried out on completely hydrated models of the HIV-1 integrase catalytic domain, one with no metal ions and another with one magnesium ion in the catalytic site. The simulations predict that the region of the active site that is missing in the published crystal structures has (at the time of this work) more secondary structure than previously thought. The flexibility of this region has been discussed with respect to the mechanistic function of the enzyme. The results of these simulations will be used as part of inhibitor design projects directed against the catalytic domain of the enzyme.

Theoretical study of inhibition of adenosine deaminase by (8R)-coformycin and (8R)-deoxycoformycin

Molecular dynamics and free energy simulations were performed to examine the binding of (8R)-deoxycoformycin and (8R)-coformycin to adenosine deaminase. The two inhibitors differ only at the 2' position of the sugar ring; the sugar moiety of conformycin is ribose, while it is deoxyribose for deoxycoformycin. The 100 ps molecular dynamics trajectories reveal that Asp 19 and His 17 interact strongly with the 5' hydroxyl group of the sugar moiety of both inhibitors and appear to play an important role in binding the sugar. The 2' and 3' groups of the sugars are near the protein-water interface and can be stabilized by either protein residues or water. The flexibility of the residues at the opening of the active site helps to explain the modest difference in binding of the two inhibitors and how substrates/inhibitors can enter an otherwise inaccessible binding site.

Binding of tacrine and 6-chlorotacrine by acetylcholinesterase

Multiconfiguration thermodynamic integration was used to determine the relative binding strength of tacrine and 6-chlorotacrine by Torpedo californica acetylcholinesterase. 6-Chlorotacrine appears to be bound stronger by 0.7+/-0.4 kcal/mol than unsubstituted tacrine when the active site triad residue His-440 is deprotonated. This result is in excellent agreement with experimental inhibition data on electric eel acetylcholinesterase. Electrostatic Poisson-Boltzmann calculations confirm that order of binding strength, resulting in deltaG of binding of -2.9 and -3.3 kcal/mol for tacrine and chlorotacrine, respectively, and suggest inhibitor binding does not occur when His-440 is charged. Our results suggest that electron density redistribution upon tacrine chlorination is mainly responsible for the increased attraction potential between pronated inhibitor molecule and adjacent aromatic groups of Phe-330 and Trp-84.

Open "back door" in a molecular dynamics simulation of acetylcholinesterase

The enzyme acetylcholinesterase generates a strong electrostatic field that can attract the cationic substrate acetylcholine to the active site. However, the long and narrow active site gorge seems inconsistent with the enzyme's high catalytic rate. A molecular dynamics simulation of acetylcholinesterase in water reveals the transient opening of a short channel, large enough to pass a water molecule, through a thin wall of the active site near tryptophan-84. This simulation suggests that substrate, products, or solvent could move through this "back door," in addition to the entrance revealed by the crystallographic structure. Electrostatic calculations show a strong field at the back door, oriented to attract the substrate and the reaction product choline and to repel the other reaction product, acetate. Analysis of the open back door conformation suggests a mutation that could seal the back door and thus test the hypothesis that thermal motion of this enzyme may open multiple routes of access to its active site.

Inversion of receptor binding preferences by mutagenesis: free energy thermodynamic integration studies on sugar binding to L-arabinose binding proteins

The Escherichia coli L-arabinose-binding protein (ABP) participates as a specific receptor in the transport of L-arabinose, D-fucose, and D-galactose through the periplasmic space. The wild-type protein binds L-arabinose about 40 times more strongly than D-fucose. A mutation of the protein at position 108 (Met-->Leu) causes a specificity change. The Met108Leu ABP slightly prefers binding of D-fucose over L-arabinose. Molecular dynamics (MD) and thermodynamic integration (TI) computer simulations were performed to study the mechanism of sugar discrimination and specificity change based on the known high-resolution X-ray structures. The specificity change was evaluated by calculating the difference in free energy of L-arabinose versus D-fucose bound to wild-type and Met108Leu ABP. The calculated free energy differences are consistent with the experimentally observed specificity of wild-type and Met108Leu ABP. The simulations indicate that the specificity change of Met108Leu ABP is accomplished mainly by reduced Lennard-Jones interactions of residue 108 with L-arabinose and improved interactions with D-fucose. In addition to MD/TI calculations on sugar binding, finite difference Poisson-Boltzmann calculations were performed to identify the most stable ionization state of buried ionizable residues in ABP.