Search

Martini 3, the latest version of the widely used Martini force field for coarse-grained molecular dynamics simulations, is a promising tool to investigate proteins in phospholipid bilayers. However, simulating other lipid environments, such as detergent micelles, presents challenges due to the absence of validated parameters for their constituent molecules. Here, we propose parameters for the micelle-forming surfactant, dodecylphosphocholine (DPC). These result in micelle assembly with aggregation numbers in agreement with experimental values. However, we identified a lack of hydrophobic interactions between transmembrane helix protein dimers and the tails of DPC molecules, preventing insertion and stabilization of the protein in the micelles. This problem was also observed for protein insertion by self-assembling 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) or dipalmitoylphosphatidylcholine (DPPC) bilayers. We propose the reduction of the non-bonded interactions between protein and water beads by 10% as a simple and effective solution to this problem that enables protein encapsulation in phospholipid micelles and bilayers without altering protein dimerization or bilayer structure.Abstract Figure

The optimization of compounds with multiple targets is a difficult multidimensional problem in the drug discovery cycle. Here, we present a systematic, multidisciplinary approach to the development of selective anti-parasitic compounds. Computational fragment-based design of novel pteridine derivatives along with iterations of crystallographic structure determination allowed for the derivation of a structure-activity relationship for multitarget inhibition. The approach yielded compounds showing apparent picomolar inhibition of T. brucei pteridine reductase 1 (PTR1), nanomolar inhibition of L. major PTR1, and selective submicromolar inhibition of parasite dihydrofolate reductase (DHFR) versus human DHFR. Moreover, by combining design for polypharmacology with a property-based on-parasite optimization, we found three compounds that exhibited micromolar EC50 values against T. brucei brucei, whilst retaining their target inhibition. Our results provide a basis for the further development of pteridine-based compounds, and we expect our multitarget approach to be generally applicable to the design and optimization of anti-infective agents.

The COVID-19 pandemic continues to pose a substantial threat to human lives and is likely to do so for years to come. Despite the availability of vaccines, searching for efficient small-molecule drugs that are widely available, including in low- and middle-income countries, is an ongoing challenge. In this work, we report the results of a community effort, the “Billion molecules against Covid-19 challenge”, to identify small-molecule inhibitors against SARS-CoV-2 or relevant human receptors. Participating teams used a wide variety of computational methods to screen a minimum of 1 billion virtual molecules against 6 protein targets. Overall, 31 teams participated, and they suggested a total of 639,024 potentially active molecules, which were subsequently ranked to find ‘consensus compounds’. The organizing team coordinated with various contract research organizations (CROs) and collaborating institutions to synthesize and test 878 compounds for activity against proteases (Nsp5, Nsp3, TMPRSS2), nucleocapsid N, RdRP (Nsp12 domain), and (alpha) spike protein S. Overall, 27 potential inhibitors were experimentally confirmed by binding-, cleavage-, and/or viral suppression assays and are presented here. All results are freely available and can be taken further downstream without IP restrictions. Overall, we show the effectiveness of computational techniques, community efforts, and communication across research fields (i.e., protein expression and crystallography, in silico modeling, synthesis and biological assays) to accelerate the early phases of drug discovery.

Accurate protein druggability predictions are important for the selection of drug targets in the early stages of drug discovery. Due to the flexible nature of proteins, the druggability of a binding pocket may vary due to conformational changes. We have therefore developed two statistical models, a logistic regression model (TRAPP-LR) and a convolutional neural network model (TRAPP-CNN), for predicting druggability and how it varies with changes in the spatial and physicochemical properties of a binding pocket. These models are integrated into TRAPP (TRAnsient Pockets in Proteins), a tool for the analysis of binding pocket variations along a protein motion trajectory. The models, which were trained on publicly available and self-augmented data sets, show equivalent or superior performance to existing methods on test sets of protein crystal structures, and have sufficient sensitivity to identify potentially druggable protein conformations in trajectories from molecular dynamics simulations. Visualization of the evidence for the decisions of the models in TRAPP facilitates identification of the factors affecting the druggability of protein binding pockets.

The binding kinetic properties of potential drugs may significantly influence their subsequent clinical efficacy. Predictions of these properties based on computer simulations provide a useful alternative to their expensive and time-demanding experimental counterparts, even at an early drug discovery stage.Herein, we perform Scaled Molecular Dynamics (ScaledMD) simulations on a set of 27 ligands of HSP90 belonging to more than 7 chemical series in order to estimate their relative residence time. We introduce two new techniques for the analysis and the classification of the simulated unbinding trajectories. The first technique, which helps in estimating the limits of the free energy well around the bound state and the second one, based on a new contact map fingerprint, allows the description and the comparison of the paths that lead to unbinding.Using these analyses, we find that ScaledMD’s relative residence time generally enables the identification of the slowest unbinders. We propose an explanation for the underestimation of the residence times of a subset of compounds and we investigate how the biasing in ScaledMD can affect the mechanistic insights that can be gained from the simulations.

Drugs that target human thymidylate synthase (hTS) are widely used in anti-cancer therapy. However, treatment with classical substrate-site-directed TS inhibitors induces its over-expression and the development of drug resistance. We thus pursued an alternative strategy that led to the discovery of TS-dimer disrupters that bind at the monomer-monomer interface and shift the dimerization equilibrium of both the recombinant and the intracellular protein toward the inactive monomers. We performed a structural, spectroscopic and kinetic investigation of the effects of these small molecules andthe best one, E7, accelerates the proteasomal degradation of hTS in cancer cells. E7 showed a superior anticancer profile to fluorouracil in a mouse model of human pancreatic and ovarian cancer. Thus, over sixty years after the discovery of the first TS prodrug inhibitor, fluorouracil, E7 breaks the link between TS inhibition and enhanced expression in response, providing a strategy to fight drug-resistant cancers.

The virtual screening of large numbers of compounds against target protein binding sites has become an integral component of drug discovery workflows. This screening is often done by computationally docking ligands into a protein binding site of interest, but this has the drawback of a large number of poses that must be evaluated to obtain accurate estimates of protein-ligand binding affinity. We here introduce a fast pre-filtering method for ligand prioritization that is based on a set of machine learning models and uses simple pose-invariant physicochemical descriptors of the ligands and the protein binding pocket. Our method, Rapid Screening with Physicochemical Descriptors + machine learning (RASPD+), is trained on PDBbind data and achieves a regression performance that is better than that of the original RASPD method and traditional scoring functions on a range of different test sets without the need for generating ligand poses. Additionally, we use RASPD+ to identify molecular features important for binding affinity and assess the ability of RASPD+ to enrich active molecules from decoys.

Pteridine reductase-1 (PTR1) is a promising drug target for the treatment of trypanosomiasis. We investigated the potential of a previously identified class of thiadiazole inhibitors of Leishmania major PTR1 for activity against Trypanosoma brucei (Tb). We solved crystal structures of several TbPTR1-inhibitor complexes to guide the structure-based design of new thiadiazole derivatives. Subsequent synthesis and enzyme- and cell-based assays confirm new, mid-micromolar inhibitors of TbPTR1 with low toxicity. In particular, compound 4m, a biphenyl-thiadiazole-2,5-diamine with IC50 = 16 μM, was able to potentiate the antitrypanosomal activity of the dihydrofolate reductase inhibitor methotrexate (MTX) with a 4.1-fold decrease of the EC50 value. In addition, the antiparasitic activity of the combination of 4m and MTX was reversed by addition of folic acid. By adopting an efficient hit discovery platform, we demonstrate, using the 2-amino-1,3,4-thiadiazole scaffold, how a promising tool for the development of anti-T. brucei agents can be obtained.

Flavonoids represent a potential source of new antitrypanosomatidic leads. Starting from a library of natural products, we combined target-based screening on pteridine reductase 1 with phenotypic screening on Trypanosoma brucei for hit identification. Flavonols were identified as hits, and a library of 16 derivatives was synthesized. Twelve compounds showed EC50 values against T. brucei below 10 μM. Four X-ray crystal structures and docking studies explained the observed structure-activity relationships. Compound 2 (3,6-dihydroxy-2-(3-hydroxyphenyl)-4H-chromen-4-one) was selected for pharmacokinetic studies. Encapsulation of compound 2 in PLGA nanoparticles or cyclodextrins resulted in lower in vitro toxicity when compared to the free compound. Combination studies with methotrexate revealed that compound 13 (3-hydroxy-6-methoxy-2-(4-methoxyphenyl)-4H-chromen-4-one) has the highest synergistic effect at concentration of 1.3 μM, 11.7-fold dose reduction index and no toxicity toward host cells. Our results provide the basis for further chemical modifications aimed at identifying novel antitrypanosomatidic agents showing higher potency toward PTR1 and increased metabolic stability.

Sequence-specific cleavage of collagen by mammalian collagenase plays a pivotal role in cell function. Collagenases are matrix metalloproteinases that cleave the peptide bond at a specific position on fibrillar collagen. The collagenase Hemopexin-like (HPX) domain has been proposed to be responsible for substrate recognition, but the mechanism by which collagenases identify the cleavage site on fibrillar collagen is not clearly understood. In this study, Brownian dynamics simulations coupled with atomic-detail and coarse-grained molecular dynamics simulations were performed to dock matrix metalloproteinase-1 (MMP-1) on a collagen IIIα1 triple helical peptide. We find that the HPX domain recognizes the collagen triple helix at a conserved R-X11-R motif C-terminal to the cleavage site to which the HPX domain of collagen is guided electrostatically. The binding of the HPX domain between the two arginine residues is energetically stabilized by hydrophobic contacts with collagen. From the simulations and analysis of the sequences and structural flexibility of collagen and collagenase, a mechanistic scheme by which MMP-1 can recognize and bind collagen for proteolysis is proposed.

Human thymidylate synthase (hTS), a target for antiproliferative drugs, is an obligate homodimer. Single-point mutations to alanine at the monomer–monomer interface may enable the identification of specific residues that delineate sites for drugs aimed at perturbing the protein–protein interactions critical for activity. We computationally identified putative hotspot residues at the interface and designed mutants to perturb the intersubunit interaction. Dimer dissociation constants measured by a FRET-based assay range from 60 nM for wild-type hTS up to about 1 mM for single-point mutants and agree with computational predictions of the effects of these mutations. Mutations that are remote from the active site retain full or partial activity, although the substrate KM values were generally higher and the dimer was less stable. The lower dimer stability of the mutants can facilitate access to the dimer interface by small molecules and thereby aid the design of inhibitors that bind at the dimer interface.Graphical abstract

Although previous experimental studies have shown the positional preference of different amino acids (AAs) to form a stable triple helical collagen motif, the structural basis for the variations in the sequence and the positional propensity has not been systematically investigated. Thus, we have here probed the origin of the structural stability offered by the 20 naturally occurring AAs to collagen by means of classical molecular dynamics (MD) simulation. Simulations were carried out on 39 collagen-like peptides employing a host-guest approach. The results show that the propensity of the different AAs to adopt collagen-like conformations depends primarily on their ϕ and ψ angle preferences. Changes in these angles upon substitution of different AAs in the X(AA) and Y(AA) positions in the canonical ((Gly-X(AA)-Y(AA))(7))(3) motif dictate the formation of interchain hydrogen bonds, solvent interactions, and puckering of neighboring imino acids and, thus, the structural stability of the collagen. The role of solvent-mediated hydrogen bonds in the stabilization of collagen has also been elucidated from the MD simulations. In addition to the conventional hydrogen bonds known to be present in collagen, a hitherto unidentified direct interchain hydrogen bond, between the X(AA) N-H group and the Hyp O-H group of the neighboring chain, was observed during the simulations. Its occupancy was ∼36% when Leu was present at the X(AA) position.

Folate analogue inhibitors of Leishmania major pteridine reductase (PTR1) are potential antiparasitic drug candidates for combined therapy with dihydrofolate reductase (DHFR) inhibitors. To identify new molecules with specificity for PTR1, we carried out a virtual screening of the Available Chemicals Directory (ACD) database to select compounds that could interact with L. major PTR1 but not with human DHFR. Through two rounds of drug discovery, we successfully identified eighteen drug-like molecules with low micromolar affinities and high in vitro specificity profiles. Their efficacy against Leishmania species was studied in cultured cells of the promastigote stage, using the compounds both alone and in combination with 1 (pyrimethamine; 5-(4-chlorophenyl)-6-ethylpyrimidine-2,4-diamine). Six compounds showed efficacy only in combination. In toxicity tests against human fibroblasts, several compounds showed low toxicity. One compound, 5c (riluzole; 6-(trifluoromethoxy)-1,3-benzothiazol-2-ylamine), a known drug approved for CNS pathologies, was active in combination and is suitable for early preclinical evaluation of its potential for label extension as a PTR1 inhibitor and antiparasitic drug candidate.

In order to study protein-inorganic surface association processes, we have developed a physics-based energy model, the ProMetCS model, which describes protein-surface interactions at the atomistic level while treating the solvent as a continuum. Here, we present an approach to modeling the interaction of a protein with an atomically flat Au(111) surface in an aqueous solvent. Protein-gold interactions are modeled as the sum of van der Waals, weak chemisorption, and electrostatic interactions, as well as the change in free energy due to partial desolvation of the protein and the metal surface upon association. This desolvation energy includes the effects of water-protein, water-surface, and water-water interactions and has been parametrized using molecular dynamics (MD) simulations of water molecules and a test atom at a gold-water interface. The proposed procedure for computing the energy terms is mostly grid-based and is therefore efficient for application to long-time simulations of protein binding processes. The approach was tested for capped amino acid residues whose potentials of mean force for binding to a gold surface were computed and compared with those obtained previously in MD simulations with water treated explicitly. Calculations show good quantitative agreement with the results from MD simulations for all but one amino acid (Trp), as well as correspondence with available experimental data on the adhesion properties of amino acids.

The factors that determine the extent to which diffusion and thermal activation processes govern electron transfer (ET) between proteins are debated. The process of ET between plastocyanin (PC) and cytochrome f (CytF) from the cyanobacterium Phormidium laminosum was initially thought to be diffusion-controlled but later was found to be under activation control (Schlarb-Ridley, B. G.; et al. Biochemistry 2005, 44, 6232). Here we describe Brownian dynamics simulations of the diffusional association of PC and CytF, from which ET rates were computed using a detailed model of ET events that was applied to all of the generated protein configurations. The proteins were modeled as rigid bodies represented in atomic detail. In addition to electrostatic forces, which were modeled as in our previous simulations of protein-protein association, the proteins interacted by a nonpolar desolvation (hydrophobic) force whose derivation is described here. The simulations yielded close to realistic residence times of transient protein-protein encounter complexes of up to tens of microseconds. The activation barrier for individual ET events derived from the simulations was positive. Whereas the electrostatic interactions between P. laminosum PC and CytF are weak, simulations for a second cyanobacterial PC-CytF pair, that from Nostoc sp. PCC 7119, revealed ET rates influenced by stronger electrostatic interactions. In both cases, the simulations imply significant contributions to ET from both diffusion and thermal activation processes.

In protein unfolding simulations, elevated temperature, significantly exceeding the melting temperature Tm, provides an important means to accelerate unfolding to a computationally accessible time range. This procedure is based on the assumption that protein thermal unfolding has Arrhenius behavior and therefore that increasing temperature does not alter the protein unfolding pathways. However, in nature, proteins can show non-Arrhenius behavior and, in practice, overly fast unfolding in high-temperature simulations can result in difficulties in identifying unfolding intermediates and distinguishing their relative stabilities. In this paper, we describe simulations of two WW domains, small protein domains that have a three-stranded β-sheet structure. Simulations were carried out at several temperatures ranging from 300 K to 500 K, starting from folded structures. The results demonstrate the temperature dependence of the unfolding pathways, showing that to obtain unfolding pathways corresponding to those observed in experiments, the elevation of the simulation temperature has to be controlled. Based on trajectory analysis, we proposed a qualitative criterion for judging when an elevated temperature is acceptable or not, namely, that the temperature must be such that the native folded state is sampled substantially before protein unfolding begins. While depending on force field parameters and protein fold complexity, this criterion can be quantified to obtain the upper bound of an "acceptable elevated temperature", which was observed to be dependent on the thermostabilities of the two WW domain proteins.

The oxidation of melatonin by the mammalian myeloperoxidase (MPO) provides protection against the damaging effects of reactive oxygen species. Indole derivatives, such as melatonin and serotonin, are also substrates of the plant horseradish peroxidase (HRP), but this enzyme exhibits remarkable differences from MPO in the specificity and reaction rates for these compounds. A structural understanding of the determinants of the reactivity of these enzymes to indole derivatives would greatly aid their exploitation for biosynthetic and drug design applications. Consequently, after validation of the docking procedure, we performed computational docking of melatonin and serotonin to structural models of the ferric and compound I and II (co I and co II, respectively) states of HRP and MPO. The substrates dock at the heme edge on the distal side, but with different orientations in the two proteins. The distal cavity is larger in MPO than in HRP; however, in MPO, the substrates make closer contacts with the heme involving ring stacking, whereas in HRP, no ring stacking is observed. The observed differences in substrate binding may contribute to the higher reaction rates and lower substrate specificity of MPO relative to those of HRP. The docking results, along with the previously measured heme-protein reduction potentials, suggest that the differentially lowered reaction rates of co II of HRP and MPO with respect to those of co I could stem from as yet undetermined conformational or electrostatic differences between the co I and co II states of MPO, which are absent in HRP.

The activity of a biological compound is dependent both on specific binding to a target receptor and its ADME (Absorption, Distribution, Metabolism, Excretion) properties. A challenge to predict biological activity is to consider both contributions simultaneously in deriving quantitative models. We present a novel approach to derive QSAR models combining similarity analysis of molecular interaction fields (MIFs) with prediction of logP and/or logD. This new classification method is applied to a set of about 100 compounds related to the auxin plant hormone. The classification based on similarity of their interaction fields is more successful for the indole than the phenoxy compounds. The classification of the phenoxy compounds is however improved by taking into account the influence of the logP and/or the logD values on biological activity. With the new combined method, the majority (8 out of 10) of the previously misclassified derivatives of phenoxy acetic acid are classified in accord with their bioassays. The recently determined crystal structure of the auxin-binding protein 1 (ABP1) enabled validation of our approach. The results of docking a few auxin related compounds with different biological activity to ABP1 correlate well with the classification based on similarity of MIFs only. Biological activity is, however, better predicted by a combined similarity of MIFs + logP/logD approach.

Specific binding of transcription factors to DNA is crucial for gene regulation. We derived models for the binding specificity of transcription factors of the nuclear receptor family to DNA using two QSAR methods: a Free-Wilson-like method and COMparative BINding Energy (COMBINE) analysis. The analysis is based on experimental data for the interaction of 20 mutant glucocorticoid receptor DNA-binding domains with 16 different response elements in a total of 320 complexes (Zilliacus, J.; Wright, A. P.; Carlstedt-Duke, J.; Nilsson, L.; Gustafsson, J. A. Proteins 1995, 21, 57-67). The predictive abilities of the models obtained by the two methods are similar. The COMBINE analysis indicates that the most important properties for determining binding specificity for this dataset are the changes upon binding of the solvation free energies of the bases that are mutated in the dataset and the electrostatic interactions of the mutated nucleotides with certain charged amino acids. Further important descriptors are the changes of solvation free energy and surface area of the side chain of the mutated residue. It is clear, however, that there are additional features important for the specificity of binding that are not included in the models, such as differences in interfacial hydration of the complexes.

Heparin−protein interactions play an important role in many steps of the immune system. Here, we evaluated the search capabilities of three widely used programsGRID, DOCK, and AutoDockfor heparin binding sites. Because of the weak surface complementarity and the high charge density of the sulfated sugar chain, the docking of heparin to its protein partners presents a challenging task for computational docking. Our protocols were tested on antithrombin and acidic and basic fibroblast growth factor, the only three proteins for which structures of their complexes with heparin are available. With all three programs, the heparin binding site for these test cases was determined correctly. We then used these protocols to predict the heparin binding site on Interleukin-8, a chemokine with a central role in the human immune response. The results indicate that His18, Lys20, Arg60, Lys64, and Arg68 in interleukin-8 bind to heparin.

NMR studies of the YTGP tetrapeptide in aqueous solution show a large structural shift on the glycine peptide proton, indicating that the peptide adopts structures in which the aromatic ring lies close to this proton in an aromatic−(i+2) amine interaction (Kemmink and Creighton, J. Mol. Biol. 1993, 234, 861). Here, molecular dynamics simulations are used to obtain details at the atomic level of the structural properties of this system. To maximize the volume of conformational space searched, 13 separate trajectories of 300 ps were calculated at 300 K, with starting structures chosen from a conformational search, protein crystal structures, and modeling. Analysis shows that together the trajectories sample all the important regions of conformational space. Structures from the trajectories were clustered into conformational states, defined as regions in coordinate space around maxima in the probability distribution. Nine conformational states with significant populations were identified in which the glycine...

The ability to reliably compute accurate protein−ligand binding affinities is crucial to understanding protein−ligand recognition and to structure-based drug design. A ligand's binding affinity is specified by its absolute binding free energy, ΔGbind, the free energy difference between the bound and unbound states. To compute accurate free energy differences by free energy perturbation (FEP), “alchemical” rather than physical processes are usually simulated by molecular dynamics simulations so as to minimize the perturbation to the system. Here, we report a novel “alchemistic” application of the FEP methodology involving a large perturbation. By mutating a ligand with 11 non-hydrogen atoms into six water molecules in the binding site of a protein, we computed a ΔGbind within 3 kJ/mol of the experimental value. This is the first successful example of the computation of ΔGbind for a protein:ligand pair with full treatment of the solvent degrees of freedom.

In order to study protein−protein and protein−DNA association, the electrostatic forces and interaction free energies for two macromolecules at different mutual orientations and separations need to be evaluated. Realistic systems typically consist of thousands of atomic charges in an environment with a nonuniform dielectric permittivity and a solvent of nonzero ionic strength. Consequently, accurate evaluations of electrostatic forces and energies for such systems are only computationally feasible for a limited number of macromolecule positions. Here, we show that, by representing each molecule by a small number of effective charges in a uniform dielectric, intermolecular electrostatic interactions can be calculated simply and with high accuracy. The effective charges are derived by fitting them to reproduce the molecular electrostatic potential calculated by numerical solution of the finite-difference linearized Poisson−Boltzmann equation. The derived charges are expected to be useful in applications suc...

The pH-dependent characteristics of a protein are determined by the pKas of its titratable residues. Continuum electrostatic approaches provide a fast means of calculating the pKas of the titratable residues in proteins immersed in aqueous ionic solution. Here, we show how optimization of the parameters used in continuum electrostatic calculations can lead to significant improvements in the accuracy of pKa predictions. Dependence on the protein dielectric constant is studied, and two classes of ionizable sites are identified: (1) a large number of mostly solvent exposed residues for which the best pKa is calculated if the protein dielectric constant is set close to that of the aqueous solvent; and (2) a small number of mostly buried residues for which the best pKa is calculated with a lower site-specific protein dielectric constant. These two classes of ionizable sites can be distinguished using a criterion based on desolvation energy. A priori determination of the optimum protein dielectric constant for...

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.

A new computational method for deducing quantitative structure-activity relationships (QSARs) using structural data from ligand-macromolecule complexes is presented. First, the ligand-macromolecule interaction energy is computed for a set of ligands using molecular mechanics calculations. Then, by selecting and scaling components of the ligand-macromolecule interaction energy that show good predictive ability, a regression equation is obtained in which activity is correlated with the interaction energies of parts of the ligands and key regions of the macromolecule. Application to the interaction of the human synovial fluid phospholipase A2 with 26 inhibitors indicates that the derived QSAR has good predictive ability and provides insight into the mechanism of enzyme inhibition. The method, which we term comparative binding energy (COMBINE) analysis, is expected to be applicable to ligand-receptor interactions in a range of contexts including rational drug design, host-guest systems, and protein engineering.

The specificity of interactions between biological macromolecules and their ligands may be partially attributed to the directional properties of hydrogen bonds. We have now extended the GRID method (Goodford, P. J. J. Med. Chem. 1985, 28, 849. Boobbyer, D. N. A.; Goodford, P. J.; McWhinnie, P. M.; Wade, R. C. J. Med. Chem. 1989, 32, 1083), of determining energetically favorable ligand binding sites on molecules of known structure, in order to improve the treatment of groups which can make multiple hydrogen bonds. In this method, the interaction energy between a probe (a small chemical group that may be part of a larger ligand) and a target molecule is calculated using an energy function which includes a hydrogen bond term which is dependent on the length of the hydrogen bond, its orientation at the hydrogen-bonding atoms, and their chemical character. The methods described in the preceding paper (Wade, R. C.; Clark, K. J.; Goodford, P. J. J. Med. Chem., preceding paper in this issue) for probes capable of making two hydrogen bonds are here extended to the following probes which have the ability to make more than two hydrogen bonds: ammonium-NH3+, amine-NH2, sp3-hybridized hydroxyl, and water. Use of the improved GRID procedure is demonstrated by the determination of the conformation of an amino acid side chain at the subunit interface in hemoglobin and of the location of water binding sites in human lysozyme.

Brownian dynamics simulations of the diffusional encounter between the glycolytic enzyme triose phosphate isomerase (TIM) and its substrate, D-glyceraldehyde phosphate (GAP), were performed. GAP was modeled hydrodynamically as two touching spheres (a “dumbbell”) using charges which reproduced the molecular dipole moment of the glyceraldehyde phosphate molecule as estimated by an ab initio molecular orbital calculation. The crystal structure of TIM was used to construct a detailed topographical and electrostatic grid on which the diffusion of the dumbbell was numerically simulated. By determining the number of diffusional encounters which resulted in GAP descending into the active sites of TIM with the appropriate orientation, the diffusioncontrolled rate constant for the reaction was estimated to be 1.7 X lo8 M-* s-l. This is in reasonable agreement with the experimentally determined diffusion-controlled rate constant of 4.8 X lo8 M-l si. By reversing the direction of the dipole moment on the GAP model, it was shown that the orientational steering of the substrate by electrostatic torques can significantly increase the reaction rate constant. This effect is in addition to the previously established translational steering of the charged substrate by electrostatic forces.

Comparative binding energy (COMBINE) analysis was conducted for 18 substrates of the haloalkane dehalogenase from Xanthobacter autotrophicusGJ10 (DhlA): 1-chlorobutane, 1-chlorohexane, dichloromethane, 1,2-dichloroethane, 1,2-dichloropropane, 2-chloroethanol, epichlorohydrine, 2-chloro- acetonitrile, 2-chloroacetamide, and their brominated analogues. The purpose of the COMBINE analysis was to identify the amino acid residues determining the substrate specificity of the haloalkane dehalogenase. This knowledge is essential for the tailoring of this enzyme for biotechnological applications. Complexes of the enzyme with these substrates were modeled and then refined by molecular mechanics energy minimization. The intermolecular enzyme-substrate energy was decomposed into residue-wise van der Waals and electrostatic contributions and complemented by surface area dependent and electrostatic desolvation terms. Partial least-squares projection to latent structures analysis was then used to establish relationships between the energy contributions and the experimental apparent dissociation constants. A model containing van der Waals and electrostatic intermolecular interaction energy contributions calculated using the AMBER force field explained 91% (73% cross-validated) of the quantitative variance in the apparent dissociation constants. A model based on van der Waals intermolecular contributions from AMBER and electrostatic interactions derived from the Poisson -Boltzmann equation explained 93% (74% cross- validated) of the quantitative variance. COMBINE models predicted correctly the change in apparent dissociation constants upon single-point mutation of DhlA for six enzyme-substrate complexes. The amino acid residues contributing most significantly to the substrate specificity of DhlA were identified; they include Asp124, Trp125, Phe164, Phe172, Trp175, Phe222, Pro223, and Leu263. These residues are suitable targets for modification by site-directed mutagenesis.

An empirical energy function designed to calculate the interaction energy of a chemical probe group, such as a carbonyl oxygen or an amine nitrogen atom, with a target molecule has been developed. This function is used to determine the sites where ligands, such as drugs, may bind to a chosen target molecule which may be a protein, a nucleic acid, a polysaccharide, or a small organic molecule. The energy function is composed of a Lennard-Jones, an electrostatic and a hydrogen-bonding term. The latter is dependent on the length and orientation of the hydrogen bond and also on the chemical nature of the hydrogen-bonding atoms. These terms have been formulated by fitting to experimental observations of hydrogen bonds in crystal structures. In the calculations, thermal motion of the hydrogen-bonding hydrogen atoms and lone-pair electrons may be taken into account. For example, in a alcoholic hydroxyl group, the hydrogen may rotate around the C-O bond at the observed tetrahedral angle. In a histidine residue, a hydrogen atom may be bonded to either of the two imidazole nitrogens and movement of this hydrogen will cause a redistribution of charge which is dependent on the nature of the probe group and the surrounding environment. The shape of some of the energy functions is demonstrated on molecules of pharmacological interest.