Protein Science Sheba protein
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

This Article
Right arrow Abstract Freely available
Right arrow Full Text (PDF)
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Right arrow Citation Map
Services
Right arrow Email this article to a friend
Right arrow Similar articles in this journal
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Likic, V. A.
Right arrow Articles by Gooley, P. R.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Likic, V. A.
Right arrow Articles by Gooley, P. R.
Social Bookmarking
Add to CiteULike   Add to Connotea   Add to Del.icio.us   Add to Digg   Add to Reddit   Add to Technorati  
What's this?
Protein Science (2003), 12:2215-2229.
Copyright © 2003 The Protein Society

Dynamics of Ca2+-saturated calmodulin D129N mutant studied by multiple molecular dynamics simulations

Vladimir A. Likic1, Emanuel E. Strehler2 and Paul R. Gooley1

1 Department of Biochemistry and Molecular Biology, Russell Grimwade School of Biochemistry, The University of Melbourne, Parkville, VIC 3052, Australia
2 Department of Biochemistry and Molecular Biology, Mayo Clinic and Foundation, Rochester, Minnesota 55905, USA

Reprint requests to: Vladimir A. Likic, Department of Biochemistry and Molecular Biology, The University of Melbourne, Parkville, VIC 3052, Australia; e-mail: vlikic{at}unimelb.edu.au; fax: 61-3-9347 7730.

(RECEIVED May 5, 2003; FINAL REVISION July 10, 2003; ACCEPTED July 11, 2003)

Article and publication are at http://www.proteinscience.org/cgi/doi/10.1110/ps.0377803.


   Abstract
TOP
 Abstract
Introduction
 Results
 Discussion
 Materials and methods
 References
 
Fifteen independent 1-nsec MD simulations of fully solvated Ca2+ saturated calmodulin (CaM) mutant D129N were performed from different initial conditions to provide a sufficient statistical basis to gauge the significance of observed dynamical properties. In all MD simulations the four Ca2+ ions remained in their binding sites, and retained a single water ligand as observed in the crystal structure. The coordination of Ca2+ ions in EF-hands I, II, and III was sevenfold. In EF-hand IV, which was perturbed by the mutation of a highly conserved Asp129, an anomalous eightfold Ca2+ coordination was observed. The Ca2+ binding loop in EF-hand II was observed to dynamically sample conformations related to the Ca2+-free form. Repeated MD simulations implicate two well-defined conformations of Ca2+ binding loop II, whereas similar effect was not observed for loops I, III, and IV. In 8 out of 15 MD simulations Ca2+ binding loop II adopted an alternative conformation in which the Thr62 >C=O group was displaced from the Ca2+ coordination by a water molecule, resulting in the Ca2+ ion ligated by two water molecules. The alternative conformation of the Ca2+ binding loop II appears related to the "closed" state involved in conformational exchange previously detected by NMR in the N-terminal domain fragment of CaM and the C-terminal domain fragment of the mutant E140Q. MD simulations suggest that conformations involved in microsecond exchange exist partially preformed on the nanosecond time scale.

Keywords: Calmodulin; CaM; dynamics; EF-hands; Ca2+ binding; multiple MD simulations


   Introduction
TOP
 Abstract
 Introduction
Results
 Discussion
 Materials and methods
 References
 
Calmodulin (CaM) is a small, acidic protein of 148 amino acid residues that acts as a principal modulator of Ca2+ signaling pathways in eukaryotic cells. CaM is a central regulatory element implicated in a number of essential biological processes, including cell proliferation, hormone secretion, transcription of certain genes, neuronal survival, and others (Cohen and Klee 1988; van Eldik and Watterson 1998). Many external stimuli increase the intracellular concentration of Ca2+ ions, which in turn, bind to CaM and cause a conformational change. Ca2+-loaded CaM is able to recognize many specific target proteins with high affinity (Crivici and Ikura 1995).

The crystal structures of Ca2+-loaded CaM show an unusual, elongated structure with two globular domains connected with a straight seven-turn {alpha}-helix (Babu et al. 1988; Wilson and Brunger 2000). Each globular domain is composed of four {alpha}-helices arranged in two EF-hand structural motifs (Ikura 1996). The EF-hand structural motif (also known as the "helix-loop-helix" motif) consists of two {alpha}-helices and the loop that contains the Ca2+ binding site (Strynadka and James 1989; Falke et al. 1994; Ikura 1996). In the Ca2+ free form, the two helices are nearly antiparallel, adopting a compact, "closed" conformation (Falke et al. 1994; Ikura 1996). Upon Ca2+ binding the two helices become nearly perpendicular to one another, exposing a patch of hydrophobic residues involved in target recognition and binding (Crivici and Ikura 1995; Ikura 1996). Ca2+ binding is therefore of central importance to biological function of CaM and other EF-hand proteins. Currently characterized EF-hand sites exhibit Ca2+ binding affinities that span six orders of magnitude (Falke et al. 1994). EF-hands perform their function within complex and interleaving Ca2+ signaling pathways, and Ca2+ binding properties of various EF-hands appear tuned to match their biological roles (Falke et al. 1994; Linse and Forsén 1995). Despite extensive physico-chemical, mutational, and structural studies, the detailed molecular determinants of Ca2+ binding properties in EF-hands remain elusive (Forsén et al. 1991; Falke et al. 1994; Linse and Forsén 1995).

It is now well understood that a pair of EF-hands forms a conserved and highly cooperative Ca2+ binding unit (Ikura 1996; Nelson et al. 2002), and asymmetry in dynamical properties of the two EF-hands has been observed in calbindin D9k (Akke et al. 1993; Mäler et al. 2000), troponin C (Li et al. 1997; Gagné et al. 1998), and CaM (Evenäs et al. 1998, 2001; Malmendal et al. 1998). NMR studies of CaM domain fragments reveal complex dynamics on picosecond and microsecond time scales, with main substates exchanging on the microsecond time scale being similar to open and closed conformations (Evenäs et al. 1997, 1999, 2001; Malmendal et al. 1998, 1999). Extensive anisotropic and anharmonic disorder was described in the crystal structure of Ca2+-CaM refined to 1.0 Å, demonstrating that the protein occupies a number of conformational substates even at cryogenic temperatures (Wilson and Brunger 2000). The recently determined structure of Ca2+-saturated CaM shows a less open conformation of the N-terminal domain in solution compared to the crystal structures, suggesting considerable intradomain plasticity in CaM (Chou et al. 2001).

MD simulations are based on empirical force fields, and provide a detailed, atomic-level model of protein dynamics (Karplus and McCammon 2002). Due to limitations in the available computational power, the time scales achievable in MD simulations of solvated proteins are currently limited to a few nanoseconds. Protein equilibrium motions have characteristic times spanning more than 10 orders of magnitude, from picoseconds to seconds or longer, and many slower motions are not properly sampled in typical MD simulations. In contrast, most water molecules within the protein solvation shell move on the subnanosecond time scale, and therefore, MD simulations are particularly useful for providing new insights into protein–water interactions (Brunne et al. 1993; Steinbach and Brooks 1993; Phillips and Pettitt 1995).

In previously reported MD simulations of solvated CaM both global structural changes and dynamics of solvating water molecules were studied (Wriggers et al. 1998; Vigil et al. 2001; Yang et al. 2001; Komeiji et al. 2002). However, a detailed analysis of the dynamics of Ca2+ binding sites has not yet been presented. In this work we report 15 independent 1-nsec MD simulations of fully solvated Ca2+-saturated CaM mutant D129N, with the focus on the dynamics and solvation of Ca2+ binding sites. MD simulations near room temperature are not ergodic, and independent MD runs yield different results (Straub et al. 1994; Auffinger et al. 1995; Caves et al. 1998). We rely on multiple independent MD runs to probe the reliability of predictions obtained from MD simulations, and also to improve the sampling of the conformational space (Elofsson and Nilsson 1993; Straub et al. 1994; Caves et al. 1998). Although all MD simulations were initialized from the fully open, crystal structure CaM conformation (Babu et al. 1988; Wilson and Brunger 2000), on the examined time scale EF-hands were observed to adopt more closed conformations consistent with the recently determined CaM structure in solution (Chou et al. 2001). The Ca2+ binding site in EF-hand II exhibited anomalous dynamics in MD simulations, related to the tendency of Ca2+ binding loop II to adopt an alternative conformation. In the alternative conformation the Thr62 >C=O group is displaced from the Ca2+ coordination sphere and the Ca2+ ion is ligated by two water molecules, which is highly reminiscent of, but different from the proposed "closed" state involved in microsecond conformational exchange detected by NMR in CaM N- and C-terminal domain fragments (Evenäs et al. 1997, 1998, 2001; Malmendal et al. 1998). The molecular mechanisms that led to the observed effects, their relationship to the experimental data, and the consequences of this observation for our understanding of dynamics and function of EF-hand Ca2+ binding sites, are discussed.


   Results
TOP
 Abstract
 Introduction
 Results
Discussion
 Materials and methods
 References
 
Global dynamics
CaM structure has an unusual global flexibility that stems from its elongated, dumbbell-like shape (Fig. 1AGo). In our MD simulations instantaneous backbone RMS deviations of up to 7.4 Å relative to the initial, crystal structure were observed (residues 5–147). However, the two globular domains remained structurally stable in all MD simulations (Table 1Go), and all four Ca2+ ions remained in their EF-hand binding sites. Figure 2Go shows the time dependence of RMS deviations in simulations that exhibited the largest and the smallest average RMS deviations, and Table 1Go details average RMS deviations for all simulations.



View larger version (24K):
[in this window]
[in a new window]
  		Figure 1. Schematic diagram of instantaneous CaM structures: (A) the crystal structure (Babu et al. 1988) used as the initial structure in all MD simulations; (B) the final structure (1 nsec) from simulation i; (C) the final structure (1 nsec) from simulation e. Ca2+ ions are shown as black spheres, and the structures from simulations i and e are best-fit superimposed to the N-terminal domain (backbone atoms of residues 5–74) of the crystal structure. The simulations i and e exhibited the global minimum and maximum in calcium virtual dihedral angle (see Fig. 3Go). Molecular graphics created with MOLSCRIPT (Kraulis 1991).

 

View this table:
[in this window]
[in a new window]
  		Table 1. Summary of parameters for multiple MD simulations of calmodulin
 


View larger version (24K):
[in this window]
[in a new window]
  		Figure 2. The time dependence of backbone RMS deviations versus the crystal structure (atoms N, C{alpha}, C') for the entire backbone (residues 5–147), the N-terminal domain (residues 5–74), and the C-terminal domain (residues 78–147). In each panel the simulations with the largest and the smallest average RMS deviations are shown (see Table 1Go). RMS deviations are shown as 2-psec averages, and data sets are labeled with the one-letter simulation code as given in Table 1Go.

 
In the crystal structure (Babu et al. 1988), the two globular lobes appear connected by a continuous seven turn {alpha}-helix (residues 65 to 92), which is both the exiting helix of EF-hand II (helix D, residues 65–75) and the entering helix of EF-hand III (helix E, residues 83–92). The central helix is slightly bent in the crystal structure (Babu et al. 1988), with the interhelical angle between the helices D and E close to 22°. Analysis of RMS deviations and interhelical angles showed that in several MD runs the central helix was interrupted within the linker region (i.e., residues 76–82). For example, in simulations d and g the D–E interhelical angle reached 63.3° and 61.8° near 829 psec and 988 psec, respectively. In simulation c, a transient unwinding of the central helix was observed: at 472 psec the structure of central helix was interrupted near residues Asp80–Glu82, and the D–E interhelical angle increased to 49.4°. After approximately 150 psec the central {alpha}-helix fully reformed and the D–E interhelical angle returned to ~20°, close to the value observed in the crystal structure.

Small-angle scattering experiments show that Ca2+ saturated CaM is more compact in solution (radius of gyration, Rg = 21.3 ± 0.2 or 21.5 ± 0.2 Å) compared to the extended conformation (Rg = 22.0 Å) found in the crystal structure (Heidorn and Trewhella 1988). Our MD simulations resulted in an overall average Rg = 21.6 ± 0.5 Å (Table 2Go; Fig. 3Go). A dihedral angle formed by four Ca2+ ions is often used to characterize relative orientation of the two CaM domains in MD studies (Pascual-Ahuir et al. 1991; Wriggers et al. 1998; Yang et al. 2001; Komeiji et al. 2002). In two MD simulations in which the electrostatic interactions were treated with a cutoff, this angle was observed to drift toward more negative values (Wriggers et al. 1998; Yang et al. 2001), while in two simulations performed with Ewald summation this angle fluctuated near the value observed in the crystal structure (Yang et al. 2001; Komeiji et al. 2002). Our MD simulations employed Ewald summation, and although in individual MD runs the Ca2+ dihedral angle was observed to drift in either direction (Table 2Go; Fig. 3Go), the overall behavior is consistent with fluctuations near the value observed in the crystal structure (Table 2Go).


View this table:
[in this window]
[in a new window]
  		Table 2. The average radius of gyration (Rg) and calcium virtual dihedral angle for 15 MD simulations b–p
 


View larger version (32K):
[in this window]
[in a new window]
  		Figure 3. Radius of gyration (upper panel) and calcium virtual dihedral angle (lower panel) for simulations that exhibited extreme values over all MD simulations. Dashed lines represent values observed in the X-ray structure (Babu et al. 1988).

 
Dynamics of N- and C-terminal domains
In all MD runs the N- and C-terminal domains remained structurally stable with average RMS deviations between 1.00 and 2.52 Å (Table 1Go; Fig. 2Go), which are values typically observed in MD simulations of globular proteins. The N-terminal domain showed consistently higher RMS deviations versus the crystal structure than the C-terminal domain (Table 1Go). Similar results were observed in previous MD simulations of Ca2+-loaded CaM: a 3-nsec MD simulation in a solvent sphere (Wriggers et al. 1998) and two 4-nsec periodic boundary MD simulations with full electrostatics (Yang et al. 2001; Komeiji et al. 2002). Because previously reported MD simulations entailed a single MD run, it was difficult to gauge the importance of this observation. Out of 15 independent MD runs reported here only in one the average RMS deviation was higher for the C-terminal domain compared to the N-terminal domain (simulation g, see Table 1Go), providing evidence that the above observation is statistically significant. This implies that in our MD simulations the N-terminal domain exhibited a tendency to drift further away from the initial, crystal structure compared to the C-terminal domain.

Each globular domain of CaM is composed of two EF-hands, which upon binding of Ca2+ change from a "closed" (the two helices are nearly antiparallel) to an "open" conformation (the two helices are nearly perpendicular), ready for target binding (Ikura 1996). Surprisingly, the structure of Ca2+-saturated CaM determined from NMR residue dipolar couplings shows that all four EF-hands adopt a more closed conformation in solution compared to the crystal structure (change of 10–26° in interhelical angles) (Chou et al. 2001). In our MD simulations EF-hands consistently adopted more closed conformations, similar to those determined for the structure in solution (Fig. 4Go). For the C-terminal domain the agreement between the interhelical angles calculated from MD simulations and the solution structure is striking (see below). For the N-terminal domain, the average interhelical angles calculated from MD simulations were mostly between the values for the initial, crystal structure (Babu et al. 1988) and the solution structure (Chou et al. 2001). The change in interhelical angles between the crystal and solution structures is significantly greater in the N-terminal domain than for the C-terminal domain (shown by arrows in Fig. 4Go). This implies that the N-terminal domain has to depart more from its crystal structure conformation to achieve the solution conformation, and provides an explanation for higher RMS deviations observed for the N-terminal domain in our MD simulations. Furthermore, Figure 4Go suggests that 1 nsec is not a sufficiently long time for N-terminal EF-hands to relax entirely into their solution conformation. For EF-hand III located in the C-terminal domain the overall average interhelical angle calculated from MD simulations is remarkably close to the value observed in the solution structure (solution structure: 73.2°, MD simulations: 73.7°). For EF-hand IV, MD simulations predict an even more closed conformation (an overall average of 71.9°) than determined in solution (78.2°). It is questionable whether this difference (6.3°) is significant, but it should be noted that the mutation D129N used in MD simulations perturbs the Ca2+ binding site in EF-hand IV, and this may also affect the corresponding interhelical angle.



View larger version (22K):
[in this window]
[in a new window]
  		Figure 4. Average interhelical angles calculated from 15 MD simulations for the N-terminal domain (upper panel) and the C-terminal domain (lower panel). The values from MD simulations are shown as points, with triangles for the first EF-hand (I and III) and asterisks for the second EF-hand (II and IV) in each domain. The solid vertical lines show interhelical angles observed in the crystal structure of Ca2+-CaM (Babu et al. 1988), while the dashed vertical lines show interhelical angles observed in the solution structure of Ca2+-CaM (Chou et al. 2001). The arrows depict the change from the crystal to solution state.

 
Dynamics of EF-hand Ca2+ binding sites
Ca2+ binding loops in canonical EF-hands consist of 12 residues consecutively numbered 1–12, and the four Ca2+ binding loops in CaM are usually labeled I–IV according to the EF-hand they belong to (Strynadka and James 1989; Falke et al. 1994). The backbone RMS deviations versus the crystal structure of individual EF-hands, {alpha}-helices, and 12 residue Ca2+ binding loops observed in our MD simulations are summarized in Table 3Go. Both EF-hands from the N-terminal domain (I and II) exhibited larger RMS deviations than EF-hands from the C-terminal domain (III and IV), which is consistent with the above analysis of RMS deviations of individual domains and interhelical angles. Overall, the Ca2+ binding loop II exhibited significantly larger RMS deviations compared to the other three Ca2+ binding loops. Apart from terminal helices (A and H), {alpha}-helix C (EF-hand II) exhibited slightly higher RMS deviations.


View this table:
[in this window]
[in a new window]
  		Table 3. Backbone RMS deviations from 15 MD simulations b–p of individual EF-hand, {alpha}-helices, and Ca2+ binding loops versus the initial X-ray structure
 
For Ca2+ binding loops I, III, and IV the observed RMS deviations versus the crystal structure were steady and without pronounced extremes (Fig. 5Go). Average RMS deviations over 15 MD runs were 0.51 ± 0.07 Å, 0.72 ± 0.08 Å, 0.50 ± 0.04 Å, and 0.56 ± 0.03 Å, for the four loops, respectively. The largest overall average observed for the Ca2+ binding loop II originated from several discrete events, such as an abrupt increase in RMS deviations observed in simulation p near 700 psec (see Fig. 5Go). A conformational exchange that involves conformations similar to open and closed states was observed in Ca2+ saturated CaM C-terminal domain mutants E140Q and E104Q (Evenäs et al. 1997, 1998, 1999, 2001), in the Mg2+-saturated N-terminal domain of CaM (Malmendal et al. 1998), and in the C-terminal domain of apo CaM (Malmendal et al. 1999). Therefore, it is of interest to examine if in our MD simulations Ca2+ binding loops showed any tendency to approach the structure of the Ca2+ free form. We examined this by using the difference between the RMS deviations calculated with respect to the Ca2+-free structure (Kuboniwa et al. 1995) and Ca2+-loaded structure (Babu et al. 1988) ({Delta}RMSD = RMSapo - RMSholo), which shows whether a given conformation is closer to the Ca2+-loaded or Ca2+-free structure in relative terms. Because MD simulations were started from the Ca2+-loaded form (RMSholo = 0), for the initial structure {Delta}RMSD = RMSapo, which is 1.72, 1.69, 1.77, and 2.34 Å for Ca2+ binding loops I, II, III, and IV, respectively. In our MD simulations Ca2+ binding loops settled near RMS deviations of 0.5–1.0 Å relative to the initial, Ca2+-loaded structure (Fig. 5Go). Therefore, if the conformation remained at the same "distance" from the Ca2+-free structure (RMSD = 1.7–2.3 Å) we expect {Delta}RMSD to fluctuate in the range of 1–2 Å. This appears to be the case for Ca2+ binding loops I, III, and IV: the global {Delta}RMSD minimum recorded over 15 MD runs was 0.74, 0.91, and 1.48 Å, respectively. On the other hand, the Ca2+ binding loop II sampled conformations with {Delta}RMSD <0.5 Å in five MD simulations (specifically, in simulations d, e, j, o, and p). In simulation p, negative {Delta}RMSD was observed near 830 psec (Fig. 6Go), implying that the Ca2+ binding loop sampled conformations more similar to the Ca2+-free form than to the initial Ca2+-loaded structure. A brief excursion of {Delta}RMSD into negative values in simulation p was the consequence of both a sharp increase in RMSholo just prior to 700 psec (shown in Fig. 5Go) and a pronounced minimum in RMSapo observed near 833.2 psec (RMSapo = 1.14 Å). This was correlated with the penetration of the second water ligand into the EF-hand II Ca2+ coordination sphere (see below). The time dependence of RMSapo and {Delta}RMSD for simulation p are shown in Figure 6Go.



View larger version (33K):
[in this window]
[in a new window]
  		Figure 5. Backbone RMS deviations versus the crystal structure of 12-residue Ca2+ binding loops that exhibited global extremes over 15 MD simulations. For each Ca2+ binding loop I–IV the simulations that exhibited a global minimum (dashed line) and a global maximum (solid line) are shown. Shown are 2-psec averages, and data sets are labeled with the one-letter code as given in Table 1Go.

 


View larger version (26K):
[in this window]
[in a new window]
  		Figure 6. Time dependence of backbone RMS deviations in simulation p calculated with the Ca2+-free CaM structure as a reference (RMSapo, upper panel), and differences in backbone RMS deviations ({Delta}RMSD = RMSapo - RMSholo) for the 12-residue Ca2+ binding loop in EF-hand II as a function of the simulation time (simulations p and o). The average coordinates of 25 structures of apo CaM determined by NMR (Kuboniwa et al. 1995) were used as the reference for the Ca2+-free form, while the crystal structure of holo CaM (Babu et al. 1988) was used as the reference for the Ca2+-loaded form. For clarity, 2-psec averages of instantaneous RMS deviations are shown. In simulation p near 830 psec the structure of the Ca2+ binding loop in EF-hand II became more similar to the Ca2+-free form than to the Ca2+-loaded form (negative {Delta}RMSD), observed concomitantly with a pronounced minimum in RMSapo. Simulation o exhibited a continuous decrease in {Delta}RMSD.

 
Average backbone mean-square fluctuations of 12-residue segments that make the four Ca2+ binding sites in the CaM mutant D129N are shown in Figure 7Go, upper panel. All four Ca2+ binding loops showed a sharp maximum in mean-square fluctuations near the residue in position 4. This residue is invariably glycine, which facilitates a sharp turn required for the proper geometry of the Ca2+ binding sites (Strynadka and James 1989). Another highly conserved glycine occurs in position 6, but this residue showed lower mobility compared to that in position 4. The most interesting feature of the data shown in Figure 7Go is that the backbone of the EF-hand II binding loop consistently exhibited higher mobility than the other three Ca2+ binding loops. The mean-square fluctuations of the four Ca2+ ions (Fig. 7Go, lower panel) are consistent with this picture. The statistics compiled from all 15 MD simulations show unambiguously that the Ca2+ binding site from EF-hand II exhibited higher mobility compared to the binding sites in EF-hand I, III, and IV.



View larger version (24K):
[in this window]
[in a new window]
  		Figure 7. Average mean-square fluctuations of the four Ca2+ binding loops (upper panel), and mean-square fluctuations of the four Ca2+ ions (lower panel). Calculated values for four 12-residue Ca2+ binding loops are depicted in the upper panel. For each residue the mean-square fluctuations were calculated an average over backbone N, C{alpha}, C' atoms, and these values were subsequently averaged over all MD simulations. The solid lines represent Ca2+ binding loops from the N-terminal domain (squares—Ca(I) site, triangles—Ca(II) site), while the dashed lines represent Ca2+ binding loops from the C-terminal domain (vertical cross—Ca(III) site, x cross—Ca(IV) site). The lower panel shows mean-square fluctuations of the four Ca2+ ions, and for each ion the values from 15 MD simulations (b–p) are shown explicitly. The labels for the four Ca2+ ions are consistent with the labels used for Ca2+ binding loops on the upper panel.

 
EF-hand binding sites: Ca2+ ligands
An overall insight into Ca2+ coordination for the four EF-hands was obtained by analyzing the ligand occupancy for each Ca2+ binding site, averaged over all trajectory coordinate frames and multiple MD runs (Fig. 8AGo). It appears that for each Ca2+ binding site the ligands observed in the CaM crystal structure persisted throughout all MD simulations. This holds both for ligands donated by the protein (leftmost bars in each panel in Fig. 8Go) and for Ca2+ ligated water molecules.




View larger version (49K):
[in this window]
[in a new window]
  		Figure 8. Cumulative Ca2+ ligand occupancy over 15 MD simulations. The plate (A) shows the N-terminal Ca2+ binding sites Ca(I) and Ca(II) (upper and lower panels), while the plate (B) shows C-terminal Ca2+ binding sites Ca(III) and Ca(IV) (upper and lower panels, respectively). In each panel, the six leftmost bars represent six ligands found in the CaM crystal structure (Babu et al. 1988). The bars in the middle section shown with dotted line represent extra ligands donated by the polypeptide chain. The bars on the right shown in dashed lines represent water ligands. For water ligands two bars are shown: One shows the occupancy of the crystal structure water molecule, and the other represents the occupancy of all other water molecules combined. This figure presents cumulative results: an atom that was persistently a Ca2+ ligand in all simulations would have the occupancy of exactly one; an atom that was a Ca2+ ligand all the time but in only one of the simulations would have the occupancy of 1/15, or 0.06. The individual ligands and contributing residues are denoted on the x-axis. Aspartate carboxylate O{delta} oxygens, as well as asparagine amide group O{delta} oxygens are denoted OD, while glutamate carboxylate O{varepsilon} oxygens are denoted OE.

 
The most striking feature of the data shown in Figure 8Go is the tendency of the Ca2+ ion from EF-hand II to acquire an extra water ligand. The contribution of an extra water molecule to the Ca2+ coordination appears concomitant with a decrease in contribution of the Thr62 main chain >C=O group. An increase in contribution of the Asp56 and Asp58 side chain carboxylates is also visible. In contrast, Ca2+ binding site I remained essentially as observed in the crystal structure in all MD simulations. In the C-terminal domain, side-chain carboxylates of Asp93 and Asp95 of EF-hand III, which are homologous to Asp56 and Asp58, showed a weak tendency to participate as Ca2+ ligands. In EF-hand IV the side chain of Asp131 (homologous to Asp58 from the site II) was a fully occupied eighth ligand across all 15 MD simulations (occupancy ~1). In our MD simulations Asp129, the first residue in the Ca2+ binding loop of EF-hand IV, was mutated to Asn. Asp in this position is absolutely conserved in all known EF-hands (Strynadka and James 1989; Falke et al. 1994), and this mutation is presumably responsible for the observed anomalous eightfold coordination. None of the Ca2+ ions from EF-hands I, III, and IV showed any propensity to bind a second water ligand.

Dynamics of EF-hand II Ca2+ binding site in individual simulations
Figure 8AGo shows that the presence of a second water ligand for the Ca2+ ion from EF-hand II is at a level of ~40%, implying that this effect was observed in more than five out of 15 MD simulations. Analysis of individual trajectories shows that the Ca2+ ion from EF-hand II acquired a second water ligand in precisely eight MD runs: d, e, f, g, i, l, n, and p. In the remaining simulations the Ca2+ ion from EF-hand II was coordinated by exactly one water molecule throughout the entire trajectory, similarly as observed for EF-hands I, III, and IV. However, in several simulations a transient displacement of the Thr62 >C=O group from the Ca2+ coordination sphere was observed, presumably lacking a properly oriented water molecule that could replace this group in Ca2+ coordination. For example, such processes were observed in simulation j, which shows conspicuously high Ca2+ mean-square fluctuations in Figure 7Go. In simulation j two unsuccessful attempts of the Thr62 >C=O group to detach from the Ca2+ ion were observed (as evidenced by a transient increase in the >C=O to Ca2+ distance to over 4.7 Å), each corresponding to a structural perturbation lasting over a hundred picoseconds.

Substantial variations in the dynamics of the protein–water interface were observed in different MD simulations, reflecting microscopic heterogeneity previously noted in MD simulations (Auffinger and Westhof 1996). This heterogeneity can be visualized by following the distance between the Thr62 carbonyl oxygen and the Ca2+ ion as a function of the simulation time, as shown in Figure 9Go. In some MD runs the displacement of the Thr62 >C=O group was a single well-defined event, which resulted in two water molecules stably bound as Ca2+ ligands (e.g., simulations f and p, see Fig. 9Go). In simulation l, the Thr62 >C=O group was temporarily displaced from the Ca2+ coordination sphere by the water molecule, which persisted as the Ca2+ ligand for about 300 psec. On the other hand, in simulation e, 10 different water molecules were observed to coordinate the Ca2+ ion as the second water ligand within 1 nsec (this type of behavior is illustrated in Fig. 9Go, simulation n). It seems that these water molecules were not able to find a stable configuration for the Ca2+ coordination: they coordinated the Ca2+ ion one at a time, with coordination times ranging from several picoseconds to several hundred picoseconds, not affecting the "permanently" bound water molecule observed in the crystal structure.



View larger version (31K):
[in this window]
[in a new window]
  		Figure 9. The distances between the carbonyl oxygen of residue Thr62 >C=O and the EF-hand II Ca2+ ion (solid line), and between the Ca2+ ion and the oxygen of extra water ligands (dots) as observed in simulations f, l, n, and p. In simulations f, l, and p one extra water ligand was observed during 1 nsec of production dynamics. In simulation n seven different extra water ligands were observed to coordinate the Ca2+ ion, with coordination times between 26 and 336 psec. The bottom panel shows trajectories of two water molecules with the longest coordination times (336 psec and 152 psec).

 
The molecular graphics showing the penetration of the second water ligand into the coordination sphere of the Ca2+ ion from EF-hand II in simulation p is shown in Figure 10Go. The changes in the distance between the Ca2+ ion and the Thr62 >C=O group (also of the incoming water ligand) are shown Figure 9Go. As noted above, the previously analyzed RMS deviations of the Ca2+ binding loop exhibited an abrupt increase near 700–800 psec that clearly correlated with the exchange of Ca2+ ligands, as shown in Figures 5Go and 6Go.



View larger version (37K):
[in this window]
[in a new window]
  		Figure 10. Trajectory of a water molecule in simulation p, which penetrated the EF-hand site II Ca2+ coordination sphere near 700 psec to become a second water ligand (see Figs. 5Go and 9Go). The structure of the Ca2+ binding loop and second water ligand positions are shown at times 652.0 psec (labeled 1), 683.0 psec (2), and 700.0 psec (3). The Ca2+ ion is gray, and the crystal structure water molecule is colored yellow. Trajectory of the second water ligand over a time period of 48.0 psec is traced with a blue line in steps of 0.2 psec. The structure of the Ca2+ binding loop at 700.0 psec is shown as a solid thick line. The position of Thr62 is labeled, and the arrow shows the displacement of Thr62 >C=O group from the Ca2+ coordination sphere. In this particular simulation the second water ligand displaced the crystal structure water molecule from its normal position. In its new position, the crystal structure water molecule appears inserted between the Ca2+ ion and the Thr62 >C=O group, and is engaged in H-bonding interactions with the Thr62 carbonyl oxygen. Molecular graphics created with MOLSCRIPT (Kraulis 1991).

 
Interactions of second water ligands with the protein
Because MD simulations provide trajectories of individual water molecules, they provide insight into questions that are often beyond the reach of conventional experiments. The following two questions are of particular interest: (1) what are the interactions that facilitate the approach of water ligands to the Ca2+ binding site; (2) what are the stabilizing interactions for the extra water ligand once it is within the Ca2+ coordination sphere? In relation to question (1), we were interested primarily in water molecules that coordinated the Ca2+ ion with reasonable stability. To separate secondary water ligands that were bound "stably" within the Ca2+ coordination sphere we chose a cutoff of 150 psec as a minimum coordination time. Although this cutoff is arbitrary, it is guided by the property of bulk phase water molecules to undergo a large rotational motion every 10 psec to find new stable H-bonding partners (Ohmine and Tanaka 1993). Thus, the cutoff of 150 psec is sufficiently large to eliminate water molecules that are not substantially stabilized within the Ca2+ coordination sphere.

In our MD simulations nine water ligands satisfied the 150 psec criterion (the Ca2+ coordination time is given in parentheses): e (229 psec), f (775 psec), i (171 psec), l (300 psec), p (305 psec). Furthermore, in each of the simulations g and n two such water molecules were found, with coordination times of 212 psec and 151 psec (simulation g), and 336 psec and 152 psec (simulation n). Trajectories of these water molecules were analyzed for the interactions with protein residues near the coordination event. Large variations in dynamic behavior reflecting microscopic heterogeneity were observed between individual trajectories. To extract some generalized conclusions we relied on a qualitative picture obtained from H-bonding interactions. For water molecules which became a stably bound second Ca2+ ligand the analysis of H-bonding patterns was initiated at the time point when the water oxygens were ~8–10 Å away from the Ca2+ ion (typically observed within 20–60 psec prior to the coordination), and was extended until the time of Ca2+ coordination. A count of H-bonding partners suggested that four protein groups (carboxylates of Asp58, Asp64, and Glu67, and the side chain amide group of Asn60) consistently "guided" the approach of water molecules which successfully penetrated the Ca2+ coordination sphere in EF-hand II.

The "stably" Ca2+-bound second water ligands were on average engaged in 1.99 H-bonds when bound to the Ca2+ ion: 1.65 hydrogen bonds were with other water molecules, and 0.34 hydrogen bonds were with protein groups. Only three protein groups were implicated in the stabilization of the second water ligand: side chain carboxylates of Asp58, Asp64, and Glu67 (H-bond scores of 0.26, 0.04, and 0.03, respectively). These groups were also involved in modulating the approach of extra water ligands to the Ca2+ binding site in EF-hand II, thus suggesting a particular spatial arrangement of negatively charged residues which modulates this event.


   Discussion
TOP
 Abstract
 Introduction
 Results
 Discussion
Materials and methods
 References
 
Global dynamics
Fifteen independent MD simulations of fully solvated Ca2+-saturated CaM mutant D129N were analyzed. Extensive global flexibility facilitated by the central {alpha}-helix was observed, in agreement with small-angle X-ray scattering experiments (Heidorn and Trewhella 1988), NMR relaxation studies (Barbato et al. 1992; Lee et al. 2000), and previous MD simulations (Wriggers et al. 1998; Vigil et al. 2001; Yang et al. 2001; Komeiji et al. 2002). The radius of gyration of CaM obtained from MD simulations was Rg = 21.6 ± 0.5 Å, in close agreement with the value obtained from small-angle X-ray scattering (SAXS) experiments (Rg = 21.3 ± 0.2 or 21.5 ± 0.2 Å ) (Heidorn and Trewhella 1988). In contrast to previously reported MD simulations (Wriggers et al. 1998; Komeiji et al. 2002), we calculated Rg as an average from 15 independent MD simulations, and the reported standard deviation reflects the dispersion of individual MD runs around the mean rather than fluctuations in any single trajectory. This approach is validated by observed large variations in Rg between individual MD runs (Table 2Go; Fig. 3Go), which is consistent with the notion that each MD run sampled different localized regions of the conformational space (Caves et al. 1998).

Internal dynamics of N- and C-terminal domains
Individual globular domains of CaM were structurally stable (Table 1Go), and all four Ca2+ ions remained in their EF-hand binding sites in all MD simulations. Recently, the structure of Ca2+-saturated CaM was determined in a liquid crystalline medium from NMR residual dipolar couplings (Chou et al. 2001). Surprisingly, the solution structure shows that the EF-hands, in particular EF-hand I and II located in the N-terminal domain, are in a more closed conformation in solution compared to the crystal structure. Our MD simulations were initiated from the fully open crystal structure conformation (Babu et al. 1988), and a strong tendency of all four EF-hands to adopt more closed conformations, similar to the solution structure (Chou et al. 2001), was observed (Fig. 4Go). The existence of more closed EF-hand conformations in solution was hypothesized even before the solution structure was available based on MD simulations carried out with a different force field (Vigil et al. 2001). Two additional points emerge from our analysis: (1) 1 nsec MD runs appear insufficiently long for the N-terminal EF-hands to relax into the solution conformation; (2) EF-hand interhelical angles calculated from different 1-nsec MD runs exhibit large variations. The latter point implies that a single nanosecond MD trajectory does not afford a reliable estimate of the average interhelical angles as predicted by the computational model. In agreement with previous MD simulations of CaM, we observed higher RMS deviations of the N-terminal domain relative to the C-terminal domain (Wriggers et al. 1998; Yang et al. 2001; Komeiji et al. 2002). This effect is explained by the greater conformational strain of the N-terminal domain EF-hands in the crystal structure conformation (Babu et al. 1988; Wilson and Brunger 2000) relative to the solution structure (Chou et al. 2001).

Dynamics of Ca2+ binding sites
In the CaM crystal structures there is precisely one water molecule that completes Ca2+ coordination in each of the Ca2+ binding sites (Babu et al. 1988; Wilson and Brunger 2000). In each of the 15 MD simulations these water molecules persisted as Ca2+ ligands throughout the entire simulated period of 1 nsec (true for all four EF-hands, Fig. 8Go), suggesting exceptional resilience of "stably" bound water molecules as Ca2+ ligands in EF-hands. In another EF-hand protein, calbindin D9k, the residence times of two calcium-coordinated water molecules were estimated to be between 5 nsec and 7 µsec by nuclear magnetic resonance dispersion (Denisov and Halle 1995).

Increased mobility of Ca2+ binding loop II relative to loops I, III, and IV was deduced from the comparison of backbone mean-square fluctuations, as well as mean-square fluctuations of bound Ca2+ ions. Backbone RMS deviations also indicated larger structural deviations of the loop II relative to the initial crystal structure. These effects were linked with the conformational change that involved a displacement of the loop II backbone carbonyl group in position 7 (Thr62 >C=O) from the Ca2+ coordination sphere by a water molecule. A similar effect was not observed in Ca2+ binding loops I, III, and IV in any of the simulations. The average distances between the Ca2+ ions and their ligands in two 4-nsec MD simulations of Ca2+-saturated CaM were reported previously (Yang et al. 2001). Interestingly, the only deviation from the crystal structure was the interruption of the interaction between the Thr62 >C=O group and the Ca2+ ion in one of the simulations, which was attributed to the incomplete treatment of electrostatic interactions (Yang et al. 2001).

The conformational "exchange" suggested by our MD simulations occurs between the conformation observed in the crystal structure (Babu et al. 1988) (conformation A), and a conformation with the Thr62 >C=O group displaced from the Ca2+ coordination by a water molecule thus resulting in the Ca2+ ion ligated with two water molecules (conformation B). The exchange between the two conformations occurs on a time scale that is significantly longer than the duration of reported MD runs, and only the initial stages of this process can be sampled within 1 nsec. The exact molecular mechanism of the Thr62 >C=O displacement by water molecules, as well as the observed dynamics, were different in each simulation. In simulations d, e, f, g, i, l, n, and p the Thr62 >C=O ligand was displaced by a water molecule (A->B barrier crossing). In simulation l, a reversible one-step displacement A->B->A was observed, with the conformation B lasting for about 200 psec (Fig. 9Go). In three other simulations (j, k, and o) a transient and reversible detachment of the Thr62 >C=O group from the Ca2+ ion was observed without the penetration of a water molecule into the Ca2+ coordination sphere; these represent unsuccessful attempts of crossing the potential energy barrier (A->B). This microscopic heterogeneity is expected to be observed if individual protein molecules are examined for short periods of time such as 1 nsec. From the standpoint of a single MD simulation the observed effect essentially represents a "rare" event. Because MD simulations are stochastic in nature, little weight can be attached to a single such observation. However, we observed this effect in eight out of 15 MD simulations, and therefore, it is unambiguously predicted by the employed computational model.

The observed conformational change is highly reminiscent of exchange processes detected in Mg2+-saturated N-terminal domain (Tr1C) of CaM (Malmendal et al. 1998), and also in Ca2+-saturated C-terminal domain (Tr2C) mutant E140Q (Evenäs et al. 1997, 1999, 2001). When Tr1C is titrated with Mg2+, a continuous broadening of backbone amide NMR signals for several residues located in the putative Ca2+ binding loop is observed (Malmendal et al. 1998). Furthermore, when both EF-hands I and II are saturated with Mg2+, chemical shifts of many resonances from loop II appear half way between the ion-free and Ca2+-loaded states (Malmendal et al. 1998). These observations suggest a conformational exchange, and it was hypothesized that loop II undergoes exchange between "lower affinity" and "higher affinity" conformations on the time scale of 100 µsec. Furthermore, it was shown previously that Ca2+ coordination by the backbone carbonyl group in position 7 causes polarization of the amido group —C(7)O—N(8)—, causing deshielding of N(8), which is a reliable indicator of metal binding by EF-hands (Biekofsky et al. 1998). Mg2+ binding to the N-terminal domain of CaM causes a change in 15N chemical shift of Ile63 (residue in position 7 of loop II) of 3.0 ppm compared to 4.9 ppm for Ca2+ binding, suggesting that the Thr62 carbonyl oxygen does not coordinate Mg2+ ion in the "lower affinity" conformation (Malmendal et al. 1998). A similar effect was observed in E140Q-Tr2C when Ca2+ binding was examined (Evenäs et al. 1997, 1999, 2001). Two sets of exclusive NOEs were observed in Ca2+ loaded states of E140Q-Tr2C, consistent with open and closed conformations (Evenäs et al. 1997). NMR relaxation data excluded a significant population of disordered species (Evenäs et al. 1999), and 15N chemical shifts of the residue in position 8 (Val136) suggested that in loop IV the backbone carbonyl group of Gln135 coordinates Ca2+ ion only in the open conformation (Evenäs et al. 2001).

The similarity between the conformational exchange between the "closed" conformation suggested by NMR (Evenäs et al. 1997, 1999, 2001; Malmendal et al. 1998) and our conformation B goes beyond the displacement of the Thr62 backbone carbonyl group. In loop II the >NH group of a highly conserved glycine in position 6 (Gly61) forms a hydrogen bond with the side chain carboxylate of Asp56 O{delta}1, which causes a downfield shift of Gly61 amide proton in NMR spectra (Malmendal et al. 1998). A continuous broadening of this signal is observed upon Mg2+ saturation of loop II, in contrast to loop I, supporting the inference about conformational exchange (Malmendal et al. 1998). In several MD simulations reported here the displacement of the Thr62 >C=O group from the Ca2+ coordination was accompanied by a loss of the interaction between the Gly61 >NH and Asp56 O{delta}1. For example, in simulation f, this interaction was lost near 200 psec, and in simulation p this interaction was lost near 700 psec (compare with plots in Fig. 9Go). Although the structures of states involved in conformational exchange reported by NMR are not known, NOEs, NMR relaxation data, as well as 1H and 15N chemical shifts suggest that conformations similar to apo and Ca2+-saturated states of CaM are involved (Evenäs et al. 1997, 2001; Malmendal et al. 1998). In our MD simulations Ca2+ binding loops generally sampled conformations near the initial Ca2+-saturated state, with backbone RMS deviations of ~0.5 Å. The relative "distance" of the observed structures to apo and Ca2+-saturated states was quantified by differences in RMS deviations relative to apo- and holo-structures ({Delta}RMSD = RMSapo - RMSholo). For loops I, III, and IV, {Delta}RMSD were steady across all 15 MD simulations in the range 1–2 Å. However, for loop II, {Delta}RMSD below 0.5 Å was observed in five MD runs (d, e, j, o, and p); negative {Delta}RMSD values were observed in simulation p, which implies that the Ca2+ binding loop sampled conformations more similar to the Ca2+-free form than to the initial, Ca2+-loaded structure (Fig. 6Go). In all five of these simulations the interaction between the Thr62 >C=O group and the Ca2+ ion was interrupted, either transiently or permanently, suggesting that the displacement of the Thr62 >C=O group from Ca2+ coordination involved a tendency of the Ca2+ binding loop II to adopt conformations more similar to the apo state.

The "closed" conformations observed in Tr1C (Malmendal et al. 1998) and E140Q-Tr2C (Evenäs et al. 2001) are likely to involve further conformational changes compared to our conformation B. For example, analysis of NMR chemical shifts suggested that in E140Q-Tr2C the "closed" state involves conformational change in the {chi}1 dihedral angle of Ile100 (Evenäs et al. 2001), which was not observed for the analogous side chain (Ile27) in our MD simulations. This suggests that the conformation B represents a substate of "closed" conformations observed by NMR, in which the Thr62 >C=O group is displaced from Ca2+ coordination but the hydrophobic cleft is still "open." In (Ca2+)2-E140Q-Tr2C, the exchange occurs between the states with open and closed hydrophobic clefts, while in (Mg2+)2-Tr1C, the exchange occurs between the two states that both have closed hydrophobic clefts (A. Malmendal, pers. commun.). The closure of the hydrophobic cleft appears to be linked to the lost Ca2+ ligation by Glu in loop position 12, and bidentate ligation by this residue persisted throughout all MD simulations.

Despite the partial similarity between the conformation B and the "closed" structure suggested by NMR studies (Evenäs et al. 1997, 1999, 2001; Malmendal et al. 1998), the exchange processes that involve these conformations occur on vastly different time scales. The conformational exchange observed by NMR in Tr1C and E140Q-Tr2C occurs on the microsecond time scale (Evenäs et al. 1997, 2001; Malmendal et al. 1998). For example, the rate constant between closed and open conformations was estimated to be (2.7 ± 1.0) x 104 sec-1 in E140Q-Tr2C (Evenäs et al. 2001). This process is too slow to be observed in our MD simulations. Let us assume a two state reaction between conformations A and B and first-order kinetics with equal rate constants for the forward and backward reaction (k = k1 = k-1), and denote the initial concentration of the conformation A with [A]0. Even for a rate constant of k = 106 sec-1, only 10-3 x [A]0 of the initial conformation A would convert into B within 1 nsec. For the conformation A to decay to 2/3 x [A]0 within 1 nsec (which would create 1/3 x [A]0 of the conformation B), the rate constant must be approximately k = 108 sec-1. This allows us to roughly estimate that the conformational exchange suggested by MD simulations occurs on the 10–100-nsec time scale. Much longer MD simulations would be required to accurately characterize the kinetics of this process. On the other hand, the fast conformational exchange suggested by our MD simulations is not observable in NMR experiments such as those used in studies of Tr1C and Tr2C (Evenäs et al. 1997, 1999, 2001; Malmendal et al. 1998). Rapid exchange on the nanosecond time scale would lead to averaging of chemical shifts, and relaxation dispersion experiments cannot resolve processes faster than about 1 µsec (Evenäs et al. 2001). 15N relaxation measurements coupled with a model-free analysis are capable of resolving rapid fluctuations on the picosecond-nanosecond time scale, and were reported for E140Q-Tr2C (Evenäs et al. 1999). However, the exchange processes suggested by our MD simulations fall in the time scale for which accurate determination of model-free order parameters is difficult (Andrec et al. 1999).

In canonical EF-hands the backbone carbonyl group in position 7 is the only main-chain group that serves as a Ca2+ ligand (Strynadka and James 1989). In both Tr1C and E140Q-Tr2C the experimental evidence suggests the displacement of the main-chain >C=O ligand in the second EF-hand (EF-hand II in Tr1C [Malmendal et al. 1998] and EF-hand IV in E140Q-Tr2C [Evenäs et al. 2001]). In several EF-hand proteins Ca2+ binding loop I is structurally more similar to loop III, and loop II is more similar to loop IV (Strynadka and James 1989). Furthermore, asymmetry in dynamic properties of paired EF-hands is well documented in CaM (Evenäs et al. 1998; Malmendal et al. 1998), troponin C (Li et al. 1997; Gagné et al. 1998), and calbindin D9k (Akke et al. 1993; Mäler et al. 2000). The origin of this asymmetry in dynamics is not understood, but is likely to play a functional role (Gagné et al. 1998). We did not observe the displacement of the backbone carbonyl group in position 7 (Gln135) in EF-hand IV. It is possible that EF-hand IV does not exhibit effects similar to EF-hand II. Alternatively, the mutation of Asp129, which is located in the first position of the binding loop IV, may alter the potential surface in such a way to eliminate this effect. Asp in the first position of the EF-hand Ca2+ binding loop is absolutely conserved in all known EF-hands (Falke et al. 1994), and analogous mutation Asp->Asn diminished Ca2+ binding in troponin C, suggesting profound effects of the mutation (Babu et al. 1992). In our MD simulations, the Ca2+ ion in EF-hand IV exhibited anomalous eightfold coordination (Fig. 8Go), which is presumably a consequence of D129N mutation. Such coordination is expected to lead to destabilized Ca2+ binding, because all known Ca2+ binding sites in proteins exhibit exclusively six- or sevenfold coordination (Strynadka and James 1989).

Proteins are known to undergo a complex, multibarrier dynamics with hierarchical ordering of conformational substates (Frauenfelder et al. 1991; Becker and Karplus 1997). Protein potential energy surfaces are rugged, characterized by many minima that cluster into "basins of attraction" (Frauenfelder et al. 1991; Becker and Karplus 1997; Caves et al. 1998), and slower processes could be understood in terms of basin-to-basin kinetics (Becker and Karplus 1997). Our MD simulations suggest that the attraction basin of Ca2+ binding loop II includes two local substates (A and B), where the substate B involves a conformation with the Thr62 >C=O group displaced from the Ca2+ coordination and the Ca2+ ion concomitantly ligated with two water molecules. The energy barrier between these two conformational substates is low, of the order of kT, and rapid exchange (estimated to involve ~100-nsec time scale) occurs at room temperature. The conformational exchange observed for Tr1C (Malmendal et al. 1998) and E140Q-Tr2C mutants (Evenäs et al. 1997, 1999, 2001) occurs on a much slower time scale, and matches the Ca2+ off-rate. This process was suggested to occur in wild-type protein, and to provide gating for Ca2+ dissociation (Evenäs et al. 2001). We hypothesize that the conformational states observed in our MD simulations represent higher level "tiers" involved in attraction basins (Frauenfelder et al. 1991; Becker and Karplus 1997) for the conformational exchange on the µsec-time scale suggested by NMR data. The resemblance of our conformation B to the "closed" conformations observed by NMR (Evenäs et al. 1997, 1999, 2001; Malmendal et al. 1998) suggests that the conformation B may function as a connectivity node on the transition pathway between basins that give rise to µsec-dynamics (Becker and Karplus 1997). The "closed" conformation state associated with µsec-dynamics observed in Tr1C and E140Q-Tr2C involves conformational changes beyond those observed in our MD simulations.


   Materials and methods
TOP
 Abstract
 Introduction
 Results
 Discussion
 Materials and methods
References
 
Definitions of secondary structure and Ca2+ sites
The eight {alpha}-helices of CaM are labeled A–H, in the order from the N terminus to the C terminus. The definition of {alpha}-helices in terms of residue numbers is as follows: 6–19 (A), 29–38 (B), 45–55 (C), 65–75 (D), 83–92 (E), 102–111 (F), 118–128 (G), and 138–146 (H). These definitions are basically the same as reported in the crystal structure (Babu et al. 1988), with only the terminal helices A and H shortened for one residue. Furthermore, the long central helix observed in the crystal structure (residues 65–92) was formally split into two helices to exclude residues 76–82, which form a flexible linker (Barbato et al. 1992; Crivici and Ikura 1995). Individual EF-hands are labeled in the standard notation as I, II, III, and IV in the order from the N terminus to the C terminus.

Preparation of MD simulations
In all 15 MD simulations (simulations b–p, Table 1Go) the initial coordinates of calmodulin were taken from the crystal structure of mammalian Ca2+-CaM (Babu et al. 1988) refined at 2.2 Å (Protein Data Bank code 3CLN [PDB] ). Four Ca2+ ions and 69 ordered water oxygen sites observed in the crystal structure were also included in the initial structure. Five residues missing in the crystal structure (1–4 and residue 148) were reconstructed in an extended conformation. All charged residues were taken in their standard states at pH = 7, which resulted in a total protein charge of -15.

The initial molecular model prior to solvation consisted of 2264 polypeptide atoms, 69 water molecules, and 4 Ca2+ ions. MD simulations were set up in a rectangular parallelepiped with dimensions 89 x 67 x 67 Å and periodic boundary conditions. The initial structure was placed in the center of a previously equilibrated box of water molecules so that the protein’s elongated dimension was oriented along the 89-Å dimension. Water molecules from the water box that overlapped with the polypeptide were deleted (crystal water was also considered to be the part of the initial structure). A water molecule was deemed overlapped if its oxygen was within a certain cutoff of any heavy atom of the initial structure. Seven different preequilibrated water boxes and five different water cutoffs were used to initialize 15 independent MD simulations, as shown in Table 1Go. In summary, a slightly different but equally plausible solvent configuration was created in each simulation to maximize the differences in initial conditions.

One molecule of CaM per volume of the water box used in these simulations corresponds to ~4 mM protein concentration. To neutralize the system a total of 15 Na+ ions were added; additional 13 Na+ ions and 13 Cl- ions were added in each simulation to mimic the solution with an ionic strength of ~100 mM. These ions were added iteratively, replacing the water molecules at the minima and maxima of the electrostatic potential calculated based on the protein atoms, four Ca2+ ions, and any previously placed Na+ and Cl- ions. Fifteen Na+ ions were added first to neutralize the system, and subsequently one Cl- ion followed by one Na+ ion was added until the total ion count was 41. Ions were not allowed to be placed closer than 6.5 Å and also not further away than 17.0 Å to any protein heavy atom; during the placement procedure any two ions were not allowed to approach one another closer than 5.5 Å. In each MD simulation the system consisted of approximately 40,000 atoms (the total number of water molecules present in each system after the placement of Na+/Cl- ions is shown in Table 1Go). Each simulation was equilibrated for 390 psec prior to running 1 nsec of the production dynamics (except for simulation b, which was equilibrated for 300 psec).

Parameters for MD simulations
All simulations were carried out with the program NAMD2, versions 22 and 23b2 (Kale et al. 1999), and CHARMM 22 force field (MacKerell et al. 1998). Water molecules were represented with the TIP3P water model, which corresponds to the water model used in previous MD simulations of CaM (Wriggers et al. 1998; Yang et al. 2001). The parameters for the calcium ion were obtained from the study of calbindin (Marchand and Roux 1998), a related EF-hand protein.

An integration step of 1 fsec was used to propagate equations of motion in the microcanonical ensemble (NVE). Nonbonded van der Waals interactions were smoothly truncated at 12.0 Å, with the switching function activated at 10.0 Å. The particle-mesh Ewald method was used to avoid truncation of electrostatic interactions. During the course of production dynamics coordinates of the entire system were saved every 0.2 psec, resulting in 5000 coordinate frames per 1-nsec simulation. The distance between the protein and its nearest image was monitored in all simulations, and was typically >20.0 Å (the distance between protein heavy atoms). In simulation h random tumbling of the protein brought th