Theory of biomolecular recognition

Specific, noncovalent binding of biomolecules can only be understood by considering structural, thermodynamic, and kinetic issues. The theoretical foundations for such analyses have been clarified in the past year. Computational techniques for both particle-based and continuum models continue to improve and to yield useful insights into an ever wider range of biomolecular systems.

Solvation studies of DMP323 and A76928 bound to HIV protease: analysis of water sites using grand canonical Monte Carlo simulations

We examine the water solvation of the complex of the inhibitors DMP323 and A76928 bound to HIV-1 protease through grand canonical Monte Carlo simulations, and demonstrate the ability of this method to reproduce crystal waters and effectively predict water positions not seen in the DMP323 or A76928 structures. The simulation method is useful for identifying structurally important waters that may not be resolved in the crystal structures. It can also be used to identify water positions around a putative drug candidate docked into a binding pocket. Knowledge of these water positions may be useful in designing drugs to utilize them as bridging groups or displace them in the binding pocket. In addition, the method should be useful in finding water sites in homology models of enzymes for which crystal structures are unavailable.

Electrostatic channeling of substrates between enzyme active sites: comparison of simulation and experiment

Recent simulation work has indicated that channeling of charged substrates between the active sites of bifunctional enzymes or bienzyme complexes can be significantly enhanced by favorable interactions with the electrostatic field of the enzymes. The results of such simulations are expressed in terms of transfer efficiencies, which describe the probability that substrate leaving the active site of the first enzyme will reach the active site of the second enzyme before escaping out into bulk solution. The experimental indicators of channeling, on the other hand, are factors such as a decrease in the transient (lag) time for appearance of the final product of the coupled enzyme reaction or a decrease in the susceptibility of the overall reaction rate to the presence of competing enzymes or competitive inhibitors. The work reported here aims to establish a connection between the transfer efficiencies obtained from simulation, with the above-mentioned experimental observables. This is accomplished by extending previously reported analytical approaches to combine the simulated transfer efficiency with the Michaelis-Menten kinetic parameters Km and Vmax of the enzymes involved; expressions are derived to allow both transient times and steady state rates to be calculated. These results are applied to the two systems that have been studied both theoretically and experimentally. In the first case, that of the bifunctional enzyme dihydrofolate reductase-thymidylate synthase (DHFR-TS), the experimentally observed decrease in transient times is found to be consistent with a transfer efficiency of >/=80%. In the second case, that of a citrate synthase-malate dehydrogenase fusion protein, a transfer efficiency of 73% is consistent with the experimental transient time measurements. Separate and independent analysis of the effects of adding the competing enzyme aspartate aminotransferase gives a transfer efficiency of 69%, in excellent agreement with the transient time results. The transfer efficiencies thus obtained from experimental results are in both cases in good agreement with those obtained from simulations that include electrostatic interactions. One important discrepancy between simulation and experiment, is however, found in the reported effects of adding a competitive inhibitor in the DHFR-TS system: qualitatively different results are expected from the theoretical analysis. A possible reason for this apparent contradiction is discussed.

pH dependence of antibody/lysozyme complexation

Association between proteins often depends on the pH and ionic strength conditions of the medium in which it takes place. This is especially true in complexation involving titratable residues at the complex interface. Continuum electrostatics methods were used to calculate the pH-dependent energetics of association of hen egg lysozyme with two closely related monoclonal antibodies raised against it and the association of these antibodies against an avian species variant. A detailed analysis of the energetic contributions reveals that even though the hallmark of association in the two complexes is the presence of conserved charged-residue interactions, the environment of these interactions significantly influences the titration behavior and concomitantly the energetics. The contributing factors include minor structural rearrangements, buried interfacial area, dielectric environment of the key titratable residues, and geometry of the residue dispositions. Modeled structures of several mutant complexes were also studied so as to further delineate the contribution of individual factors to the titration behavior.

Kinase conformations: a computational study of the effect of ligand binding

Protein function is often controlled by ligand-induced conformational transitions. Yet, in spite of the increasing number of three-dimensional crystal structures of proteins in different conformations, not much is known about the driving forces of these transitions. As an initial step toward exploring the conformational and energetic landscape of protein kinases by computational methods, intramolecular energies and hydration free energies were calculated for different conformations of the catalytic domain of cAMP-dependent protein kinase (cAPK) with a continuum (Poisson) model for the electrostatics. Three protein kinase crystal structures for ternary complexes of cAPK with the peptide inhibitor PKI(5-24) and ATP or AMP-PNP were modeled into idealized intermediate and open conformations. Concordant with experimental observation, we find that the binding of PKI(5-24) is more effective in stabilizing the closed and intermediate forms of cAPK than ATP. PKI(5-24) seems to drive the final closure of the active site cleft from intermediate to closed state because ATP does not distinguish between these two states. Binding of PKI(5-24) and ATP is energetically additive.

Molecular properties of amphotericin B membrane channel: a molecular dynamics simulation

Amphotericin B is a powerful but toxic antifungal antibiotic that is used to treat systemic infections. It forms ionic membrane channels in fungal cells. These antibiotic/sterol channels are responsible for the leakage of ions, which causes cell destruction. The detailed molecular properties and structure of amphotericin B channels are still unknown. In the current study, two molecular dynamic simulations were performed of a particular model of amphotericin B/cholesterol channel. The water and phospholipid environment were included in our simulations, and the results obtained were compared with available experimental data. It was found that it is mainly the hydrogen bonding interactions that keep the channel stable in its open form. Our study also revealed the important role of the intermolecular interactions among the hydroxyl, amino, and carboxyl groups of the channel-forming molecules; in particular, some hydroxyl groups stand out as new "hot spots" that are potentially useful for chemotherapeutic investigations. Our results also help to clarify why certain antibiotic derivatives, with a blocked amino group, are less active. We present a hypothesis for the role of membrane lipids and cholesterol in the channel.

Electrostatic influence on the kinetics of ligand binding to acetylcholinesterase. Distinctions between active center ligands and fasciculin

To explore the role that surface and active center charges play in electrostatic attraction of ligands to the active center gorge of acetylcholinesterase (AChE), and the influence of charge on the reactive orientation of the ligand, we have studied the kinetics of association of cationic and neutral ligands with the active center and peripheral site of AChE. Electrostatic influences were reduced by sequential mutations of six surface anionic residues outside of the active center gorge (Glu-84, Glu-91, Asp-280, Asp-283, Glu-292, and Asp-372) and three residues within the active center gorge (Asp-74 at the rim and Glu-202 and Glu-450 at the base). The peripheral site ligand, fasciculin 2 (FAS2), a peptide of 6.5 kDa with a net charge of +4, shows a marked enhancement of rate of association with reduction in ionic strength, and this ionic strength dependence can be markedly reduced by progressive neutralization of surface and active center gorge anionic residues. By contrast, neutralization of surface residues only has a modest influence on the rate of cationic m-trimethylammoniotrifluoroacetophenone (TFK+) association with the active serine, whereas neutralization of residues in the active center gorge has a marked influence on the rate but with little change in the ionic strength dependence. Brownian dynamics calculations for approach of a small cationic ligand to the entrance of the gorge show the influence of individual charges to be in quantitative accord with that found for the surface residues. Anionic residues in the gorge may help to orient the ligand for reaction or to trap the ligand. Bound FAS2 on AChE not only reduces the rate of TFK+ reaction with the active center but inverts the ionic strength dependence for the cationic TFK+ association with AChE. Hence it appears that TFK+ must traverse an electrostatic barrier at the gorge entry imparted by the bound FAS2 with its net charge of +4.

Modeling the electrophoresis of lysozyme. II. Inclusion of ion relaxation

In this work, boundary element methods are used to model the electrophoretic mobility of lysozyme over the pH range 2-6. The model treats the protein as a rigid body of arbitrary shape and charge distribution derived from the crystal structure. Extending earlier studies, the present work treats the equilibrium electrostatic potential at the level of the full Poisson-Boltzmann (PB) equation and accounts for ion relaxation. This is achieved by solving simultaneously the Poisson, ion transport, and Navier-Stokes equations by an iterative boundary element procedure. Treating the equilibrium electrostatics at the level of the full rather than the linear PB equation, but leaving relaxation out, does improve agreement between experimental and simulated mobilities, including ion relaxation improves it even more. The effects of nonlinear electrostatics and ion relaxation are greatest at low pH, where the net charge on lysozyme is greatest. In the absence of relaxation, a linear dependence of mobility and average polyion surface potential, (lambda zero)s, is observed, and the mobility is well described by the equation [formula: see text] where epsilon 0 is the dielectric constant of the solvent, and eta is the solvent viscosity. This breaks down, however, when ion relaxation is included and the mobility is less than predicted by the above equation. Whether or not ion relaxation is included, the mobility is found to be fairly insensitive to the charge distribution within the lysozyme model or the internal dielectric constant.

Simulation of electrostatic and hydrodynamic properties of Serratia endonuclease

We analyze the electrostatic and hydrodynamic properties of a nuclease from the pathogenic gram-negative bacterium Serratia marcescens using finite-difference Poisson-Boltzmann methods for electrostatic calculations and a bead-model approach for diffusion coefficient calculations. Electrostatic properties are analyzed for the enzyme in monomeric and dimeric forms and also in the context of DNA binding by the nuclease. Our preliminary results show that binding of a double-stranded DNA dodecamer by nuclease causes an overall shift in the charge of the protein by approximately three units of elementary charge per monomer, resulting in a positively charged protein at physiologic pH. In these calculations, the free enzyme was found to have a negative (-1 e) charge per monomer at pH 7. The most dramatic shift in pKa involves His 89 whose pKa increases by three pH units upon DNA binding. This shift leads to a protonated residue at pH 7, in contrast to the unprotonated form in the free enzyme. DNA binding also leads to a decrease in the energetic distances between the most stable protonation states of the enzyme. Dimerization has no significant effect on the electrostatic properties of each of the monomers for both free enzyme and that bound to DNA. Results of hydrodynamic calculations are consistent with the dimeric form of the enzyme in solution. The computed translational diffusion coefficient for the dimer model of the enzyme is in very good agreement with measurements from light scattering experiments. Preliminary electrooptical calculations indicate that the dimer should possess a large dipole moment (approximately 600 Debye units) as well as substantial optical anisotropy (limiting reduced linear electric dichroism of about 0.3). Therefore, this system may serve as a good model for investigation of electric and hydrodynamic properties by relaxation electrooptical experiments.

Electrostatic effects in homeodomain-DNA interactions

We report here an investigation of the role of electrostatics in homeodomain-DNA interactions using techniques based around the use of the Poisson-Boltzmann equation. In the present case such a study is of particular interest, since in contrast to other proteins previously studied with this method, the homeodomain is a small, highly charged protein that forms extensive ion pairs upon binding DNA. We have investigated the salt dependence of the binding constant for specific association and for a variety of models for non-specific association. The results indicate that, in line with the models proposed by Manning and Record, the entropy of counterion release accounts for a significant fraction of the salt dependence of the binding free energy, though this is perhaps due to fortuitous cancellation of other contributing terms. The thermodynamic effects of a number of specific homeodomain mutants were also investigated, and partly rationalized in terms of favorable electrostatic interactions in the major goove of DNA. Investigation of the temperature-dependence of the free energy of association indicates that the electrostatic contributions become increasingly favorable as the temperature rises. For this particular system, however, there appears to be no significant electrostatic contribution to the delta(delta C(p)) of association. Finally, an analysis of the free energy of interaction when the homeodomain is moved ca one Debye length from the DNA suggests that pure electrostatic forces are able to steer the homeodomain into a partially correct orientation for binding to the DNA.

The statistical-thermodynamic basis for computation of binding affinities: a critical review

Although the statistical thermodynamics of noncovalent binding has been considered in a number of theoretical papers, few methods of computing binding affinities are derived explicitly from this underlying theory. This has contributed to uncertainty and controversy in certain areas. This article therefore reviews and extends the connections of some important computational methods with the underlying statistical thermodynamics. A derivation of the standard free energy of binding forms the basis of this review. This derivation should be useful in formulating novel computational methods for predicting binding affinities. It also permits several important points to be established. For example, it is found that the double-annihilation method of computing binding energy does not yield the standard free energy of binding, but can be modified to yield this quantity. The derivation also makes it possible to define clearly the changes in translational, rotational, configurational, and solvent entropy upon binding. It is argued that molecular mass has a negligible effect upon the standard free energy of binding for biomolecular systems, and that the cratic entropy defined by Gurney is not a useful concept. In addition, the use of continuum models of the solvent in binding calculations is reviewed, and a formalism is presented for incorporating a limited number of solvent molecules explicitly.

Enzyme-inhibitor association thermodynamics: explicit and continuum solvent studies

Studying the thermodynamics of biochemical association reactions at the microscopic level requires efficient sampling of the configurations of the reactants and solvent as a function of the reaction pathways. In most cases, the associating ligand and receptor have complementary interlocking shapes. Upon association, loosely connected or disconnected solvent cavities at and around the binding site are formed. Disconnected solvent regions lead to severe statistical sampling problems when simulations are performed with explicit solvent. It was recently proposed that, when such limitations are encountered, they might be overcome by the use of the grand canonical ensemble. Here we investigate one such case and report the association free energy profile (potential of mean force) between trypsin and benzamidine along a chosen reaction coordinate as calculated using the grand canonical Monte Carlo method. The free energy profile is also calculated for a continuum solvent model using the Poisson equation, and the results are compared to the explicit water simulations. The comparison shows that the continuum solvent approach is surprisingly successful in reproducing the explicit solvent simulation results. The Monte Carlo results are analyzed in detail with respect to solvation structure. In the binding site channel there are waters bridging the carbonyl oxygen groups of Asp189 with the NH2 groups of benzamidine, which are displaced upon inhibitor binding. A similar solvent-bridging configuration has been seen in the crystal structure of trypsin complexed with bovine pancreatic trypsin inhibitor. The predicted locations of other internal waters are in very good agreement with the positions found in the crystal structures, which supports the accuracy of the simulations.

Prediction of titration properties of structures of a protein derived from molecular dynamics trajectories

This paper explores the dependence of the molecular dynamics (MD) trajectory of a protein molecule on the titration state assigned to the molecule. Four 100-ps MD trajectories of bovine pancreatic trypsin inhibitor (BPTI) were generated, starting from two different structures, each of which was held in two different charge states. The two starting structures were the X-ray crystal structure and one of the solution structures determined by NMR, and the charge states differed only in the ionization state of N terminus. Although it is evident that the MD simulations were too short to sample fully the equilibrium distribution of structures in each case, standard Poisson-Boltzmann titration state analysis of the resulting configurations shows general agreement between the overall titration behavior of the protein and the charge state assumed during MD simulation: at pH 7, the total net charge of the protein resulting from the titration analysis is consistently lower for the protein with the N terminus assumed to be neutral than for the protein with the N terminus assumed to be charged. For most of the ionizable residues, the differences in the calculated pKaS among the four trajectories are statistically negligible and remain in good agreement with the data obtained by crystal structure titration and by experiment. The exceptions include the N terminus, which responds directly to the change of its imposed charge; the C terminus, which in the NMR structure interacts strongly with the former; and a few other residues (Arg 1, Glu 7, Tyr 35, and Arg 42) whose pKaS reflect the initial structure and the limited trajectory lengths. This study illustrates the importance of the careful assignment of protonation states at the start of MD simulations and points to the need for simulation methods that allow for the variation of the protonation state in the calculation of equilibrium properties.

Structure-based drug design: computational advances

Structure-based computational methods continue to enhance progress in the discovery and refinement of therapeutic agents. Several such methods and their applications are described. These include molecular visualization and molecular modeling, docking, fragment methods, 3-D database techniques, and free-energy perturbation. Related issues that are discussed include the use of simplified potential energy functions and the determination of the positions of tightly bound waters. Strengths and weaknesses of the various methods are described.

Evidence for electrostatic channeling in a fusion protein of malate dehydrogenase and citrate synthase

Brownian dynamics simulations were performed to investigate a possible role for electrostatic channeling in transferring substrate between two of the enzymes of the citric acid cycle. The diffusion of oxaloacetate from one of the active sites of malate dehydrogenase (MDH) to the active sites of citrate synthase (CS) was simulated in the presence and absence of electrostatic forces using a modeled structure for a MDH-CS fusion protein. In the absence of electrostatic forces, fewer than 1% of substrate molecules leaving the MDH active site are transferred to CS. When electrostatic forces are present at zero ionic strength however, around 45% of substrate molecules are successfully channeled. As expected for an electrostatic mechanism of transfer, increasing the ionic strength in the simulations reduces the calculated transfer efficiency. Even at 150 mM however, the inclusion of electrostatic forces results in an increase in transfer efficiency of more than 1 order of magnitude. The simulations therefore provide evidence for the involvement of electrostatic channeling in guiding substrate transfer between two of the enzymes of the citric acid cycle. Similar effects may operate between other members of the citric acid metabolon.

Electrostatic channeling in the bifunctional enzyme dihydrofolate reductase-thymidylate synthase

The bifunctional enzyme dihydrofolate reductase-thymidylate synthase (DHFR-TS) carries out two distinct reactions, with the dihydrofolate produced by the TS-catalyzed reaction acting as the substrate for the DHFR-catalyzed reaction. Brownian dynamics simulation techniques were used to investigate the possible role of electrostatics in determining efficient channeling of the substrate, by explicitly simulating substrate diffusion between the two active sites. With a substrate charge of -2, almost all (> 95%) substrate molecules leaving the TS active site reached the DHFR active site at zero ionic strength. Under the same conditions, but in the absence of electrostatic effects, successful channeling was reduced to only around 6%: electrostatic effects therefore appear essential to explain the efficient channeling observed experimentally. The importance of substrate charge, the relative contributions of specific basic residues in the protein, the role played by the second monomer of the dimer in channeling and the effects of changing ionic strength were all investigated. Simulations performed for substrate transfer in the opposite direction suggest that channeling in DHFR-TS is not strongly directional and that the role of electrostatics is perhaps more one of restricting diffusion of the substrate than one of actively guiding it from the TS to the DHFR active site. The results demonstrate that electrostatic channeling can be a highly efficient means of transferring charged substrates between active sites in solvent-exposed environments.

Acetylcholinesterase: role of the enzyme's charge distribution in steering charged ligands toward the active site

The electrostatic steering of charged ligands toward the active site of Torpedo californica acetylcholinesterase is investigated by Brownian dynamics simulations of wild type enzyme and several mutated forms, in which some normally charged residues are neutralized. The simulations reveal that the total ligand influx through a surface of 42 A radius centered in the enzyme monomer and separated from the protein surface by 1-14 A is not significantly influenced by electrostatic interactions. Electrostatic effects are visible for encounters with a surface of 32 A radius, which is partially hidden inside the protein, but mostly within the solvent. A clear accumulation of encounter events for that sphere is observed in the area directly above the entrance to the active site gorge. In this area, the encounter events are increased by 40% compared to the case of a neutral ligand. However, the differences among the encounter rates for the various mutants considered here are not pronounced, all rate constants being within +/- 10% of the average value. The enzyme charge distribution becomes more important as the charged ligand moves toward the bottom of the gorge, where the active site is located. We show that neither the enzyme's total charge, nor its dipole moment, fully account for the electrostatic steering of ligand to the active site. Higher moments of the enzyme's charge distribution are also important. However, for a series of mutations for which the direction of the enzyme dipole moment is constant within a few degrees, one observes a gradual decrease in the diffusional encounter rate constant with the number of neutralized residues. On the other hand, for other mutants that change the direction of the dipole moment from that of the wild type, the calculated encounter rate constants can be very close to that of the wild type. The present work yields two new insights to the kinetics of acetylcholinesterase. First, evolution appears to have built a redundant electrostatic steering capability into this important enzyme through the overall distribution of its thousands of partially charged atoms. And second, roughly half of the rate enhancement due to electrostatics arises from steering of the substrate outside the enzyme; the other half of the rate enhancement arises from improved trapping of the substrate after it has entered the gorge. The computational results reproduce qualitatively, and help to rationalize, many surprising experimental results obtained recently for human acetylcholinesterase.

The determinants of pKas in proteins

Although validation studies show that theoretical models for predicting the pKas of ionizable groups in proteins are increasingly accurate, a number of important questions remain: (1) What factors limit the accuracy of current models? (2) How can conformational flexibility of proteins best be accounted for? (3) Will use of solution structures in the calculations, rather than crystal structures, improve the accuracy of the computed pKas? and (4) Why does accurate prediction of protein pKas seem to require that a high dielectric constant be assigned to the protein interior? This paper addresses these and related issues. Among the conclusions are the following: (1) computed pKas averaged over NMR structure sets are more accurate than those based upon single crystal structures; (2) use of atomic parameters optimized to reproduce hydration energies of small molecules improves agreement with experiment when a low protein dielectric constant is assumed; (3) despite use of NMR structures and optimized atomic parameters, pKas computed with a protein dielectric constant of 20 are more accurate than those computed with a low protein dielectric constant; (4) the pKa shifts in ribonuclease A that result from phosphate binding are reproduced reasonably well by calculations; (5) the substantial pKa shifts observed in turkey ovomucoid third domain result largely from interactions among ionized groups; and (6) both experimental data and calculations indicate that proteins tend to lower the pKas of Asp side chains but have little overall effect upon the pKas of other ionizable groups.

Study of global motions in proteins by weighted masses molecular dynamics: adenylate kinase as a test case

The weighted masses molecular dynamics (WMMD) technique is applied to the protein adenylate kinase. A novel set of restraints has been developed to allow the use of this technique with proteins. The WMMD simulation is successful in predicting the flexibility of the two mobile domains of the protein. The end product of the simulation is similar to the known open and AMP bound forms of the enzyme. The biological relevance of the restraints used and potential methods of improving the technique are discussed.

Orientational steering in enzyme-substrate association: ionic strength dependence of hydrodynamic torque effects

The effect of hydrodynamic torques on the association rate constants for enzyme-ligand complexation is investigated by Brownian dynamics simulations. Our hydrodynamic models of the enzyme and ligand are composed of spherical elements with friction forces acting at their centers. A quantitative measure of hydrodynamic torque orientational effects is introduced by choosing, as a reference system, an enzyme-ligand model with the same average hydrodynamic interactions but without orientational dependence. Our simple models show a 15% increase in the rate constant caused by hydrodynamic torques at physiological ionic strength. For more realistic hydrodynamic models, which are not computationally feasible at present, this effect is probably higher. The most important finding of this work is that hydrodynamic complementarity in shape (i.e. like the fitting together of pieces of a puzzle) is most effective for interactions between molecules at physiological ionic strength.