INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSION REFERENCES

Calcium is widely used in biological systems as a regulator of physiological function. Eukaryotic intracellular resting levels of free Ca²⁺ are maintained at a concentration of ~10⁷ M via membrane-bound Ca²⁺-ATP pumps. Ca²⁺ influx responses increase the cytosolic Ca²⁺ concentration to ~10⁵ M as a result of a specific signal stimulus, as in neuron stimulation of muscle cells or hormonal stimulation of a cell through the binding of specific receptors. Regulatory Ca²⁺-binding proteins that direct physiological processes through Ca²⁺-induced conformational changes have been designed by nature such that their Ca²⁺ affinity is sufficient to bind Ca²⁺ during an influx, but not at resting levels (for review see da Silva and Reinach, 1991; Kawasaki and Kretsinger, 1994; Chazin, 1995; and Ikura, 1996).

The EF-hand family is a large class of Ca²⁺-binding proteins that contain homologous Ca²⁺-binding sites within a characteristic helix-loop-helix motif (da Silva and Reinach, 1991; Kawasaki and Kretsinger, 1994; Falke et al., 1994; Chazin, 1995). EF-hands are generally found back-to-back in anti-parallel pairs with -sheet-like hydrogen bonding occurring between the loops of the coupled sites. These coupled sites are often found to have cooperative metal ion binding, as in the case of calmodulin. The two parvalbumin sites do form an anti-parallel pair, but they do not exhibit cooperative binding. However, EF-hand loops are all strongly linked to context and do not exhibit the same binding properties when separated from their protein environments. The EF-hand binding loop is traditionally defined as the 12 sequential residues starting with the first coordinating residue of the loop and ending with the last coordinating residue, although the last three residues of this sequence are actually the first three residues of the trailing helix. Six of the 12 loop residues coordinate the bound metal ion, and all coordinating ligands are oxygens (Fig. 1). There is a highly conserved glutamate at loop position 12 that changes from a bidentate ligand in the sevenfold ligation of Ca²⁺ to a monodentate ligand in the sixfold coordination of Mg²⁺. The position 12 glutamate was found to be 92% conserved in an analysis of 567 EF-hand sequences performed by Falke et al. (1994). The remaining 8% of the EF-hand sequences analyzed that did not contain glutamate at position 12 all contained aspartate. Our search for EF-hands with aspartate at position 12 turned up only varying species of myosin regulatory light chain (RLC); the remainder of the well-known regulatory Ca²⁺-binding proteins, such as troponin C and calmodulin, contain glutamate at position 12 in all of their Ca²⁺-binding sites. It is interesting to note that, in the published crystal structure of myosin RLC, the binding site containing aspartate at position 12 is conjectured to be Mg²⁺ specific under physiological conditions (Houdusse and Cohen, 1996). It has previously been postulated that this last coordinating residue of the binding loop plays a crucial role in determining the metal ion affinities and specificities of individual EF-hands (Strynadka and James, 1991; Declercq et al., 1991; da Silva et al., 1995; Houdusse and Cohen, 1996).

View larger version (36K): [in this window] [in a new window]   FIGURE 1   Coordinating residues of an EF-hand site. The CD site of carp parvalbumin is used here to illustrate the coordinating positions of a typical EF-hand site. The residues at loop positions 1, 3, 5, 7, 9, and 12 coordinate the bound metal ion. The residue at loop position 1 is almost always an aspartate. Positions 3 and 5 are most commonly aspartates also, although here, in the CD site of parvalbumin, there is a serine at position 5. The backbone carbonyl oxygen of the residue at position 7 coordinates the metal ion, so there is a great deal of variability in amino acid identity at this position. When the residue at position 9 does not have a sidechain oxygen to donate to the metal ion, a water molecule will be recruited into the binding site at this position. This is the case in the EF site of parvalbumin (not shown). The glutamate at position 12 is the residue that transforms the coordination sphere from sevenfold Ca2+ coordination (as illustrated), to sixfold Mg2+ coordination, through switching from bidentate ligation to monodentate (see Falke et al., 1994 for review of amino acid conservation in EF-hands).

The responses to changes in Ca²⁺ levels exhibited by individual members of the EF-hand family are quite diverse, and the metal-ion binding parameters play a principal role in determining these responses. Ca²⁺ and Mg²⁺ are the most physiologically relevant metal ions that interact with EF-hand proteins. Both Mg²⁺ and Ca²⁺ are small, closed-shell, spherical metal ions, more likely to form ionic bonds than covalent. There are subtle differences between the two divalent cations that the EF-hand binding site can exploit to select between the two. Mg²⁺ favors sixfold, octahedral coordination, whereas Ca²⁺ is most commonly found coordinated by seven or eight ligands. Similarly, Ca²⁺-ligand bond distances are typically 2.3-2.6 Å, whereas Mg²⁺-ligand bond lengths average ~2.0-2.1 Å in length (Martin, 1990). Finally, as a result of its greater surface charge density, Mg²⁺ has a 10³ times slower desolvation rate than Ca²⁺ (Falke et al., 1994).

Members of the EF-hand family include parvalbumin, and many prominent regulatory proteins such as calmodulin, troponin C, calpain, myosin regulatory light chains, and calcineurin. EF-hand proteins are involved in a wide variety of physiological processes, including signaling, cell cycle regulation, second messenger production, muscle contraction, and vision. At least 200 different EF-hand proteins have already been putatively identified in the human genome (http://genome.wustl.edu/est/hmm_results.html). Ca²⁺-binding proteins have been implicated in many serious disease states such as Alzheimer's disease (Vito et al., 1996), and other age-related cognitive defects (Krzywkowski et al., 1996), diseases of the nervous system (Zimmer et al., 1995), leukemia (Calabretta et al., 1985), and various forms of cancer (Blum and Berchtold, 1994; Polans et al., 1995). The numerous disease states that have been associated with EF-hand proteins emphasize the immense physiological importance of this family of Ca²⁺-binding proteins. A better understanding of the mechanisms by which these proteins distinguish one metal ion from another would be beneficial for effectively addressing these disease states with drug-targeting strategies (Chazin, 1995).

Parvalbumin is a small, ~12-kDa protein with two EF-hand type Ca²⁺-binding sites. It is commonly divided into three domains, each containing a helix-loop-helix motif. The domains AB, CD, and EF are named according to the two flanking helices. Domain AB contains a two amino-acid deletion in the loop region and consequently does not bind metal ions (Kretsinger and Nockolds, 1973); domain CD contains the N-terminal binding site, and domain EF, the C-terminal site.

Parvalbumin has been chosen here as a model chiefly because of its close structural relationship to larger, regulatory proteins. The two EF-hand binding sites present in wild-type parvalbumin are homologous to the EF-hand binding sites found in troponin C, calmodulin, and many other regulatory members of the EF-hand family (Chazin, 1995). These regulatory proteins, because of their increased structural and functional complexity, are more difficult to deal with than parvalbumin, both computationally and experimentally.

Parvalbumin lends itself well to molecular dynamics simulations because it is small and relatively stable. An important factor affecting the reliability of molecular dynamics (MD) simulations is the accuracy of the starting model (Melton et al., 2000). Declercq et al. (1999) have recently published a 0.91-Å resolution crystal structure of  parvalbumin from pike. This structure provides ideal starting coordinates for our computational modeling.

The simulations in this article were devised to expound theories derived from a crystallographic analysis of a carp parvalbumin mutation, E101D, whereby the glutamate at EF loop position 12 has been mutated to aspartate (Cates et al., 1999). As a result of this mutation, the EF site binds Ca²⁺ with sixfold coordination, instead of the characteristic wild-type sevenfold coordination. Out of 108 total residues, 84 are identical between pike and carp parvalbumin (Table 1), constituting a 78% identity between the two proteins. The coordinating residues in the pike versus carp parvalbumin EF-hands are identical, with the exception that the amino acid that offers its backbone carbonyl oxygen to the bound cation in the second binding site changes from lysine in wild-type carp parvalbumin to methionine in wild-type pike parvalbumin. The identity of the sidechain does not affect metal ion coordination at this position; moreover, these sidechains are seen to be protruding into the solvent at quite similar orientations when the pike and carp crystal structures are superimposed. When the first eight residues of the dynamic N-terminus are removed from both proteins, the structural root mean square deviation (RMSD) between the superimposed backbones of these two parvalbumins is 0.5 Å. Therefore, it was determined that the pike parvalbumin model system would yield results that were readily comparable to experimental results in carp parvalbumin.

MD simulations can be used advantageously in conjunction with experimental data to investigate proposed mechanisms of metal ion selectivity in parvalbumin. The simulations allow a visualization of the chain of events leading up to an experimental observation, and this approach often provides insights that cannot be immediately derived from experimental results. Moreover, the simulations approximate the energy fluctuations that correlate with specific biochemical events and thereby provide us with a useful quantitative analysis, although it is important to understand the limitations of the methodology when interpreting these results. The accuracy of the simulation is limited in that the electronic configurations are represented by point charges, instead of electronic distribution probabilities. Also, though much progress has been made regarding solvent modeling, the ability to comprehensively represent the varied properties of ordinary water is still developing. However, for systems where protein structure, nonbonded van der Waals interactions, and electrostatic interactions are the prevailing influences on the inherent functional properties, MD simulations often provide reasonable representations. Although there is much room for improvement in the degree of accuracy with which MD calculations represent proteins, it is nonetheless true that these methods, despite the approximations and truncations used, are frequently able to correctly predict or explain the behavior of proteins.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSION REFERENCES

Molecular dynamics starting model

The parvalbumin starting coordinates for each of the simulations originated from the 0.91-Å atomic resolution crystal structure of pike parvalbumin published by Declercq et al., 1999. These coordinates were obtained from the Protein Data Bank (PDB), accession code 2PVB. The hydrogen atoms in the model were built in X-PLOR (Brunger, 1992), and CHARMm22 parameters and topology were used (MacKerell et al., 1998) except for the residues where we have specifically mentioned deriving modified parameters.

Model solvation

Parvalbumin was solvated in a cubic periodic boundary cell with edge lengths equal to 61.8 Å to generate an explicitly modeled bulk solvent consisting of 7736 H₂O molecules, 21 sodium ions, and 19 chloride ions. The protein was solvated by placing the pike parvalbumin coordinate set at the center of the water box and subtracting all water molecules within 2.3 Å of any protein atom. Randomly selected water molecules were replaced with sodium and chloride ions, added in amounts calculated to offset the charged amino acids present in the protein (Ibragimova and Wade, 1998). The sodium and chloride ions were then subjected to simulated annealing while holding all protein heavy atoms rigid, thereby allowing the sodium and chloride ions to achieve minimized positions without the added expense of in-depth electrostatic field calculations.

Model parameterization

Representing a Ca²⁺-binding site in MD simulations required analysis of the suitability of the important model parameters. We chose to use the CHARMm parameter set (MacKerell et al., 1998), but some of the parameters were not appropriate for a metal-containing environment. In our initial equilibration simulation, a glutamate in the CD site, glutamate 58, bound the metal ion in a bidentate fashion throughout the simulation. This resulted in eightfold Ca²⁺ coordination, because, normally, only the glutamate at the last coordinating position of the binding loop, glutamate 61, binds the metal ion with both oxygens. This finding prompted an evaluation of the Lennard-Jones van der Waals parameters for Ca²⁺ in the model, because the behavior of glutamate 58 was indicative of a larger metal ion than Ca²⁺.

As we parameterized the coordinating residues in the two parvalbumin binding sites, we experimented with varying van der Waals values within the range of values found in CHARMm and CNS (Brunger et al., 1998) parameter files. The final Lennard-Jones parameter values that best represented wild-type coordinating behavior in our model system were 2.4 Å for Ca²⁺ and 1.9 Å for Mg²⁺. This was in good agreement with the Leonard-Jones parameters for the Ca²⁺ ion calculated by Marchand and Roux (1998) for MD studies of calbindin D9k. Marchand et al. (1998) derived a final value for the van der Waals parameter for Ca²⁺ of 2.436 Å by determining the values that would reproduce the experimental free energy of hydration for the Ca²⁺ ion.

After several trial simulations, it also became clear that the serine at position 5 of the parvalbumin CD binding site consistently tended to drift away from the metal ion. The carbonyl oxygens that coordinate bound metal ions at position 7 of both sites also drifted away from the ion during simulations. Because the CHARMm charge distribution allotted to the serine sidechain and the backbone carbonyl was not calculated in the close proximity of a divalent cation, it was reasonable to assume that the oxygens at these specific positions in our model should be more polarized. Force fields are expressed in classical molecular dynamics programs in such a way as to make them nonpolarizable, except through manual re-parameterization of specific atoms in the model. To compute a more representative charge distribution, restricted Hartree-Fock single-point energy calculations were performed in Gaussian (Frisch et al., 1995) on a model system consisting of the coordinating sidechains of the binding loop and the calcium ion. The ratio of the resulting charge distribution was used to derive a more polarized charge distribution for the serine sidechain at position 5 in the CD site, Serine 55, and the carbonyl groups at position 7 in both sites.

CHARMm uses the TIP3P 3-point intermolecular potential function model for water (Jorgensen et al., 1983). Other representations, such as the TIP4P or TIP5P (Jorgensen et al., 1983; Mahoney and Jorgensen, 2000) and extended SPC/E (Berendsen et al., 1987) models, better mimic the properties of solvent water in many cases. However, the computational expense of using these models is greater than that incurred using the TIP3P model. Another water model cannot be substituted for use with the CHARMm parameters without leading to inconsistencies that could imbalance the intermolecular parameterization, because the CHARMm protein parameters were determined specifically for the TIP3P solvent representation (MacKerell et al., 1998). However, in our EF site simulations, where a water molecule is recruited as a Ca²⁺ ligand, the TIP3P water parameters were not adequate. Exchange of the EF site H₂O molecule with solvent water molecules might be expected, but, instead, the coordinated H₂O drifted away from the metal ion, and was not replaced. Therefore, it was necessary to tether the water to the ion by creating a new residue and specifying a water-metal ion bond length to correctly model the behavior of the water molecule that coordinates the Ca²⁺ in the EF site. This was deemed acceptable because the simulations are not designed to test the coordinating behavior of the water, instead they are designed for the purpose of testing the coordinating carboxylate sidechains, particularly the glutamate at position 12. The bond length was constrained to 2.37 Å for the water oxygen-Ca²⁺ bond distance, and 2.15 Å for the water oxygen-Mg²⁺ bond distance. Both of these distances were derived from analysis of the crystal structures of parvalbumins.

Simulations

The NAMD molecular dynamics simulation software (Nelson et al., 1996) was used for the simulations of parvalbumin. The cutoff distance for the van der Waals pairs was 8.5 Å. A smooth switching function was used at a switching distance of 8 Å. Full electrostatics were used in all of the NAMD parvalbumin simulations using particle-mesh-based methods for fast Ewald summation (Darden et al., 1997) to compute the electrostatic potential energy. The coupling constant used in conjunction with the Berendsen temperature-coupling algorithm was 0.40. The Berendsen pressure target was 1 bar, the compressibility was 0.000049 bar ¹, the relaxation time was 500 fs, and the number of time steps between application of pressure scaling was 12. The ShakeH algorithm was used to fix the bond between each hydrogen and its mother atom to the ideal bond length specified in the parameter file, and the timestep was 2 fs.

The steered MD applied in the Force 0 simulations used a spring constant value (k) of 414 pN/Å, a pulling velocity of 0.5 Å/ps, and specified the initial coordinates of the Ca²⁺ as the vector defining the direction of movement for the reference position. The duration of the simulation was 3000 fs, or 1500 timesteps, which would allow a maximum movement of 1.5 Å toward the Ca²⁺. The atom constrained to the moving reference position was the -carbon (C) of aspartate 101. (For further discussion of steered MD in NAMD, see Bhandarkar et al., 1999 and Lu and Schulten, 1999.)

The constraints in the Force 1 simulations were applied through the NAMD harmonic constraints parameters. The PDB file representing the coordinates to which the aspartate 101 sidechain was constrained was generated in the program O (Jones et al., 1991), by manually forcing the aspartate to be close enough to the Ca²⁺ ion to obtain bidentate coordination. The duration of the Force 1 simulations was 6000 fs, or 3000 timesteps.

The distance constraints in the Force 2 simulations were user-supplied forcing restraints, applied using a free energy perturbation script that required that the initial 94 C to 101 C be decreased by 1.5 Å over 100-ps simulation time. The potential energy function was changed using a coupling parameter () which was first incremented from a value of 0 to 1, then the simulation equilibrated at  = 1, and finally the value was reduced from  = 1 back to 0 again, following the slow growth protocol illustrated in the NAMD User's Guide (Bhandarkar et al., 1999).

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSION REFERENCES

The Alchemy simulations

The purpose of the first simulation, called Alchemy, was to investigate whether our classical MD model could reproduce the sevenfold coordination of Ca²⁺ and the transition to the sixfold coordination of Mg²⁺. This was a prerequisite for modeling amino acid substitutions in the wild-type protein. The ability to reproduce known parvalbumin metal ion binding behavior was particularly important in the EF site because the E101D mutation is located in this site.

Depiction of the transition from sevenfold Ca²⁺ coordination to sixfold Mg²⁺ coordination was accomplished by beginning the simulation with van der Waals parameters that represent Ca²⁺ in the binding sites, then reducing the radius of the bound cations during the simulation until the radius became representative of Mg²⁺. We were also interested in observing whether this transition was reversible as described by our model, because we know empirically that it is reversible in wild-type parvalbumin. Finally, it was of interest to see whether the glutamate at position 12, the last coordinating residue of the loop, was correctly positioned in the simulations to provide bidentate coordination of Ca²⁺. There are several carboxylate sidechains coordinating the Ca²⁺ ion in both binding sites, and the corresponding carboxylate oxygens are all described with equivalent electrostatic charge distributions as in the CHARMm parameter files. Therefore, with the nonpolarizable representation of our system, none of the carboxylate oxygens should be preferred by the metal ion over another on the basis of charge. If, in our simulations, the glutamate at position 12 were consistently and reversibly predicted to be the source of bidentate ligation of Ca²⁺, this would provide strong evidence that properties of the structural framework of EF-hand binding sites impose this role on the residue at position 12.

The protocol for the Alchemy simulation began with 50 ps of simulation, and a van der Waals radius for the two bound metal ions representative of Ca²⁺ (see Methods for parameterization of the Lennard-Jones  values). Then the radius was reduced by 0.1-Å increments over 20 ps in 5 equal steps. When the radius reached a van der Waals value for the bound cations representative of Mg²⁺, this was simulated for 50 ps. To analyze the reversibility of the simulation, the ionic van der Waals radius was increased back to the Ca²⁺ value over 20 ps, in an equivalent, but reciprocal manner. The Ca²⁺-bound parvalbumin was modeled at the last stage for an additional 50 ps of simulation, making the duration of the entire simulation 190 ps.

The Alchemy simulation was stable, and the maximum RMSD from the starting structure was 1.8 Å. This RMSD is within the acceptable range for MD simulations where the system is being perturbed, as it is here by the incremental shifts in the radii of the bound metal ions.

There are coordinating aspartate sidechains in the CD site at positions 1 and 3 of the binding loop, and in the EF site at positions 1, 3, and 5. All coordinating aspartates exhibited the same metal ion-binding behavior in the simulation that they exhibit in experimentally determined crystal structures (Fig. 2).

View larger version (47K): [in this window] [in a new window]   FIGURE 2   The Ca-O distances for the coordinating aspartates of the CD loop of the Alchemy simulation are presented in these graphs as examples of wild-type aspartate behavior in our simulations. In both residue 51 (A and B) and 53 (C and D), one of the oxygens of the aspartate carboxylate sidechain coordinated the Ca2+ ion, then moved in to an even shorter Mg2+-oxygen bond distance during the Mg2+ phase of the simulation. Then, finally, that same oxygen moved back out again to an appropriate Ca2+-oxygen distance in the final stage. Also, in both of these aspartates, the other oxygen remained too far away to coordinate the metal ion during all phases of the simulation. The wild-type aspartates in each binding site in all of the simulations exhibited this kind of customary metal ion binding behavior.

As described in Methods, the model parameterization required that the charge distribution be modified for the serine at position 5 of the CD site and the backbone carbonyl group that coordinates Ca²⁺ at position 7 of both loops. This modification was made to more accurately represent the polarization of the coordinating oxygens in the presence of a divalent cation. After adjustment, the backbone carbonyl oxygen in our simulations emulated the same Mg²⁺ and Ca²⁺ coordination that has been observed in crystal structures. Serine 55 at loop position 5 stayed within coordination distance of the metal ion when it was represented as Ca²⁺; however, when the ionic radius of the bound cation decreased to represent Mg²⁺, the serine had a propensity to drift away (Fig. 3). Interestingly, when the radius increased again, the serine moved back in to coordinate the Ca²⁺ ion. This movement away from the Mg²⁺ ion by the serine at position 5 might be an artifact of the imprecision of the charge distribution representation for serine in the MD simulation. In contrast, it might be suggestive of another way EF-hands discriminate between Mg²⁺ and Ca²⁺ binding.

View larger version (25K): [in this window] [in a new window]   FIGURE 3   Serine 55 at loop position 5 of the CD loop coordinates Ca2+ well, but moves away from Mg2+. (A) In this Alchemy simulation, the beginning radius is representative of Mg2+, the transition to Ca2+ is complete by frame 700 (70 ps, simulation time), and at frame 1400 the metal ion becomes Mg2+ again. The serine is coordinating the Ca2+ in a stable manner, but tends to dissociate from Mg2+. (B) This simulation begins with the Ca2+ ion bound, transitions to Mg2+ at frame 700, and reverts back to Ca2+ by frame 1400. Here, serine is clearly not coordinating the metal ion when its identity is Mg2+.

There are a variety of EF-hand proteins that contain serine in one of the binding loops (Nelson and Chazin, 1992). The structural information available about EF-hands that include serine at one of the coordinating loop positions implies that serine may play a role in metal ion selectivity. First of all, in the Mg²⁺-loaded pike parvalbumin structure published by Declercq et al. (1991), the CD site with the serine at position 5 still retains Ca²⁺; it is only the EF site that coordinates Mg²⁺. The calcium persisted in this site despite the fact that the protein solution used for crystallization was dialyzed successively against 1.5 M MgSO₄, until a Ca²⁺ concentration corresponding to 0.53 Ca atom per protein molecule was determined by atomic absorption spectrometry. There has not been a doubly loaded parvalbumin/Mg²⁺ structure published to date, to illustrate Mg²⁺ binding in the CD site. Moreover, NMR evidence (Blancuzzi et al., 1993) indicates that titration of apo-parvalbumin with Mg²⁺ results in only one singly loaded Mg²⁺ intermediate before formation of the doubly loaded parvalbumin/Mg²⁺ complex. This singly loaded intermediate has Mg²⁺ in the EF site and nothing bound in the CD site. The intermediate with Mg²⁺ in the CD site and nothing bound in the EF site has not been observed in these NMR Mg²⁺ titrations. Furthermore, the titration of fully loaded parvalbumin/Mg²⁺ with Ca²⁺ exhibits only a single intermediate species, in which Ca²⁺ is bound in the CD site and Mg²⁺ is bound in the EF site. The opposite intermediate, Mg²⁺ in the CD site and Ca²⁺ in the EF site, is never observed.

It should be mentioned, though, that a theoretical study has been performed by Allouche et al. (1999) wherein the Ca²⁺ in the CD site of this pike parvalbumin crystal structure has been transformed computationally into Mg²⁺. In the Allouche et al. theoretical model, the serine at position 5 coordinates the Mg²⁺ ion adequately. The obtained value of roughly 10³ for the ratio of Ca²⁺ and Mg²⁺ affinity constants in the CD site was in good agreement with experimental observations. However, many constraints were applied in the free energy perturbation model that would make it difficult for any of the ligands to move very far away from the bound metal ion. The metal ion and all heavy atoms lying more than 9 Å away from the metal were held fixed, and all atoms of any type lying more than 11 Å from the metal ion were also fixed. This makes it impossible for the bulk of the protein just outside of the binding site to be flexible and therefore makes it unlikely that the loop could accommodate much movement away from the metal ion, even if one of the coordinating residues were inclined to do so.

In an effort to investigate the behavior of serine in a Mg²⁺-bound EF-hand, we looked for Mg²⁺-loaded structures in the PDB. There were no Mg²⁺/EF-hand complexes for EF-hand sites that contain serine at position 5, similar to parvalbumin. However PDB files for both calbindin D9k, PDB accession code 5ICB (Andersson et al., 1997), and the myosin RLC, PDB accession code 1WDC (Houdusse and Cohen, 1996), were analyzed. These structures contain a serine at loop position 9 in the binding sites containing Mg²⁺. This serine does not directly coordinate the Mg²⁺ ion in either structure; instead, a water molecule is recruited as a ligand (Fig. 4). Other EF-hand proteins have serine residues at key loop positions, but their structures are either unknown or they are available only for the Ca²⁺-bound species. For example, serine can be found at position 5 in the second binding sites of certain species of calmodulin and in most species of skeletal muscle troponin C. However, these sites in calmodulin and troponin C are Ca²⁺ specific at physiological levels of Ca²⁺ and Mg²⁺.

View larger version (20K): [in this window] [in a new window]   FIGURE 4   Serine in the EF-hand sites of Mg2+-bound Calbindin D9k and Myosin RLC. (A) Serine 62 in Calbindin D9k is at position 9 of the binding loop, but in this Mg2+-bound EF-hand site, loop position 9 is too far out of the coordination sphere to bind metal ions. Instead, a water molecule is recruited in both Mg2+ and Ca2+ coordination (Andersson et al., 1997). (B) Serine 36 at position 9 in the Mg2+-bound site of myosin RLC also does not directly coordinate the Mg2+ ion. It does, however, coordinate the water molecule that substitutes for a sidechain ligand at position 9. This site is believed to be a Mg2+ site in vivo (Houdusse and Cohen, 1996).

Blumenschein and Reinach (2000) have introduced a serine-for-aspartate substitution at position 5 in the myosin RLC, which lowers the Mg²⁺ affinity by a factor of 33 in comparison to wild-type. However, in conjunction with a glutamate-for-aspartate substitution at position 9 (causing both of these coordinating residues to mimic the correlating positions in the parvalbumin CD site) a double mutant results that only lowers the Mg²⁺ affinity by a factor of 2 in comparison to wild-type RLC. It can be seen in parvalbumin/Ca²⁺ complexes (PDB accession codes 5CPV [Swain et al., 1989] and 2PVB [Declercq et al., 1999]) that the sidechains of the serine at position 5 and the glutamate at position 9 are within hydrogen-bonding distance. By reason of the smaller coordination sphere defined by the structural framework of the RLC binding loop, Blumenschein and Reinach propose that this hydrogen bond can occur in their site only when Mg²⁺ is bound, and suggest that this interaction stabilizes the binding site and is putatively the source of the improvement in Mg²⁺ affinity in the D5S/D9E double mutant. An analysis of the interaction between serine 55 and glutamate 59 of the parvalbumin CD site in the Alchemy simulation shows that the serine O and the noncoordinated glutamate O₂ are close enough to hydrogen bond only during the Ca²⁺-bound phase of the simulation.

It would require further investigation to state definitively that serine plays a role in metal ion selectivity through stabilizing the binding site region when ions are bound that allow it to make favorable hydrogen bonds. It does seem worth exploring the possibility, though, that EF-hand proteins might exploit such a property to confer a distinction between Mg²⁺- and Ca²⁺-binding parameters in certain binding sites.

In the Alchemy simulations, the transition from sevenfold Ca²⁺ coordination to sixfold Mg²⁺ coordination is observed in both binding sites. Furthermore, this transition is reversible in the simulations and is a function of the glutamate at position 12. This glutamate exhibits bidentate coordination of Ca²⁺ and monodentate coordination of Mg²⁺ in the simulation, just as it has been observed to do in crystal structures (Fig. 5, A-D).

View larger version (38K): [in this window] [in a new window]   FIGURE 5   The glutamate at position 12 is responsible for the transition in metal ion binding geometry in the Alchemy simulation. (A) In this frame chosen from the Mg2+ phase of the simulation, the glutamate's coordinating sidechain oxygen is 1.99 Å from the Mg2+ ion. The noncoordinating oxygen is 2.97 Å from the Mg2+. (B) In this frame from the Ca2+ phase, the uppermost oxygen is 2.24 Å from the metal ion, and the other oxygen is 2.58 Å away. (C) This graph of the distance between the position 12 glutamate coordinating sidechain oxygen, O1, and the metal ion during the Alchemy simulation shows that it stays bound to the metal ion during all three phases of the simulation (Ca2+/Mg2+/Ca2+). (D) The second sidechain oxygen, O2, stays bound to Ca2+ in the first and last stages of this simulation, but moves out of coordination distance in the Mg2+ stage of the simulation. This result accurately mimics the bidentate coordination of Ca2+ and the monodentate coordination of Mg2+ that is observed in parvalbumin crystal structures.

What is noteworthy is that glutamate 59, the coordinating residue at position 9 of the CD loop, does not bind either metal ion in a bidentate fashion. It is known from crystal structures that this glutamate is a monodentate ligand. The charge distributions on all of the carboxylate sidechains in the CD loop are equivalent; therefore, there is no inherent electrostatic preference for one of the carboxylate sidechains over another. The aspartates in the binding site might be less likely to attain bidentate coordination because they are shorter and less flexible. However, the only apparent reason for the model to select the glutamate at position 12 for bidentate coordination of Ca²⁺ over the glutamate at position 9 would be because the geometry of the binding loop prescribes that choice. This observation is a strong indication that the crucial role the residue at position 12 plays in transforming the coordinating sphere through reducing the number of ligands by one when Mg²⁺ is bound, and decreasing the radius of the coordination sphere by extending one oxygen toward the Mg²⁺ ion, is not solely attributable to the length and flexibility of the glutamate sidechain. Instead, it seems this role is at least partly assigned to the residue at position 12 as a function of the structural configuration of the loop.

There are no nearby steric hindrances that would prevent both position 12 sidechain oxygens from binding Ca²⁺. Also, the other coordinating sidechains are cushioned within the binding loop, where nearby residues can easily make intramolecular hydrogen bonds with the noncoordinating carboxylate oxygens. In contrast, the position 12 glutamate is positioned at the outermost edge of the loop and thus is very solvent exposed. This leaves both of the carboxylate oxygens very available to the bound metal ion. Hence, it is likely that these circumstances conferred on the residue at position 12 by the structural composition of the loop act in combination with the length and flexibility of the glutamate to prescribe its role as the transitional factor in Ca²⁺ coordination versus Mg²⁺.

The Aspartate simulation

In the parvalbumin E101D/Ca²⁺ crystal structure previously reported (Cates et al., 1999), when the residue at EF loop position 12 has been mutated from glutamate to aspartate, it only binds Ca²⁺ in a monodentate fashion, resulting in sixfold Ca²⁺ coordination. It is highly unusual for Ca²⁺ to experience sixfold coordination, and it has never been reported in EF-hand binding sites before the E101D mutant. Indeed, sixfold coordination is characteristic of EF-hand Mg²⁺ binding. The E101D substitution affects binding kinetics such that the Ca²⁺ affinity decreases 100-fold and Mg²⁺ affinity increases 10-fold. Additionally, the Ca²⁺ off-rate in the E101D site increases from 3/s to 630/s. This is consistent with the fact that all Ca²⁺-specific, regulatory EF-sites, such as those found in calmodulin and the two N-terminal sites in troponin C, contain glutamate at position 12, whereas the myosin RLC domain I EF-hand site is putatively Mg²⁺ specific (Houdusse and Cohen, 1996), and this site contains aspartate at position 12.

The purpose of the Aspartate simulation was to investigate whether our model would correlate the crystallographic results. There were two motivations for performing this simulation. First, it was of interest to visualize the movement of the F helix in toward the metal ion; the end result of this movement is seen in the parvalbumin E101D crystal structures. Second, if the simulation could accurately correlate the results seen experimentally, this would indicate that our model is a good system in which to perform further simulations designed to test the flexibility of the loop through applying external forces.

It was noteworthy to observe the movement of the F helix, the helix C terminal to the binding loop, because it was observed to move into the loop ~1 Å in the parvalbumin E101D/Mg²⁺ crystal structure, and the F helix moved in ~0.7 Å in the parvalbumin E101D/Ca²⁺ structure. In fact, it seems that there is a tug-of-war between the attraction of the carboxylate sidechain oxygens for Ca²⁺ and their requirement to drag the anchoring protein along with them as they move in toward the metal ion; and it is likely that this tug-of-war is a strong determining factor defining the range of coordination spheres that can be achieved in a particular EF-hand site.

The protocol for the Aspartate simulation began with the replacement of the glutamate at coordinating position 12 (residue no. 101) with aspartate in the initial coordinate set. A simulation was run for 25,000 steps, representing 50 ps of simulation time, with a van der Waals radius for the bound metal ions that was representative of Ca²⁺. The RMS deviation from the starting structure never increased above 1.4 Å for the Aspartate simulation. The unaltered CD site coordinated Ca²⁺ with normal wild-type pentagonal bipyramid Ca²⁺-binding geometry.

Because the substitution of aspartate for glutamate was performed by altering only the sidechain coordinates, at the start of the Aspartate simulation, the mainchain of aspartate 101 was in the same position usually occupied by the mainchain of glutamate 101. For that reason, the two aspartate 101 carboxylate oxygens were both out of coordination range of the Ca²⁺ ion at the beginning of the simulation. In the starting coordinate file the two aspartate sidechain carboxylate oxygens, O₁ and O₂, are 4.97 Å and 2.94 Å from the Ca²⁺, respectively. After only 80 fs, they moved in to a Ca-O distance of 2.20 and 3.87 Å (Fig. 6, A and B). At this point, the closest oxygen is definitely coordinating the metal ion, and the other oxygen is not. Throughout the remainder of the simulation, the aspartate at position 12 displays monodentate coordination of the Ca²⁺ ion (Fig. 6, C and D). Just as in the crystal structure, the aspartate in our simulation cannot move in enough to allow Ca²⁺ coordination with both oxygens.

View larger version (37K): [in this window] [in a new window]   FIGURE 6   Aspartate 101 moves in to coordinate the Ca2+ within the first 80 fs of simulation. In the Aspartate simulation, the starting Ca-O distances for the two sidechain oxygens of Aspartate 101 were 2.94 Å for the closest, and 4.97 Å for the most distant. (A) After 40 fs of simulation, the Ca-O distances were 4.72 and 2.74 Å. (B) By the end of the first 80 fs of simulation, the closest sidechain oxygen was coordinating the Ca2+ at 2.2 Å. Aspartate 101 never achieved bidentate coordination, however, as the second sidechain oxygen remained about 3.5 Å from the Ca2+ throughout the simulation. (C) The coordinating sidechain oxygen has moved in to bind the Ca2+ by the second frame of the simulation (80 fs). It stays ~2.1 Å from the Ca2+ throughout the simulation. (D) This graph, illustrating the Ca-O distance for the noncoordinating oxygen, shows that it stays ~3.5 Å from the Ca2+ ion throughout most of the simulation.

The F helix shows obvious movement into the binding site toward the metal ion during the Aspartate simulation. A comparison of the distances of the F helix backbone nitrogens from the Ca²⁺ ion reveals that the entire helix moves inward toward the binding site during the first 20 ps of the simulation (Table 2).

The Force simulations

It is shown in the parvalbumin E101D x-ray structure (Cates, 1999), and in the Aspartate simulation, that an aspartate at the last coordinating position in the parvalbumin EF site is unable to offer bidentate coordination to a Ca²⁺ ion in the same way the wild-type glutamate does. It is difficult to distinguish between effects caused by the difference in the flexibility of an aspartate in comparison to a glutamate, and effects due to constraints placed on the radius of the coordination sphere by the backbone structure of the binding loop. In the first case, the aspartate might not be able to achieve bidentate coordination because, lacking the added flexibility bestowed by the extra ₃ torsion angle of the glutamate, it cannot position both of its two oxygens in a favorable orientation with respect to the rest of the oxygens in the binding loop. In the second case, the aspartate may fail to provide bidentate coordination merely because it is not long enough, and the backbone of the binding loop is not sufficiently pliant to move in far enough for the aspartate to reach the Ca²⁺ with both oxygens. Therefore, the Force simulations were designed to test the plasticity of the binding loop through applying external forces. All of the Force simulations used the final coordinates of the Aspartate simulation as the initial coordinate set, wherein the loop backbone has moved in toward the binding site about 1 Å in comparison to the wild-type loop. We wanted to see if an external force could compel the backbone to move in closer and allow the aspartate to coordinate Ca²⁺ with both oxygens; and if it could, we were interested in the effect this movement would have on the loop conformation.

Our first indication that at least part of the reason that the aspartate does not achieve bidentate Ca²⁺ coordination is attributable to resistance from the loop backbone was in the Force 0 simulations. These simulations used NAMD steered MD (Nelson, 1996) to apply a harmonic restraining force to the aspartate 101 C coupled to a moving reference position pulling in the direction of the bound Ca²⁺. With none of the backbone atoms fixed, the result of this simulation was merely to drag the entire protein through the solvent in the direction of the moving reference position. The loop backbone was rigid enough that movement of the 101 C does not merely cause the backbone to move further into the binding site, but instead displaces the whole protein relatively intact. The same force was applied in additional simulations with several of the backbone loop atoms fixed on the opposite side of the loop from residue 101 (residue 93 C, and carbonyl C; residue 94 backbone N, C, and carbonyl C; and residue 95 backbone N, and C), to prevent a simple displacement of the protein in the solvent. The external force was insufficient to place the coordinating oxygens of the aspartate sidechain close enough to the Ca²⁺ for bidentate coordination. Instead, a harmonic movement of the sidechain into and back out of the binding loop resulted, but the sidechain carboxylate oxygen that remained farthest from the Ca²⁺ never moved in closer than 2.8 Å to the ion (Fig. 7). The RMSD from the initial coordinates was always below 1 Å for the Force 0 simulations.

View larger version (17K): [in this window] [in a new window]   FIGURE 7   These graphs represent the Force 0 simulation where several of the backbone atoms were held fixed on the opposite side of the loop from residue 101. (A) The noncoordinating sidechain oxygen of aspartate 101 never comes closer than 2.8 Å to the Ca2+ ion in the Force 0 simulations. (B) The coordinating sidechain oxygen remains ~2.1 Å from the Ca2+.

A similar result was observed in the Force 1 simulations where harmonic constraints were used to drive the sidechain atoms of aspartate 101 toward coordinates that would place both of the carboxylate oxygens within binding distance of the Ca²⁺ ion. However, even in the presence of such constraints, the second sidechain oxygen never gets closer to the bound ion than 2.7 Å (Fig. 8). This Ca-O distance of 2.7 Å is slightly closer than that achieved in the Force 0 simulations, but no significant rearrangement of the loop backbone is required to obtain this distance (Fig. 9 A). Although bidentate conformation was not achieved during 6 ps of simulation (3000 timesteps), additional trials showed that increased simulation time did not help. In fact, during the 6-ps simulation, the noncoordinating sidechain oxygen moved closest to the ion during the first 1.5 ps of simulation (Fig. 8 A). This result was yet another indication that the loop backbone strongly resists moving into the binding site any farther than the movement seen in the Aspartate simulation. The Force 1 maximum RMSD from initial coordinates is less than 1 Å in each simulation.

View larger version (27K): [in this window] [in a new window]   FIGURE 8   In the Force 1 simulation, harmonic constraints were used to drive the sidechain atoms of aspartate 101 toward coordinates that would place both of the carboxylate oxygens within binding distance of the Ca2+ ion. Only one oxygen (B) coordinates the Ca2+ ion, however. The second sidechain oxygen (A) never gets closer to the bound ion than 2.7 Å.

View larger version (30K): [in this window] [in a new window]   FIGURE 9   The simulations were all performed using the entire protein as the model; however, only the backbone and coordinating sidechains from residues 92 to 101 of the EF loop are shown in these stereo illustrations for clarity. Superpositions were performed by setting the Ca2+ ion as the origin and applying the rotation matrix calculated by the Kabsch algorithm [Kabsch, 1978] to the remaining atoms. (a) The progression of the Force 1 simulation is shown in this stereo superposition of frames, where the initial frame is shown in green, the 400-fs frame is shown in yellow, and the frame at 1240 fs is blue. The noncoordinating oxygen of aspartate 101 comes as close as 2.85 Å at 1240 fs, but the loop backbone does not need to rearrange significantly to allow this. Aspartate 101 never achieves bidentate coordination of the Ca2+, despite the fact that constraints have been placed on the sidechain heavy atoms to bring both sidechain oxygens within coordinating distance of the ion. (b) The progression of the aspartate 101 sidechain oxygens toward the Ca2+ ion in the Force 2 simulation is shown here. The initial frame is shown in green, the 5-ps frame is shown in red, the frame at 152 ps is yellow, and the frame at 157 ps is blue. Application of constraints reducing the distance between the C atoms of residues 94 and 101 by 1.5 Å necessitates considerable constriction of the loop backbone. The backbone nitrogen of residue 100 at 157 ps is ~1.5 Å closer to the Ca2+ ion than the same atom at 0 ps. The -carbon of residue 92 moves ~2 Å closer to Ca2+ ion between 0 and 157 ps. The second sidechain carboxylate oxygen of aspartate 92 appears to swing out of the plane of the page as the loop constricts, but this orientation aligns better with the orientation of aspartate 92 in the crystal structure than does the simulation's starting orientation, as is shown in part c. (c) Frames representing the shortest Ca-O distances of the Force 1 (green) and Force 2 (blue) simulations are superimposed here on the PVEF-E101D/Ca2+ crystal structure coordinates (red). The Force 2 simulation, which applies distance constraints that pull the residue 94 and 101 -carbons 1.5 Å closer together, shows some significantly different backbone conformation from the other two pictured. There are no intramolecular clashes or breaks in bonds in the Force 2 model, however, implying that the aspartate is physically able to achieve bidentate coordination of the Ca2+ ion, it is just not energetically favorable for the loop to constrict to the degree necessary for bidentate coordination to occur.

In a final simulation, entitled Force 2, the backbone was successfully compelled to move close enough to the Ca²⁺ ion to allow bidentate coordination by aspartate 101 (Fig. 10). This was achieved by applying distance constraints between the C atoms of residues 94 and 101, two residues that lie directly across the binding loop from each other. The constraints required that the distance between the two atoms be reduced by 1.5 Å. The change in the total energy of the model as a result of forcing bidentate Ca²⁺ coordination was too small to distinguish from the normal energy fluctuations. Each force-field energy expression was also individually evaluated over time, but again changes in the various energies are too small to distinguish from normal energy fluctuations. The RMSD from the initial structure increased to ~2.2 Å, but this was expected given the dramatic nature of the distance constraints applied. The Force 2 simulations were stable and the protein remained folded and well-behaved throughout.

View larger version (30K): [in this window] [in a new window]   FIGURE 10   The Force 2 simulation successfully achieved bidentate Ca2+ coordination by aspartate 101 by applying distance constraints between the C atoms of residues 94 and 101, two residues that lie directly across the binding loop from each other. The total simulation time was 240 ps (2400 frames), where the external force was applied slowly during the first 120 ps, the simulation was allowed to equilibrate for another 20 ps, then the force was slowly released during the last 100 ps of simulation. (A) The noncoordinating aspartate 101 sidechain oxygen binds Ca2+ near the end of the equilibration period (frame 1400). (B) The coordinating oxygen shows typical coordination distance throughout the Force 2 simulation.

The loop backbone constricts rather dramatically as a result of the Force 2 distance constraints (Fig. 9, B and C), and aspartate 101 successfully achieves bidentate coordination of the Ca²⁺ ion. There are no significant intramolecular clashes or breaks in coordination bonds in the Force 2 model and no rearrangement of nearby sidechains is necessary to allow the loop constriction, so it evidently is not direct steric clashes that prevent bidentate coordination. The aspartate is shown to be flexible enough to orient itself properly for bidentate coordination, it appears that the simulations are not able to distinguish the force that prevents the loop from constricting to the degree necessary for bidentate coordination to occur.

	CONCLUSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSION REFERENCES

These MD simulations were used to investigate some of the putative mechanisms the EF-hand binding loop uses to determine metal ion affinity and specificity. The first set of simulations, the Alchemy simulations, accurately and reversibly predicted that the glutamate at position 12 is the source of bidentate ligation of Ca²⁺. This outcome not only verified that our model was appropriate for performing additional simulations, it provided evidence that EF-hand loop structure, to some extent, imposes this role in selectivity on the residue at position 12. This result was construed from the fact that the electronic properties of glutamate and aspartate are described by identical electrostatic charge distributions in the model system, thereby discounting the idea that the glutamate at position 12 may be providing a slightly different electrostatic environment than the other carboxylate sidechains. Additionally, the CD binding site contains an additional glutamate, the coordinating residue at position 9, that could offer bidentate ligation of Ca²⁺ if the greater length and flexibility of a glutamate residue over an aspartate were all that was required. It is likely that the residue at position 12 is assigned the bidentate to monodentate transition because 1), the glutamate oxygens are in a favorable geometric position with respect to the other coordinating oxygens in the binding loop for bidentate coordination and 2), the last coordinating residue of the loop is positioned at the outermost edge of the loop, with only water molecules available for hydrogen bonding, and thus, the sidechain carboxylate oxygens are very available to the bound metal ion. In contrast, the coordinating sidechains at loop positions 1, 3, 5, and 9 are in a more crowded region, protected within the binding loop where they can easily make intramolecular hydrogen bonds with nearby residues.

An additional and unexpected piece of information that came from the Alchemy simulation was the implication that a coordinating serine in an EF-hand binding site may play a role in metal ion selectivity through its need to satisfy sidechain hydrogen bonding requirements. Our model displays a trend for the serine at position 5 of the CD binding site to move away from the Mg²⁺ ion during MD simulations. It is possible that this result is an artifact of the imprecise charge distribution representation for serine in the MD, and further investigation would be required to discern definitively whether certain EF-hands might exploit serine to discriminate between Mg²⁺ and Ca²⁺.

The Aspartate simulation produced results that correlate well with the experimental result that an E101D substitution at EF loop position 12 resulted in monodentate Ca²⁺ coordination. The simulation was intended not only to observe whether our model would predict sixfold, octahedral binding geometry in agreement with the E101D crystal structure, but also to scrutinize the loop backbone movement that would be required for the aspartate in our model to offer bidentate ligation. It is our postulate that the tug-of-war between the attraction of the carboxylate sidechain oxygens for Ca²⁺, and their requirement to drag the loop backbone along with them as they move in toward the metal ion, to a large degree defines the range of coordination spheres that can be achieved in a particular EF-hand site.

This postulate is supported by the results of the Force simulations where it was seen that the loop backbone is quite resistant to the degree of constriction required to allow an aspartate at position 101 to obtain bidentate Ca²⁺ coordination. However, this constriction can be imposed on the loop through distance constraints without significant rearrangement of nearby sidechains and without generating intramolecular clashes or bond breakages. The Force simulations illustrate that the aspartate is capable of attaining a suitable orientation for bidentate coordination, thus implying that it is the inherent rigidity of the loop that prevents bidentate coordination in the parvalbumin E101D mutant.

The EF-hand binding site is extraordinary in its ability to discriminate between two small cations, similar in charge, size, and electronic configuration. This study has shown that these binding sites exploit subtle, but crucial, differences in the properties of all the pertinent components involved in metal ion binding. These properties include those dictated by the structural configuration of the binding loop, properties of the amino acid sidechain moieties involved in metal ion coordination, and properties of the metal ions themselves.