Protein Engineering, Vol. 15, No. 9, 739-743,
September 2002
© 2002 Oxford University Press
Free-energy contributions to complex formation between botulinum neurotoxin type B and synaptobrevin fragment
Department of Cell Biology and Biochemistry, USAMRIID,1425 Porter Street, Frederick, MD 21702, USA
| Abstract |
|---|
|
|
|---|
Free-energy terms that contribute to complex formation between the catalytic domain of botulinum neurotoxin type B (BoNT/B-LC) and a 36-residue synaptobrevin fragment were estimated by using a combination of microscopic simulations and continuum methods. The complex for a non-hydrolyzed substrate was calculated by optimizing an energy function applied to the X-ray co-crystal structure of BoNT/B-LC bound with reaction products from a cleaved synaptobrevin peptide, refined to high crystallographic thermal factors. The estimated absolute binding affinity of the simulation structure is in good qualitative agreement with the experimental free energy of Michaelis complex formation, given the approximations of the model calculations. The simulation structure revealed significant complex stabilization from the hydrophobic effect, while the electrostatic cost of releasing water molecules from the interface determined to be highly unfavorable. By partitioning the total electrostatic and hydrophobic terms into residue free-energy contributions, a binding-affinity signature for synaptobrevin was developed from the optimized conformation. The results demonstrate the effect of substrate length on complex formation and identify a peripheral high-affinity binding site near the N-terminal region that might initiate cooperative activation responsible for the large minimal substrate length requirement. The so-called SNARE motif is observed to contribute negligible free energy of binding.
Keywords: continuum models/electrostatics/hydration/hydrophobic interaction
| Introduction |
|---|
|
|
|---|
Botulinum neurotoxins are zinc proteases that function by cleaving SNARE proteins involved in neuronal synaptic vesicle fusion. Produced by Clostridium botulinum, this family of toxins is comprised of seven immunologically distinct types, referred to as BoNT/AG, and is related to tetanus toxin (TeNT), produced by Clostridium tetani. Recently, X-ray crystallographic structures have been reported (Hanson and Stevens, 2000
35 residues (Shone et al., 1993Here, we provide a model for complex formation between a 36-residue non-hydrolyzed SB2 fragment bound to the BoNT/B-LC. Our model is developed from the co-crystal structure by generating an ensemble of alternative conformations extracted from a molecular dynamics (MD) simulated-annealing simulation. Free-energy terms that contribute to binding affinity are calculated by applying continuum methods for modeling interactions of the complex and solvent. Numerical solutions to the PoissonBoltzmann (PB) equation are used to estimate electrostatic terms and a molecular-surface area relationship is applied for the hydrophobic term. We compare the modeled complex with the co-crystal structure and discuss which residues of SB2 contribute significantly to binding affinity and their possible role in substrate specificity by means of a regulatory site.
| Computational methods |
|---|
|
|
|---|
The molecular association between the BoNT/B-LC (denoted P) and the synaptobrevin fragment (
) can be described as a two-step process (see, e.g., Olson, 2001
![]() | (1) |
![]() | (2) |
conformers to the functional conformers P* and
* and the second step describes the binding of the two molecules in their functional state. The Gibbs free energy of binding is obtained by combining Equations 1 and 2
![]() | (3) |
Gconf is the free-energy shift upon conformational changes and
Gstatic is the non-relaxation free energy. The term
Gstatic is partitioned into independent energetic contributions, where
W represents the change in the potential of mean force for P * and L* due to solutesolvent interactions,
Gint is the interaction free energy between the bound molecules and T
S is the change in non-electrostatic entropy determined at temperature T. The entropic change is defined as a combination of rotational, translational and main-chain flexibility contributions:
![]() | (4) |
Calculation of
W given in Equation 3
is obtained from the expression
![]() | (5) |
Gcav is the change in free energy required to form the solute-sized cavities in the solvent,
Gs,vdW is the free-energy change in van der Waals (vdW) interactions of inserting solute molecules into the cavities and
Gs,ele is the change in the free energy of solvent polarization. The cavitation term models the hydrophobic effect and is proportional to the change in solvent-exposed surface area upon association with a constant surface tension,
:
![]() | (6) |
is set at 69 cal/mol.Å2 (Jackson and Sternberg, 1994
The contribution
Gint contains polar and non-polar intermolecular interaction terms between P* and
*:
![]() | (7) |
Gm,vdW is the vdW contribution and
Gm,ele the electrostatic interaction component. The term
Gm,vdW combined with
Gs,vdW of Equation 5
The electrostatic free energies given in Equations 5 and 6![]()
can be conveniently calculated from a thermodynamic process of charging and uncharging P* and L * on complex formation (Gilson and Honig, 1988
; Muegge et al., 1997
):
![]() | (8) |
and
are the spatially dependent charge density and electrical potential, respectively, in volume element dr. The first term represents the loss of solutesolvent interaction energy through the partial desolvation of the electrostatically charged P* on binding the uncharged
* (denoted
*'), and similarly, the second term corresponds to the desolvation of
* on binding P*'; the third term is the electrostatic interaction energy between P* and
* in the solvent dielectric.
Coordinates for the complex between the SB2 products and BoNT/B-LC were taken from the PDB file 1F83. Placement of the hydrogen atoms was modeled by InsightII (Molecular Simulations, v97.0), with ionized residues selected for neutral pH. Modeling of the non-hydrolyzed substrate was built from the co-crystal structure by reconnecting the peptide linkage at the scissile bond (Q76F77) and optimizing Equation 3
calculated for two models: (1) local residues 7283 of the SB2 peptide or (2) the complete SB2 peptide and surrounding BoNT/B-LC residues within an 8 Å region (163 residues) of the interface. Optimization was performed by a MD simulated-annealing (MD-SA) approach. All MD simulations treated Coulombic interactions by a cell multipole method (Ding et al., 1992
) with a constant dielectric of
= 1 and the vdW interactions by a group-based approach with a cutoff radius of 12.0 Å. The force field and atomic charges were taken from AMBER (Weiner et al., 1986
). Constraints were applied via the RATTLE algorithm (Andersen, 1983
) to bond lengths of both molecules plus the solvent during dynamics runs. The integration time step was set at 2.0 fs.
MD-SA calculations were initiated with 200 cycles of energy minimization using a conjugate gradient method, followed by a MD simulation for 500 steps at a temperature of 1000 K. This was then followed by an annealing protocol for a temperature range of 2501000 K, employing an interval of 25 K from the initial high temperature. Each simulation was completed for 1000 steps and started from the previous temperature by using the final calculated structure. Active-site regions of both models were immersed in a 9 Å layer of explicit solvent water, modeled by the TIP3P potential (Jorgensen et al., 1983
). Atomic parameters for the interfacial zinc ion were estimated at a charge of +2 and radius of 1.4 Å. During all MD-SA calculations, the zinc moiety and one of crystallographic waters near the ion were held fixed at their crystallographic position. To keep the solvent layer from escaping at high temperatures, the outermost shell of the layer was harmonically tethered at the initial modeled position with a force constant of 50 kcal/Å. The final ensemble of 300 structures from five independent MD-SA simulations was subjected to 1000 cycles of minimization.
The lowest energy conformation from annealing was further refined by MD simulation at a temperature of 298 K. The annealing simulation waters were replaced with a new set of waters and were minimized for 1000 cycles. This was followed by 500 steps of MD equilibration; additional energy minimization of 1000 cycles, second equilibration for 500 MD steps and production MD run for 50 ps. Conformations from the final MD trajectory were saved every 500 steps.
To evaluate Equation 8
, we applied a finite-difference PB method implemented in the program DelPhi (Gilson and Honig, 1988
) for each conformation (
500) extracted from the MD run. The electrostatic potentials for each molecule were calculated by using the solvent-accessible surfaces to define regions of low dielectric medium embedded in a high dielectric solvent water of ionic strength set at 0.145 M. DelPhi calculations were carried out on a cubic grid of 1693 grid points. The protein dielectric constant was set at either
p = 2 or 4 and a dielectric of 80 was used to model bulk water. The zinc ion and the nearby water were retained explicitly as part of the protein, while other waters were deleted. Boundary conditions were set at full Coulombic for all calculations. Dummy atoms were implemented for the unbound structures in retaining an identical scale and position on the grid used with the complex configuration.
Cavitation energies were determined by applying a molecular-surface algorithm (Connolly, 1981
) with the solvent probe radius set at 1.4 Å. For entropic contributions, estimates were applied for rotational and translational terms (Karplus and Janin, 1999
). The loss of main-chain conformational entropy was determined from an empirical scale that includes side-chain effects (Pal and Chakrabarti, 1999
).
The final free energy for complex formation was taken from the lowest energy conformer sampled by the MD simulation. Because of the conformational entropic term, the free energy for this conformer represents a macroscopic state of the native basin and can be partitioned on a per-residue basis for modeling individual contributions to binding. Dipolar fluctuations are modeled by
p
2. An alternative method is the calculation of an average free energy over the simulation trajectory, which accounts for statistical fluctuations of the ensemble. This latter approach can be compared with the more rigorous linear response method of extracting
Gs,ele from explicit treatment of solvent in all-atom simulations (Muegge et al., 1998
; Olson, 2001
).
| Results and discussion |
|---|
|
|
|---|
Summarized in Table I
Gbind determined for the three modeled SB2BoNT/B-LC complexes: (1) the cleaved SB2 fragment evaluated directly from the co-crystal structure with no calculated optimization; (2) SB2 modeled as a non-hydrolyzed substrate by local optimization at the scissile bond; and (3) a global conformational search of the modeled non-hydrolyzed substrate and surrounding interfacial BoNT/B-LC region. For comparison purposes, the experimental binding free energy for the 36-mer SB2 was estimated from the expression
Gexpt = -RTlnKm, where Km is the reported Michaelis constant (Shone et al., 1993
|
Before we discuss the results of the calculations, a few comments are needed regarding the meaning of the dielectric constant and its application in continuous approaches. The parameter
p in PB models is a scaling factor that represents all of the electrostatic contributions that are not treated explicitly, rather than a true physical dielectric constant (King and Warshel, 1991
p set at 2 models at the implicit level the induced protein dipoles, while the dipolar reorganization from conformational changes upon binding is to be taken into account explicitly. With this approach, the determination of
Gconf for the structural transition P P* is required and must have a value less than the free energy of unfolding. Reorganization free energies can be evaluated from all-atom simulations (Muegge et al., 1998
Gconf
5 kcal/mol for modeling the small structural changes in BoNT/B-LC due to binding SB2 (Hanson and Stevens, 2000
Gconf by applying
p = 4, which embeds reorganization into the 
Gs,ele and
Gm,ele terms. This continuum treatment is more common for modeling the effects of amino acid substitutions on complex stabilization, rather than the absolute free energy of binding (see, e.g., Olson, 1998
* conformational change, the reorganization penalty should be insignificant. This is justified from the observation that SB2 binds in a random coil conformation and contains little periodic structure in its free solution state (Hazzard et al., 1999
The first model listed in Table I
calculated with
p = 2 clearly shows a
Gbind unfavorable for complex formation. A generalization observed of results from continuum model calculations is that the hydrophobic effect provides the driving force underlying favorable complex formation (namely,
Gcav < 0), while the net electrostatic effect opposes binding (
Gele > 0), yet confers conformational specificity. Here, for the crystal structure of SB2 products, the total electrostatic contribution is significantly unfavorable at
246 kcal/mol and entirely offsets the favorable hydrophobic term calculated at
-186 kcal/mol. The end result is that the calculated
Gcav +
Gele > 0 and thus the complex is unstable. As noted further below, this unfavorable balance arises from the high positional uncertainty in the SB2 products located in the binding groove of BoNT/B-LC.
Keeping with the dielectric treatment of
p = 2, the model for the non-hydrolyzed substrate obtained by reconstructing the cleaved peptide bond and optimizing by a local conformational search fails to improve upon the crystal structure. In contrast, global optimization of the same model yields the correct relationship between specific and non-specific binding terms. The highest scoring conformation found by the global search lowers the desolvation penalty of both BoNT/B-LC and SB2 from the two alternative models. Moreover, the optimized complex yields a significant gain in stabilization from electrostatic intermolecular interactions determined by
Gint, however, at a cost of diminishing the hydrophobic contribution estimated from
Gcav. The final calculated hydrophobic term is, nevertheless, extraordinarily favorable at -173 kcal/mol and indicates a large interfacial surface area from binding SB2. Located at the interface, interactions for Zn2+ show a minimal contribution in the binding affinity of the substrate, which is consistent with a functional role.
The favorable sum of
Gcav and
Gele from optimizing the conformation is partially offset by the entropic loss in main-chain freedom of the SB2 fragment upon binding, calculated at -T
Smc
41 kcal/mol. Theoretical estimates of the rotational and translational loss for proteinprotein assemblies suggest values in the range 715 kcal/mol at room temperature (Karplus and Janin, 1999
; and references cited therein). We selected a value for -T
Srot,trans of
10 kcal/mol (see, e.g., Froloff et al., 1997
). Given these approximations,
Gbind is estimated at
-10 kcal/mol, while
Gexpt is
-5 kcal/mol (Shone and Roberts, 1994
). In comparison with other reported continuum model estimates of absolute binding affinities (Jackson and Sternberg, 1995; Froloff et al., 1997
; Olson, 1998
), the calculations outlined here are reasonable for the SB2BoNT/B-LC complex and suggest that the model captures the correct physics of binding.
For the dielectric model of
p = 4, the complex with either products or the substrate calculated with local conformational optimization produced a better result for
Gbind than the equivalent models obtained with a smaller dielectric constant. Further reconciliation of the calculated values with
Gexpt requires additional scaling of
p. A disadvantage of larger
p models is that chargecharge interactions and solvent polarization are significantly coarse-grained, yielding much weaker physical interpretation of the electrostatic effects of binding.
Figure 1
compares the simulation structure and the X-ray crystal structure of SB2 products. Plotted is the coordinate root-mean-square deviation (r.m.s.d.) between the two structures as a function of SB2 residues and their estimated crystallographic B-factors. Several observations can be made from the results. First, no clear correlation exists between the r.m.s.d. and the B-factors reported by Hanson and Stevens (Hanson and Stevens, 2000
). The simulation structure finds significant r.m.s.d. near residues Q58 and D65, yet their B values are not drastically different from those of the other residues. The second observation is that the simulation model and its conformational deviation from the products can easily fit within the positional envelope outlined by the high B-factors. In other words, rather than occupancy of cleaved products in the crystal structure (Hanson and Stevens, 2000
), a non-hydrolyzed substrate is equally plausible. The r.m.s.d. for the cleaved peptide bond between Q76 and F77 is no larger than the average value for SB2 and includes a different side-chain rotamer for phenylalanine to that observed in the crystal structure. Other stereochemical differences are similarly found between the two structures.
|
Having established a good qualitative free-energy model for complex formation, we now ask the question of which residues of SB2 contribute the most to binding BoNT/B-LC. Using the optimized complex with
p = 2, the absolute free energy of binding is decomposed on a per-residue level into individual energetic contributions. The modeling scheme used to derive this binding-affinity signature is that previously developed and applied (Olson and Cuff, 1999
Gres, which includes the net sum of desolvation, electrostatic intermolecular interaction and the hydrophobic contribution. Also reported is the accumulation of
Gres as a function of residue count, where S1 designates residues 5376 located in the N-terminal region of the scissile bond and S2 residues 7788 in the COOH region.
|
The calculations show that Q76 interacts as a critical residue in stabilizing the orientation of the scissile bond, contributing significant electrostatic free energy of binding. Determination of Q76 as a hot residue in complex formation is consistent with binding studies of substrate analogues (Shone and Roberts, 1994
The accumulation of
Gres illustrated in Figure 2c
establishes the connection between substrate length and free energy of complex formation. Residues downstream from Q76 along S1 yield only a moderate net gain in binding affinity, until reaching D64. Other than perhaps the latter residue and L63, the so-called SNARE motif of sequence 6371 (Montecucco and Schiavo, 1995
) contributes an overall weak binding signature and is not likely to lead significantly to SB2 specificity. In contrast, the sharp gradient in
Gres for the region 5864 reveals considerable free-energy contributions in achieving the Michaelis complex. Calculated distal to the scissile bond, this high-affinity binding region provides substrate length specificity by a combination of electrostatic (Q58, S61 and D64) and hydrophobic (L60 and L63) interaction terms. The region contacts two surface loops on BoNT/B-LC, loops 50 (residues 5179) and loop 200 (residues 200219). Upon complex formation, the loops undergo a disorderedordered transition observed in the B-factors (Hanson and Stevens, 2000
). A possible implication of the transition is allosteric activation of BoNT/B-LC. Buttress for this idea is a set of elegant experiments on TeNT byCornille et al.(1997)
, who reported a cooperative mechanism based on kinetic data. What is presented here is the identification of a peripheral binding site on SB2 combined with its free-energy determinant of loop stabilization in BoNT/B-LC.
For the S2 region, overall contributions to binding affinity are less than S1, yet two residues, K83 and Y88, provide significant free energies for complex stabilization. Among the ionic residues of SB2, the charge interactions of K83 are the most energetically favorable. This residue is mostly variable among the substrates (exceptions are the same proteolytic site of synaptobrevin for TeNT and syntaxin for BoNT/C) and its interacting residues on BoNT/B-LC, Tyr53 and Asp170, are variable among the toxins. The side chain of SB2 Y88 provides the strongest cavitation term of any residue, while F77 is an important contributor to binding. Both residues are variable among substrates and the latter residue has been shown in substrate analogues to play a key role (Shone and Roberts, 1994
).
Combining S1 and S2 regions shows a free-energy minimum at
32 residues, which qualitatively reproduces the experimental observation of 35 residues as an SB2 length requirement (Shone and Roberts, 1994
). Each residue contribution is not required to be favorable, but rather binding takes place by a conformational search for the net minimum free energy of complex formation. The most noticeable unfavorable contribution is from R56, which shows a large desolvation penalty. The residue D68 displays electrostatic strain from weak complementarity with the electrical potential of BoNT/B-LC.
The free-energy function proposed in this paper and its optimization by means of an ensemble of conformers taken from a MD simulation trajectory show a promising method for suggesting alternative structures where uncertainty is unusually high in the crystallographic determination. Although the method cannot address the issue of poorly defined electron density maps of the substrate products (Rupp and Segelke, 2001
), the orientation of SB2 in the active site suggested by Hanson and Stevens (Hanson and Stevens, 2000
) appears to be a good starting point for computational refinement. Our model for SB2 complex formation yields new insight and should prove useful in gauging binding free energies of inhibitors found by either database searching or de novo design methods.
| Notes |
|---|
1 To whom correspondence should be addressed. Email: molson{at}ncifcrf.gov
2 Present address: Department of Biology, Virginia Polytechnic Institute and State University, Blacksburg, VA 24061, USA ![]()
| Acknowledgments |
|---|
We thank Drs C.Millard, F.Lebeda and J.Schmidt for many valuable discussions. This work was supported by a generous grant of computer time from the Advanced Biomedical Computing Center of the Frederick Cancer Research and Development Center, National Cancer Institute.
| References |
|---|
|
|
|---|
Andersen,H.C. (1983) J. Chem. Phys., 52, 2434.
Connolly,M.L. (1981) Molecular Surface Program 429. Quantum Chemistry Program Exchange, Indiana University, Bloomington, IN.
Cornille,F., Martin,L., Lenoir,C. Cussac,D., Roques,B.P. and Fournie-Zaluski,M-C. (1997) J. Biol. Chem., 272, 34593464.
Ding,H.Q., Karasawa,N. and Goddard,W.A.,III (1992) J. Chem. Phys., 97, 43094315.[CrossRef]
Froloff,N., Windemuth,A. and Honig,B. (1997) Protein Sci., 6, 12931301.[Web of Science][Medline]
Gilson,M.K. and Honig,B. (1988) Proteins, 4, 718.[CrossRef][Web of Science][Medline]
Hanson,M.A. and Stevens,R.C. (2000) Nature Struct. Biol., 7, 687692.[CrossRef][Web of Science][Medline]
Hazzard,J., Sudhof,T.C. and Rizo,J. (1999) J. Biomol. NMR, 14, 203207.[CrossRef][Web of Science][Medline]
Jackson,R.M. and Sternberg,M.J.E. (1994) Protein Eng., 7, 371383.
Jorgensen,W.L., Chandrasekhar,J., Madura,J.D., Impey,R.W. and Klein,M.L. (1983) J. Chem. Phys., 79, 926935.[CrossRef]
Karplus,M. and Janin,J. (1999) Protein Eng., 12, 185186.
King,G.F. and Warshel,A. (1991) J. Chem. Phys., 95, 43664377.[CrossRef]
Montecucco,C. and Schiavo,G. (1995) Q. Rev. Biophys., 28, 423472.[Web of Science][Medline]
Muegge,I., Tao,H. and Warshel,A. (1997) Protein Eng., 10, 13631372.
Muegge,I., Schweins,T. and Warshel,A. (1998) Proteins, 30, 407423.[CrossRef][Web of Science][Medline]
Nicholls,A., Sharp,K.A. and Honig,B. (1991) Proteins, 11, 281296.[CrossRef][Web of Science][Medline]
Olson,M.A. (1998) Biophys. Chem., 75, 115128.
Olson,M.A. (2001) Biophys. J., 81, 18411853.[Web of Science][Medline]
Olson,M.A. and Cuff,L. (1999) Biophys. J., 76, 2839.[Web of Science][Medline]
Pal,D. and Chakrabarti,P. (1999) Proteins, 36, 332339.[CrossRef][Web of Science][Medline]
Rupp,B. and Segelke,B. (2001) Nature Struct. Biol., 8, 663664.
Schiavo,G., Malizio,C., Trimble,W.S. Polerino de Laureto,P., Milan,G., Sugiyama,H., Johnson,E. and Montecucco,C. (1994) J. Biol. Chem., 269, 2021320216.
Schiavo,G., Shone,C.C., Bennett,M.K., Scheller,R.H. and Montecucco,C. (1995) J. Biol. Chem., 270, 1056610570.
Schmidt,J.J. and Bostian,K.A. (1997) J. Protein Chem., 16, 1926.[CrossRef][Web of Science][Medline]
Shone,C.C. and Roberts,A.K. (1994) Eur. J. Biochem., 225, 26370[Web of Science][Medline]
Shone,C.C., Quinn,C.P., Wait,R., Hallis,B., Fooks,S. and Hambleton,P. (1993) Eur. J. Biochem., 217, 965971.[Web of Science][Medline]
Soleihac,J.M., Cornille,F., Martin,L., Lenoir,C., Fournié-Saluski,M.C. and Roques,B.P. (1996) Anal. Biochem., 241, 120127.[CrossRef][Web of Science][Medline]
Warshel,A. and Åqvist,A. (1991) Annu. Rev. Biophys. Biophys. Chem., 20, 267298.[CrossRef][Web of Science][Medline]
Weiner,S.J., Kollman,P.A., Nguyen,D.T. and Case,D.A. (1986) J. Comput. Chem., 7, 230252.[CrossRef][Web of Science]
Received December 23, 2001; revised May 8, 2002; accepted May 31, 2002.
![]()
CiteULike
Connotea
Del.icio.us What's this?
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||










4; 4 < yellow