Protein Engineering, Vol. 13, No. 1, 21-26,
January 2000
© 2000 Oxford University Press
Effect of the reaction field electrostatic term on the molecular dynamics simulation of the activation domain of procarboxypeptidase B
Departament de Quimica Analítica, Universitat de Barcelona, Martí i Franqués, 1-11, 08028, Barcelona and 1 Institut de Biologia Fonamental and Departament de Bioquímica, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain
 |
Abstract |
|---|
 Top Abstract IntroductionMaterials and methodsResultsDiscussionReferences |
|---|
Molecular dynamics simulations of the activation domain of porcine procarboxypeptidase B (ADBp) were performed in order to examine the effects of the inclusion of a reaction field (RF) term into the calculation of electrostatics forces for highly charged proteins. Two simulations were performed with the GROMOS96 package, studying the influence of counterions on the final results. Comparison with previous results without the inclusion of the RF term (Martí-Renom,M.A., Mas,J.M., Oliva,B., Querol,E. and Avilés,F.X., Protein Engng, 1998, 11, 101–110) shows that the structure is well maintained when the RF term is included. Moreover, the analysis of the trajectories shows that simulations of solvated highly-charged proteins are sensitive to the presence of counterions, the secondary structures being more stable when their charges are neutralized.
Keywords: electrostatics/molecular dynamics/procarboxypeptidase B/reaction field/salt effects
|
Introduction |
|---|
|
TopAbstractIntroductionMaterials and methodsResultsDiscussionReferences |
|---|
The conformational stability of a protein results from a large array of local and non-local interactions, electrostatic interactions being one of the most important contributions to the latter. The occurrence of stabilizing or destabilizing electrostatic effects in proteins can be tested by determining the effects of salt concentration upon protein stability. This is based on the common assumption that a high salt concentration screens electrostatic interactions (Kohn et al., 1997
).
While experimental knowledge is essential for the understanding of the effects of counterions on the structure and dynamic properties of proteins in solution, theoretical studies and computer simulations, such as molecular dynamics (MD), complement the experimental data. One of the most demanding systems for such studies is the highly charged globular proteins lacking disulfide bridges. A previous study examined the influence of the counterions and volume on the simulated dynamics of the activation domain of procarboxypeptidase B (ADBp), a highly charged globular protein (Martí-Renom et al., 1998). The enlargement of the water box stabilizes the system and a similar trend was also observed when counterions were included in the MD simulation. However, the overall structure of ADBp was not maintained during the simulations owing to unsatisfactory treatment of electrostatic forces, and all the secondary structures of ADP were distorted.
The most important factor in the calculation of forces in MD simulations is that of the electrostatic interactions. Besides, the calculation of the electrostatic forces takes most of the CPU time. In systems with a large number of atoms, the computational load can be reduced by the use of rapid but approximate methods. Two areas in which accuracy is exchanged for speed are long-range interactions and electronic polarizability (Gilson, 1995). Electrostatic interactions are long-range forces and truncation, i.e. the use of cut-offs significantly affects results (Schreiber and Steinhauer, 1992). Moreover, long MD simulations of proteins without electrostatic cut-offs yield trajectories that resemble known crystallographic structures more than similar simulations with cut-offs.
Improvements in MD studies of highly charged systems that reduce CPU time have been proposed (Berendsen, 1993; Luty et al., 1994; Saito, 1994; Cheatham et al., 1995). Some of them separate the force of each charge into short- and long-range parts, and use special methods for dealing with the latter. Among them, the inclusion in the calculation of electrostatic forces of a reaction field (RF) term based on the Poison–Boltzman approach has been used in the study of polar and ionic systems (Tironi et al., 1995). Here we report the effect of the inclusion of an RF term in the treatment of the electrostatic contribution for the study of ADBp dynamics in solution.
We used as a model the activation domain of porcine procarboxypeptidase B, which corresponds to the 75 residue N-terminal globular part of the pro-segment. The domain is attached to the carboxypeptidase B moiety by a 15 residue-connecting segment. It consists of a four-stranded antiparallel ß-sheet with two 
-helices and one 310-helix packed over its external face in an antiparallel /antiparallel ß topology (Figure 1). The structure has two internal salt bridges (4Glu–49Arg, 49Arg–30Asp) and there are no disulfide bridges in this domain. The activation domain has charged residues and its stability in solution depends on the ionic composition of the medium (Conejo-Lara et al., 1991).
|
|
Materials and methods |
|---|
|
TopAbstractIntroductionMaterials and methodsResultsDiscussionReferences |
|---|
The latest release of the GROMOS package (van Gunsteren et al., 1996
) was used to perform solvated simulations using SPC/E water molecules. The crystallographic coordinates of ADBp, obtained at 2.3 Å resolution (Coll et al., 1991), were used to seed the MD simulations (see the technical aspects in Table I).
|
The initial X-ray structure was minimized with 100 steps of steepest descent (Levitt and Lifson, 1969
) without SHAKE (Ryckaert et al., 1997) and 100 steps with SHAKE, all in a vacuum. The minimized protein was solvated by a rectangular box of water molecules with the PROBOX program (van Gunsteren and Berendsen, 1977). The solvent was energetically minimized using 500 steps of steepest descent and positional constraints over the protein structure. In simulation 1, which is representative of a non-counterbalanced protein charged system, the minimized solvated structure was the starting point for an MD simulation run under periodic boundary conditions at constant temperature and pressure and neutral pH.
For simulation 2 the charges of ADBp were neutralized by seven Cl– ions and 18 Na+ ions. Therefore, the water molecules closest to the charges in the minimized solvated structure were replaced by the counterions and the system was again energetically minimized using 100 steps of steepest descent. This minimized structure (protein, solvent and ions) was the starting point for the MD simulation in the same conditions of temperature, pressure and pH. The time step was 0.002 ps.
The Coulomb force acting on a charge qi, at the center of the cut-off sphere of radius Rrf, due to a charge qj, for solute and solvent sites, including the PB generalized reaction field contribution is
|
|
0)–1 is 138.9354 kJmol–1e–2nm and 1 = 1. Crf is the coefficient governing the size of the reaction field force
 |
2 is the relative dielectric permittivity of the electrostatic continuum (2 = 54) and 
is an inverse Debye screening length outside the reaction field cut-off radius Rrf ( = 0 nm–1). The first and second terms in eqn 1 refer to the short- and long-range forces, respectively. The last term is constant and has been added to make the interaction zero at the reaction-field cut-off distance Rrf. A charge group pairlist was built in the first step of the simulation and updated every 5 steps. A cut-off radius of 8 Å was chosen to select nearest-neighbor charge groups and a cut-off radius of 12 Å for the long range electrostatics. The continuum RF term was calculated for the region beyond 12 Å (Rrf).
Four physical properties of the system were analyzed throughout the simulations (root mean square deviation, temperature B-factors, radius of gyration and hydrogen bonding net). The root mean square deviation (r.m.s.d.) was analyzed in order to characterize the conformational changes of the protein. MD B-factors derived from the matrix of the atomic fluctuations were calculated using the equation B-factor = 82(r.m.s.d.)2/3. The analysis of the changes in the radius of gyration provided an additional measure of the global changes in the protein structure. Finally, the hydrogen bonding net was used to characterize the stability and distortions of the secondary structure.
|
Results |
|---|
|
TopAbstractIntroductionMaterials and methodsResultsDiscussionReferences |
|---|
Structural properties of the protein system
Figure 2 shows a time series of r.m.s.d. values for the whole protein, -helices and ß-sheet C atoms for simulations 1 and 2 (in the absence and presence of counterions, respectively). R.m.s.d. results for the whole protein (Figure 2a) show that both simulations reached structural equilibrium. Simulation 1 shows r.m.s.d. values near 3 Å from 600 to 2000 ps. The whole of simulation 2 between 400 and 2000 ps showed the best agreement with respect to X-ray experimental structure, with r.m.s.d. values near 2 Å. Figure 2b and c show the r.m.s.d. values for the -helix 1 and -helix 2. For both simulations, the values are lower than 1 Å. Finally, Figure 2d shows the r.m.s.d. values for the ß-sheet. As in the case of the whole protein, the whole of simulation 2 shows the lowest r.m.s.d. series from 400 to 2000 ps. R.m.s.d. values for simulation 1 were near 1.3 Å and near 1 Å for simulation 2.
|
Table II
shows some of the structural parameters of the X-ray structure and of the averaged structures for simulations 1 and 2. R.m.s.d. values for C atoms of acidic and basic groups are higher in simulation 1 than in simulation 2, reflecting the influence of counterions on their structural stability. Only 310-helix r.m.s.d. values are slightly worse in simulation 2 than in simulation 1. R.m.s.d. values of heavy atoms show a similar trend to that described for r.m.s.d. values of C atoms.
|
The radius of gyration time series for both simulations are shown in Figure 3
and the mean values are also included in Table II
. The radius of gyration through the dynamics gives information on the overall shape behavior of the protein. Radii of gyration larger than that in the original imply an expansion of the structure. The graph shows that the whole protein reaches a structural equilibrium at 11.2 Å, a value which is near the value obtained for the X-ray structure (11.0 Å). From these results, it can be deduced that the packing is maintained during both simulations.
|
Figure 4
shows temperature B-factors for the side chains as a function of residue number. The results were obtained from an average structure of the last 100 ps of each simulation. The B-factor values are quite similar for both simulations and the highest values are located at the N-terminus, the C-terminus and the region connecting ß-strand 2 and ß-strand 3. The presence of counterions decreases the value of the fluctuations in some charged residues. This effect was particularly observed in the range of residues from 50–70 in -helix 2 (Glu53, Asp54, Glu63 and Glu66). On the other hand, simulation 2 shows higher B-factor values for the residues located in ß-strand 4.
|
Protein–protein hydrogen bonds
Table III shows the backbone hydrogen bonds found in the X-ray structure and the frequency of their occurrence in the simulations (>40%). In general, most of the hydrogen bonds present initially in the X-ray structure were maintained throughout the simulations. Hence, hydrogen bonds in ß-sheet 1, ß-sheet 3 and ß-sheet 4 structures were maintained in both simulations and the same happened with hydrogen bonds in -helix 1 and -helix 2. In contrast, ß-strand 2 lost most of its original hydrogen bonds in both simulations, as was the case for the 310-helix.
|
|
Discussion |
|---|
|
TopAbstractIntroductionMaterials and methodsResultsDiscussionReferences |
|---|
Protein structures, dynamics and stability can be simulated on computers under conditions similar to those observed experimentally (Honig and Yang, 1995
). The ideal simulation would be that in which the protein concentration is of the same order of magnitude as in the experimental study. However, such a simulation would need the inclusion of a large number of water molecules, increasing notably the CPU time. In a previous study of the ADBp system with the GROMOS87 package (Martí-Renom et al., 1998), it was shown that an increase in the box dimensions stabilizes the protein structure in the MD simulations, both in the presence and in the absence of counterions. The box dimensions optimized in that study were applied here.
From the results obtained when the RF term was not included in the calculation of electrostatic forces, it was concluded that the whole structure of ADBp was not well maintained for all simulations (Martí-Renom et al., 1998). This was probably due to the fact that ADBp is a highly charged protein with a clear excess of negative charges, so a large number of counterions and an appropriate treatment of the electrostatics was considered necessary in a later stage to compensate the system. The present work has been carried out in similar conditions to those described in the previous study but including the RF term in order to assess the importance of this correction on the final results.
In the present work, the structure of ADBp was maintained for the simulations carried out with and without counterions in the system. R.m.s.d. values for both simulations show that the protein reached an equilibrium in a short time. Moreover, r.m.s.d. values for both simulations in the equilibrium were acceptable, particularly in simulation 2 (r.m.s.d. near 2 Å). A similar pattern was also observed for the individual secondary structures. On the whole, -helices are the most stable part of the entire protein, particularly in the presence of counterions. This is probably due to the fact that -helix 1 is the most charged secondary structure of ADBp and the presence of counterions is needed to equilibrate the electrostatic forces. This trend was also observed when this system was studied without the inclusion of the RF term into the calculation of electrostatic forces (Martí-Renom et al., 1998). ß-Sheets were also maintained in the simulations showing r.m.s.d. values under 2 Å. As in the case for -helices, the inclusion of counterions also improved the stability of their structure. Surprisingly, this was not observed for the 310-helix, where the largest fluctuations were observed when counterions were added.
A partial unfolding of ADBp during simulations without the RF term is observed when water molecules interact with the protein core (Martí-Renom et al., 1998), competing with the original protein–protein hydrogen bonds. In the present study, water penetration into the protein core was not observed, a fact which is reflected in the radius of gyration, which was equilibrated along the simulations around 11 Å.
Figure 5 shows the averaged structures calculated for simulation 1 (in the time range from 600 to 2000 ps) and for simulation 2 (from 400 to 2000 ps). The averaged structures are similar to the initial X-ray structure shown in Figure 1a and they do not show serious unfolding. Only the structure obtained in simulation 1 is slightly different from the X-ray structure; the main difference is the small displacement of the -helices from a parallel arrangement in the X-ray structure to a slightly crossed arrangement in the average structure from simulation 1. This small displacement is reflected in the loss of the hydrogen bond in the X-ray structure between 66Glu and 19Glu (-helix 1) in simulation 1. A new hydrogen bond appeared between 26Thr (-helix 1) and 65Asn, which was not present in the X-ray structure. On the other hand, simulation 2 shows a similar spatial arrangement of -helices and hydrogen bond pattern to that observed in the X-ray structure. Hence, the hydrogen bond between 66Glu and 19Glu was maintained in the simulation when counterions were present. This is probably due to the fact that the charges on both Glu were screened by ions.
|
Both simulations show the loss of hydrogen bonds in the ß-strand 2 and 310-helix moieties. This is also reflected in the high r.m.s.d. values for this part of the protein (see Table II
). An explanation for this behavior could be that this region is the part of the activation domain which interacts with the rest of the protein, which was not simulated in this study. The lack of the rest of the protein prevents the stabilization of this region of the activation domain. On the other hand, hydrogen bonds in -helices were maintained which reflects the fact that these structures are stable throughout the simulations.
Ionic interactions are the main focus of a study of highly charged biomolecules. In proteins, ionic interactions are extremely complex because they can exist on the fully exposed surface of the protein, fully buried in the interior of the protein or in a partially buried environment with varying degrees of hydrophobicity of the surrounding residues. In addition, buffer conditions (salts and pH) can have dramatic effects on the contribution of ion pairs to protein stability (Goto et al., 1990; Kohn et al., 1997).
The stability of the ADBp conformation in simulation 1, carried out without counterions, can be related to the presence of hydrophobic residues (Leu23, Ile29, Val58, Leu62 and Leu67) that form a hydrophobic face on the -helices of the protein. Theoretical calculations (Hendsch and Tidor, 1994; Sindelar et al., 1998) and experimental studies (Wimley et al., 1996) suggested that hydrophobic interactions are more stabilizing than salt bridges in protein folding. That is, salt bridges in proteins are generally destabilizing relative to replacement with hydrophobic groups of the same size and shape. Kohn and co-workers also studied the role of surface-accessible ion pairs in protein stability by determining the effect of added salt on the stability of two-stranded -helical coiled-coils (Zhou et al., 1992; Kohn et al., 1997). Moreover, there are also charged residues in -helices, which have been predicted to be involved in interhelical ion pairs and thereby to contribute to coiled-coil stability (Talbot and Hodges, 1982).
In addition to the hydrophobic interactions, the higher stability of ADBp in simulation 2 is due to the presence of counterions. Here, it has been shown that the presence of counterions stabilizes ADBp, reaching structural equilibrium in the protein 200 ps before that in simulation 1. The higher stability of simulation 2 is also reflected in the stabilization of the internal hydrogen bonds. This stabilization is correlated with the decrease in fluctuation of the charged side chain when a counterion neutralizes the residue charge, as shown in the B-factor analyses. Table II shows that the higher r.m.s.d. values obtained in simulation 1 are those related to charged groups, the acidic groups being the most affected. The results of this study illustrate the influence of salt concentration on the stability of proteins. In addition to the ionic strength effect of general charge screening, specific counterion-binding to charged residues on the protein as well as effects of the salt on the solvent structure and properties influence protein stability (Goto et al., 1990).
MD simulations of the ADBp system yield a globally stable system with small yet significant internal readjustments. The results confirmed the quality of the standard GROMOS96 force field and the validity of including the RF term in the calculation of electrostatics forces for highly charged biomolecules. The present study complements the work of Roxström et al. (1998) on MD simulations of the DNA–protein complex. Those results show that the GROMOS87 force field can provide reliable results in the study of such highly-charged systems with some corrections for the appropriate representation of coion (ion atmosphere) effects on the counterions, and corrections for hydrophobicity default in the water to solute parameter.
|
Acknowledgments |
|---|
The support by C4 (CEPBA-CESCA) is gratefully acknowledged. This work was supported by grants BIO97-0511 and BIO98-0362 from the CICYT (Ministerio de Educación y Ciencia, Spain).
|
Notes |
|---|
2 To whom correspondence should be addressed; fx.aviles{at}blues.uab.es
|
References |
|---|
|
TopAbstractIntroductionMaterials and methodsResultsDiscussionReferences |
|---|
Berendsen,H.J.C. (1993) In Van Gunsteren,W.F., Weiner,P.K. and Wilkinson,A.J. (eds), Computer Simulation of Biomolecular Systems. Theoretical and Experimental Applications, Vol. 2. Leiden, ESCOM, pp. 161–181.
Cheatham,I.T.E., Miller,J.L., Fox,T., Darden,T.A. and Kollman,P.A. (1995) J. Am. Chem. Soc., 117, 4193–4194.
Coll,M., Guasch,A., Avilés,F.X. and Huber,R. (1991) EMBO J., 10, 1–9.[Web of Science][Medline]
Conejo-Lara,F., Sánchez-Ruiz,J.M., Mateo,P.L., Burgos,F.J., Vendrell,J. and Avilés,F.X. (1991) Eur. J. Biochem., 200, 663–670.[Web of Science][Medline]
Gilson,M.K. (1995) Curr. Opin. Struct. Biol., 5, 216–223.[Web of Science][Medline]
Goto,Y., Takahashi,N. and Fink,A.L. (1990) Biochemistry, 29, 3480–3488.[Medline]
Hendsch,Z.S. and Tidor,B. (1994) Protein Sci., 3, 211–226.[Web of Science][Medline]
Honig,B. and Yang,A.W. (1995) Adv. Protein Chem., 46, 27–58.[Web of Science][Medline]
Kohn,W.D., Kay,C.M. and Hodges,R.S. (1997) J. Mol. Biol., 267, 1039–1052.[Web of Science][Medline]
Levitt,M. and Lifson,S. (1969) J. Mol. Biol., 46, 269–279.[Web of Science][Medline]
Luty,B.A., Davis,M.E., Tironi,I.G. and Van Gunsteren,W.F. (1994) Mol. Simul., 14, 11–20.
Martí-Renom,M.A., Mas,J.M., Oliva,B., Querol,E. and Avilés,F.X. (1998) Protein Engng, 11, 101–110.
Roxström,G., Velázquez,I., Paulino,M. and Tapia,O. (1998) J. Phys. Chem. B, 102, 1828–1832.
Ryckaert,J.P., Ciccotti,G. and Berendsen,H.J.C. (1977) J. Comp. Chem., 23, 327–341.
Saito,M. (1994) J. Chem. Phys., 101, 4055–4061.
Screiber,H. and Steinhauer,O. (1992) Biochemistry, 31, 5856–5860.[Medline]
Sindelar,C.V., Hendsch,Z.S. and Tidor,B. (1998) Protein Sci., 7, 1898–1914.[Web of Science][Medline]
Talbot,J.A. and Hodges,R.S. (1982) Acc. Chem. Res., 15, 224–230.
Tironi,I.G., Sperb,R., Smith,P.E. and Van Gunsteren,W.F. (1995) J. Chem. Phys., 102, 5451–5459.
van Gunsteren,W.F. and Berendsen,H.J.C. (1977) Mol. Phys., 34, 1311–1327.[Web of Science]
van Gunsteren,W.F., Billeter,S.R., Eising,A.A., Hünenberger,P.H., Krüger,P., Mark,A.E., Scott,W.R.P. and Tironi,I.G. (1996) In Groningen Molecular Simulation (GROMOS) Library Manual. BIOMOS B.V., Zürich.
Wimley,W.C., Gawrisch,K., Creamer,T.P. and White,S.H. (1996) Proc. Natl Acad. Sci. USA, 93, 2985–2990.
Zhou,N.E., Kay,C.M. and Hodges,R.S. (1992) Biochemistry, 31, 5739–5746.[Medline]
Received July 28, 1999; revised November 10, 1999; accepted November 16, 1999.
CiteULike 
Connotea 
Del.icio.us What's this?
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||