Discovery of HIV-1 integrase inhibitors through a novel combination of ligand and structure-based drug design.

Over the past 10 years, classical computer-aided molecular design methods have not been frequently applied for the discovery of novel HIV-1 integrase (IN) inhibitors, due to the intrinsic challenges that this enzyme presents. Therefore, a novel approach that combines the chemical information of known integrase inhibitors with the enzyme's detailed 3D structure in a stepwise fashion is proposed: (I) use of a pharmacophore model (PM), which takes into account in a weighted fashion the chemical features of known ligands, in analogous manner to the to search the Maybridge and the NCI 3D databases; (II) drug-likeness optimization; (III) virtual high-throughput screening of the hits matching the PM query against 1QS4 wild-type IN structure using different Docking/Scoring combinations; (IV) visual inspection and selection of the hits in function of: binding free energies; binding mode type within the active site; retrieval among the best 20% hits in more than 6 Docking/Scoring protocols at the same time. This approach aims at a rational selection of new potential HIV-1 integrase inhibitors.

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.

A model for enzyme-substrate interaction in alanine racemase.

We report on a theoretical model for the complex of the enzyme alanine racemase with its natural substrate (L-alanine) and cofactor (pyridoxal 5'-phosphate). Electrostatic potentials were calculated and ionization states were predicted for all of the ionizable groups in alanine racemase. Some rather unusual charge states were predicted for certain residues. Tyr265' has an unusually low predicted pK(a) of 7.9 and at pH 7.0 has a predicted average charge of -0.37, meaning that 37% of the Tyr265' residues in an ensemble of enzyme molecules are in the phenolate form. At pH 8-9, the majority of Tyr265' side groups will be in the phenolate form. This lends support to the experimental evidence that Tyr265' is the catalytic base involved in the conversion of L-alanine to D-alanine. Residues Lys39 and Lys129 have predicted average charges of +0.91 and +0.14, respectively, at pH 7.0. Lys39 is believed to be the catalytic base for the conversion of D-alanine to L-alanine, and the present results show that, at least some of the time, it is in the unprotonated amine form and thus able to act as a base. Cys311', which is located very close to the active site, has an unusually low predicted pK(a) of 5.8 and at pH 7.0 has a predicted average charge of -0.72. The very low predicted charge for Lys129 is consistent with experimental evidence that it is carbamylated, since an unprotonated amine group is available to act as a Lewis base and form the carbamate with CO(2). Repeating the pK(a) calculations on the enzyme with Lys129 in carbamylated form predicts trends similar to those of the uncarbamylated enzyme. It appears that the enzyme has the ability to stabilize negative charge in the region of the active site. Implications for selective inhibitor design are discussed.

Brownian dynamics simulations of interaction between scorpion toxin Lq2 and potassium ion channel.

The association of the scorpion toxin Lq2 and a potassium ion (K(+)) channel has been studied using the Brownian dynamics (BD) simulation method. All of the 22 available structures of Lq2 in the Brookhaven Protein Data Bank (PDB) determined by NMR were considered during the simulation, which indicated that the conformation of Lq2 affects the binding between the two proteins significantly. Among the 22 structures of Lq2, only 4 structures dock in the binding site of the K(+) channel with a high probability and favorable electrostatic interactions. From the 4 candidates of the Lq2-K(+) channel binding models, we identified a good three-dimensional model of Lq2-K(+) channel complex through triplet contact analysis, electrostatic interaction energy estimation by BD simulation and structural refinement by molecular mechanics. Lq2 locates around the extracellular mouth of the K(+) channel and contacts the K(+) channel using its beta-sheet rather than its alpha-helix. Lys27, a conserved amino acid in the scorpion toxins, plugs the pore of the K(+) channel and forms three hydrogen bonds with the conserved residues Tyr78(A-C) and two hydrophobic contacts with Gly79 of the K(+) channel. In addition, eight hydrogen-bonds are formed between residues Arg25, Cys28, Lys31, Arg34 and Tyr36 of Lq2 and residues Pro55, Tyr78, Gly79, Asp80, and Tyr82 of K(+) channel. Many of them are formed by side chains of residues of Lq2 and backbone atoms of the K(+) channel. Thirteen hydrophobic contacts exist between residues Met29, Asn30, Lys31 and Tyr36 of Lq2 and residues Pro55, Ala58, Gly79, Asp80 and Tyr82 of the K(+) channel. These favorable interactions stabilize the association between the two proteins. These observations are in good agreement with the experimental results and can explain the binding phenomena between scorpion toxins and K(+) channels at the level of molecular structure. The consistency between the BD simulation and the experimental data indicates that our three-dimensional model of Lq2-K(+) channel complex is reasonable and can be used in further biological studies such as rational design of blocking agents of K(+) channels and mutagenesis in both toxins and K(+) channels.

Comparative molecular field analysis (coMFA) study of epothilones-tubulin depolymerization inhibitors: pharmacophore development using 3D QSAR methods.

A three-dimensional quantitative structure-activity relationship (3D QSAR) study has been carried out on epothilones based on comparative molecular field analyses (CoMFA) using a large data set of epothilone analogs, which are potent inhibitors of tubulin depolymerization. Microtubules, which are polymers of the a/beta-tubulin heterodimer, need to dissociate in order to form the mitotic spindle, a structure required for cell division. A rational pharmacophore searching method using 3D QSAR procedures was carried out and the results for the epothilones are described herein. One-hundred and sixty-six epothilone analogs and their depolymerization inhibition properties with tubulin were used as a training set. Over a thousand molecular field energies were generated and applied to generate the descriptors of QSAR equations. Using a genetic function algorithm (GFA) method, combined with a least square approach, multiple QSAR models were considered during the search for pharmacophore elements. Each GFA run resulted in 100 QSAR models, which were ranked according to their lack of fit (LOF) scores, with a total of 40 GFA runs having been performed. The 40 best QSAR equations from each run had adequate fitted correlation coefficients (R from 0.813 to 0.863) and were of sufficient statistical significance (F value from 7.2 to 10.9). The pharmacophore elements for epothilones were studied by investigating the hit frequency of descriptors (i.e. the sampling probabilities of grid points from the GFA studies) from the set of the 4,000 top scoring QSAR equations. By comparing the frequency with which each grid point appeared in the QSAR equations, three candidate regions in the epothilones were proposed to be pharmacophore elements. Two of them are completely compatible with the recent model proposed by Ojima et al. [Proc. Natl. Acad. Sci. USA, 96 (1999) 4256], however, one is quite different and is necessary to accurately predict the activities of all 166 epothilone molecules used in our training set. Finally, by visualizing the 35 most probable grid points, it was found that changes related to the C6, C7, C8, C 12, S20, and C21 atoms of the epothilones were highly correlated to their activity.

Metabolic complications of HIV and AIDs.

Numerous anecdotal reports, clinical observations, and newly published studies provide evidence of the growing interest and concern about body shape changes and metabolic complications seen in HIV infection. At first believed to be a single complication caused by the protease inhibitor class of antiretroviral therapy, the focus of research is shifting to distinct syndromes with multiple causes. In this unfolding story, peripheral fat atrophy, central fat accumulation, dyslipidemia, and glucose disregulation characterize the more commonly recognized syndromes. Osteoporosis and osteopenia have been recently observed. While the etiologies await discovery, the long-term consequences of these metabolic changes demand the expertise of clinicians not formerly considered "front-line" in HIV/AIDS treatment, such as orthopaedic nurses.

Effects of pore mutations and permeant ion concentration on the spontaneous gating activity of OmpC porin.

Porins are trimers of beta-barrels that form channels for ions and other hydrophilic solutes in the outer membrane of Gram-negative bacteria. The X-ray structures of OmpF and PhoE show that each monomeric pore is constricted by an extracellular loop that folds into the channel vestibule, a motif that is highly conserved among bacterial porins. Electrostatic calculations have suggested that the distribution of ionizable groups at the constriction zone (or eyelet) may establish an intrinsic transverse electrostatic field across the pore, that is perpendicular to the pore axis. In order to study the role that electrostatic interactions between pore residues may have in porin function, we used spontaneous mutants and engineered site-directed mutants that have an altered charge distribution at the eyelet and compared their electrophysiological behavior with that of wild-type OmpC. We found that some mutations lead to changes in the spontaneous gating activity of OmpC porin channels. Changes in the concentration of permeant ions also altered this activity. These results suggest that the ionic interactions that exist between charged residues at the constriction zone of porin may play a role in the transitions between the channel's closed and open states.

Ionization state and molecular docking studies for the macrophage migration inhibitory factor: the role of lysine 32 in the catalytic mechanism.

The macrophage migration inhibitory factor (MIF) is a cytokine that is structurally similar to certain isomerases and for which multiple immune and catalytic roles have been proposed. Different catalytic activities have been reported for MIF, yet the exact mechanism by which MIF acts is not completely known. As a tautomerase, the enzyme uses a general acid-base mechanism of proton transfer in which the amino-terminal proline has been shown to function as the catalytic base. We report the results of molecular docking simulations of macrophage migration inhibitory factor with three substrates, D-dopachrome, L-dopachrome methyl ester and p-(hydroxyphenyl)pyruvate. Electrostatic pK(a) predictions were also performed for the free and complexed forms of the enzyme. The predicted binding mode of p-(hydroxyphenyl)pyruvate is in agreement with the recently published X-ray structure. A model for the binding mode of D-dopachrome and L-dopachrome methyl ester to MIF is proposed which offers insights into the catalytic mechanism of D-dopachrome tautomerase activity of MIF. The proposed catalytic mechanism is further supported by the pK(a) predictions, which suggest that residue Lys32 acts as the general acid for the enzymatic catalysis of D-dopachrome.

Developing a dynamic pharmacophore model for HIV-1 integrase.

We present the first receptor-based pharmacophore model for HIV-1 integrase. The development of "dynamic" pharmacophore models is a new method that accounts for the inherent flexibility of the active site and aims to reduce the entropic penalties associated with binding a ligand. Furthermore, this new drug discovery method overcomes the limitation of an incomplete crystal structure of the target protein. A molecular dynamics (MD) simulation describes the flexibility of the uncomplexed protein. Many conformational models of the protein are saved from the MD simulations and used in a series of multi-unit search for interacting conformers (MUSIC) simulations. MUSIC is a multiple-copy minimization method, available in the BOSS program; it is used to determine binding regions for probe molecules containing functional groups that complement the active site. All protein conformations from the MD are overlaid, and conserved binding regions for the probe molecules are identified. Those conserved binding regions define the dynamic pharmacophore model. Here, the dynamic model is compared to known inhibitors of the integrase as well as a three-point, ligand-based pharmacophore model from the literature. Also, a "static" pharmacophore model was determined in the standard fashion, using a single crystal structure. Inhibitors thought to bind in the active site of HIV-1 integrase fit the dynamic model but not the static model. Finally, we have identified a set of compounds from the Available Chemicals Directory that fit the dynamic pharmacophore model, and experimental testing of the compounds has confirmed several new inhibitors.

Investigations on human immunodeficiency virus type 1 integrase/DNA binding interactions via molecular dynamics and electrostatics calculations.

The complete three-dimensional structure of the active site region of the human immunodeficiency virus type 1 (HIV-1) integrase (IN) is not unambiguously known. This region includes a flexible loop comprising residues 141-148 and the N-terminal portion of the helix alpha-4, which contains E152, the third catalytic residue, and Y143, which plays a secondary role in catalysis. Relatively high B-factors exist for most of the residues in the aforementioned region. The HIV-1 IN belongs to the polynucleotidyl transferase superfamily, whose members have been proposed to use two divalent metal ions for catalysis. Although only the position of the first metal ion has been determined crystallographically for the HIV-1 IN, we recently have proposed a binding site for the second metal ion. Based on this information, we have performed two 500-psec molecular dynamics simulations of the catalytic domain of the HIV-1 IN containing two Mg(2)+ ions. In one of the simulations, we included a dianionic phosphate group (HPO(4)(2)-) in the active site to mimic a portion of the DNA backbone of a substrate for the integration reaction. Electrostatics calculations and ionization state predictions were carried out on representative structures taken from the molecular dynamics simulations. Different conformational behaviors of the enzyme were observed, depending upon whether two Mg(2)+ ions were bound or two Mg(2)+ ions plus phosphate. The electrostatic calculations performed on the dynamical structures provide a further refinement about which regions of the catalytic domain of the HIV-1 IN may be involved in the DNA binding.

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 modeling of the intrastrand guanine-guanine DNA adducts produced by cisplatin and oxaliplatin.

Intrastrand DNA adducts formed by cisplatin and oxaliplatin were modeled with molecular mechanics minimization and restrained molecular dynamics simulations in a comparative study. A reasonable set of force field parameters for the Pt atom were refined by using the available cisplatinated DNA crystal structure as a guide. This crystal structure was also used as the starting structure for the simulations. Analysis of the resulting structures indicated that the covalent effects of oxaliplatin coordination on DNA structure were very similar to those of cisplatin. The most prominent difference between the two structures resulted from the presence of the 1,2-diaminocyclohexane ring in the oxaliplatin adduct. The modeling indicated that this ring protrudes directly outward into, and fills much of, the narrowed major groove of the bound DNA, forming a markedly altered and less polar major groove in the area of the adduct. The differences in the structure of the adducts produced by cisplatin and oxaliplatin are consistent with the observation that they are differentially recognized by the DNA mismatch repair system.

Poisson-Boltzmann model studies of molecular electrostatic properties of the cAMP-dependent protein kinase.

Protonation equilibria of residues important in the catalytic mechanism of a protein kinase were analyzed on the basis of the Poisson-Boltzmann electrostatic model along with a cluster-based treatment of the multiple titration state problem. Calculations were based upon crystallographic structures of the mammalian cAMP-dependent protein kinase, one representing the so called closed form of the enzyme and the other representing an open conformation. It was predicted that at pH 7 the preferred form of the phosphate group at the catalytically essential threonine 197 (P-Thr197) in the closed form is dianionic, whereas in the open form a monoanionic ionization state is preferred. This dianionic state of P-Thr197, in the closed form, is stabilized by interactions with ionizable residues His87, Arg165, and Lys189. Our calculations predict that the hydroxyl of the Ser residue in the peptide substrate is very difficult to ionize, both in the closed and open structures of the complex. Also, the supposed catalytic base, Asp166, does not seem to have a pK(a) appropriate to remove the hydroxyl group proton of the peptide substrate. However, when Ser of the peptide substrate is forced to remain ionized, the predicted pK(a) of Asp166 increases strongly, which suggests that the Asp residue is a likely candidate to attract the proton if the Ser residue becomes deprotonated, possibly during some structural change preceding formation of the transition state. Finally, in accord with suggestions made on the basis of the pH-dependence of kinase kinetics, our calculations predict that Glu230 and His87 are the residues responsible for the molecular pK(a) values of 6.2 and 8.5, observed in the experiment.

Docking of 4-oxalocrotonate tautomerase substrates: implications for the catalytic mechanism.

The enzyme 4-oxalocrotonate tautomerase catalyzes the ketonization of dienols, which after further processing become intermediates in the Krebs cycle. The enzyme uses a general acid-base mechanism for proton transfer: the amino-terminal proline has been shown to function as the catalytic base and Arg39 has been implicated as the catalytic acid. We report the results of molecular docking simulations of 4-oxalocrotonate tautomerase with two substrates, 2-hydroxymuconate and 5-carboxymethyl-2-hydroxymuconate. pKa calculations are also performed for the free enzyme. The predicted binding mode of 2-hydroxymuconate is in agreement with experimental data. A model for the binding mode of 5-carboxymethyl-2-hydroxymuconate is proposed which explains the lower catalytic efficiency of the enzyme toward this substrate. The pKa predictions and docking simulations support residue Arg39 as the general acid for the enzyme catalysis.

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.

Calculation of the pKa values for the ligands and side chains of Escherichia coli D-alanine:D-alanine ligase.

Poisson-Boltzmann electrostatics methods have been used to calculate the pKa shifts for the ligands and titratable side chains of D-alanine:D-alanine ligase of the ddlb gene of Escherichia coli (DdlB). The focus of this study is to determine the ionization state of the second D-alanine (D-Ala2) in the active site of DdlB. The pKa of the amine is shifted over 5 pKa units more alkaline in the protein, clearly implying that D-Ala2 is bound to DdlB in its zwitterionic state and not in the free-base form as had been previously suggested. Comparisons are made to the depsipeptide ligase from the vancomycin-resistance cascade, VanA. It is suggested that VanA has different enzymatic properties due to a change in binding specificity rather than altered catalytic behavior and that the specificity of binding D-lactate over D-Ala2 may arise from the difference in ionization characteristics of the ligands.

Rapid binding of a cationic active site inhibitor to wild type and mutant mouse acetylcholinesterase: Brownian dynamics simulation including diffusion in the active site gorge.

It is known that anionic surface residues play a role in the long-range electrostatic attraction between acetylcholinesterase and cationic ligands. In our current investigation, we show that anionic residues also play an important role in the behavior of the ligand within the active site gorge of acetylcholinesterase. Negatively charged residues near the gorge opening not only attract positively charged ligands from solution to the enzyme, but can also restrict the motion of the ligand once it is inside of the gorge. We use Brownian dynamics techniques to calculate the rate constant kon, for wild type and mutant acetylcholinesterase with a positively charged ligand. These calculations are performed by allowing the ligand to diffuse within the active site gorge. This is an extension of previously reported work in which a ligand was allowed to diffuse only to the enzyme surface. By setting the reaction criteria for the ligand closer to the active site, better agreement with experimental data is obtained. Although a number of residues influence the movement of the ligand within the gorge, Asp74 is shown to play a particularly important role in this function. Asp74 traps the ligand within the gorge, and in this way helps to ensure a reaction.

Correlation between rate of enzyme-substrate diffusional encounter and average Boltzmann factor around active site.

The utility of the average Boltzmann factor around the active site of an enzyme as the predictor of the electrostatic enhancement of the substrate binding rate is tested on a set of data on wild-type acetylcholinesterase and 18 charge mutants recently obtained by Brownian dynamics simulations. A good correlation between the average Boltzmann factors and the substrate binding rate constants is found. The effects of single charge mutations on both the Boltzmann factor and the substrate binding rate constant are modest, i.e., < 5 fold increase or decrease. This is consistent with the experimental results of Shafferman et al. but does not support their suggestion that the overall rate of the catalytic reaction is not limited by the diffusional encounter of acetylcholinesterase and its substrate.

Brownian and essential dynamics studies of the HIV-1 integrase catalytic domain.

The three-dimensional structure of the active site region of the enzyme HIV-1 integrase is not unambiguously known. This region includes a flexible peptide loop that cannot be well resolved in crystallographic determinations. Here we present two different computational approaches with different levels of resolution and on different time-scales to understand this flexibility and to analyze the dynamics of this part of the protein. We have used molecular dynamics simulations with an atomic model to simulate the region in a realistic and reliable way for 1 ns. It is found that parts of the loop wind up after 300 ps to extend an existing helix. This indicates that the helix is longer than in the earlier crystal structures that were used as basis for this study. Very recent crystal data confirms this finding, underlining the predictive value of accurate MD simulations. Essential dynamics analysis of the MD trajectory yields an anharmonic motion of this loop. We have supplemented the MD data with a much lower resolution Brownian dynamics simulation of 600 ns length. It provides ideas about the slow-motion dynamics of the loop. It is found that the loop explores a conformational space much larger than in the MD trajectory, leading to a "gating"-like motion with respect to the active site.