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

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

Binding of tacrine and 6-chlorotacrine by acetylcholinesterase

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

Determination of the pKa values of titratable groups of an antigen-antibody complex, HyHEL-5-hen egg lysozyme

The titration behavior of the ionizable residues of the HyHEL-5-hen egg lysozyme complex and its individual components has been studied using continuum electrostatic calculations. Several residues of HyHEL-5 had pKa values shifted away from model values for isolated residues by more than three pH units. Shifts away from the model values were smaller for the residues of hen egg lysozyme. A moderate variation in the pKa values of the titratable groups was observed upon increase of the ionic strength from 0 to 100 mM, amounting to 1-2 pH units in most cases. Under physiological conditions, the net charge of HyHEL-5 was opposite that for hen egg lysozyme. Several residues, including those involved in the Arg-Glu salt bridges that have been proposed to be important in antibody-antigen binding, had pKa values that were changed significantly upon binding. The main titration event upon antibody-antigen binding appears to be loss of a proton from residue GluH50 of the Fv molecule. The limitations of our calculation methods and the role they might play in the design of antibodies for use in assays, sensors and separations are discussed.

Free energy simulations of the HyHEL-10/HEL antibody-antigen complex

Free energy simulations are reported for the N31LD mutation, both in the HyHEL-10-HEL antibody-lysozyme complex and in the unliganded antibody, using the thermodynamic-cycle perturbation method. The present study suggests that the mutation would change the free energy of binding of the complex by -5.6 kcal/mol (unrestrained free energy simulations), by -0.5 kcal/mol (free energy simulations with a restrained backbone) and by 1.8 kcal/mol (Poisson-Boltzmann calculations, which also use a restrained geometry model). A detailed structural analysis helps in estimating the contributions from various residues and regions of the system. Enhanced recognition of HEL by the mutant HyHEL-10 would arise from the combination of thermodynamically more favorable conformational changes of the CDR loops upon association and subsequent charge pairing with Lys96 in the antigen.

Electrostatic and hydrodynamic orientational steering effects in enzyme-substrate association

Diffusional encounters between a dumbbell model of a cleft enzyme and a dumbbell model of an elongated ligand are simulated by Brownian dynamics. The simulations take into account electrostatic and hydrodynamic interactions between the molecules. It is shown that the primary effect of inclusion of hydrodynamic interactions into the simulation is an overall decrease in the rate constant. Hydrodynamic orientational effects are of modest size for the systems considered here. They are manifested when changes in the rate constants for diffusional encounters favored by hydrodynamic interactions are compared with those favored by electrostatic interactions as functions of the overall strength of electrostatic interactions. The electrostatic interactions modify the hydrodynamic torques by modifying the drift velocity of the substrate toward the enzyme. We conclude that simulations referring only to electrostatic interactions between an enzyme and its ligand may yield rate constants that are somewhat (e.g., 20%) too high, but provide realistic descriptions of the orientational steering effects in the enzyme-ligand encounters.

Simulation of charge-mutant acetylcholinesterases

A recent experimental study of human acetylcholinesterase has shown that the mutation of surface acidic residues has little effect on the rate constant for hydrolysis of acetylthiocholine. It was concluded, on this basis, that the reaction is not diffusion controlled and that electrostatic steering plays only a minor role in determining the rate. Here we examine this issue through Brownian dynamics simulations on Torpedo californica acetylcholinesterase in which the surface acidic residues homologous with those mutated in the human enzyme are artificially neutralized. The computed effects of the mutations on the rate constants reproduce quite well the modest effects of the mutations upon the measured encounter rates. Nonetheless, the electrostatic field of the enzyme is found to increase the rate constants by about an order of magnitude in both the wild type and the mutants. We therefore conclude that the mutation experiments do not disprove that electrostatic steering substantially affects the catalytic rate of acetylcholinesterase.

Secondary structure prediction of the H5 pore of potassium channels

The 'H5' segment located between the putative fifth and sixth transmembrane helices is the most highly conserved region in voltage-gated potassium channels and it is believed to constitute a major part of the ion conduction path (pore). Here we present a two-step procedure, comprising secondary structure prediction and hydrophobic moment profiling, to predict the structure of this important region. Combined results from the application of the procedure to the H5 region of four voltage-gated and five other K+ channel sequences lead to the prediction of a beta-strand-turn-beta-strand structure for H5. The reasons for the application of these soluble protein methods to parts of membrane proteins are: (i) that pore-lining residues are accessible to water and (ii) that a large enough database of high-resolution membrane protein structures does not yet exist. The results are compared with experimental results, in particular spectroscopic studies of two peptides based on the H5 sequence of SHAKER potassium channel. The procedure developed here may be applicable to water-accessible regions of other membrane proteins.

Conservative and nonconservative mutations in proteins: anomalous mutations in a transport receptor analyzed by free energy and quantum chemical calculations

Experimental studies on a bacterial sulfate receptor have indicated anomalous relative binding affinities for the mutations Ser130-->Cys,Ser130-->Gly, and Ser130-->Ala. The loss of affinity for sulfate in the former mutation was previously attributed to a greater steric effect on the part of the Cys side chain relative to the Ser side chain, whereas the relatively small loss of binding affinity for the latter two mutations was attributed to the loss of a single hydrogen bond. In this report we present quantum chemical and statistical thermodynamic studies of these mutations. Qualitative results from these studies indicate that for the Ser130-->Cys mutation the large decrease in binding affinity is in part caused by steric effects, but also significantly by the differential work required to polarize the Cys thiol group relative to the Ser hydroxyl group. The Gly mutant cobinds a water molecule in the same location as the Ser side chain resulting in a relatively small decrease in binding affinity. Results for the Ala mutant are in disagreement with experimental results but are likely to be limited by insufficient sampling of configuration space due to physical constraints applied during the simulation.

Acetylcholinesterase: diffusional encounter rate constants for dumbbell models of ligand

For some enzymes, virtually every substrate molecule that encounters the entrance to the active site proceeds to reaction, at low substrate concentrations. Such diffusion-limited enzymes display high apparent bimolecular rate constants ((kcat/KM)), which depend strongly upon solvent viscosity. Some experimental studies provide evidence that acetylcholinesterase falls into this category. Interestingly, the asymmetric charge distribution of acetylcholinesterase, apparent from the crystallographic structure, suggests that its electrostatic field accelerates the encounter of its cationic substrate, acetylcholine, with the entrance to the active site. Here we report simulations of the diffusion of substrate in the electrostatic field of acetylcholinesterase. We find that the field indeed guides the substrate to the mouth of the active site. The computed encounter rate constants depend upon the particular relative geometries of substrate and enzyme that are considered to represent successful encounters. With loose reaction criteria, the computed rates exceed those measured experimentally, but the rate constants vary appropriately with ionic strength. Although more restrictive reaction criteria lower the computed rates, they also lead to unrealistic variation of the rate constants with ionic strength. That these simulations do not agree well with experiment suggests that the simple diffusion model is incomplete. Structural fluctuations in the enzyme or events after the encounter may well contribute to rate limitation.

Prediction of pH-dependent properties of proteins

We describe what may be the most accurate approach currently available for the calculation of the pKas of ionizable groups in proteins. The accuracy is assessed by comparison of computed pKas with 60 measured pKas in a total of seven proteins. The overall root-mean-square error is 0.89 pKa units. Linear regression analysis of computed versus measured pKas yields a slope of 0.95, y-intercept of -0.02 and a correlation coefficient of 0.96. The proposed approach also picks out many of the shifted pKas of groups in enzyme active sites and special salt bridges. However, it does yield several over-shifted pKas and tends to underestimate pKa shifts which result from desolvation effects. We examine the ability of the new approach to reproduce the dependence of protein stability upon pH, using the ionization polynomial formalism. Overall features of the stability curves are reproduced, but the quantitative agreement is not particularly good. The reasons for the disagreement may have to do both with insufficient accuracy in the theory and with uncertainty in the nature of the unfolded state of proteins. The methodology described here is based upon finite difference solutions of the Poisson-Boltzmann equation. Its success depends upon the use of the rather high protein dielectric constant of 20. However, theoretical considerations and the fact that pKa shifts which result from desolvation are underestimated here imply that the dielectric constant of the protein interior actually is lower than 20. We suggest that the high protein dielectric constant improves the overall agreement with experiment because it accounts approximately for phenomena which tend to mitigate pKa shifts and which are not specifically included in the model. These include conformational relaxation and specific ion-binding. Future models based upon a low protein dielectric constant and treating such phenomena explicitly might yield improved agreement with experiment.

Combined conformational search and finite-difference Poisson-Boltzmann approach for flexible docking. Application to an operator mutation in the lambda repressor-operator complex

The N-terminal domain of the phage lambda repressor binds as a dimer to its palindromic DNA operator sequence. In addition to a helix-turn-helix DNA recognition motif, the first six amino acids of the phage lambda repressor form a flexible peptide segment which wraps around DNA. Site-directed mutagenesis studies have shown that amino acid replacements or partial removal of the arm structure, or changes in the DNA sequence contacting the N-terminal arm, can lower the repressor-operator binding affinity by several orders of magnitude. The finite-difference Poisson-Boltzmann approach in combination with a conformational search procedure was used to study energetic contributions of the lambda arm to repressor-operator recognition based on the high resolution X-ray structure. It allows for the local relaxation of the structure upon changing the DNA sequence in the lambda arm binding region. A simplified potential energy function including torsional, truncated Lennard-Jones and approximate electrostatic terms is used in the initial step to screen out energetically unfavorable structures. The electrostatic energy of selected conformations is subsequently calculated more accurately using the finite-difference Poisson-Boltzmann approach. The method was applied to study the effect of a C-->T mutation at position 6 of the consensus half-site of the operator. This base-pair contacts Lys4 which is part of the arm segment. Keeping only the Lys4 side-chain mobile and with the wild-type DNA operator sequence, several conformations close to the X-ray structure were identified as those with lowest energy. In the case of the DNA mutation, lowest energy conformations differed significantly from those selected for the wild-type sequence. These initial calculations indicate that the approach might be a useful tool to estimate conformational and energetic effects upon mutagenesis of protein-DNA complexes.

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

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

Simulation of enzyme-substrate encounter with gated active sites

We describe a brownian dynamics simulation method that allows investigation of the effects of receptor flexibility on ligand binding rates. The method is applied to the encounter of substrate, glyceraldehyde 3-phosphate, with triose phosphate isomerase, a diffusion-controlled enzyme with flexible peptide loops at its active sites. The simulations show that while the electrostatic field surrounding the enzyme steers the substrate into its active sites, the flexible loops appear to have little influence on the substrate binding rate. The dynamics of the loops may therefore have been optimized during evolution to minimize their interference with the substrate's access to the active sites. The calculated and experimental rate constants are in good agreement.

Molecular modeling methods. Basic techniques and challenging problems

An overview is presented of computer modeling and simulation methods that play an increasing role in drug design: quantum chemical methods, molecular mechanics, molecular dynamics and Brownian dynamics. The application of molecular dynamics for the prediction of thermodynamic properties like free energy differences and binding constants is discussed. The Brownian dynamics method is presented in connection with the calculation of effective electrostatic forces using the Poisson-Boltzmann equation, which allows one to sample ligand-binding geometries and to predict the kinetics of diffusion-limited enzyme reactions. New techniques that have recently been extensively developed, such as the global energy minimization and quantum-classical dynamics methods, are also introduced. The molecular modeling methods are illustrated with selected examples.

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

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

Acetylcholinesterase: electrostatic steering increases the rate of ligand binding

Brownian dynamics simulations have been used to calculate the diffusion-controlled rate constants for the binding of a positively charged ligand to several models of acetylcholinesterase (AChE). The crystal structure was used to define the detailed topography and the active sites of the dimeric enzyme. The electric field around AChE was then computed by solving the Poisson equation for different charge distributions in the enzyme at zero ionic strength. These fields were used in turn to calculate the forces on the diffusing ligand. Significant increases in the rate constant resulted in going from a model with no charges to one with the net charges concentrated at the centers of the monomers and then to a model with a realistic distribution of charges throughout the enzyme. The results show that electrostatic steering of ligands contributes to the high rate constants that are observed experimentally for AChE.

Gating of the active site of triose phosphate isomerase: Brownian dynamics simulations of flexible peptide loops in the enzyme

The enzyme triose phosphate isomerase has flexible peptide loops at its active sites. The loops close over these sites upon substrate binding, suggesting that the dynamics of the loops could be of mechanistic and kinetic importance. To investigate these issues, the loop motions in the dimeric enzyme were simulated by Brownian dynamics. The two loops, one on each monomer, were represented by linear chains of appropriately parameterized spheres, each sphere corresponding to an amino acid residue. The loops moved in the electrostatic field of the rest of the enzyme, which was held rigid in its crystallographically observed conformation. In the absence of substrate, the loops exhibited gating of the active site with a period of about 1 ns and occupied "closed" conformations for about half of the time. As the period of gating is much shorter than the enzyme-substrate relaxation time, the motion of the loops does not reduce the rate constant for the approach of substrate from its simple diffusion-controlled value. This suggests that the flexible loops may have evolved to create the appropriate environment for catalysis while, at the same time, minimizing the kinetic penalty for gating the active site.

Poisson-Boltzmann analysis of the lambda repressor-operator interaction

A theoretical study of the ion atmosphere contribution to the binding free energy of the lambda repressor-operator complex is presented. The finite-difference form of the Poisson-Boltzmann equation was solved to calculate the electrostatic interaction energy of the amino-terminal domain of the lambda repressor with a 9 or 45 base pair oligonucleotide. Calculations were performed at various distances between repressor and operator as well as at different salt concentrations to determine ion atmosphere contributions to the total electrostatic interaction. Details in the distribution of charges on DNA and protein atoms had a strong influence on the calculated total interaction energies. In contrast, the calculated salt contributions are relatively insensitive to changes in the details of the charge distribution. The results indicate that the ion atmosphere contribution favors association at all protein-DNA distances studied. The theoretical number of ions released upon repressor-operator binding appears to be in reasonable agreement with experimental data.

Binding of an antiviral agent to a sensitive and a resistant human rhinovirus. Computer simulation studies with sampling of amino acid side-chain conformation. I. Mapping the rotamers of residue 188 of viral protein 1

The mutation of valine 188 to leucine in the viral protein 1 of human rhinovirus 14 renders the virus resistant to certain antiviral compounds. Thermodynamic-cycle perturbation theory provides a means of calculating the difference in the binding free energies of an antiviral compound to the wild-type virus and to the mutant virus. In calculating the relevant free-energy differences in molecular dynamics simulations, it is important to sample the multiple rotational isomers of residue 188 correctly. In general, these rotamers will not be fully sampled during a single molecular dynamics simulation. However, the contributions of all the rotamers to the free-energy differences associated with mutation of residue 188 may be considered explicitly once they have been identified and their relative free energies determined. Therefore, we describe here the mapping of the rotamers of residue 188 by steric-bump search and energy minimization techniques, and by the computation of potentials of mean force (p.m.f.s.) using umbrella sampling. The usefulness, validity and efficiency of these methods of examining rotameric states is discussed. Adiabatic mapping by energy minimization was found to be unreliable for this residue due to the small magnitude of its interactions with the surrounding protein atoms. Ambiguities in the adiabatic maps were resolved by computing p.m.f.s. The p.m.f. for valine 188 in the unliganded wild-type virus shows a minimum corresponding to the crystallographically observed conformation of valine 188. The p.m.f.s. for valine 188 in the liganded virus and for leucine 188 in the unliganded mutant virus suggest that the experimentally observed conformations may be interpreted as averages of a number of conformations corresponding to those at the minima in the p.m.f.s. The calculations suggest also that the conformation of leucine 188 may change when the ligand binds. The use of the calculated p.m.f.s. to compute the difference in the free energy of binding of an antiviral compound to the wild-type and mutant rhinoviruses is described in the accompanying article.

Binding of an antiviral agent to a sensitive and a resistant human rhinovirus. Computer simulation studies with sampling of amino acid side-chain conformations. II. Calculation of free-energy differences by thermodynamic integration

Thermodynamic-cycle perturbation theory and molecular dynamics simulations were used to calculate the difference in the free energy of binding of the antiviral compound WIN53338 to the wild-type human rhinovirus 14 and to a drug-resistant mutant of the virus in which valine 188 of the viral protein 1 is mutated to leucine. Because of the difficulty of achieving adequate sampling of all of the rotational isomers of amino acid side-chains in molecular dynamics simulations, an explicit treatment of the effects of the existence of multiple rotational isomers of residue 188 on the calculated free energies was used. The rotamers of residue 188 were first mapped by steric and energetic techniques as described in the accompanying article. Thermodynamic integration was then carried out during simulations of the virus, both with and without the antiviral compound bound, by mutating residue 188 while restraining its side-chain to one conformation. The contributions of the other rotamers of residue 188 to the free-energy changes for this mutation were then added to those calculated by thermodynamic integration as correction factors. Binding of WIN53338 to the wild-type virus was calculated to be favored over binding to the mutant virus by 1.7(+/- 3.0) kcal/mol. This is consistent with experimental data which, if differences in activity are assumed to be due to differences in binding, indicate that the binding affinity of WIN53338 for the wild-type virus is at least 0.15 to 1.7 kcal/mol greater than for the mutant virus. Thermodynamic integration was also performed in the conventional manner without restraints and was found to give less accurate results.

Anti-insulin antibody structure and conformation. I. Molecular modeling and mechanics of an insulin antibody

A knowledge-based three-dimensional model of an anti-insulin antibody, 125, was constructed using the structures of conserved residues found in other known crystallographic immunoglobulins. Molecular modeling and mechanics were done with the 125 amino acid sequences using QUANTA and CHARMm on a Silicon Graphics 4D70GT workstation. A minimal model was made by scaffolding using crystallography coordinates of the antibody HyHEL-5, because it had the highest amino acid sequence homology with 125 (84% light chain, 65% heavy chain). The three hypervariable loop turns that are longer in 125 than in HyHEL-5 (L1, L3, and H3) were modeled separately and incorporated into the HyHEL-5 structure; then other amino acid substitutions were made and torsions optimized. The 125 model maintains all the structural attributes of an antibody and the structures conserved in known antibodies. Although there are many polar amino acids (especially serines) in this site, the overall van der Waals surface shape is determined by positions of aromatic side chains. Based on this model, it is suggested that hydrogen bonding may be key in the interaction between the human insulin A chain loop antigenic epitope and 125.

Anti-insulin antibody structure and conformation. II. Molecular dynamics with explicit solvent

Molecular dynamics at 300 K was used as a conformation searching tool to analyze a knowledge-based structure prediction of an anti-insulin antibody. Solvation effects were modeled by packing water molecules around the antigen binding loops. Some loops underwent backbone and side-chain conformational changes during the 95-ps equilibration, and most of these new, lower potential energy conformations were stable during the subsequent 200-ps simulation. Alterations to the model include changes in the intraloop, main-chain hydrogen bonding network of loop H3, and adjustments of Tyr and Lys side chains of H3 induced by hydrogen bonding to water molecules. The structures observed during molecular dynamics support the conclusion of the previous paper that hydrogen bonding will play the dominant role in antibody-insulin recognition. Determination of the structure of the antibody by x-ray crystallography is currently being pursued to provide an experimental test of these results. The simulation appears to improve the model, but longer simulations at higher temperatures should be performed.