Department of Cell Biology and Biochemistry, USAMRIID,1425 Porter Street, Frederick, MD 21702, USA
![]() |
Abstract |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
Keywords: continuum models/electrostatics/hydration/hydrophobic interaction
![]() |
Introduction |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
Here, 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 |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
![]() | (1) |
![]() | (2) |
![]() | (3) |
![]() | (4) |
Calculation of W given in Equation 3
is obtained from the expression
![]() | (5) |
![]() | (6) |
The contribution Gint contains polar and non-polar intermolecular interaction terms between P* and
*:
![]() | (7) |
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) |
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 |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
|
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 1compares 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.
|
|
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 |
---|
2 Present address: Department of Biology, Virginia Polytechnic Institute and State University, Blacksburg, VA 24061, USA
![]() |
Acknowledgments |
---|
![]() |
References |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
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][ISI]
Froloff,N., Windemuth,A. and Honig,B. (1997) Protein Sci., 6, 12931301.
Gilson,M.K. and Honig,B. (1988) Proteins, 4, 718.[ISI][Medline]
Hanson,M.A. and Stevens,R.C. (2000) Nature Struct. Biol., 7, 687692.[CrossRef][ISI][Medline]
Hazzard,J., Sudhof,T.C. and Rizo,J. (1999) J. Biomol. NMR, 14, 203207.[CrossRef][ISI][Medline]
Jackson,R.M. and Sternberg,M.J.E. (1994) Protein Eng., 7, 371383.[Abstract]
Jorgensen,W.L., Chandrasekhar,J., Madura,J.D., Impey,R.W. and Klein,M.L. (1983) J. Chem. Phys., 79, 926935.[CrossRef][ISI]
Karplus,M. and Janin,J. (1999) Protein Eng., 12, 185186.
King,G.F. and Warshel,A. (1991) J. Chem. Phys., 95, 43664377.[CrossRef][ISI]
Montecucco,C. and Schiavo,G. (1995) Q. Rev. Biophys., 28, 423472.[ISI][Medline]
Muegge,I., Tao,H. and Warshel,A. (1997) Protein Eng., 10, 13631372.[Abstract]
Muegge,I., Schweins,T. and Warshel,A. (1998) Proteins, 30, 407423.[CrossRef][Medline]
Nicholls,A., Sharp,K.A. and Honig,B. (1991) Proteins, 11, 281296.[ISI][Medline]
Olson,M.A. (1998) Biophys. Chem., 75, 115128.[CrossRef][ISI]
Olson,M.A. (2001) Biophys. J., 81, 18411853.
Olson,M.A. and Cuff,L. (1999) Biophys. J., 76, 2839.
Pal,D. and Chakrabarti,P. (1999) Proteins, 36, 332339.[CrossRef][ISI][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.[ISI][Medline]
Shone,C.C. and Roberts,A.K. (1994) Eur. J. Biochem., 225, 26370[Abstract]
Shone,C.C., Quinn,C.P., Wait,R., Hallis,B., Fooks,S. and Hambleton,P. (1993) Eur. J. Biochem., 217, 965971.[Abstract]
Soleihac,J.M., Cornille,F., Martin,L., Lenoir,C., Fournié-Saluski,M.C. and Roques,B.P. (1996) Anal. Biochem., 241, 120127.[CrossRef][ISI][Medline]
Warshel,A. and Åqvist,A. (1991) Annu. Rev. Biophys. Biophys. Chem., 20, 267298.[CrossRef][ISI][Medline]
Weiner,S.J., Kollman,P.A., Nguyen,D.T. and Case,D.A. (1986) J. Comput. Chem., 7, 230252.[ISI]
Received December 23, 2001; revised May 8, 2002; accepted May 31, 2002.