The Structure of Neuronal Calcium Sensor-1 in Solution Revealed by Molecular Dynamics Simulations

Neuronal calcium sensor-1 (NCS-1) is a protein able to trigger signal transduction processes by binding a large number of substrates and re-shaping its structure depending on the environmental conditions. The X-ray crystal structure of the unmyristoilated NCS-1 shows a large solvent-exposed hydrophobic crevice (HC); this HC is partially occupied by the C-terminal tail and thus elusive to the surrounding solvent. We studied the native state of NCS-1 by performing room temperature molecular dynamics (MD) simulations starting from the crystal and the solution structures. We observe relaxation to a state independent of the initial structure, in which the C-terminal tail occupies the HC. We suggest that the C-terminal tail shields the HC binding pocket and modulates the affinity of NCS-1 for its natural targets. By analyzing the topology and nature of the inter-residue potential energy, we provide a compelling description of the interaction network that determines the three-dimensional organization of NCS-1.

The structure of NCS-1 is characterized by helix-loop-helix motifs, called EF-hand motifs, which can bind Ca 2z ions. There are 4 EF-hand motifs (Figure 1), namely EF1, EF2, EF3 and EF4, but EF-1 is not able to bind Ca 2z ions [3,4]. The N-terminal glycine is myristoylated in vivo through the action of Nmiristoyltransferase. It is widely accepted that the miristoyl chain confers the capability to the NCS family proteins to link to the membrane through the ''myristoyl switching'' mechanism, [3,5,14,15] which, eventually, can be modulated by changing the Ca 2z concentration [16]. In the opposite position of the EFmotifs there is a central and wide hydrophobic crevice (HC), that can interact with various physiological targets activating specific biological processes [2,6,7,[16][17][18].
Recently, Heidarsson et al. [19] determined the NMR structure of the non-myristoylated calcium-bound NCS-1. The solution structure shows that the N-domain is more flexible than the C-domain, and in the absence of a binding partner the hydrophobic C-terminal residues partially occupy the HC.
The NMR experiments, therefore, suggest that the C-terminal segment of NCS-1 may dock into the hydrophobic crevice and act as a ligand-mimic in the absence of an interaction partner, conferring conformational stability to the NCS-1 structure. The structure proposed by Heidarsson et al. is free from polyethylene glycol (PEG) molecules and not constrained by crystal packing interactions. Note, however, that the resolution of the C-terminal region is rather poor, compromised by the small number of NMR signals available in this region [19].
In the crystal structure of NCS-1 [4], the C-terminal tail is instead outside the crevice, suggesting a poor affinity between the C-terminal tail and the HC. In this structure the HC is occupied by PEG molecules, mimicking the position of the C-terminal tail observed in the NMR structure and likely inhibiting the binding of the C-terminal tail in the HC.
To investigate the origin of the differences between the solution and crystal structures and to clarify the biological role of the Cterminal, long molecular dynamics simulations in explicit solvent were performed for the non-myristoylated Ca 2z -bound NCS-1. Specifically, we carried out two simulations: (i) one starting from the NMR structure and (ii) the other starting from the protein crystal structure modeled without the PEG molecules docked in the HC.
During both simulations the C-terminal tail docked into the HC in a very similar manner. This, in turn, suggests that the Cterminal can dock into the HC acting as an auto-inhibitor, by blocking substrate access to the HC and conferring stability to the protein by shielding the HC from the environment.
To find the critical interactions that stabilize the protein structure, we evaluated the inter-residue interaction potential energy maps including all pairs of residues. The analysis of the inter-residue maps highlights the presence of several salt bridges that interconnect the secondary structures, thus modulating the stability of the NCS-1 tertiary structure [20][21][22]. Based on this discovery, we suggest a way to modify the stability of the NCS-1 by mutating few crucial residues.

Substantial differences exist between X-ray and NMR structures
The crystal and solution structures are shown in Figure 1. We carried out a quantitative comparison of the two structures by evaluating total and partial root mean square deviations (RMSDs). This analysis reveals striking differences between the two structures.
The main difference concerns the conformation of the L3 segment. In the X-ray structure, ( Figure 1A) L3 is in helical conformation and outside the HC. In the NMR structure, ( Figure 1B) L3 is partially docked into the HC. Appreciable differences also involve the relative orientation of the ''external'' H1, H2, H8 and H9 helices. The RMSD evaluated for the backbone of residues 11 to 174, the protein core (PC), of the NMR structure with respect to the crystallographic structure is 4.4 Å , which reflects substantial divergence. Partial RMSDs for the loops, a-helices and b-strands between the crystal and solution structures are given in the supplementary information (see Table S1).
Both structures are affected by measurement conditions and may deviate from the inherent protein conformation. On one hand, the X-ray static structure is affected by the presence of glycol molecules [4] and by crystal packing. On the other hand, the NMR solution structure has a poor resolution in the C-terminal segment, because of a small number of NOE signals [19].

Molecular dynamics simulations converge to a unique state
The MD technique has been assessed as a reliable tool to refine experimental structures [23][24][25][26]. Thus, two long MD simulations at the temperature of 310 K were performed starting from the Xray structure (MD-XR) and the NMR structure (MD-NMR).
The structural adjustments of the protein during the simulations were monitored by evaluating the RMSD of backbone residues 11 to 174, i.e. the PC, with respect to the original X-ray crystallographic structure. The RMSD values calculated along the MD-XR and MD-NMR simulations are reported in Figure 2 (gray and orange curves, respectively). The RMSD during the MD-XR simulation reaches a plateau after 20 ns and remains stable until the end of the simulation (250 ns). The RMSD does not exceed 3 Å and fluctuates roughly between 2 Å and 2.5 Å . These are evidences that the PC did not undergo substantial conformational rearrangements during the MD-XR simulation (see Figure S1 and Table S2 for more details).
Nevertheless, during the first 20 ns of the MD-XR simulation L3 docked into the HC and gained a placement similar to that in the NMR structure. In Figure 3, the distance between the center of mass (COM) of L3 and the COM of the whole protein is plotted: the plot shows a major deviation of L3 from the whole protein at the beginning, with an abrupt transition around 20 ns. This behavior, along with structure visualization at various stages of the simulation (inset of Figure 3), manifests that L3 docked into the HC at the early stage of the simulation and remained trapped into the HC for the entire duration afterwards; the initial a-helix conformation of L3 was lost during the settling of the segment into  In the first 100 ns of the MD-NMR simulation, the RMSD values are larger than 5 Å : this means that during the early stage of the MD-NMR simulation the protein core remarkably differs from the crystallographic structure. Yet, we observe that L3 entered the HC in this early stage and remained trapped into the crevice for the remainder of the time evolution, as in the MD-XR simulation. Between 100 ns and 300 ns, the RMSD values decrease and reach a local basin between 5 Å and 4 Å . After 300 ns the RMSD values settle a plateau around 2 Å . The RMSD trends attest that the solution structure during the MD-NMR simulation clearly underwent major conformational rearrangements (Table S3) and eventually converged to the crystal structure.
In fact, the RMSD between the average structures evaluated over the last 100 ns of each simulation is 1.9 Å , much smaller than the initial 4.4 Å ( Figure S2). Representative structures were extracted from the trajectories by a clustering algorithm applied during the final 100 ns of each simulation (see Materials and Methods). The representative structures of the most populated clusters (also called in the following most representative structures) from MD-XR and MD-NMR are visualized in Figure 4. The two structures show a considerable overlap, which is consistent with the low protein-core RMSD value of 2.1 Å between them. Furthermore, Figure 5 shows that the shapes of the HC are equal in these two representative structures, and that L3 is almost totally accommodated in the groove in both cases.
The RMSD values between any two representative structures of all the clusters obtained from the two trajectories span between 1.9 and 2.8 Å (Table S4), demonstrating that protein conformations in the last 100 ns of both simulations belong to the same state. 80% of the NOE restraints evaluated in the last 100 ns of both simulations was respected, showing a good agreement between experimental and simulation data. More interesting is the evolution of the NOEs violation during the MD-NMR simulation. At the early stage of the MD-NMR simulation, namely between 50 and 150 ns the average violation is 2.2 Å . In the last 100 ns of the MD-NMR simulation, the average violation decreases to 1.9 Å . During the MD-XR the average violation is constant, 1.5 Å . These NOE trends are an additional indication that during the MD-NMR simulation the protein core of the NCS-1 relaxed toward the state of the protein core of the MD-XR.
To inspect in more detail the evolution of the protein conformation during the MD-NMR simulation, the RMSDs for significant protein regions were calculated (Table 1). Such partial RMSD values were calculated for each separate region between  the most representative structure of the MD-NMR trajectory and the corresponding initial structure. The RMSD data reported in Table 1 clearly show that the buried a-helices H2, H3, H4, H5, H6 and H7, as well as the b-strands involved in the EF motifs, undergo minor rearrangements (small RMSD). As expected, the random coils L3 and L2 have the highest RMSD values. The high RMSD value for L3 is consistent with the docking of L3 into the HC at the early stage of the time evolution, as noted above. The high RMSD value for L2 is due to the high mobility of this segment. The large RMSD value of 7.3 Å for H1 suggests that this segment is very mobile. The H8 and H9 helices also have a fairly large RMSD value of 5 Å . Indeed, these helices underwent an appreciable conformational rearrangement during the dynamics. In particular, during the MD-NMR simulation H9 was ''wedged'' between H8 and H6, further closing the edge of the HC, whereas H8 was reoriented to extract the L2 segment from the HC ( Figure 6).
A test simulation, carried out at 350 K starting from the last structure of the MD-NMR simulation, shows that the L3 segment remains firmly docked into the HC and the PC does not undergo any major conformational changes ( Figure S3). These evidences confirm that during the MD-NMR simulation the system converged to a stable state.
To gain insights into the importance of specific residues in stabilizing the NCS-1 structure, pairwise residue-residue potential interaction energies were calculated for both simulations during the respective final 100 ns and then averaged (see Material and Methods). The result of this procedure for each trajectory is the inter-residue energy map (IEM), which reflects the basic interactions that stabilize the native structure of a proteins [27,28].
The computed IEMs are illustrated in Figure 7. The interresidue energy values range between 2120 kcal/mol and +20 kcal/mol. The strongest attractive interactions are shown in darkblue (negative values), whereas the most unfavorable interactions (positive values) are shown in red. The yellow and light-green zones span an energy interval between 22 and 0 kcal/mol and are mainly due to weak dispersion and hydrophobic interactions. The green zone spans the interval between 215 and 22 kcal/mol: this range is typical of hydrogen bonds. The cyan and light-blue zones   Partial RMSD evaluated after structure alignment over the backbone atoms of protein core residues. RMSD values for the helices H1 to H9, the b-strands b1 to b6, and the loops L1, L2, and L3 between the most representative structure of the MD-NMR trajectory and the corresponding initial structure. b2 and b5 are sub-structures of the L1 and L2 segments respectively. doi:10.1371/journal.pone.0074383.t001 span the interval between 230 and 215 kcal/mol and highlight strong polar interactions or weak ionic interactions. Dark points reflect ionic interactions. The green points that form straight segments on the two sides of the diagonal mainly reflect the 1-4 hydrogen bonds that define the secondary structure of the ahelices. The lack of interactions along the very diagonal is due to the loops (L1, L2 and L3) and to the presence of the calcium binding EF motifs. The calcium ions are coordinated in the binding pocket by negatively charged residues. Because the calcium ions were not considered in the potential energy calculation, only the repulsive interactions between negatively charged residues are reflected in the IEMs. The two maps in Figure 7 highlight the same interactions and same topology, thus yielding another proof that the two simulations fell into the same stable conformation.
The ovals in Figure 7 marks the interactions that mainly define the topology of the IEMs. The N1 and N2 zones refer to interresidue interactions between residues of the N-domain. More specifically, the N1 zone embodies the interactions between H1 and H2 (H1-H2), whereas N2 expresses the H3-H5, H1-H5 and H2-H5 interactions. The small circle inside the N2 oval includes the b1{b3 interaction. The NC zone refers to interactions between the N-domain and C-domain that involve the H5-H6, H4-H6 and H4-H7 helix pairs. The C1 and C2 zones refer to the inter-residue interactions between residues of the C-domain. The C1 zone highlights the H6-H9 and H7-H8 interactions and the internal circle in C1 pertains to the b4{b6 interaction. The C2 zone contains the H9-H8 interaction. The rectangular C3 zone highlights inter-residue interactions between residues of L3 and residues of the PC that shape the HC, and show that the L3 segment interacts with the NCS-1 PC mainly through hydrophobic interactions. In this structural portion the Asp187 residue could interacts with Lys100 or Arg94 than with other residues, eventually anchoring the very terminal tail of the L3 segment to the HC.
Among the various point interactions, we focus our analysis on the strongest and conserved interactions in both IMEs, i.e., the salt bridges (see also the filtered maps in Figures S4 and S5). The conserved salt bridges with energy lower than 230 kcal/mol are marked with diamonds in Figure 7. Let us start our comments from intra-helix salt bridges located in H1 and H6. In H1 the salt bridges are Glu14:Arg18, Glu15:Arg18 and Glu15:Lys19, whereas in the H6 they are Arg102:Glu99 and Arg102:Asp98. It is worth nothing that Asp98 of H6, however, can also interact with the Lys174 residue of H9 to accomplish an inter-helix salt bridge between H6 and H9. The stabilization of the inter-versus intrahelix salt bridge involving Asp98 (i.e., Asp98:Lys174 versus Asp98:Arg102) relies, therefore, in the possibility of the Arg98 to interact with Arg102 or Lys174. In this sense, we propose that a mutation that consists in replacing Arg102 with an uncharged residue could affect the mobility of H9 and L3, compromising the docking process of L3 into the HC. In fact, it is known that the substitution of Arg102 with glutamine affects the structural and dynamical features of the C-terminal tail, providing a molecular counterpart to the hypothesis that the Arg102Gln mutation is implicated in the autism disease [10,19,29].
Other salt bridges exist between residues that belong to distinct protein segments. The conserved Arg151:Glu141 and Arg151:-Glu142 bridges, for example, connect L2 to H8. The Ar-g118:Asp150 bridge connects H7 to H8 and the Glu171:Lys158 bridge connects H8 to H9. The Asp44:Arg79 bridge stabilizes the b1-b3 interaction, whereas the Glu26:Arg94 bridge connects H2 to H5. The Glu26 and Arg94 residues interact in an ''end-on'' manner [21], thus suggesting a key role in the conformational specificity of the N-domain. The Lys63:Asp123 and Lys63:Asp126 bridges involve one residue of the N-Domain and two residues of the C-domain: these two salt bridges can coexist or exist exclusively; they start from Lys63 and, bifurcating towards Asp123/126, are able to connect the H4 helix of the N-domain to the H7 helix of the C-domain.

Discussion
The biological function of the NCS protein family is related to the ability of interacting with different target proteins and regulating their action. NCS proteins bind target ligands trough a conserved hydrophobic binding pocket, which can differ in shape and size to regulate the target specificity [1,3,17,18,29]. However, to guarantee the target specificity, each member of the NCS family seems to use multiple regulatory mechanisms. The biological function of recoverin for example is strongly related to the myristoyl switch mechanism induced by the Ca 2z concentration. Weiergrä ber et al. [30] showed that truncating the last 12 basic residues of the C-terminus altered the Ca 2z affinity of recoverin, demonstrating that C-terminal tail acts as a ''internal modulator of Ca 2z sensitivity''. In particular, they observed that in truncated recoverin (Rc 2{190 , PDB code 2HET), the Cterminal tail is docked in the hydrophobic pocket, resembling the binding mode of L3 to the HC in NCS-1. On the contrary, experimental evidences indicate that the biological function of NCS-1 is not strictly regulated by the Ca 2z concentration and/or myristoyl switch mechanism as in recoverin [1,30], corroborating the hypothesis that the hydrophobic nature of the C-terminal tail of NCS-1 is the main regulator of the NCS-1 specificity (in fact, the C-terminal tail of the NCS-1 is hydrophobic, whereas the Cterminal tail of recoverin is hydrophilic [17]).
The molecular dynamics simulations presented here assess a consensus structural state, which is reached from either the NMR solution structure or the X-ray structure, if the PEG molecules are removed from the HC in the latter. The consensus structure has a core similar to the X-ray structure, while the location of the Cterminal inside the HC is akin to the NMR solution structure: we infer that in the experimental crystal phase the consensus state is inhibited by the presence of PEG molecules that occupy the HC. Our results suggest that the L3 location is sensitive to the environmental condition and in the absence of a substrate it has the role of stabilizing the protein fold by segregating the HC from water. These findings are in agreement with NMR experimental evidences [19] and confirm the hypothesis that the C-terminal residues occupy the hydrophobic crevice as a ligand mimic. Our simulations also grant a dynamical description of the docking of L3 in the HC: L3 and HC interact mainly by hydrophobic interactions, eventually stabilized by salt-bridge interactions. The competition between the intra-molecular ligand mimic L3 and the proper inter-molecular protein target, therefore, could establish a dynamical mechanism to modulate the target specificity of NCS-1. Moreover, the affinity of L3 towards the HC could facilitate the extrusion of the myristoyl from the HC during the membrane binding and prevent the binding of the myristoyl chain to the HC.
The analysis of inter-residue interaction maps has evidenced the presence of various salt bridges that connect different secondary structure elements, thus contributing to determine the protein tertiary structure. Assigning a biological function to an individual salt bridge is a controversial issue. However, it is a widespread idea that the stability of the salt bridges is sensitive to the environment and that they have an important role in stabilizing/destabilizing protein structures [20][21][22]31,32]. Therefore, the identified salt bridges define residues whose mutation may lead to profound effects on the NCS-1 stability and biological function. In particular, we have remarked the salt-bridge network involving Arg98, Arg102 and Lys174, which could affect the mobility of H9 and L3, providing a possible rationalization on how the Arg102Gln mutation could be connected to the autism disease [10].
Two different experimental structures were used as the starting points for two independent MD simulations of the NCS-1 protein.
One simulation started from the crystal structure [4], chain B of the PDB 1G8I, after removing all crystallographic waters and residual chain glycols located in the HC. Another simulation started from the first of the NMR structures deposited by Heidarsson et al. [19], PDB code 2LCP. The simulations were conducted using periodic boundary conditions (PBC). The water layer separating the protein from any of its periodic images was at least 14 Å thick in each simulation box. To guarantee neutrality 3 Na z ions were placed randomly at a distance larger than 5 Å from the protein. The bonds between hydrogen and heavy atoms were constrained with SHAKE [37]. The r-RESPA multiple time step method [38] was employed with 2 fs for bonded potentials, 2 fs for the short-range part of non-bonded potentials, and 4 fs for the long-range part of the non-bonded potentials [33]. The long-range part of the electrostatic interaction was treated with the Particle-Mesh-Ewald (PME) method [39]. The distance cut off for nonbonded interactions was set to 12 Å , and a switching function was applied to smooth interactions between 10 and 12 Å .
All simulations were conducted in the NPT ensemble. The temperature was set to 310 K and regulated via a Langevin thermostat [40]; the pressure was set to 1 atm and regulated via an isotropic Langevin piston manostat [41].
Each simulation was preceded by a preparation phase consisting of both geometry optimization and constrained MD. 2000 steps of conjugate gradient geometry optimization were initially conducted with harmonic restraints on the protein heavy atoms using a force constant of 1 kcal/(mol : Å 2 ). After energy minimization, the system was simulated for 1 ns: in the first 0.5 ns harmonic restraints were applied on the protein heavy atoms with a force constant of 1 kcal/(molÅ 2 ); in the last 0.5 ns the force constant was scaled to 0.1 kcal/(mol : Å 2 ). Subsequently 25000 steps of minimization without restraints were performed. The final structures were then used as the distinct initial conditions for unrestrained molecular dynamics.
The MD-XR simulation that started from the X-ray structure was 250 ns long. The MD-NMR simulation that started from the NMR structure was 525 ns long. Two further test simulations starting from the crystal structure and from the last structure of the MD-NMR simulation were conducted at 350 K; both simulations were carried out for 50 ns.

Analysis
Root mean square deviations (RMSD) were evaluated after the structure alignment of the backbone residues 11 to 174, the protein core.
Cluster analysis was performed with WORDOM [42], by comparing the RMSD values for the PC calculated at every 50 ps over the last 100 ns of the MD simulation. The cut off was set equal to 2 Å and clustering was accomplished by using a hierarchical algorithm.
Nuclear overhauser enhancemnt (NOE) interproton distances were calculated from the simulation and compared to the experimental data [19]. Due to the presence of multiple equivalent proton definitions (e.g., interproton distances between two methyl groups), an effective interproton distance was evaluated by summing each distance between pairs of proton weighted by the sixth power; interproton distances were evaluted for each structure and averaged over time by the sixth power [43]. A NOE distance was considered violated if it exceeds the upper distance limit increased by 0.6 Å ; each distance violation (Vij) was defined as difference between the NOE distance and the upper limit, whereas the average violation vVw was evaluated by averaging all the Vij values over the number of total violated NOEs.
Inter-residue interaction potential energies for each pair of residues a and b were evaluated. The long-range part of the electrostatic interactions was neglected and the interaction energy between residues with indexes a and b was calculated when the condition (bwaz3) was satisfied. The calcium ions were not included in the calculation. To account for the fluctuations of the protein, the pair interaction energies were calculated every 50 ps  Figure 1, is reported along the axes. Interactions between L3 and the PC are highlighted by a rectangle. Small circles in N2 and C1 highlight the b1-b3 and b4-b6 interactions. The conserved salt bridges with energy lower than 230 kcal/mol are marked with diamonds. doi:10.1371/journal.pone.0074383.g007 over the last 100 ns of the MD simulation and averaged over the time. Figure S1 Comparison between the most representative MD-XR structure (light color) and the crystallographic structure 1G8I (dark color). The color code that distinguishes the various protein segments is the same as used in the manuscript. The most representative structure pertains to the final 100 ns of the MD-XR simulations. Alignment was performed over backbone atoms of residues 11 to 174. The aligned structures do not show particular structural differences in the orientation of the ahelices. The most significant difference concerns the location of the L3 segment. L3 is external to the HC in the crystal structure (used as the starting point of the MD-XR dynamics), whereas it is docked in the HC in the most representative dynamical structure. See Table S2 for more details.  Table S2 Partial RMSDs evaluted between the representative structure of the most populated cluster of the MD-XR and the crystallographic structure. The partial RMSDs clearly show that the initial crystallographic structure and the representative structure of the most populated cluster obtained from MD-XR simulation do not differ. Alignment was performed over backbone atoms of residues 11 to 174. (TIFF)

Supporting Information
Table S3 Partial RMSDs in Å evaluated between the most representative MD-NMR structure and experimental NMR structures. The labels S1-S6 mark the first six structures contained in the pdb file 2LCP. The partial RMSDs clearly show that the most representative MD-NMR structure differs substantially from the pdb structures, not only in the mobile loops but also in some helices and strands. This is true not only for the structure S1 that we used as a starting point for the MD-NMR run, but also for other deposited structures, meaning that our results and conclusions do not depend on the arbitrary initial condition of the MD-NMR simulation. Alignment was performed over backbone atoms of residues 11 to 174. (TIFF) Table S4 Backbone RMSDs between representative structures from the MD-XR and MD-NMR trajectories. Cluster analysis was performed in the last 100 ns of the MD-NMR and MD-XR trajectories: 10 clusters were found for MD-NMR and 3 clusters were found for MD-XR. Each cluster contains similar structures and is represented by one such structure. RMSD values are in Å . The order of the clusters is related to the population, namely the number of snapshots that it contains: cluster ''1'' is the most populated. The values in this Table clearly show a high similarity between the representative structures of both trajectories during the final 100 ns when the consensus status has been attained. Alignment was performed over backbone atoms of residues 11 to 174. (TIFF)

Author Contributions
Conceived and designed the experiments: LB SC RDF EP. Performed the experiments: LB. Analyzed the data: LB SC RDF EP. Contributed reagents/materials/analysis tools: LB SC RDF EP. Wrote the paper: LB RDF EP.