1 Faculty of Chemistry, University of Gdansk, ul. Sobieskiego 18, 80-952, Gdansk and 2 Department of Public Health, University School of Physical Education, ul. Wiejska 1, 80-336, Gdansk, Poland
![]() |
Abstract |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
Keywords: constrained simulated annealing/counterions/molecular dynamics/papain/protein stability
![]() |
Introduction |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
Thus, the electrostatic fields in and around biological macromolecules appear to be of central importance for both structure and function. The development of computational methods for the accurate determination of protein structure and dynamics is of both fundamental and applied importance (Zheng et al., 1993). However, the modeling and calculation of these fields presents a serious challenge (Gilson et al., 1987
). While molecular simulations have proved to be very successful, there are still many questions to be answered concerning the various approximations that are commonly used in order to make the calculations computationally practical (Smith and Pettitt, 1991
). Among molecular dynamics (MD) simulations, one of the challenges is the incorporation of the counterions into protein simulations (Ibragimova and Wade, 1998
; Marti-Renom et al., 1998
). While experimental knowledge is essential to understand the effects of ions on the structure and dynamic properties of proteins in solution, theoretical studies and computer simulations are necessary approaches to supplement the experimental data (Marti-Renom et al., 1998
).
The problem of computationally representing charged biological macromolecules can be clearly exemplified in the treatment given to DNA whose conformational stability is modulated by the surrounding environment and is dependent upon the neutralization of its charged groups. Although it is still not clear if counterions allow a much better reproduction of the properties of DNA, it seems that the addition of counterions improves behavior of the system in MD simulations (Cheatham and Kollman, 1996; Tapia and Velazquez, 1997
; Flatters and Lavery, 1998
; Pastor et al., 1999
). Also, studies were conducted to investigate counterions in simulations of not only DNA alone, but in complexes with polypeptide in proteinnucleic acid recognition research (Bishop et al., 1997
; Kosztin et al., 1997
; Roxstrom et al., 1998a
,b
). However, in the case of proteins, the value of counterion simulations has not been studied to such a great extent as for DNA. For proteins, charge-balancing counterions are positioned close to charged groups in the regions with the most favorable electrostatic potential for binding on the basis of either geometric or energetic criteria (Ibragimova and Wade, 1998
). The addition of ions can strongly alter the thermodynamic and physical properties of proteins in solution. This is experimentally evident in the process of denaturation and salting proteins in and out of solution. The conformation of proteins is sensitive to the composition of the solution environment (Marlow et al., 1993
).
Here, we investigate sensitivity to the placement and modeling of counterions in protein simulations. It is hoped that the latter would help in analyzing the factors which allow us to simulate a more physically correct and biologically realistic behavior. For the simulation in the present study, papain (EC 3.4.22.2) from Carica papaya was chosen. It is a globular 212-residue protein, with a net charge of +10e, three disulfide bridges, and has been investigated several times for charge distribution and related studies (Wada and Nakamura, 1981; Nakamura and Wada, 1985
; Wada et al., 1985
).
![]() |
Computational methods |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
Modeling procedures
Starting co-ordinates for all heavy atoms of two papain structures were obtained from the X-ray data [Protein Data Bank (PDB) files 1pe6 and 1pad for systems A and B, respectively]. The ligands were removed and the C25 sulfhydryl was left unmodified. Both systems were subjected to two different procedures denoted as Procedure I and II, respectively (Table I). The details of Procedure I are described in our previous article (Czaplewski et al., 1999
). Two variants have been generated for this procedure via choosing a different seed for the random number generator, giving eight systems (two different structures, two different seeds and the presence of the counterions). In Procedure II, the crude systems were initially minimized in vacuo to remove close van der Waals contacts. This (and all the following) minimization(s) consisted of 20 000 steps (after 1000 steps the minimizer was switched from the steepest descent to the conjugate gradient mode). Afterwards, the system was heated in constrained simulated annealing (CSA) in vacuo from 0 to 1000 K for 1 ps, then, during thermal stabilization of 10 ps duration, snapshots were generated every 1 ps. From the two sets of 10 snapshots, three snapshots of each set (Nos 5, 8, 9 and 1, 6, 10 for systems A and B, respectively) were randomly selected for further simulations. The snapshots were subsequently cooled down for 17 ps (Table II
; the temperature relaxation times for particular snapshots were assigned individually to assure the slowest cooling possible). During the CSA general constraints (positional: all C
atoms in papain; all peptide bonds; and all improper dihedrals to retain proper chirality) were used. Subsequently, the systems were diluted in TIP3P water (Jorgensen et al., 1983
), minimized, heated and equilibrated (5 ps of heating, 12 ps for water equilibration, and finally 10 ps for equilibration of the whole system; 27 ps in total) as described in Table III
. The 12 ps time for water equilibration was selected, because the solvent equilibration phase should be sufficiently extensive to allow the solvent to completely readjust to the potential field of the solute. For MD this implies that the length of this phase should be longer than the relaxation time of the solventthe time taken for a molecule to lose any memory of its original orientation, which for water is about 10 ps (Leach, 1996
). Thus, 12 ps should be enough for water surrounding a protein, which seems to be confirmed in our study by monitoring the physical parameters of a box (data not shown). All simulations were performed with and without adding counterions. Before surrounding the protein by water molecules, 10 Cl- charge-balancing counterions were added in positions with favorable ionresidue interactions to neutralize charges on the protein surface using the AMBER program LEAP, which adds counterions in a shell around a molecule using a Coulombic potential on a grid. Grid resolution was 1 Å. Thus, in total, 20 systems were generated for simulations: eight for Procedure I (A1-1, A1-1-I, A1-2, A1-2-I, B1-1, B1-1-I, B1-2, B1-2-I) and 12 for Procedure II (A2-5, A2-5-I, A2-8, A2-8-I, A2-9, A2-9-I, B2-1, B2-1-I, B2-6, B2-6-I, B2-10, B2-10-I). The first number denotes Procedure number, the second number denotes two different seeds in Procedure I and snapshot number in Procedure II, and I denotes inclusion of the counterions in the simulation.
|
|
|
Potentially equilibrated systems were subjected to MD productive runs. Simulations of complexes in solution were performed under periodic boundary conditions in a closed, isothermal, isobaric (NTP) ensemble. Throughout the simulation the solute and solvent were coupled to a constant-temperature (T = 300 K) heat bath and a constant-pressure (P = 1 bar) bath (Berendsen et al., 1984). All hydrogen-containing bonds were constrained using the SHAKE algorithm (Ryckaert et al., 1977
) for holonomic constraints with a relative geometric tolerance for co-ordinate resetting of 0.0005 Å, allowing a time step of 1 fs. The Leapfrog version (Hockney, 1970
) of the Verlet algorithm (Verlet, 1967
) was employed to integrate the equations of motion. A double residue-based cutoff distance of 10/14 Å was used for non-bonded interactions. The TIP3P model was used for water molecules. A typical rectangular box size was 78x66x60 Å. Approximately 8200 TIP3P water molecules were in the box, i.e. 28 000 atoms in total. Co-ordinates were saved every 1000 steps. A set of 230 and 300 ps MD runs was carried out during Procedures I and II, respectively.
Parametric mean square approximation
To quantitatively compare time-averaged C atoms mobilities as a function of sequence in dependence on simulation conditions, the minimum of
, where
is a function of one variable (a), was sought:
![]() | (1) |
Parameter a was introduced because average C mobilities differ from X-ray temperature factors by approximately one order of magnitude. Parameter a was calculated to minimize
:
![]() | (2) |
Statistical analysis
Statistical calculations were performed using the STATISTICA 5.0 program (StatSoft, Inc., 1995). One-way analysis of variance (ANOVA) of repeated measurements (two levels of within-subjects factor: ions versus non-ions) was used to determine the differences in MD simulation behavior with and without counterions. Normality of the distribution of the data was assessed by a KolmogorowSmirnov test with Lilliefors probabilities (Lilliefors, 1967) and ShapiroWilks analysis (Shapiro et al., 1968
) with Royston's algorithm (Royston, 1982
). Homogeneity of variance was evaluated by means of Levene's test. Results of normality and homogeneity of variance tests allowed the usage of ANOVA statistics (data not shown). Statistical significance was defined at the p < 0.05 level.
![]() |
Results and discussion |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
Experimentally, it has been shown elsewhere that the increase in the concentration of added salts leads to significant variations in the stability of proteins, usually to stabilization of non-chaotropic salts, probably owing to changes in the apparent hydrophobic effect and charge screening (Kohn et al., 1997). In the theoretical studies, improved protein stability has been observed several times when adding counterions (York et al., 1993
; Yelle et al., 1995
) which screen charged side chains, preventing unnatural electrostatic interactions between neighboring parts of the macromolecule (Ibragimova et al., 1998). A similar conclusion was reached from the studies of the highly charged activation domain of procarboxypeptidase B, which suggests that the use of counterions significantly contribute to the maintenance of native secondary structure (Marti-Renom et al., 1998
). Ibragimova and Wade (Ibragimova and Wade, 1998
) investigated counterion effects in simulations of the WW domain of YAP (Yes kinase-associated protein). The small size and weak stability of this domain make it very sensitive to conditions in MD simulations, and thus, in the opinion of the authors, they constitute a good system on which different protocols and models could be tested. The usual procedure of only placing charge-balancing counterions near charged groups on the protein surface did not result in a stable protein fold, nor did simulation without any counterions. For maintaining the stability of the WW domain, explicit modeling of full ionic strength (at 0.2 M) was necessary.
As it may seem quite natural that adding ions would help to stabilize highly charged polypeptides or proteins with marginal stability, the aim of our study was to investigate if such a procedure influences a typical protein of the type used widely in MD research, it not being a sophisticated example. Thus, we have used a different approach in choosing the model to simulate. Papain is a typical 212-residue globular protein and, when well equilibrated, is very stable during long MD simulations giving relaxed stable trajectories (data not shown). Four systems (two different starting structures and two different modeling procedures) in two variants each (with and without the counterions) were generated and investigated to maximize different simulation conditions and, in that way, to eliminate (at least partly) the modeling setup.
As MD simulations are chaotic and critically dependent on the starting conditions (Braxenthaler et al., 1997), to strengthen the statistical inference 10 sets of different runs (with a total time of about 5.5 ns) have been performed with differences in seeds for the random number generator, investigated systems and equilibration procedures. Figure 1
shows the root mean square deviations (RMSD) for all simulated systems. Subjected to Procedure I non-neutralized systems failed to equilibrate, the RMSD values tended to increase until the very end of the MD run reaching 5 Å in the B1-1 simulation. The influence of added counterions on these systems is somewhat ambiguous. In system A ions efficiently helped the protein to equilibrate, drastically reducing the overall RMSD (Table IV
). Surprisingly, the system B reacts on the same procedure in a different fashion: regardless of net charge, the system fluctuations rapidly increased throughout the simulation, although in the case of neutralized simulations they were much smaller in size (Table IV
). The equilibration protocol used in Procedure I apparently gives no steady-state phase (in terms of RMSD). However, it was used successfully in many MD simulations performed on similar systems (see for example: Kazmierkiewicz et al., 1997
; Czaplewski et al., 1999
). It should be stressed that although there is an insufficiency in the applied simulation protocol, it still gives some insight into the problem of the addition of counterions.
|
|
It seems reasonable that simulations with counterions should behave better; they render more accurately biophysical conditions (solvent composition) in living organisms. When performing simulations without counterions, strong electrostatic interactions between charged side chains could disturb the structure of the protein during the simulation.
The secondary structure was examined to characterize conformational changes of the papain during simulations. The results of parameter fitting (see Computational methods section) of the calculated structures to the crystallographic data, giving quantitative measure of the total charge influence on the stability of various elements of the secondary structure, are shown in Table V. The scatter of different responses of various systems' secondary structures to counterions addition is clearly seen. It is hardly possible to find any order in these responses: some (6 of 10) systems apparently are stabilized by ions in the regions of both regular (
-helices, ß-sheets and turns) and irregular (mostly loops) secondary structures (A1-1, A1-2, B1-1, B1-2, A2-5 and B2-6 simulations), the A2-8 system seems to be stabilized by counterions in all residues except those forming
-helices, whereas neutralization driven improvement of A2-9 and B2-1 systems is clear only in mobile (irregular) residues. The B2-10 system is stabilized only in terms of turn regions. It seems that there is no significant superiority of counterions usage on the basis of the analysis of the four pairs of MD runs presented. It is rather hard to evaluate counterions effects on the local fluctuations of molecule mobility because this strongly depends on the details of specific simulations that were performed.
|
![]() |
Notes |
---|
![]() |
Acknowledgments |
---|
![]() |
References |
---|
![]() ![]() ![]() ![]() ![]() ![]() |
---|
Bishop,T.C., Kosztin,D. and Schulten,K. (1997) Biophys. J., 72, 20562067.[Abstract]
Bodkin,M.J. and Goodfellow,J.M. (1995) Protein Sci., 4, 603612.
Braxenthaler,M., Unger,R., Auerbach,D., Given,J.A. and Moult,J. (1997) Proteins, 29, 417425.[ISI][Medline]
Case,D.A. et al. (1997) Amber 5.0. University of California, San Francisco.
Cheatham,T.E. and Kollman,P.A. (1996) J. Mol. Biol., 259, 434444.[ISI][Medline]
Cornell,W.D., Cieplak,P., Bayly,C.I., Gould,I.R., Merz,K.M.,Jr, Ferguson,D.M., Spellmeyer,D.C., Fox,T., Caldwell,J.W. and Kollman,P.A. (1995) J. Am. Chem. Soc., 117, 51795197.[ISI]
Czaplewski,C., Grzonka,Z., Jaskólski,M., Kasprzykowski,F., Kozak,M., Politowska,E. and Ciarkowski,J. (1999) Biochim. Biophys. Acta, 1431, 290305.[ISI][Medline]
Flatters,D. and Lavery,R. (1998) Biophys. J., 75, 372381.
Gilson,M.K., Sharp,K.A. and Honig,B.H. (1987) J. Comput. Chem., 9, 327335.[ISI]
Hockney,R.W. (1970) Methods Comput. Phys., 9, 136211.
Honig,B. and Yang,A.-S. (1995) Adv. Protein Chem., 46, 2758.[ISI][Medline]
Ibragimova,G.T. and Wade,R.C. (1998) Biophys. J., 74, 29062911.
Jorgensen,W.L., Chandreskhar,J., Madura,J.D., Imprey,R.W. and Klein,M.L. (1983) J. Chem. Phys., 79, 926935.[ISI]
Jug,K. and Gerwens,H. (1998) J. Phys. Chem. B, 102, 52175227.[ISI]
Karplus,M. and McCammon,J.A. (1983) Annu. Rev. Biochem., 52, 263300.[ISI][Medline]
Kazmierkiewicz,R., Czaplewski,C., Ciarkowski,J. (1997) Acta Biochim. Pol., 44, 453466.[ISI][Medline]
Kohn,W.D., Kay,C.M. and Hodges,R.S. (1997) J. Mol. Biol., 267, 10391052.[ISI][Medline]
Kosztin,D., Bishop,T.C. and Schulten,K. (1997) Biophys. J., 73, 557570.[Abstract]
Leach,A.R. (1996) Molecular Modelling. Principles and Applications. Longman Singapore Publishers (Pte) Ltd, Singapore, p. 326.
Lilliefors,H. (1967) J. Am. Stat. Assoc., 64, 399402.
Marlow,G.E., Perkyns,J.S. and Pettitt,B.M. (1993) Chem. Rev., 93, 25032521.[ISI]
Marti-Renom,M.A., Mas,J.M., Oliva,B., Querol,E. and Aviles,F.X. (1998) Protein Eng., 11, 881890.[Abstract]
Nakamura,H. and Wada,A. (1985) J. Phys. Soc. Jpn, 54, 40474052.[ISI]
Nolting,B., Golbik,R., Neira,J.L., Soler-Gonzalez,A.S., Schreiberg,G. and Ferscht,A. (1997) Proc. Natl Acad. Sci. USA, 94, 826830.
Pastor,N., MacKerell,A.D. and Weinstein,H. (1999) J. Biomol. Struct. Dynam., 16, 787810.[ISI][Medline]
Petsko,G.A. and Ringe,D. (1984) Annu. Rev. Biophys. Bioeng., 13, 331371.[ISI][Medline]
Roxstrom,G., Velazquez,I., Paulino,M. and Tapia,O. (1998a) J. Phys. Chem. B., 102, 18281832.[ISI]
Roxstrom,G., Velazquez,I., Paulino,M. and Tapia,O. (1998b) J. Biomol. Struct. Dynam., 16, 301312.[ISI][Medline]
Royston,J.P. (1982) Appl. Stat., 31, 115124.[ISI]
Ryckaert,J.P., Ciccotti,G. and Berendsen,H.J.C. (1977) J. Comput. Phys., 23, 327341.[ISI]
Shapiro,S.S., Wilk,M.B. and Chen,H.J. (1968) J. Am. Stat. Assoc., 63, 13431372.[ISI]
Smith,P.E. and Pettitt,B.M. (1991) J. Chem. Phys., 95, 84308441.[ISI]
StatSoft, Inc. (1995) STATISTICA for Windows (Computer program manual). StatSoft, Inc., Tusla, OK.
Tapia,O. and Velazquez,I. (1997) J. Am. Chem. Soc., 119, 59345938.[ISI]
Verlet,L. (1967) Phys. Rev., 159, 98103.[ISI]
Wada,A. and Nakamura,H. (1981) Nature, 293, 757758.[ISI][Medline]
Wada,A., Nakamura,H. and Sakamoto,T. (1985) J. Phys. Soc. Jpn, 54, 40424046.[ISI]
Yelle,R.B., Park,N.-S. and Ichiye,T. (1995) Proteins, 22, 154167.[ISI][Medline]
York,D.M., Darden,T.A., Pedersen,L.G. and Anderson,M.W. (1993) Biochemistry, 32, 14431453.[ISI][Medline]
Zheng,Q., Rosenfeld,R., Vajda S. and DeLisi,C. (1993) Protein Sci., 2, 12421248.
Received June 20, 2001; revised July 6, 2001; accepted July 10, 2001.