Adjusting force distributions in functional site of scorpion toxin BMK M1 by cooperative effect of disulfide bonds

BBRC Biochemical and Biophysical Research Communications 351 (2006) 1037–1042 www.elsevier.com/locate/ybbrc

Adjusting force distributions in functional site of scorpion toxin BMK M1 by cooperative effect of disulfide bonds Majid Taghdir, Hossein Naderi-Manesh

*

Department of Biophysics, Faculty of Science, Tarbiat Modarres University, P.O. Box 14115-175, Tehran, Iran Received 26 October 2006 Available online 7 November 2006

Abstract We decided to investigate the influence of the presence of disulfide bonds on the pattern of force distribution in the functional site of scorpion toxin BMK M1 in its functional state. Therefore, a series of short time molecular dynamics (MD) simulations were performed on this toxin in the native state and disulfide bond broken states. The comparison of disulfide bond broken states with the native state showed that the electrostatic potential energy of important functional residues in the reverse turn and C-terminal regions were modulated by the cooperative effect of all disulfide bonds in the molecule. Furthermore, our results revealed that disulfide bonds also play a cooperative role in modulating (1) the amplitude of the fluctuations of the functional segments and (2) the correlation of motions between important functional residue pairs in this toxin. Therefore, we can conclude that the disulfide bonds have cooperation to adjust the pattern of force distribution in the functional site of this toxin in its functional state.  2006 Elsevier Inc. All rights reserved. Keywords: Adjusting force distribution; Cooperative effect; Disulfide bonds; Functional state; BMK M1

Disulfide bonds are the only common covalent crosslinks in proteins. Experimental and theoretical studies so far have revealed the contribution of these covalent bonds to folding, conformational stability, and catalytic activity in proteins. Disulfide bonds as well as non-covalent interactions overcome the loss of conformational entropy associated with folding that destabilizes native conformation [1]. On the other hand, according to the two current models, these cross-links enhance stability mainly through denatured state effects [2–4]. Since disulfide bonds seem to play different roles in different states, one may ask ‘‘what is the role of disulfide bonds in the functional state of proteins?’’ There is some experimental evidence that some disulfide bonds are more important for tuning function than for structural stability and folding [5]. An attractive example is the disulfide bonds in scorpion toxins; it has been reported that in the scorpion toxin scaffold, modifications of conserved and interior cysteine residues could per*

Corresponding author. Fax: +9821 88009730. E-mail address: [email protected] (H. Naderi-Manesh).

0006-291X/$ - see front matter  2006 Elsevier Inc. All rights reserved. doi:10.1016/j.bbrc.2006.10.156

mit modulation of function, without significantly affecting folding efficiency and structure [6]. In addition, disulfide bonds seem to influence the patterns of correlated and concerted motions, which are expected to have a key role as modulators of the functionally important motions [7,8]. Beside these covalent bonds, electrostatic interactions as long-range forces control important aspects of structure and function in proteins. These forces play a crucial role in protein interactions with ligands or other proteins. Since these forces guide binding process, the calculation of electrostatic potentials and the investigation of the factors that influence these forces is a prolific field of research. Electrostatics, in particular, operates at large distances and can enhance or impede bimolecular collision rates [9–11]. Since the pattern of internal motions, correlation of motions, and electrostatic force distributions in proteins seem to play an important role in molecular recognition and function, we investigated the influence of the presence of disulfide bonds on these patterns in the functional state of scorpion toxin BMK M1. In experimental studies, this effect was usually submerged by their more rugged role in

1038

M. Taghdir, H. Naderi-Manesh / Biochemical and Biophysical Research Communications 351 (2006) 1037–1042

folding and stability. A short timescale of MD simulations does not allow significant structural deviations from the native state. Therefore, we can use short time MD simulations to investigate this aspect of disulfide bond efficacy without being influenced by large-scale structural changes. Representative scorpion toxin BMK M1 with 64 residue cross-linked by four disulfide bonds from the scorpion Buthus martensii karsch is an good model for this study (Fig. 1). Since scorpion toxins are able to adopt stable and biologically active structures, constitute a priori interesting candidates as core structure for protein engineering and design [12]. These neurotoxins are known to interact specifically with the voltage-dependent sodium channel and affect sodium conductance in various excitable tissues, and thus serve as important pharmacological tools for the study of excitability and sodium channel structure. BMK M1 belongs to a-like toxins targeting sodium channel and it has been reported that electrostatic interactions play a crucial role in the function of this toxin [13]. This miniprotein is composed of a dense core of secondary structure elements, including an a-helix and a three-stranded antiparallel b-sheet. Three b-turns connect these secondary structures. There are also four loops including Lys8-Cys12, Trp38-Asn44, Arg58-His64 residues and a long loop spanning Val13-Arg18 between the first loop and a-helix. In fact, the segment Lys8-Cys12 is a reverse turn that together with C-terminal loop Arg58-His64 forms functional site in this toxin, called site RC. Three disulfide bridges formed by residues Cys16-Cys36, Cys22-Cys46, and Cys26-Cys48 stabilize the special babb motif in this toxin. The fourth disulfide bond formed by Cys12-Cys63 stabilizes the loop Arg58-His64 in site RC [14]. Elimination of the disulfide bonds by site-directed mutagenesis in this molecule showed that Cys22-Cys46 and Cys26- Cys48 disulfide bonds are essential for the general folding, Cys16-Cys36 bond is a crucial structural element for stabilizing the general fold, and the Cys12-Cys63 bond is essenlong-loop

16 β-turn II

36 63

46

22 strand II

12 reverse turn

alpha-helix 26

48 C-terminal loop strand III strand I

β-turn I

β-turn III

Fig. 1. Structure of scorpion toxin BMK M1.

tial for the toxic activity and the binding property with the Na+ channel [15]. Results of toxicity tests and binding studies have revealed that site RC is a functional site in this toxin. The C-terminal basic residues Arg58, Lys62, and His64 together with Lys8 in the reverse turn are critical for bioactivity of the molecule [13]. So the pattern of force distribution in this site should be important for function in this toxin. In this study, we performed a series of short MD simulation times (2-ns) on BMK M1 in the native state, individual disulfide bond broken states, and also all disulfide bonds broken state. We investigated the influence of the presence of disulfide bonds on the patterns of internal motions, correlation of motions, and electrostatic force distributions in the functional site of this toxin. Methods Computational procedure. The initial coordinate of BMK M1 was obtained from the Protein Data Bank (PDB) with entry code 1SN1 [14]. The same conformations with the reduced form of cysteines were used as the initial structure of the disulfide bond broken states in MD simulations. For all simulations, the 6 Na+ and 8 Cl charge-balancing counterions were added in order to neutralize charges on the protein surface. Both initial setup and dynamics runs were carried out by Amber 8 [16]. ˚ . The protein All calculations were performed with a cutoff distance of 10 A was solvated with TIP3P model of water [17] in an octahedron box with a ˚ distance between the box edge and the closest portion of the minimum 8 A molecule. First, the system was energy minimized for 500 steps of steepest descent minimization to remove close van der Waals contacts and to allow formation of hydrogen bonds between water molecules in the periodic box and the protein. Temperature was increased from 200 to 300 K during a 100-ps MD simulation in the canonical ensemble (NVT). A 100-ps MD simulation in the isobaric–isothermal ensemble (NPT) was carried out to equilibrate the system in constant pressure. The density was stabilized around 1.02 g/cm3 during equilibration phase in constant pressure. Temperature and pressure were controlled by applying a weak coupling method [18], with temperature and pressure relaxation times sT = 1.0-ps and sp = 1.0-ps, respectively. MD simulations were performed in the NPT ensemble for 2-ns during production phase. Time step of 2-fs was used for all simulations and X–H bonds were constrained using the SHAKE algorithm [19]. Translational center of the mass motions was removed every1000 steps and coordinates were saved every 0.4-ps for analysis. Analysis procedure. The Ptraj module of AMBER was used to extract Root-mean-square deviations (RMSD) and root-mean-square fluctuations (RMSF) data from trajectories. One of the most frequently used measures to assess the stability of a MD simulation over the course of time is the RMSD between experimental coordinate and the generated structures in the trajectory, as follows: rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 X exp: RMSD ¼ ðri  rgen: ð1Þ Þ2 i N and rgen: denote the Cartesian coordinates for atom i in the where rexp: i i experimental and generated structures, respectively. The dynamical properties of Ca atoms have been reported to contain sufficient information to investigate the most important motions in proteins [20]. Therefore, RMSF of Ca atoms was used to investigate structural flexibility. Prior to calculating the RMSF, we removed the overall translational and rotational motions by superimposing backbone of each snapshot structure onto the one in the starting structure of the trajectory, using the least-squares fitting method.

M. Taghdir, H. Naderi-Manesh / Biochemical and Biophysical Research Communications 351 (2006) 1037–1042

trajectories of 2-ns simulations is displayed in Fig. 2 for all cases. All simulations showed good stability and no significant structural deviations from the native state during 2-ns. Total energy changes during 100-ps equilibration at constant pressure as function of time were also used to check simulation stability. The system was equilibrated and showed stable behavior after 20-ps in all cases (Fig. 2). The kinetic and potential energy changes in this equilibration phase also showed stable behavior in all simulations (not shown). Temperature and density also became stabilized around 300 K and 1.02 g/cm3, respectively, during this period.

The dynamic behaviors of the molecules in the MD simulations have been analyzed by using the dynamical cross-correlation matrix (DCCM) to yield information about possible correlated and anticorrelated motions. The correlation of motions between residues is indicated by the magnitude of the corresponding correlation coefficient between their Ca atoms. The cross-correlation coefficient for the displacement of each pair of Ca atoms i and j is given by: hDri  Drj i C ij ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hDr2i ihDr2j i

ð2Þ

where Dri is the displacement from the mean position of the ith atom and the Ææ represents the time average over the whole trajectory. Cross-correlation coefficients range from a value 1 (completely anticorrelated motions) to a value +1 (completely correlated motions). These coefficients do not bear any information about the magnitude of the motions, which can be small local oscillations as well as large-scale collective motions. Cross-correlation coefficients will reflect correlation displacements along a straight line. In other words, two atoms moving exactly in phase and with the same period, but in perpendicular lines, will have a cross-correlation of zero. The energy analysis module of AMBER was applied to recalculate and analyze electrostatic potential energy values for separate parts of the system. All the coordinates produced in the production phase were considered in the calculations.

Electrostatic force distributions The charged residues in the functional site of BMK M1 play a critical role in binding this molecule to sodium channel. Therefore, the mean values of electrostatic potential energy for charged residues were obtained from the trajectories of the simulations and the energy changes more than 20% compared with the native state were considered to be notable. In Table 1, we have presented mean energy difference for the charged residues in the disulfide bond broken states compared with the native state. Most changes occur in the first residue of N-terminal segment, the reverse turn, b-turn III, and the last residues of C-terminal loop in the molecule. The residues Lys8 in

Results and discussion Stability of trajectories The time evolution of the Ca RMSD values of the coordinates with respect to the experimental structure in the Native

5

4

4

3

3

2

2

1

1

0

0

5

5

C16-C36

C12-C63

C22-C46

0

RMSD values (A )

5

1039

4

4

3

3

2

2

1

1

0

0

5

C26-C48

5

4

4

3

3

2

2

1

1

0

0

400

800

1200

1600

2000

0

All bonds

0

400

800

1200

1600

2000

time (ps) ˚ ) of the structures in the trajectories over simulation time (ps) with respect to the experimental structure for the studied states. Fig. 2. RMSD values (A Total energy changes during 100-ps equilibration in constant pressure for each state are also shown (without showing time and energy scales).

1040

M. Taghdir, H. Naderi-Manesh / Biochemical and Biophysical Research Communications 351 (2006) 1037–1042

Table 1 Mean energy difference for the charged residues in the disulfide bond broken states compared with the native state Residue (Location)

C12–C63

C16–C36

C22–C46

C26–C48

All bonds

ARG2 (N-terminal) ASP3 (b-strand I) LYS8 (reverse turn) GLU15 (long loop) ARG18 (long loop) GLU20 (a-helix) ASP24 (a-helix) LYS28 (a-helix) LYS32 (strand II) LYS41 (b-turn II) GLU50 (strand III) ASP53 (b-turn III) ARG58 (C-terminal loop) LYS62 (C-terminal loop) HIS64 (C-terminal loop)

28.13 — — — — — — — — 20.89 — 45.07 17.33 29.17 52.18

22.95 — 48.40 — — — — — — — — 49.88 18.60 40.25 —

27.43 — 26.28 — — — — 32.57 — 14.06 — 54.71 — — —

45.60 38.66 42.23 — 28.95 — — — — — — — — 26.56 82.88

37.16 — 37.95 — — — — — — — — 28.28 23.10 40.56 32.20

Only the differences more than 20% are reported. The positive and the negative values denote increasing and decreasing of the energy compared with the native state, respectively. Units: kcal/mol.

the reverse turn and Arg58, Lys62, and His64 in C-terminal loop are ‘‘hot spot’’ residues in function. In Cys12-Cys63 disulfide bond broken state, the potential energy of functional residues except for Lys8 is changed considerably. By elimination of this bond, the energy of Lys8 is not changed considerably, while the potential energy of Arg58 and His64 is decreased and for Lys62 it is increased in this state. So this disulfide bond seems to play an important role in modulating force distribution in the functional segment in the C-terminal region of the molecule. In Cys16-Cys36 bond broken state, Lys8 shows a drastic decrease in potential energy in comparison with the native state, but the energy of Arg58 and Lys62 is increased. The energy of His64 is not changed considerably in this state. In Cys22Cys46 disulfide bond broken state, the energy of Lys8 is decreased, but the other functional residues show no or few changes. The presence of this bond in the molecule seems to be important for adjusting force distribution in the reverse turn, without considerable influence on this distribution in the second functional segment in C-terminal region. In Cys26-Cys48 disulfide bond broken state, Lys8 and His64 show a drastic decrease in energy, while the energy of Lys62 is increased and Arg58 shows no considerable change. This observation can show the importance of the presence of this bond for modulating force distribution in both of the functional segments. When all disulfide bonds are broken, the energy changes in functional residues compared with the native state show a similar pattern with the ones in the other disulfide bond broken states. These observations result the following: (1) The potential energy of most functional residues is changed considerably in all cases compared with the native state. Therefore, the presence of all disulfide bonds seems to be necessary for modulating energy in the functional site of the molecule. (2) The presence of Cys12-Cys63 disulfide bond seems to be important for modulating the distribution of electrostatic forces in the functional segment in

C-terminal region without considerable influence on these forces in the reverse turn, while Cys22-Cys46 modulate force distribution in the reverse turn without influence on these forces in C-terminal region. (3) The potential energy of functional residues is influenced by most disulfide bond breakages in the same way. (4) The presence of each disulfide bond in the structure seems to influence on electrostatic force distribution in global scale. These results suggest that in this toxin, disulfide bonds play a cooperative role to adjust electrostatic force distributions in the functional segments in distant regions. Structural mobility and correlation of motions Comparison of internal motions in the native state with the disulfide bond broken states shows that most changes occur in the reverse turn, b-turn II, and C-terminal loop (Fig. 3). In Cys12-Cys63 disulfide bond broken state, the flexibility of the b-turn II and C-terminal loop increases drastically, while the motions of the reverse turn are not changed in comparison with the native state. In Cys16-Cys36 disulfide bond broken state, the flexibility of the reverse turn is decreased and the other regions show the same flexibility compared with the native state. In Cys22-Cys46 disulfide bond broken state, the fluctuations of the reverse turn decrease considerably, while the b-turn II shows a minor increase in its motions compared with the native state and the flexibility of the functional segment in the C-terminal region is not changed. For Cys26-Cys48 disulfide bond broken state, the flexibility of the reverse turn and b-turn II shows a considerable decrease and increase, respectively, in comparison with the native state, while the fluctuation of the C-terminal region is not changed. Finally in all disulfide bonds broken state, the molecule shows a minor decrease in the flexibility of the reverse turn and a minor increase in its motions in the other regions.

M. Taghdir, H. Naderi-Manesh / Biochemical and Biophysical Research Communications 351 (2006) 1037–1042 4

4

Native

3

3

2

2

1

1

1

8

15

22

29

36

43

50

4

o

C12-C63

0

0

RMSF values (A )

1041

57

1

64

8

15

22

29

36

43

50

57

4

C16-C36

3

3

2

2

1

1

0

64

C22-C46

0 1

8

15

22

29

36

43

50

57

1

64

8

15

22

29

36

43

50

57

4

4

All bonds

C26-C48

3

3

2

2

1

1

0

64

0 1

8

15

22

29

36

43

50

57

64

1

8

15

22

29

36

43

50

57

64

Fig. 3. Atomic positional fluctuations of Ca atoms in disulfide bond broken states (heavy line) compared with the native state (thin line).

These findings result in the following: (1) The flexibility of the reverse turn is independent of the presence of Cys12-Cys63 disulfide bond, while by elimination of this bond the fluctuation of the second functional segment in the C-terminal region is increased drastically. (2) The fluctuation of the reverse turn shows the same pattern of changes in three other individual disulfide bond broken states. Accordingly, the disulfide bonds Cys16Cys36, Cys22-Cys46, and Cys26-Cys48 seem to have a cooperation to modulate the flexibility of this functionally important region. (3) Comparison of Cys12-Cys63 disulfide bond broken state with all disulfide bonds broken state is noticeable; the flexibility of C-terminal region in the absence of Cys12-Cys63 disulfide bond and in the presence of the other disulfide bonds is changed considerably, while when all disulfide bonds are broken, it shows a minor difference compared with the native state. This reveals the

influence of the presence of disulfide bonds on the pattern of force distribution in distant regions in the molecule. Consequently, these results suggest that disulfide bonds have a cooperation to adjust force distributions in the distant regions that are important for function in this toxin. In Table 2, we have also presented correlations between the motions of functionally important residue pairs for different states. The correlation coefficients above 0.25 and below 0.25 are considered to indicate correlated and anticorrelated motions, respectively, and the values from 0.15 to 0.25 and 0.15 to 0.25 are also considered to indicate weak correlated and weak anticorrelated motions, respectively. The residue Lys8 has anticorrelated motions with Arg58 in the native state that change to weak correlated motions in Cys12-Cys63 and Cys16-Cys36 disulfide bond broken states and change to weak anticorrelated motions in all

Table 2 Correlations of motions between the functionally important residue pairs in the functional site Residue pairs

Native state

C12–C63

C16–C36

C22–C46

C26–C48

All bonds

LYS8-ARG58 LYS8-LYS62 LYS8-HIS64 ARG58-LYS62 ARG58-HIS64

A A — WA A

WC A A A A

WC A — A A

— A WA — A

— WA WC WC —

WA A A C —

C, correlated motions; A, anticorrelated motions; WC, weak correlated motions; WA, weak anticorrelated motions.

1042

M. Taghdir, H. Naderi-Manesh / Biochemical and Biophysical Research Communications 351 (2006) 1037–1042

disulfide bonds broken state. The motions of these two residues are uncorrelated in Cys22-Cys46 and Cys26-Cys48 disulfide bond broken states. The residues Lys8 and Lys62 show anticorrelated motions in all cases except for Cys26-Cys48 disulfide bond broken state that change to weak anticorrelated motions in this state. The motions of the residues Lys8 and His64 are uncorrelated in the native and Cys16-Cys46 disulfide bond broken states. These residues show anticorrelated motions in Cys12-Cys63, Cys22Cys46 and all disulfide bonds broken states that change to weak correlated motions in Cys26-Cys48 disulfide bond broken state. The residues Arg58 and Lys62 show anticorrelated motions in the native state, Cys12-Cys63 and Cys16-Cys36 disulfide bond broken states. These residues have weak and strong correlated motions in Cys26-Cys48 and all disulfide bonds broken states, respectively, and there is no correlation between the motions of these two residues in Cys22-Cys46 disulfide bond broken state. Finally, the motions of the residues Arg58 and His64 are anticorrelated in the native state, Cys12-Cys63, Cys16-Cys36, and Cys22-Cys46 disulfide bond broken states, while these motions show no correlation in Cys26-Cys48 and all disulfide bonds broken states. These findings result: (1) in all states, the correlations of motions between most functional residues are changed compared with the native state. (2) By the elimination of Cys26-Cys48 disulfide bond, the correlations of motions between all functional residue pairs are changed in the same way, while in Cys12-Cys63 bond broken state these correlations show no or few changes for most residue pairs. This reveals that the presence of the distant disulfide bond Cys26-Cys48, which is critical in folding process, seems to play a more important role in modulating the correlations of motions between important functional residues in the functional site of this toxin. Although the alterations of the correlation of motions in this range may not have a significant effect on function, once again our results suggest that disulfide bonds play a cooperative role to adjust force distributions in distant regions, which are important for function in this toxin. Conclusion It seems that disulfide bonds have a functionally important role in adjusting force distribution pattern in native state of proteins. The pattern of this force distribution in different regions of molecule can be reflected in amplitude of fluctuations, correlation of motions, and potential energy values of residues in these regions. Our results showed that each individual disulfide bond in scorpion toxin BMK M1 seems to play a crucial role in modulating electrostatic potential energy of functional residues, the amplitude of fluctuations of the functional segments, and the correlation of motions between the functional residue pairs. On the other hand, the presence of all disulfide bonds seems to be necessary for adjusting force distribution in the functional site of the molecule. Finally, these results sug-

gest that the pattern of force distribution in the functional site of this toxin is adjusted by the cooperative effect of all disulfide bonds. References [1] T.A. Klink, K.J. Woycechowsky, K.M. Taylor, R.T. Raines, Contribution of disulfide bonds to the conformational stability and catalytic activity of ribonuclease A, Eur. J. Biochem. 267 (2000) 566–572. [2] S.H. Lin, Y. Konishi, M.E. Denton, H.A. Scheraga, Influence of an extrinsic cross-link on the folding pathway of ribonuclease A. Conformational and thermodynamic analysis of cross-linked (lysine7-lysine41)-ribonuclease A, Biochemistry 23 (1984) 5504–5512. [3] C. N Pace, G.R. Grimsley, J.A. Thomson, B.J. Barnett, Conformational stability and activity of ribonuclease T1 with zero, one, and two intact disulfide bonds, J. Biol. Chem. 263 (1988) 11820–11825. [4] A.J. Doig, D.H. Williams, Is the hydrophobic effect stabilizing or destabilizing in proteins: the contribution of disulfide bonds to protein stability, J. Mol. Biol. 217 (1991) 389–398. [5] K.A. Dill, Dominant forces in protein folding, Biochemistry 29 (1990) 7133–7155. [6] Q. Zhu et al., Role of disulfide bonds in folding and activity of leiurotoxin I: just two disulfides surface, Biochemistry 41 (2002) 11488–11494. [7] B. Rizzuti, L. Sportelli, R. Guzzi, Evidence of reduced flexibility in disulfide bridge-depleted azurin: a molecular dynamics simulation study, Biophys. Chem. 94 (2001) 107–120. [8] M.E. Moghaddam, H. Naderi-Manesh, Role of disulfide bonds in modulating internal motions of proteins to tune their function: molecular dynamics simulation of scorpion toxin Lqh III, Proteins 63 (2006) 188–196. [9] A. Warshel, J. Aqvist, Electrostatic energy and macromolecular function, Annu. Rev. Biophys. Biophys. Chem. 19 (1991) 267–298. [10] S.H. Northrup, J.O. Boles, J.C.L. Reynolds, Electrostatic effects in the brownian dynamics of association and orientation of heme proteins, J. Phys. Chem. 91 (1987) 5991–5998. [11] E. Getzoff et al., Faster superoxide dismutase mutants designed by enhancing electrostatic guidance, Nature 358 (1992) 347–351. [12] C. Vita, C. Roumestand, F. Toma, A. Menez, Scorpion toxins as natural scaffolds for protein engineering, Proc. Natl. Acad. Sci. USA 92 (1995) 6404–6408. [13] C.G. Wang et al., Exploration of the functional site of a scorpion alpha-like toxin by site-directed mutagenesis, Biochemistry 42 (2003) 4699–4708. [14] X.L. He et al., Crystal structures of two alpha-like scorpion toxins: non-proline cis peptide bonds and implications for new binding site selectivity on the sodium channel, J. Mol. Biol. 292 (1999) 125–135. [15] Y.M. Sun et al., Roles of disulfide bridges in scorpion toxin BmK M1 analyzed by mutagenesis, J. Pept. Res. 60 (2002) 247–256. [16] D.A Case, T.A. Darden, T.E. Cheatham III, C.L. Simmerling, J. Wang, R.E. Duke, R. Luo, K.M. Merz, B. Wang, D.A. Pearlman, M. Crowley, S. Brozell, V. Tsui, H. Gohlke, J. Mongan, V. Hornak, G. Cui, P. Beroza, C. Schafmeister, J.W. Caldwell, W.S. Ross, P.A. Kollman, AMBER version 8. San Francisco: University of California. 2000. [17] W.L. Jorgensen, J. Chandrasekhar, J. Madura, M.L. Klein, Comparison of simple potential functions for simulating liquid water, J. Chem. Phys. 79 (1983) 926–935. [18] H.J.C. Berendsen et al., Molecular dynamics with coupling to an external bath, J. Chem. Phys. 81 (1984) 3684–3690. [19] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes, J. Comput. Phys. 23 (1977) 327– 341. [20] A. Amadei, A.B. Linssen, H.J. Berendsen, Essential dynamics of proteins, Proteins 17 (1993) 412–425.