Journal of Molecular Biology
Volume 291, Issue 1, 6 September 1999, Pages 149-162
Journal home page for Journal of Molecular Biology

Regular article
Computer simulation of protein-protein association kinetics: acetylcholinesterase-fasciculin1

https://doi.org/10.1006/jmbi.1999.2919Get rights and content

Abstract

Computer simulations were performed to investigate the role of electrostatic interactions in promoting fast association of acetylcholinesterase with its peptidic inhibitor, the neurotoxin fasciculin. The encounter of the two macromolecules was simulated with the technique of Brownian dynamics (BD), using atomically detailed structures, and association rate constants were calculated for the wild-type and a number of mutant proteins. In a first set of simulations, the ordering of the experimental rate constants for the mutant proteins was correctly reproduced, although the absolute values of the rate constants were overestimated by a factor of around 30. Rigorous calculations of the full electrostatic interaction energy between the two proteins indicate that this overestimation of association rates results at least in part from approximations made in the description of interaction energetics in the BD simulations. In particular, the initial BD simulations neglect the unfavourable electrostatic desolvation effects that result from the exclusion of high dielectric solvent that accompanies the approach of the two low dielectric proteins. This electrostatic desolvation component is so large that the overall contribution of electrostatics to the binding energy of the complex is unlikely to be strongly favourable. Nevertheless, electrostatic interactions are still responsible for increased association rates, because even if they are unfavourable in the fully formed complex, they are still favourable at intermediate protein-protein separation distances. It therefore appears possible for electrostatic interactions to promote the kinetics of binding even if they do not make a strongly favourable contribution to the thermodynamics of binding. When an approximate description of these electrostatic desolvation effects is included in a second set of BD simulations, the relative ordering of the mutant proteins is again correctly reproduced, but now association rate constants that are much closer in magnitude to the experimental values are obtained. Inclusion of electrostatic desolvation effects also improves reproduction of the experimental ionic strength dependence of the wild-type association rate.

Introduction

Macromolecular association is a phenomenon of fundamental importance in biology, and much experimental and theoretical work has been aimed at understanding the factors governing both its thermodynamic and kinetic aspects. In terms of kinetics, it has been shown that the bimolecular association rate constants for the formation of dimeric complexes are typically in the range 0.5 × 106to 5.0 × 106M−1s−1 (Northrup & Erickson, 1992), although much higher values of around 109M−1s−1have been reported for several systems Eltis et al 1991, Schreiber and Fersht 1996, Stone et al 1986, Radic et al 1997. The latter values are at, or close to, the diffusion limit for uniformly reactive molecules. Systems exhibiting such high association rates are typically subject to strong electrostatic steering forces, an aspect generally indicated by a strong dependence of the association rate constant on the ionic strength of the solution. An understanding of the role of electrostatic interactions (and other factors) in promoting efficient protein association is therefore of interest from both academic and protein-design viewpoints. In the present work we study protein association kinetics using computer simulations: we describe the use of Brownian dynamics (BD) methods to simulate the association of acetylcholinesterase and fasciculin, two proteins whose association is known to be both diffusion-limited and dramatically accelerated by electrostatic interactions (Radic et al., 1997).

The use of BD simulation techniques for calculating association or encounter rate constants of enzymes with their (small molecule) substrates has been well documented Wade 1996, Madura et al 1994. The method involves explicit simulation of the diffusion of substrate from bulk solution to the enzyme’s active site. Throughout the course of the simulation, the substrate is subjected to random forces that model collisions with solvent molecules, and electrostatic forces that model the interaction of the substrate’s charge with the electrostatic field generated by the distribution of charges in the enzyme. Typically, thousands of simulations are run in order to obtain statistically robust estimates of the probability that substrate will successfully bind; once this probability is obtained, it can be used to calculate the bimolecular rate constant for binding. The BD method has been successfully applied to the diffusion-limited enzymes superoxide dismutase (Getzoff et al., 1992), triose phosphate isomerase Wade et al 1994, Wade et al 1998, and acetylcholinesterase Radic et al 1997, Tara et al 1998, allowing the sometimes surprising effects of ionic strength and residue mutations on reaction rate constants to be rationalised. More recently, the same technique has been used to simulate the diffusion of substrates between the active sites of enzymes constituting multi-enzyme complexes; such simulations have provided strong evidence that channeling of substrates is often mediated by electrostatic interactions Elcock and McCammon 1996b, Elcock et al 1996, Elcock et al 1997.

The reported applications of the BD method in simulating the association of two proteins (rather than a protein and a small molecule) are fewer (for a review, see Gabdoulline & Wade, 1998), a point which is perhaps not surprising given the significant step up in complexity involved Castro et al 1998, Gabdoulline and Wade 1997, Kozack et al 1995, Nambi et al 1991, Northrup et al 1993. A demonstration of the potential power of BD simulation methods in this area was given by the study reported by Northrup & Erickson (1992)which showed why the kinetics of even non-electrostatically controlled protein-protein association is surprisingly rapid. Using simplified spherical models of proteins, it was demonstrated that fast association results because when two proteins encounter each other, they typically remain close together for a time period similar to the rotational correlation times of the individual proteins. As a result, during the course of any given encounter, substantial rotational motion of both proteins occurs, thus increasing the likelihood that the two proteins will adopt a suitable orientation for proper binding. Northrup and co-workers have subsequently moved beyond the use of these simple models, and have used atomically detailed structures in a study of electron transfer rates in cytochrome systems, obtaining good qualitative agreement with the experimental effects of a number of residue mutations Northrup et al 1993, Castro et al 1998.

Recently, Gabdoulline & Wade (1997)have employed a similar approach in a study of the association kinetics of barnase and its inhibitor barstar. Experimental work on this system had shown that the association rate constant exhibits an ionic strength dependence strongly indicative of electrostatic acceleration (the rate constant extrapolates to ∼1 × 1010M−1s−1at zero ionic strength), and a viscosity dependence indicative of diffusion control (Schreiber & Fersht, 1996). Additionally, single-residue (and one double-residue) mutations had been investigated, and rate constants covering the rather wide range of 3 × 107M−1s−1to 2 × 109M−1s−1in 50 mM NaCl were obtained. This wealth of experimental information therefore provided a good test system for assessing the usefulness of BD methods for simulating (and ultimately predicting) protein-protein association kinetics. On the whole, the results obtained by Gabdoulline & Wade (1997)were excellent: the relative rate constants of the various mutants were correctly reproduced, as was the ionic strength dependence of the wild-type proteins’ association rate constant. Two aspects of the simulations were deemed to be particularly important for the good agreements obtained and so are of special relevance for the present work. Firstly, the authors used “effective charges” to model the charge distribution in each protein. Effective charges for each protein are derived so that when immersed in a uniform dielectric medium they accurately reproduce the electrostatic potential obtained when the same protein is treated (more rigorously) as a low-dielectric macromolecule immersed in a high dielectric solvent. This approach contrasts with the more usual test-charge approximation used in other BD simulations of protein-protein encounters, and is considerably more accurate, especially so at intermediate to long-range separation distances (Gabdoulline & Wade, 1996). Secondly, a number of different methods for defining successful formation of the final barnase-barstar complex were investigated: a suitable choice of complexation or “reaction” criteria is essential if the simulations are to adequately reproduce experimental association rates. The best agreement with experiment was obtained when the complex was considered to be successfully formed only when any two of the intermolecular hydrogen bond contacts present in the crystal structure were satisfied to within a distance of 6.25 Å. Alternative definitions, such as requiring only one or three contacts to be formed, or using RMS structural deviations from the crystal structure as a guide, were tried, but all were less successful in reproducing the experimental data on the effects of protein mutation than the two-contact criterion.

The success of the BD method in the barnase-barstar system is encouraging, but it remains to be seen whether the same approach can be applied to other protein systems with equal success. In particular, it is important to demonstrate whether the two-contact reaction criterion is likely to be of more general use, since a transferable (system-independent) criterion is an obvious requirement if the method is to be used to predict association rate constants for systems for which no experimental kinetic data are available. We therefore investigate the application of identical simulation methodology to the association of the enzyme acetylcholinesterase (AChE) with fasciculin (Fas), since this is another system for which excellent experimental data exist (Radic et al., 1997). AChE is well known for its role in terminating nerve impulses by hydrolysing the neurotransmitter acetylcholine. Given its physiological role, it is probably not surprising to find that kinetics studies indicate that AChE is a diffusion-limited enzyme (Quinn, 1987), i.e. its catalytic machinery has been refined to such an extent that the binding of the substrate to the active site is the rate-limiting step in the catalytic cycle. The distribution of charged residues on the surface of the enzyme has been shown to be critically important for accelerating this encounter step: acidic residues are strategically positioned around the entrance to the 20 Å-long gorge that leads to the active site, in such a way as to attract and trap the positively charged substrate acetylcholine (Radic et al., 1997). BD simulations have already been successful in rationalising the effects of mutations on the association rate constants for the enzyme with small molecule substrates (see Tara et al., 1998, and references therein).

Fasciculins are 61-residue peptides isolated from snake venom, which bind and inhibit AChE with very high affinity. The crystal structure of the AChE-Fas complex (Bourne et al., 1995) reveals that inhibition of the enzyme results because the peptide binds at the mouth of the gorge leading to the active site, and largely (although not completely) occludes the entrance, thus making it difficult for substrate to reach the active site. Of more relevance to the present work, however, is the fact that kinetics studies have shown that the binding of Fas to AChE is both diffusion-limited and strongly modulated by electrostatic interactions (Radic et al., 1997). The binding kinetics of a number of mutants have been characterised, including one in which a total of six charged residues thought to be important for promoting association were substituted by neutral residues, and a range of association rate constants from ∼3 × 106M−1s−1to ∼8 × 108M−1s−1(when extrapolated to zero ionic strength) was obtained.

We describe here the use of a BD method (Gabdoulline & Wade, 1998) to reproduce the effects of these mutations, as well as the effects of changing ionic strength on the association rate constants. As will be shown, the method is again successful in reproducing the relative rate constants of wild-type and mutant proteins, although the absolute values of the association rate constants are overestimated. The method is less successful in reproducing the ionic strength dependence of the wild-type association rate constants. As we show, both the overestimation of absolute rate constants, and the poor reproduction of ionic strength effects appear to be due to approximations made in modelling the interaction energetics of the two proteins. By comparing with the results of more complete and rigorous electrostatic calculations, it is shown that at close separation distances the energetic cost of desolvating the proteins, which is neglected in the first set of BD simulations, becomes increasingly important and therefore ought not to be omitted. When an approximate but computationally inexpensive approach to estimating these desolvation effects is included in a second set of BD simulations the ionic strength dependence is correctly reproduced, and lower absolute values for the association rate constants are obtained.

Section snippets

Comparison of mutant rate constants using the effective charge approximation

Figure 1 shows the ratio of the simulated rate constants to the experimental rate constants at 50 mM ionic strength, for the wild-type AChE and seven mutants, obtained using a two-contact association criterion and plotted as functions of the required contact distance. The overall shapes of the plots are similar to those seen in the previous study of barnase and barstar (Gabdoulline & Wade, 1997): the calculated rate constants decrease steadily as the criteria used to determine successful

Discussion

According to the complete electrostatic binding energy calculations shown in Figure 4, the overall contribution of electrostatics to the binding of AChE and Fas is unfavourable, because although favourable coulombic interactions are formed in the final complex, they are not sufficient to offset the energetic cost of desolvating the two proteins. The conclusion that the overall effect of electrostatics is to oppose binding is open to dispute, since the magnitude of the calculated results will

Conclusion

The present work has shown that as with previous studies, it is possible for computer simulation methods to reproduce the experimental effects of residue mutations on the association rate constants of two proteins even when, as in the present case, the mutant proteins contain as many as six altered residues Northrup et al 1993, Gabdoulline and Wade 1997. Particularly important is the finding that the two-contact criterion used successfully in the barnase-barstar system also appears to work well

Structures

Structures for AChE and Fas were taken from the crystal structure of the mouse AChE-Fas complex solved to 3.2 Å resolution by Bourne et al. (1995)(PDB code 1mah; see Figure 9). AChE residues 1–3 and 258–264 are not visible in the crystal and were added manually using QUANTA (Molecular Simulations Inc., San Diego, CA) and energy minimized using the molecular simulation program CHARMM (Brooks et al., 1983); since they are some distance (at least 25 Å) from the interface with Fas, completely

Acknowledgements

A.H.E. thanks Dr David A. Sept, Dr Zoran Radic and Professor Palmer A. Taylor for valuable discussions. This work was supported in part by grants (to J.A.M) from the National Institutes of Health, the National Science Foundation and the National Science Foundation MetaCenter Program.

References (39)

  • A.-S. Yang et al.

    Free energy determinants of secondary structure formation. 2. Anti-parallel beta-sheets

    J. Mol. Biol.

    (1995)
  • B.R. Brooks et al.

    CHARMMa program for macromolecular energy, minimization and dynamics calculations

    J. Comp. Chem

    (1983)
  • G. Castro et al.

    Dynamics of protein-protein dockingcytochrome c and cytochrome c peroxidase revisited

    J. Biomol. Struct. Dynam.

    (1998)
  • T.E. Creighton

    Proteins

    (1993)
  • A.H. Elcock et al.

    The low dielectric interior of proteins is sufficient to cause major structural changes in DNA on association

    J. Am. Chem. Soc.

    (1996)
  • A.H. Elcock et al.

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

    Biochemistry

    (1996)
  • A.H. Elcock et al.

    Electrostatic channeling of substrates between enzyme active sitescomparison of simulation and experiment

    Biochemistry

    (1997)
  • L.D. Eltis et al.

    Reduction of horse heart ferricytochrome-c by bovine liver ferrocytochrome-b5 - experimental and theoretical analysis

    Biochemistry

    (1991)
  • D.L. Ermak et al.

    Brownian dynamics with hydrodynamic interactions

    J. Chem. Phys.

    (1978)
  • Cited by (169)

    • Computation of FRAP recovery times for linker histone – chromatin binding on the basis of Brownian dynamics simulations

      2020, Biochimica et Biophysica Acta - General Subjects
      Citation Excerpt :

      Thus, not only the specific binding interactions, but also the crowding effects inside the nucleus should be accounted for in comparing dilute solution in-vitro and in-silico studies of molecular interactions with in-vivo or in-cellulo studies. Brownian dynamics (BD) simulation is a valuable tool for studying the diffusional binding kinetics of biomolecules and has been previously used to predict kon values and the structures of diffusional encounter complexes of a wide range of macromolecular complexes [12–19]. Here, by performing BD simulations of gH1.0 - nucleosome diffusional association, kon and koff values were calculated for diffusional encounter complex formation by wild-type (WT) and mutant gH1.0.

    • Brownian dynamics simulations of biomolecular diffusional association processes

      2023, Wiley Interdisciplinary Reviews: Computational Molecular Science
    View all citing articles on Scopus
    1

    Edited by B. Honig

    View full text