 |
|
 |
 | We have recently changed our website. Looking for the previous site? |
Abstracts of Articles in 2003
To request any journal reprints, please contact us.- The Relaxed Complex Method: Accommodating Receptor Flexibility for Drug Design with an Improved Scoring Scheme
- Modeling HIV-1 Integrase Complexes Based on their Hydrodynamic Properties
- Electrostatic Interactions
- Protein Flexibility and Computer-aided Drug Design
- On the Evaluation and Optimisation of Protein X-ray Structures for pKa Calculations
- Studying the Affinity and Kinetics of Molecular Association with Molecular Dynamics Simulation
- Protein Simulation and Drug Design
- The pH Dependence of Stability of the Activation Helix and the Catalytic site of GART
- Brownian Dynamics Simulation of Helix-Capping Motifs
- Computational Analysis of the Interactions between the Angiogenesis Inhibitor PD173074 and Fibroblast Growth Factor Receptor 1
- Finite Element Simulations of Acetylcholine Diffusion in Neuromuscular Junctions
- Influence of Structural Fluctuation on Enzyme Reaction Energy Barriers in Combined Quantum Mechanical/Molecular Mechanical Studies
- A Computational Model of Binding Thermodynamics: The Design of Cyclin-dependent Kinase 2 Inhibitors
- Nanosecond Dynamics of the Mouse Acetylcholinesterase Cys69-Cys96 Omega Loop
- From Model Complexes to Metalloprotein Inhibition: A Synergistic Approach to Structure-based Drug Discovery
- Calculating pKa Values in Enzyme Active Sites
- Brownian Dynamics Simulations of Ion Atmospheres around Polyalanine and B-DNA: Effects of Biomolecular Dielectric
- The Physical Basis of Microtubule Structure and Stability
- The Dynamics of Ligand Barrier Crossing inside the Acetylcholinesterase Gorge
- Asymmetric Structural Motions of the Homomeric α7 Nicotinic Receptor Ligand Binding Domain Revealed by Molecular Dynamics Simulation
- Studying the roles of W86, E202, and Y337 in binding of acetylcholine to acetylcholinesterase using a combined molecular dynamics and multiple docking approach
The Relaxed Complex Method: Accommodating Receptor Flexibility for Drug Design with an Improved Scoring SchemeBiopolymers, Vol. 68, Issue 1, pp. 47-62 (2003, Kollman memorial issue) [PubMed 12579579]
An extension of the new computational methodology for drug design, the "relaxed complex" method (http://dx.doi.org/10.1021/ja0260162" target="_blank" class="ref">J.-H. Lin, A.L. Perryman, J. Schames, J.A. McCammon, J. Amer. Chem. Soc. 124, 2002, 5632-5633), that accommodates receptor flexibility is described. This "relaxed complex" method recognizes that ligand may bind to conformations that occur only rarely in the dynamics of the receptor. We have shown that the ligand-enzyme binding modes are very sensitive to the enzyme conformations, and our approach is capable of finding the best ligand enzyme complexes. Rapid docking serves as an efficient initial filtering method to screen a myraid of docking modes to a limited set, and it is then followed by more accurate scoring with the MM/PBSA approach to find the best ligand-receptor complexes. The MM/PBSA scorings consistently indicate that the calculated binding modes which are most similar to those observed in the X-ray crystallographic complexes are the ones with the lowest free energies.
An extension of the new computational methodology for drug design, the "relaxed complex" method (http://dx.doi.org/10.1021/ja0260162" target="_blank" class="ref">J.-H. Lin, A.L. Perryman, J. Schames, J.A. McCammon, J. Amer. Chem. Soc. 124, 2002, 5632-5633), that accommodates receptor flexibility is described. This "relaxed complex" method recognizes that ligand may bind to conformations that occur only rarely in the dynamics of the receptor. We have shown that the ligand-enzyme binding modes are very sensitive to the enzyme conformations, and our approach is capable of finding the best ligand enzyme complexes. Rapid docking serves as an efficient initial filtering method to screen a myraid of docking modes to a limited set, and it is then followed by more accurate scoring with the MM/PBSA approach to find the best ligand-receptor complexes. The MM/PBSA scorings consistently indicate that the calculated binding modes which are most similar to those observed in the X-ray crystallographic complexes are the ones with the lowest free energies.
Modeling HIV-1 Integrase Complexes Based on their Hydrodynamic PropertiesBiopolymers, Vol. 68, Issue 1, pp. 110-120 (2003, Kollman memorial issue) [PubMed 12579583]
We present a model structure of a candidate tetramer for HIV-1 integrase. The model was built in 3 steps using data from fluorescence anisotropy, structures of the individual integrase domains, cross-linking data, and other biochemical data. First, the structure of the full-length integrase monomer was modeled using the individual domain structures and the hydrodynamic properties of the full-length protein that was recently measured by fluorescence depolarization. We calculated the rotational correlation times for different arrangements of three integrase domains, revealing that only structures with close proximity among the domains satisfied the experimental data. The orientations of the domains were constrained by iterative tests against the data on cross-linking and fooprinting in integrase-DNA complexes. Second, the structure of an integrase dimer was obtained by joining the model monomers in accordance with the available dimeric crystal structures of the catalytic core. The hydrodynamic properties of the dimer were in agreement with the experimental values. Third, the active sites of the two model dimers were placed in agreement with the spacing between the sites of integration on target DNA as well as the integrase-DNA cross-linking data, resulting in 2-fold symmetry of a tetrameric complex. The model is consistent with the experimental data indicating that the F185K substitution, which is found in the model at tetramerization interface, selectively disrupts correct complex formation in vitro and HIV replication in vivo. Our model of the integrase tetramer bound to DNA may help to design anti-integrase inhibitors.
We present a model structure of a candidate tetramer for HIV-1 integrase. The model was built in 3 steps using data from fluorescence anisotropy, structures of the individual integrase domains, cross-linking data, and other biochemical data. First, the structure of the full-length integrase monomer was modeled using the individual domain structures and the hydrodynamic properties of the full-length protein that was recently measured by fluorescence depolarization. We calculated the rotational correlation times for different arrangements of three integrase domains, revealing that only structures with close proximity among the domains satisfied the experimental data. The orientations of the domains were constrained by iterative tests against the data on cross-linking and fooprinting in integrase-DNA complexes. Second, the structure of an integrase dimer was obtained by joining the model monomers in accordance with the available dimeric crystal structures of the catalytic core. The hydrodynamic properties of the dimer were in agreement with the experimental values. Third, the active sites of the two model dimers were placed in agreement with the spacing between the sites of integration on target DNA as well as the integrase-DNA cross-linking data, resulting in 2-fold symmetry of a tetrameric complex. The model is consistent with the experimental data indicating that the F185K substitution, which is found in the model at tetramerization interface, selectively disrupts correct complex formation in vitro and HIV replication in vivo. Our model of the integrase tetramer bound to DNA may help to design anti-integrase inhibitors.
Electrostatic InteractionsIn "Structural Bioinformatics," H. Weissig and P.E. Bourne, Eds., John Wiley & Sons, New York, pp. 427-440 (2003) [PubMed 12647398]
An understanding of electrostatic interactions is essential for the full development of structural bioinformatics. The structures of proteins and other biopolymers are being determined at an increasing rate through structural genomics and other efforts. Specific linkages of these biopolymers in cellular pathways or supramolecular assemblages are being detected by genetic and other experimental efforts. To integrate this information in physical models for drug discovery and other applications requires the ability to evaluate the energetic interactions within and among biopolymers. Among the various components of molecular energetics, the electrostatic interactions are of special importance due to the long range of these interactions and the substantial charges of typical components of biopolymers. Indeed, electrostatics can be used to help assign biopolymers such as proteins to functional families, since particular kinds of ligand binding sites may be indicated by the spatial distribution of the charges in the proteins. In what follows, we provide a brief overview of the roles of electrostatics in biopolymers and supramolecular assemblages, and then outline some of the methods that have been developed for analyzing electrostatic interactions.
An understanding of electrostatic interactions is essential for the full development of structural bioinformatics. The structures of proteins and other biopolymers are being determined at an increasing rate through structural genomics and other efforts. Specific linkages of these biopolymers in cellular pathways or supramolecular assemblages are being detected by genetic and other experimental efforts. To integrate this information in physical models for drug discovery and other applications requires the ability to evaluate the energetic interactions within and among biopolymers. Among the various components of molecular energetics, the electrostatic interactions are of special importance due to the long range of these interactions and the substantial charges of typical components of biopolymers. Indeed, electrostatics can be used to help assign biopolymers such as proteins to functional families, since particular kinds of ligand binding sites may be indicated by the spatial distribution of the charges in the proteins. In what follows, we provide a brief overview of the roles of electrostatics in biopolymers and supramolecular assemblages, and then outline some of the methods that have been developed for analyzing electrostatic interactions.
Protein Flexibility and Computer-aided Drug DesignAnnual Review of Pharmacology and Toxicology, Vol. 43, No. 1, pp. 31-45 (2003) [PubMed 12142469]
Although computational techniques are increasingly being used in computer-aided drug design, the effects due to protein flexibility are still ignored in many applications. This review revisits rigorous statistical mechanical methods for predicting binding affinity, discusses some recent developments for improving their speed and reliability, and examines faster approximate models for facilitating virtual screening and lead optimization.
Although computational techniques are increasingly being used in computer-aided drug design, the effects due to protein flexibility are still ignored in many applications. This review revisits rigorous statistical mechanical methods for predicting binding affinity, discusses some recent developments for improving their speed and reliability, and examines faster approximate models for facilitating virtual screening and lead optimization.
On the Evaluation and Optimisation of Protein X-ray Structures for pKa CalculationsProtein Science, Vol. 12, Issue 2, pp. 313-326 (2003) [PubMed 12538895]
The calculation of the physical properties of a protein from its X-ray structure is of importance in virtually every aspect of modern biology. Although computational algorithms have been developed for calculating everything from the dynamics of a protein to its binding specificity, only limited information is available on the ability of these methods to give accurate results when used with a particular X-ray structure. We examine the ability of a pKa calculation algorithm to predict the proton-donating residue in the catalytic mechanism of Hen Egg White Lysozyme. We examine the correlation between the ability of the pKa calculation method to get the correct result and the overall characteristics of 41 X-ray structures such as crystallisation conditions, resolution and the output of structure validation software. We furthermore examine the ability of Energy Minimisations, Molecular Dynamics simulations and structure-perturbation methods to optimise the X-ray structures such that these give correct results with the pKa calculation algorithm. We propose a set of criteria for identifying the proton donor in a catalytic mechanism, and demonstrate that the application of these criteria give highly accurate prediction results when using unmodified X-ray structures. More specifically we are able to successfully identify the proton donor in 85% of the X-ray structures when excluding structures with crystal contacts near the active site. Neither the use of the overall chacteristics of the X-ray structures nor the optimisation of the structure by EM, MDL or other methods improves the results of the pKa calculation algorithm. We discuss these results and their implications for the design of structure-based energy calculation algorithms in general.
The calculation of the physical properties of a protein from its X-ray structure is of importance in virtually every aspect of modern biology. Although computational algorithms have been developed for calculating everything from the dynamics of a protein to its binding specificity, only limited information is available on the ability of these methods to give accurate results when used with a particular X-ray structure. We examine the ability of a pKa calculation algorithm to predict the proton-donating residue in the catalytic mechanism of Hen Egg White Lysozyme. We examine the correlation between the ability of the pKa calculation method to get the correct result and the overall characteristics of 41 X-ray structures such as crystallisation conditions, resolution and the output of structure validation software. We furthermore examine the ability of Energy Minimisations, Molecular Dynamics simulations and structure-perturbation methods to optimise the X-ray structures such that these give correct results with the pKa calculation algorithm. We propose a set of criteria for identifying the proton donor in a catalytic mechanism, and demonstrate that the application of these criteria give highly accurate prediction results when using unmodified X-ray structures. More specifically we are able to successfully identify the proton donor in 85% of the X-ray structures when excluding structures with crystal contacts near the active site. Neither the use of the overall chacteristics of the X-ray structures nor the optimisation of the structure by EM, MDL or other methods improves the results of the pKa calculation algorithm. We discuss these results and their implications for the design of structure-based energy calculation algorithms in general.
Studying the Affinity and Kinetics of Molecular Association with Molecular Dynamics SimulationJournal of Chemical Physics, Vol. 118, pp. 1821-1827 (2003)
Given a long molecular dynamics trajectory which consists of hundreds of association and dissociation events, the theoretical formulas to calculate the affinity and dissociation rate constant are presented. The derivation is based on the survival function of the associated complex, and it emphasizes the nearest neighbor distribution function. The applicability of this brute-force approach is demonstrated by the simulation of methane association in water.
Given a long molecular dynamics trajectory which consists of hundreds of association and dissociation events, the theoretical formulas to calculate the affinity and dissociation rate constant are presented. The derivation is based on the survival function of the associated complex, and it emphasizes the nearest neighbor distribution function. The applicability of this brute-force approach is demonstrated by the simulation of methane association in water.
Protein Simulation and Drug DesignAdvances in Protein Chemistry, Vol. 66, pp. 87-121 (2003)
The rapid advance of computer technology, force-field developments, and simulation methods and algorithms has gradually made force-field-based methods more easily and reliably useful in computer-aided drug design. In some applications, simple molecular mechanics geometry optimization and energy calculations can already suggest compounds that are worthwhile to make and/or screen. Such calculations can be improved by using modern continuum solvent models in improving the description of solvation effects. Proper account of protein flexibility, via molecular dynamics simulations, for example, can further enhance the predictive reliability of simulation models. Combined with rigorous statistical mechanical theories, these simulations can help predict binding free energy (Tembe and McCammon, 1984; Wong and McCammon, 1986b; Wong and McCammon, 1986a; Mezei and Beveridge, 1986; Bash et al., 1987; Jorgenson, 1989; Beveridge and DiCapua, 1989). While these simulations were limited to several tens of picoseconds in the early stages, they had moved a step forward beyond the fixed-conformation picture that dominated many scientists' thinking at the time. Recent developments in simulation methods and computer hardware have made it possible to extend the time scale of these simulations so that more reliable predictions can be made, and that results can be obtained more quickly. Methods that improve conformational sampling should further faciliate these simulations. Various approximate methods with different degrees of sophistication have also been introduced to faciliate the modeling of molecular recognition. In practice drug design, it should be beneficial to use a hierarchical approach in which simple but fast methods are first used for preliminary screening of real and virtual compound libraries and for suggesting potentially useful drug candidates to make/screen, followed by further evaluations by models with increasing levels of sophistication. In this article, we will review various approximate and rigorous force-field-based methods that can fit into such a hierarchical approach.
The rapid advance of computer technology, force-field developments, and simulation methods and algorithms has gradually made force-field-based methods more easily and reliably useful in computer-aided drug design. In some applications, simple molecular mechanics geometry optimization and energy calculations can already suggest compounds that are worthwhile to make and/or screen. Such calculations can be improved by using modern continuum solvent models in improving the description of solvation effects. Proper account of protein flexibility, via molecular dynamics simulations, for example, can further enhance the predictive reliability of simulation models. Combined with rigorous statistical mechanical theories, these simulations can help predict binding free energy (Tembe and McCammon, 1984; Wong and McCammon, 1986b; Wong and McCammon, 1986a; Mezei and Beveridge, 1986; Bash et al., 1987; Jorgenson, 1989; Beveridge and DiCapua, 1989). While these simulations were limited to several tens of picoseconds in the early stages, they had moved a step forward beyond the fixed-conformation picture that dominated many scientists' thinking at the time. Recent developments in simulation methods and computer hardware have made it possible to extend the time scale of these simulations so that more reliable predictions can be made, and that results can be obtained more quickly. Methods that improve conformational sampling should further faciliate these simulations. Various approximate methods with different degrees of sophistication have also been introduced to faciliate the modeling of molecular recognition. In practice drug design, it should be beneficial to use a hierarchical approach in which simple but fast methods are first used for preliminary screening of real and virtual compound libraries and for suggesting potentially useful drug candidates to make/screen, followed by further evaluations by models with increasing levels of sophistication. In this article, we will review various approximate and rigorous force-field-based methods that can fit into such a hierarchical approach.
The pH Dependence of Stability of the Activation Helix and the Catalytic site of GARTBiophysical Chemistry, Vol. 105, Issues 2-3, pp. 279-291 (2003) [PubMed 14499900]
We have predicted the free energy of unfolding for the pH-dependent helix-coil transition of the activtion helix of GART using continuum electrostatic calculations and structural modeling. We have assigned the contributions of each element of secondary structure and of each ionizable residue, within and in the vicinity of the activation helix, to the stability of several fragments of GART that participate in the formation of the catalytic site. We demonstrate that the interaction of His121-His132 contributes 2.2 kcal/mole to the ionization free energy between pH 0 and ~6. We also show that the ionization state of a network of five histidines, His108, His119, His121, His132, and His137, and two aspartic acids Asp141 and Asp144, contributes ~12 kcal/mol to the stability of the catalytic site of GART, out of a total stability of 16 kcal/mole of the whole enzyme. These interactions are important for the formation of the catalytic site of GART.
We have predicted the free energy of unfolding for the pH-dependent helix-coil transition of the activtion helix of GART using continuum electrostatic calculations and structural modeling. We have assigned the contributions of each element of secondary structure and of each ionizable residue, within and in the vicinity of the activation helix, to the stability of several fragments of GART that participate in the formation of the catalytic site. We demonstrate that the interaction of His121-His132 contributes 2.2 kcal/mole to the ionization free energy between pH 0 and ~6. We also show that the ionization state of a network of five histidines, His108, His119, His121, His132, and His137, and two aspartic acids Asp141 and Asp144, contributes ~12 kcal/mol to the stability of the catalytic site of GART, out of a total stability of 16 kcal/mole of the whole enzyme. These interactions are important for the formation of the catalytic site of GART.
Brownian Dynamics Simulation of Helix-Capping MotifsBiopolymers, Vol. 70, Issue 2, pp. 252-259 (2003) [PubMed 14517913]
Helix-capping motifs are believed to play an important role in stabalizing α-helices and defining helix start and stop signals. We performed microsecond scale Brownian dynamics simulations to study ten XAAD sequences with x=(A,E,I,L,N,Q,S,T,V,Y), to examine their propensity to form helix capping motifs and correlate these reulsts with those obtained from analyzing a structural database of proteins. For the widely studied capping box motif S**D, where * can be any amino acid residue, the simulations suggested that one of the two hydrogen bonds proposed earlier as a stabilizing factor might not be as important. On the other hand, side chain interactions between the capping residue and the third residue downstream on the polypeptide chain might also play a role in stabilizing this motif. These results are consistent with the explicit-solvent molecular dynamics simulations of two capping box motifs found in the proteins BPTI and α-dendrotoxin. Principal component analysis of the SAAD trajectory showed that the first three principal components, after those corresponding to translational-rotational motion were removed, accounted for more than half of the conformational fluctuations. The first component separated the conformational space into two parts with the all-helical conformation and the capping box motif lying largely in one part. the second component, on the other hand, could be used to describe conformational transitions between the all-helical form and the capping box motif.
Helix-capping motifs are believed to play an important role in stabalizing α-helices and defining helix start and stop signals. We performed microsecond scale Brownian dynamics simulations to study ten XAAD sequences with x=(A,E,I,L,N,Q,S,T,V,Y), to examine their propensity to form helix capping motifs and correlate these reulsts with those obtained from analyzing a structural database of proteins. For the widely studied capping box motif S**D, where * can be any amino acid residue, the simulations suggested that one of the two hydrogen bonds proposed earlier as a stabilizing factor might not be as important. On the other hand, side chain interactions between the capping residue and the third residue downstream on the polypeptide chain might also play a role in stabilizing this motif. These results are consistent with the explicit-solvent molecular dynamics simulations of two capping box motifs found in the proteins BPTI and α-dendrotoxin. Principal component analysis of the SAAD trajectory showed that the first three principal components, after those corresponding to translational-rotational motion were removed, accounted for more than half of the conformational fluctuations. The first component separated the conformational space into two parts with the all-helical conformation and the capping box motif lying largely in one part. the second component, on the other hand, could be used to describe conformational transitions between the all-helical form and the capping box motif.
Computational Analysis of the Interactions between the Angiogenesis Inhibitor PD173074 and Fibroblast Growth Factor Receptor 1Journal of Theoretical and Computational Chemistry, Vol. 2, No. 1, pp. 43-56 (2003)
We have carried out computational sensitivity analysis to analyze the interactions between the inhibitor PD173074 and FGFR1 in order to identify the determinants of their recognition and generate insights into further refining the inhibitor. The analysis has identified the parts of the inhibitor that are already useful for binding, e.g., the part that recognizes the linker connecting the N-terminal and C-terminal lobes of the kinase domain. These parts are profitably kept during a lead optimization process. The analysis has also pointed out regions of the inhibitors that may be useful to modify to improve its binding affinity, e.g., the dimethoxyphenyl ring. Comparative structural analysis of the binding pocket of almost 400 protein kinases also suggests that modifying the dimethoxyphenyl moiety might improve selective binding. Selectivity may be achieved not only by introducing groups to the 3 and 5 positions but also to the 1 and 6 positions. Replacing the tertiary amines by hydrocarbon might also improve binding affinity.
We have carried out computational sensitivity analysis to analyze the interactions between the inhibitor PD173074 and FGFR1 in order to identify the determinants of their recognition and generate insights into further refining the inhibitor. The analysis has identified the parts of the inhibitor that are already useful for binding, e.g., the part that recognizes the linker connecting the N-terminal and C-terminal lobes of the kinase domain. These parts are profitably kept during a lead optimization process. The analysis has also pointed out regions of the inhibitors that may be useful to modify to improve its binding affinity, e.g., the dimethoxyphenyl ring. Comparative structural analysis of the binding pocket of almost 400 protein kinases also suggests that modifying the dimethoxyphenyl moiety might improve selective binding. Selectivity may be achieved not only by introducing groups to the 3 and 5 positions but also to the 1 and 6 positions. Replacing the tertiary amines by hydrocarbon might also improve binding affinity.
Finite Element Simulations of Acetylcholine Diffusion in Neuromuscular JunctionsBiophysical Journal, Vol. 84, No. 4, pp. 2234-2241 (2003) [PubMed 12668432]
A robust infrastructure for solving time-dependent diffusion using the finite element package FEtk has been developed to simulate synaptic transmission in a neuromuscular junction with realistic postsynaptic folds. Simplified rectilinear synapse models serve as benchmarks in initial numerical studies of how variations in geometry and kinetics relate to endplate currents associated with fast-twitch, slow-twitch, and dystrophic muscles. The flexibility and scalability of FEtk affords increasingly realistic and complex models that can be formed in concert with expanding experimental understanding from electron microscopy. Ultimately, such models may provide useful insight on the functional implications of controlled changes in processes, suggesting therapies for neuromuscular diseases.
A robust infrastructure for solving time-dependent diffusion using the finite element package FEtk has been developed to simulate synaptic transmission in a neuromuscular junction with realistic postsynaptic folds. Simplified rectilinear synapse models serve as benchmarks in initial numerical studies of how variations in geometry and kinetics relate to endplate currents associated with fast-twitch, slow-twitch, and dystrophic muscles. The flexibility and scalability of FEtk affords increasingly realistic and complex models that can be formed in concert with expanding experimental understanding from electron microscopy. Ultimately, such models may provide useful insight on the functional implications of controlled changes in processes, suggesting therapies for neuromuscular diseases.
Influence of Structural Fluctuation on Enzyme Reaction Energy Barriers in Combined Quantum Mechanical/Molecular Mechanical StudiesJournal of Physical Chemistry B, Vol. 107, No. 18, pp. 4459-4463 (2003)
To account for protein dynamics and to investigate the effect of different conformations on the enzyme reaction energy barrier, we have studied the initial step of the acylation reaction catalyzed by acetylcholinesterase (AChE) with a multiple QM/MM reaction path approach. The approach consists of two main components: generating enzyme-substrate conformations with classical molecular dynamics simulation, and mapping out the minimum reaction energy path for each conformational snapshot with combined quantum mechanical/molecular mechanical (QM/MM) calculations. It is found that enzyme-substrate conformation fluctuations lead to significant differences in the calculated reaction energy barrier; however, the qualitative picture of the role of the catalytic triad and oxyanion hole in AChE catalysis is very consistent. Our results emphasize the importance of employing multiple starting structures in the QM/MM study of enzyme reactions, and indicate that structural fluctuation is an integral part of the enzyme reaction process.
To account for protein dynamics and to investigate the effect of different conformations on the enzyme reaction energy barrier, we have studied the initial step of the acylation reaction catalyzed by acetylcholinesterase (AChE) with a multiple QM/MM reaction path approach. The approach consists of two main components: generating enzyme-substrate conformations with classical molecular dynamics simulation, and mapping out the minimum reaction energy path for each conformational snapshot with combined quantum mechanical/molecular mechanical (QM/MM) calculations. It is found that enzyme-substrate conformation fluctuations lead to significant differences in the calculated reaction energy barrier; however, the qualitative picture of the role of the catalytic triad and oxyanion hole in AChE catalysis is very consistent. Our results emphasize the importance of employing multiple starting structures in the QM/MM study of enzyme reactions, and indicate that structural fluctuation is an integral part of the enzyme reaction process.
A Computational Model of Binding Thermodynamics: The Design of Cyclin-dependent Kinase 2 InhibitorsJournal of Medicinal Chemistry, Vol. 46, No. 15, pp. 3314-3325 (2003) [PubMed 12852762]
The cyclin-dependent protein kinases are important targets in drug discovery because of their role in cell cycle regulation. In this computational study, we have applied a continuum solvent model to study the interactions between cyclin-dependent kinases 2 (CDK2) and analogs of the clinically tested anticancer agent, flavopiridol. The continuum solvent model uses Coulomb's law to account for direct electrostatic interactions, solves the Poisson equation to obtain the electrostatic contributions to solvation energy, and calculates scaled solvent accessible surface area to account for hydrophobic interactions. The computed free energy of binding gauges the strength of protein-ligand interactions. Our model was first validated through a binding study on a number of flavopiridol derivatives to CKD2, and its ability to identify potent inhibitors was observed. The model was then used to aid in the design of novel CDK2 inhibitors with the aid of a computational sensitivity analysis. Some of these hypothetical structures could be significantly more potent than the lead compound, flavopiridol. We applied two approaches to gaining insights into designing selective inhibitors. One relied on the comparative analysis of the binding pocket for several hundred protein kinases to identify the parts of a lead compound whose modifications might lead to selective compounds. The other was based on building and using homology for energy calculations. The homology models appear to be able to classify ligand potency into groups but cannot yet give reliable quantitative results.
The cyclin-dependent protein kinases are important targets in drug discovery because of their role in cell cycle regulation. In this computational study, we have applied a continuum solvent model to study the interactions between cyclin-dependent kinases 2 (CDK2) and analogs of the clinically tested anticancer agent, flavopiridol. The continuum solvent model uses Coulomb's law to account for direct electrostatic interactions, solves the Poisson equation to obtain the electrostatic contributions to solvation energy, and calculates scaled solvent accessible surface area to account for hydrophobic interactions. The computed free energy of binding gauges the strength of protein-ligand interactions. Our model was first validated through a binding study on a number of flavopiridol derivatives to CKD2, and its ability to identify potent inhibitors was observed. The model was then used to aid in the design of novel CDK2 inhibitors with the aid of a computational sensitivity analysis. Some of these hypothetical structures could be significantly more potent than the lead compound, flavopiridol. We applied two approaches to gaining insights into designing selective inhibitors. One relied on the comparative analysis of the binding pocket for several hundred protein kinases to identify the parts of a lead compound whose modifications might lead to selective compounds. The other was based on building and using homology for energy calculations. The homology models appear to be able to classify ligand potency into groups but cannot yet give reliable quantitative results.
Nanosecond Dynamics of the Mouse Acetylcholinesterase Cys69-Cys96 Omega LoopJournal of Biological Chemistry, Vol. 278, Issue 33, pp. 30905-30911 (2003) [PubMed 12759360]
The paradox of high substrate turnover occurring within the confines of a deep, narrow gorge through which acetylcholine must traverse to reach the catalytic site of acetylcholinesterase has suggested the existence of transient gorge enlargements that would enhance substrate accessibility. To establish a foundation for the experimental study of transient fluctuations in structure, site-directed labeling in conjunction with time-resolved fluorescence anisotropy were utilized to assess the possible involvement of the omega loop, a segment that forms the outer wall of the gorge. Specifically, the flexibility of three residues (L76C, E81C, and E84C) in the Cys69-Cys96 omega loop and one residue (Y124C) across the gorge from the omega loop were studied in the absence and presence of two inhibitors of different size, fasciculin and huperzine. Additionally, to validate the approach molecular dynamics was employed to simulate anisotropy decay of the side chains. The results show that the omega loop residues are significantly more mobile than the non-loop residue facing the interior of the gorge. Moreover, fasciculin, which binds at the mouth of the gorge, well removed from the active site, decreases the mobility of IAEDANS reporter groups attached to L76C and Y124C, but increases the mobility of the reporter groups attached to E81C and E84C. Huperzine, which binds at the base of active site gorge, has no effect on the mobility of reporter groups attached to L76C and Y124C, but increases the mobility of reporter groups attached to E81C and E84C. Besides showing that fluctuations of the omega-loop residues are not tightly coupled, the results indicate that residues in the omega loop exhibit distinctive conformational fluctuations and therefore are likely to contribute to transient gorge enlargements in the non-liganded enzyme.
The paradox of high substrate turnover occurring within the confines of a deep, narrow gorge through which acetylcholine must traverse to reach the catalytic site of acetylcholinesterase has suggested the existence of transient gorge enlargements that would enhance substrate accessibility. To establish a foundation for the experimental study of transient fluctuations in structure, site-directed labeling in conjunction with time-resolved fluorescence anisotropy were utilized to assess the possible involvement of the omega loop, a segment that forms the outer wall of the gorge. Specifically, the flexibility of three residues (L76C, E81C, and E84C) in the Cys69-Cys96 omega loop and one residue (Y124C) across the gorge from the omega loop were studied in the absence and presence of two inhibitors of different size, fasciculin and huperzine. Additionally, to validate the approach molecular dynamics was employed to simulate anisotropy decay of the side chains. The results show that the omega loop residues are significantly more mobile than the non-loop residue facing the interior of the gorge. Moreover, fasciculin, which binds at the mouth of the gorge, well removed from the active site, decreases the mobility of IAEDANS reporter groups attached to L76C and Y124C, but increases the mobility of the reporter groups attached to E81C and E84C. Huperzine, which binds at the base of active site gorge, has no effect on the mobility of reporter groups attached to L76C and Y124C, but increases the mobility of reporter groups attached to E81C and E84C. Besides showing that fluctuations of the omega-loop residues are not tightly coupled, the results indicate that residues in the omega loop exhibit distinctive conformational fluctuations and therefore are likely to contribute to transient gorge enlargements in the non-liganded enzyme.
From Model Complexes to Metalloprotein Inhibition: A Synergistic Approach to Structure-based Drug DiscoveryAngewandte Chemie International Edition, Vol. 42, Issue 32, pp. 3772-3774 (2003) [PubMed 12923840]
Matrix metalloproteinases (MMPs) are zinc-containing hydrolytic enzymes that are involved in restructuring the connective tissue. The design of effective inhibitors of matrix metalloproteinases is a significant goal in chemotherapeutic development, due to the correlation of MMP activity with a variety of illnesses including cancer, arthritis, and inflammatory disease. However, the design of inhibitors for MMPs and other metalloproteins is limited by the ability to predict the interaction of a given inhibitor with the active site metal ion. In most cases, the elucidation of a protein structure with the bound inhibitor, by using X-ray diffraction or NMR spectroscopic methods, is necessary for revealing the metal-inhibitor interactions. Innovative approaches to increasing the efficiency and speed of the drug discovery process may provide attractive, alternative routes to identifying new drug candidates. For example, the structure-activity-relationship by nuclear magnetic resonance (SAR by NMR) approach has been used to reveal probable binding modes of metal chelators in order to facilitate the development of improved MMP inhibitors. Despite being a very effective approach, the use of SAR by NMR still requires substantial amounts of 15N-labeled metalloprotein to determine the interactions of potential inhibitor fragments with the macromolecule. In addition, it is likely that SAR by NMR, in terms of metalloprotein drug design, would be limited to metalloenzymes that contain diamagnetic metal centers. Another approach for metalloprotein drug design is to reproduce the drug-metalloenzyme interactions using small molecule models that can be readily characterized. Recently we have shown that synthetic inorganic model complexes can be used to determine metal-ligand interactions relevant to MMP inhibitors.
Matrix metalloproteinases (MMPs) are zinc-containing hydrolytic enzymes that are involved in restructuring the connective tissue. The design of effective inhibitors of matrix metalloproteinases is a significant goal in chemotherapeutic development, due to the correlation of MMP activity with a variety of illnesses including cancer, arthritis, and inflammatory disease. However, the design of inhibitors for MMPs and other metalloproteins is limited by the ability to predict the interaction of a given inhibitor with the active site metal ion. In most cases, the elucidation of a protein structure with the bound inhibitor, by using X-ray diffraction or NMR spectroscopic methods, is necessary for revealing the metal-inhibitor interactions. Innovative approaches to increasing the efficiency and speed of the drug discovery process may provide attractive, alternative routes to identifying new drug candidates. For example, the structure-activity-relationship by nuclear magnetic resonance (SAR by NMR) approach has been used to reveal probable binding modes of metal chelators in order to facilitate the development of improved MMP inhibitors. Despite being a very effective approach, the use of SAR by NMR still requires substantial amounts of 15N-labeled metalloprotein to determine the interactions of potential inhibitor fragments with the macromolecule. In addition, it is likely that SAR by NMR, in terms of metalloprotein drug design, would be limited to metalloenzymes that contain diamagnetic metal centers. Another approach for metalloprotein drug design is to reproduce the drug-metalloenzyme interactions using small molecule models that can be readily characterized. Recently we have shown that synthetic inorganic model complexes can be used to determine metal-ligand interactions relevant to MMP inhibitors.
Calculating pKa Values in Enzyme Active SitesProtein Science, Vol. 12, Issue 9, pp. 1894-1901 (2003) [PubMed 12930989]
The ionization properties of the active-site residues in enzymes are of considerable interest in the study of the catalytic mechanisms of enzymes. Knowledge of these ionization constants (pKa values) often allows the researcher to identify the proton donor and the catalytic nucleophile in the reaction mechanism of the enzyme. Estimates of protein residue pKa values can be obtained by applying pKa calculation algorithms to protein X-ray structures. We show that pKa values accurate enough for identifying the proton donor in an enzyme active site can be calculated by considering in detail only the active-site residues and their immediate electrostatic interaction partners, thus allowing for a large decrease in calculation time. More specifically we omit the calculation of site-site interaction energies, and the calculation of desolvation and background interaction energies for a large number of pairs of titratable groups. The method presented here is well suited to be applied on a genomic scale, and can be implemented in most pKa calculation algorithms to give significant reductions in calculation time with little or no impact on the accuracy of the results. The work presented here has implications for understanding of enzymes in general and for the design of novel biocatalysts.
The ionization properties of the active-site residues in enzymes are of considerable interest in the study of the catalytic mechanisms of enzymes. Knowledge of these ionization constants (pKa values) often allows the researcher to identify the proton donor and the catalytic nucleophile in the reaction mechanism of the enzyme. Estimates of protein residue pKa values can be obtained by applying pKa calculation algorithms to protein X-ray structures. We show that pKa values accurate enough for identifying the proton donor in an enzyme active site can be calculated by considering in detail only the active-site residues and their immediate electrostatic interaction partners, thus allowing for a large decrease in calculation time. More specifically we omit the calculation of site-site interaction energies, and the calculation of desolvation and background interaction energies for a large number of pairs of titratable groups. The method presented here is well suited to be applied on a genomic scale, and can be implemented in most pKa calculation algorithms to give significant reductions in calculation time with little or no impact on the accuracy of the results. The work presented here has implications for understanding of enzymes in general and for the design of novel biocatalysts.
Brownian Dynamics Simulations of Ion Atmospheres around Polyalanine and B-DNA: Effects of Biomolecular DielectricBiopolymers, Vol. 70, Issue 3, pp. 391-402 (2003) [PubMed 14579311]
We have extended an earlier Brownian dynamics simulation algorithm for simulating the structural dynamics of ions around biomolecules to accommodate dielectric inhomogeneity. The electrostatic environment of a biomolecule immersed in water was obtained by numerically solving the Poisson equation with the biomolecule treated as a low dielectric region and the solvent treated as a high dielectric region. Instead of using the mean-field type approximations of ion interactions as in the Poisson-Boltzmann model, the ions were treated explicitly by allowing them to evolve dynamically under the electrostatic field of the biomolecule. This model thus accounts for ion-ion correlations and the finite-size effects of the ions. For a 13-residue α-helical polyalanine and a 12-bp B-form DNA, we found that the choice of the dielectric constant of the biomolecule has much larger effects on the mean ionic structure around the biomolecule than on the fluctuational and dynamical properties of the ions surrounding the biomolecule.
We have extended an earlier Brownian dynamics simulation algorithm for simulating the structural dynamics of ions around biomolecules to accommodate dielectric inhomogeneity. The electrostatic environment of a biomolecule immersed in water was obtained by numerically solving the Poisson equation with the biomolecule treated as a low dielectric region and the solvent treated as a high dielectric region. Instead of using the mean-field type approximations of ion interactions as in the Poisson-Boltzmann model, the ions were treated explicitly by allowing them to evolve dynamically under the electrostatic field of the biomolecule. This model thus accounts for ion-ion correlations and the finite-size effects of the ions. For a 13-residue α-helical polyalanine and a 12-bp B-form DNA, we found that the choice of the dielectric constant of the biomolecule has much larger effects on the mean ionic structure around the biomolecule than on the fluctuational and dynamical properties of the ions surrounding the biomolecule.
The Physical Basis of Microtubule Structure and StabilityProtein Science, Vol. 12, Issue 10, pp. 2257-2261 (2003) [PubMed 14500883]
Microtubules are cylindrical polymers found in every eukaryotic cell. They have a unique helical structure that has implications at both the cellular level, in terms of the functions they perform, and multi-cellular level, such as determining the left-right asymmetry in plants. Through the combination of an atomically-detailed model for a microtubule and new large scale computational techniques for computing electrostatic interactions, we are able to explain the observed microtubule structure. Based on the lateral interactions between protofilaments, we have determined that a helix with a subunit rise of 8-9 Å is the most favourable configuration. Further, we find that these lateral bonds are significantly weaker than the longitudinal bonds along protofilaments. This explains observations of microtubule disassembly and may serve as another step toward understanding the basis for dynamic instability.
Microtubules are cylindrical polymers found in every eukaryotic cell. They have a unique helical structure that has implications at both the cellular level, in terms of the functions they perform, and multi-cellular level, such as determining the left-right asymmetry in plants. Through the combination of an atomically-detailed model for a microtubule and new large scale computational techniques for computing electrostatic interactions, we are able to explain the observed microtubule structure. Based on the lateral interactions between protofilaments, we have determined that a helix with a subunit rise of 8-9 Å is the most favourable configuration. Further, we find that these lateral bonds are significantly weaker than the longitudinal bonds along protofilaments. This explains observations of microtubule disassembly and may serve as another step toward understanding the basis for dynamic instability.
The Dynamics of Ligand Barrier Crossing inside the Acetylcholinesterase GorgeBiophysical Journal, Vol. 85, No. 4, pp. 2267-2272 (2003) [PubMed 14507691]
The dynamics of ligand movement through the constricted region of the acetylcholinesterase gorge is important in understanding how the ligand gains access to and is released from the active site of the enzyme. Molecular dynamics simulations of the simple ligand, tetramethylammonium, crossing this bottleneck region are conducted using umbrella potential sampling and activated flux techniques. The low potential of mean force obtained is consistent with the fast reaction rate of acetylcholinesterase observed experimentally. From the results of the activated dynamics simulations, local conformational fluctuations of the gorge residues and larger scale collective motions of the protein are found to correlate highly with the ligand crossing.
The dynamics of ligand movement through the constricted region of the acetylcholinesterase gorge is important in understanding how the ligand gains access to and is released from the active site of the enzyme. Molecular dynamics simulations of the simple ligand, tetramethylammonium, crossing this bottleneck region are conducted using umbrella potential sampling and activated flux techniques. The low potential of mean force obtained is consistent with the fast reaction rate of acetylcholinesterase observed experimentally. From the results of the activated dynamics simulations, local conformational fluctuations of the gorge residues and larger scale collective motions of the protein are found to correlate highly with the ligand crossing.
Asymmetric Structural Motions of the Homomeric α7 Nicotinic Receptor Ligand Binding Domain Revealed by Molecular Dynamics SimulationBiophysical Journal, Vol. 85, No. 5, pp. 3007-3018 (2003) [PubMed 14581202]
A homology model of the ligand binding domain of the α7 nicotinic receptor is constructed based on the AChBP crystal structure. This structure is refined in a 10 ns molecular dynamics simulation. The model structure proves fairly resilient, with no significant changes at the secondary or tertiary structural levels. The hypothesis that the AChBP template is in the activated or desensitized state and the absence of a bound agonist in the simulation suggests that the structure may also be relaxing from this state to the activatable state. Candidate motions that take place involve not only the side chains of residues lining the binding sites, but also the subunit positions that determine the overall shape of the receptor. In particular, two non-adjacent subunits move outward while their partners counter-clockwise to them move inward, leading to a marginally wider interface between themselves and an overall asymmetric structure. This in turn affects the binding sites, producing two that are more open and characterized by distinct side chain conformations of W54 and L118, although motions of the side chains of all residues in every binding site still contribute to a reduction in binding site size, especially the outward motion of W148, which hinders ACh binding. The Cys loop at the membrane interface also displays some flexibility. While the short simulation time scale is unlikely to sample adequately all the conformational states, the pattern of observed motions suggests how ligand binding may correlate with larger scale subunit motions that would connect with the transmembrane region that controls the passage of ions. Further, the shape of the asymmetry with binding sites of differing affinity for ACh, characteristic of other nicotinic receptors, may be a natural property of the relaxed, activatable state of α7.
A homology model of the ligand binding domain of the α7 nicotinic receptor is constructed based on the AChBP crystal structure. This structure is refined in a 10 ns molecular dynamics simulation. The model structure proves fairly resilient, with no significant changes at the secondary or tertiary structural levels. The hypothesis that the AChBP template is in the activated or desensitized state and the absence of a bound agonist in the simulation suggests that the structure may also be relaxing from this state to the activatable state. Candidate motions that take place involve not only the side chains of residues lining the binding sites, but also the subunit positions that determine the overall shape of the receptor. In particular, two non-adjacent subunits move outward while their partners counter-clockwise to them move inward, leading to a marginally wider interface between themselves and an overall asymmetric structure. This in turn affects the binding sites, producing two that are more open and characterized by distinct side chain conformations of W54 and L118, although motions of the side chains of all residues in every binding site still contribute to a reduction in binding site size, especially the outward motion of W148, which hinders ACh binding. The Cys loop at the membrane interface also displays some flexibility. While the short simulation time scale is unlikely to sample adequately all the conformational states, the pattern of observed motions suggests how ligand binding may correlate with larger scale subunit motions that would connect with the transmembrane region that controls the passage of ions. Further, the shape of the asymmetry with binding sites of differing affinity for ACh, characteristic of other nicotinic receptors, may be a natural property of the relaxed, activatable state of α7.
Studying the roles of W86, E202, and Y337 in binding of acetylcholine to acetylcholinesterase using a combined molecular dynamics and multiple docking approachProtein Science, Vol. 12, Issue 12, pp. 2675-2684 (2003) [PubMed 14627729]
A combined molecular dynamics simulation and multiple ligand docking approach is applied to study the roles of the anionic subsite residues (W86, E202, Y337) in the binding of acetylcholine (ACh) to acetylcholinesterase (AChE). We find that E202 stabilizes docking of ACh via electrostatic interactions. However, we find no significant electrostatic contribution from the aromatic residues. Docking energies of ACh to mutant AChE show a more pronounced effect due to size/shape complementarity. Mutating to smaller residues results in poorer binding, both in terms of docking energy and statistical docking probability. Besides separating out electrostatics by turning off the partial charges from each residue and comparing to the native, the mutations in this study are W86F, W86A, E202D, E202Q, E202A, Y337F, and Y337A. We also find that all perturbations result in a significant reduction in binding of extended ACh in the catalytically productive orientation. This effect is primarily due to a small shift in preferred position of the quaternary tail.
A combined molecular dynamics simulation and multiple ligand docking approach is applied to study the roles of the anionic subsite residues (W86, E202, Y337) in the binding of acetylcholine (ACh) to acetylcholinesterase (AChE). We find that E202 stabilizes docking of ACh via electrostatic interactions. However, we find no significant electrostatic contribution from the aromatic residues. Docking energies of ACh to mutant AChE show a more pronounced effect due to size/shape complementarity. Mutating to smaller residues results in poorer binding, both in terms of docking energy and statistical docking probability. Besides separating out electrostatics by turning off the partial charges from each residue and comparing to the native, the mutations in this study are W86F, W86A, E202D, E202Q, E202A, Y337F, and Y337A. We also find that all perturbations result in a significant reduction in binding of extended ACh in the catalytically productive orientation. This effect is primarily due to a small shift in preferred position of the quaternary tail.
