Calmodulin (CaM) is a calcium sensing protein that regulates the function of a large number of proteins, thus playing a crucial part in many cell signaling pathways. CaM has the ability to bind more than 300 different target peptides in a Ca2+-dependent manner, mainly through the exposure of hydrophobic residues. How CaM can bind a large number of targets while retaining some selectivity is a fascinating open question. Here, we explore the mechanism of CaM selective promiscuity for selected target proteins. Analyzing enhanced sampling molecular dynamics simulations of Ca2+-bound and Ca2+-free CaM via spectral clustering has allowed us to identify distinct conformational states, characterized by interhelical angles, secondary structure determinants and the solvent exposure of specific residues. We searched for indicators of conformational selection by mapping solvent exposure of residues in these conformational states to contacts in structures of CaM/target peptide complexes. We thereby identified CaM states involved in various binding classes arranged along a depth binding gradient. Binding Ca2+ modifies the accessible hydrophobic surface of the two lobes and allows for deeper binding. Apo CaM indeed shows shallow binding involving predominantly polar and charged residues. Furthermore, binding to the C-terminal lobe of CaM appears selective and involves specific conformational states that can facilitate deep binding to target proteins, while binding to the N-terminal lobe appears to happen through a more flexible mechanism. Thus the long-ranged electrostatic interactions of the charged residues of the N-terminal lobe of CaM may initiate binding, while the short-ranged interactions of hydrophobic residues in the C-terminal lobe of CaM may account for selectivity. This work furthers our understanding of the mechanism of CaM binding and selectivity to different target proteins and paves the way towards a comprehensive model of CaM selectivity.
Calmodulin (CaM), Fig 1a and 1b, is a promiscuous Ca2+-sensing protein that plays part in many physiologically important cellular processes [1]. Its two lobes, connected by a flexible linker, have one beta sheet and two EF-hand motifs each, Fig 1a–1c. The EF-hand binds a Ca2+ ion which induces tertiary structure rearrangements of lobe helices, exposing hydrophobic residues [2]. This allows CaM to bind and regulate a myriad of target proteins such as ion channels, kinases and G-protein coupled receptors. The Ca2+-signaling and olfactory transduction pathways (Fig 1d and 1e) are two examples of cell-signaling pathways where CaM is involved. In the Ca2+-signaling pathway, CaM activates and regulates the myosin light chain kinase IV (MLCK) and calcineurin (CaN), among others [3]. Through this, CaM plays part in regulating a variety of biological processes including metabolism, proliferation and learning. In olfactory transduction, CaM’s role is instead to inhibit the CNG channel and activate the Ca2+/CaM-dependent protein kinase (CAMK) which inhibits the adenylate cyclase 3 (ADCY3). CaM’s promiscuity can be mainly ascribed to its ubiquity and interactions to specific target proteins depend partly on local cell environments such as availability of target proteins. Moreover, high affinity binding to target proteins is partly due to the flexible linker that allows wrapping around target proteins but also stems from the structural properties of the two lobes, thus linked to the pockets formed by their hydrophobic interfaces [4]. Here, we study the interactions between CaM and various target proteins, including proteins that are involved in the Ca2+-signaling and olfactory transduction pathways. Ca2+-independent CaM regulation, with Ca2+-free CaM (apo CaM) regulating the target protein may occur, although not as frequently as Ca2+-dependent CaM regulation. The IQ binding motif, named after the first two conserved residues of the motif, is often associated to apo-CaM binding while 1-5-10 or 1-5-8-14 motifs, named after the position of conserved hydrophobic residues, are associated with Ca2+-bound CaM (holo CaM) binding. [5–7] These motifs are, however, lacking the complexity to completely distinguish holo-CaM and apo-CaM binding, as well as different types of binding. An example is seen in Fig 2, where apo and holo CaM C-terminal domain (C-CaM) expose different hydrophobic interfaces, but both bind the voltage-gated sodium channel NaV1. 2 at the same IQ binding motif. Holo CaM often interacts with a higher affinity, but in some cases, such as the IQ motif of Nav1. 2, apo-CaM binding is more favorable than holo-CaM binding [8,9]. Both CaM dynamics and motif-dependent target-protein binding have been extensively studied via both experiments [2,8, 10–25] and simulations. Indeed, all-atom molecular dynamics (MD) simulations allow for a detailed description with molecular insights. Early MD simulations of the holo CaM conformational ensemble have shown the holo-CaM linker to be flexible [26–28] and indicated structural destabilization when removing Ca2+ from one site in the N-term domain [29]. Moreover, CaM target-protein binding was approached by studying the dynamics as well as the thermodynamics of CaM or CaM complexes [4,26,30], with insights into binding of specific target proteins using both regular MD and metadynamics [30]. However, little is known about the conformations of the binding interfaces and different binding modes of the two lobes, as well as their relation to protein-unbound CaM conformational ensemble. Part of this was addressed through the analysis of multiple aggregated MD simulations of C-CaM using Markov state modeling (MSM), in which target protein binding was proposed to occur by C-CaM adopting the bound conformation before binding to the protein [31]. However, the analysis did not cover the differences of CaM N-terminal domain (N-CaM) and C-CaM dynamics, mechanisms and binding modes, nor did the study explore the possibility of different binding mechanisms linked to different binding modes. Binding a ligand can be conceptualized using two frameworks: conformational selection and induced fit. In pure conformational selection, the apo protein adopts a holo-like state before binding [32]. In pure induced fit, the ligand binds in a typical apo state that is not ideal, which causes subsequent rearrangements before reaching a high-affinity holo state [33]. In reality, binding likely involves a combination of the two mechanisms. However, spectroscopy experiments, as well as extensive MD simulations, may shed light on which mechanism is dominating. If the apo protein samples the holo-like state, conformational selection is typically assumed to be dominating, otherwise induced fit is assumed [34]. Here, we analyze thermally enhanced MD simulations of calmodulin with different Ca2+-occupancy and use spectral clustering to elucidate calmodulin selective promiscuity. We search for indicators of conformational selection by mapping solvent exposure of residues from sampled states (clusters) to contacts of already existing structures of CaM/peptide complexes. Moreover, we gain knowledge about the characteristics behind different binding modes of the two lobes, as well as the difference between holo and apo binding modes.
For this project, we considered four different binding states: holo and apo calmodulin, as well as Ca2+ bound only in the N-term (N-holo) and Ca2+ bound only in the C-term (C-holo), Table 1. The simulations of N-holo, C-holo and holo CaM used structure 3CLN [35], while the apo simulations used structure 1LKJ [36]. N-holo was generated by removing Ca2+ from C-CaM and C-holo by removing Ca2+ from N-CaM from the holo structure. The systems were built using Charmm-gui [37,38], where the protein was solvated in a box of about 21000 water molecules. The systems were then ionized with 0. 15 M NaCl. Charmm36 was chosen as force field [39], and TIP3P [40] as water model. The modified parameters of Charmm27 force field from Marinelli and Faraldo-Gomez [41] were used to parameterize Ca2+ ions. 5000 steps of minimization were carried out, followed by a 50 ps NVT ensemble equilibration with strong harmonic restraints on the protein atoms. The box was then scaled, relaxing pressure with Berendsen barostat [42], while gradually releasing the position restraints for 350 ps. The MD parameters used in these simulations are extensively described elsewhere [43]. Calmodulin was simulated in an NPT ensemble with a 1 atm pressure and 2 fs time step. The short-ranged electrostatic interactions were modeled with a 1. 2 nm cutoff where the switching function started at 1. 0 nm, and the long-ranged ones with PME [44]. We used a Nose-Hoover thermostat [45], an isotropic Parinello-Rahman barostat [46], and constrained hydrogen bonds with LINCS [47]. The MD simulations were run at a constant temperature of 303. 5 K. In addition to regular MD simulations, temperature enhanced simulations were performed; temperature replica exchange [48] (T-REMD) and replica exchange solute tempering [49–51] (REST). In T-REMD, parallel simulations of independent replicas are run at different temperatures. A random walk in temperature space is achieved by employing the Metropolis criterion periodically to accept coordinate exchanges between neighboring replicas. Conformations obtained at higher temperatures are propagated down to the lower temperatures through exchanges. Because barriers are more easily passed at higher temperatures, the efficiency is increased when the free energy landscape is rugged. However, the number of replicas needed to span a certain temperature interval scales with system size [52]. For this reason, T-REMD may not always be efficient. To alleviate this, REST only modifies the hamiltonian of the system for the solute (protein), and not the solvent [51]. However, the relative efficiencies of regular MD, T-REMD and REST depend on the ruggedness of the free energy landscape [43]. We used 25 replicas in both T-REMD and REST. The replicas of T-REMD spanned a temperature range of 299. 13-326. 09 K, while the REST replicas were simulated at temperatures between 300. 0-545. 0 K. The temperature ranges for REST and T-REMD were chosen using the “Temperature generator for REMD simulations” [53], considering only the protein for REST. Exchanges between neighboring replicas were attempted every 2 ps, where half of the replicas were involved in each attempt. The REST simulations were performed with the Plumed 2. 3b [54] plug-in patched with Gromacs version 5. 1. 2 [55], where the charge of the atoms in the hot region were scaled, as well as the interactions between the two regions and the proper dihedral angles [56]. Analysis was performed using the replica at the lowest temperature. and was carried out on the protein heavy atoms. The first four residues in the apo structure were removed, because those are missing in the 3CLN structure. The apo CaM ensemble is more diffusive and generally well sampled by regular MD, while holo CaM free energy landscape is more hierarchical, and thus sampled more efficiently by temperature enhanced methods [43]. For this reason, more MD data is used in the apo analysis while more REST data is used for holo, Table 1. The three methods may sample different regions with varying efficiency, which is why data from all three simulations is used. Fig 3 illustrates the steps carried out to analyze the MD simulation trajectories. In a first step, the protein was divided into quasi-rigid domains using SPECTRUS [57]. This procedure exploits fluctuations between residues to determine which parts of the protein move together. The clustering and post-processing described hereafter were done on each quasi-rigid domain. To further reduce dimensionality and complexity of the dataset while preserving the most important features, we used spectral clustering [58]. The advantage of spectral clustering compared to regular clustering techniques like k-means is three-fold. First, it is able to accurately assign points to non-convex clusters. Second, non-linear dimensionality reduction is intrinsic to the algorithm. For high-dimensional data, such as MD trajectories, non-linear dimensionality reduction improves clustering by circumventing the curse of dimensionality where the sparsity of the data increases with increased dimensionality [59]. Third, the number of clusters is the same as the number of dimensions onto which the points are projected. This feature becomes advantageous as the number of free parameters is reduced. In spectral clustering, the data manifold is first approximated by a graph with adjacency, or similarity, matrix A. It contains the local relationships between points and is constructed given the matrix of distances d. The distances dij are passed through a radial basis function, or Gaussian kernel, with scaling parameter σ, yielding the graph edge-weights A i j = e - d i j 2 2 σ 2, (1) where dij is the dissimilarity between conformation xi and xj. The geodesic distance between two points is the distance between these points along the manifold, the shortest distance on the graph. The size of the scaling parameter, σ, influences the accuracy of geodesic distances, and should not be too large nor too small. A too large σ would yield short-cuts, thus causing non-convex clusters to be poorly defined. A too small σ, on the other hand, would result in a disconnected graph. Here, σ is the average nearest neighbor distance between configurations. The random walk matrix, related to the Laplacian [58], is then constructed L = D - 1 / 2 A D - 1 / 2. (2) The degree matrix, D, is diagonal with Dii being the degree of node i. The first k eigenvectors are computed and normalized per row to obtain points projected onto the k-sphere. These points are clustered into k clusters using k-means with centers projected onto the sphere. The choice of number of clusters is guided by the maximum eigengap, the difference between two consecutive eigenvalues. The representative structure of a cluster is chosen as the structure with smallest RMSD with respect to the other structures in the cluster. A cluster with all structures represents the conformational heterogeneity of one state. Here, in practice, the dissimilarity between conformations, dij, is measured as the distance between contact maps of inverse inter-atomic distances (iiad-cmap). This general metric is effective for all proteins and, unlike root-mean-square-deviation (RMSD), does not rely on structural alignment. The inverse distances make larger distances small so that far-away motions are cancelled out, without requiring a cutoff. The frame-to-frame distance matrix is compiled by computing the distances between each frame iiad-cmap. Each state, or cluster of frames, was analyzed to provide statistical information about its specific characteristics and molecular features. The features included interhelical angles, secondary structure and importance of states for target protein binding.
We sampled conformations of apo and holo CaM using MD, T-REMD and REST, Table 1. The analysis began by identifying two quasi-rigid domains, namely N-CaM and C-CaM, with SPECTRUS [57]. The identified domains are shown as red and blue CaM ribbons in Fig 3. The clustering and post-processing were therefore carried out on these two domains separately. S1–S4 Figs show the representative structures of holo and apo CaM obtained after clustering on C-CaM and N-CaM. Together with these states, experimental structures of similar conformations are shown for comparison. The molecular features were extracted from the conformational states by computing the interhelical angles and secondary structures. Fig 4 shows the interhelical angles of each state (cluster) from the holo simulations, while Fig 5 shows the interhelical angles for each state of the apo simulations. Each dot is a projection of a conformation in interhelical angle space, colored according to the state (cluster) it belongs to. On top of this, experimentally obtained structures are plotted as squares. The white squares denote structures where CaM is not bound to target peptides, Table 3, while pink squares correspond to structures with CaM bound to a target peptide, Tables 4 and 5. The black circles are mean values of CaM interhelical angles inferred from NMR data [2,24], used for comparison. Furthermore, to allow a comparison between interhelical angles in the holo and apo ensemble, the sampled apo conformations are shown as a gray shadow in Fig 4 (holo) and the sampled holo conformations are shown as a gray shadow in Fig 5 (apo). A simple sanity check, S5–S10 Figs, showed that most states were sampled by all three simulation methods, with a few exceptions. Examples of such are states that are separated by a high free energy barrier and therefore not sampled by the regular MD simulations within these timescales, or states that could not be sampled by the replica exchange simulations because of the relatively short time scales simulated. Similar conclusions can be drawn from the kinetic analysis of MD trajectories (S11 Fig). A change in interhelical angles due to Ca2+-binding can be seen in Fig 4 (holo) and Fig 5 (apo), in agreement with previous work [2,24,26] (see S1 Text for more extensive description and details). We note also that the conformational ensemble generated by MD simulations is broader than the ones observed experimentally for both protein-bound and unbound CaM. In agreement with NMR studies of holo CaM [60], the protein-unbound CaM adopts the conformations of CaM bound to specific target proteins, which is seen as an overlap between the pink squares and the MD data (colored dots) in Fig 4, thus hinting at a possible conformational selection mechanism of target protein binding. Fig 6 displays the difference in secondary structure frequency for each residue compared to holo. Red dots denote the helical, black the strand and blue the coil content. Positive values indicate a gain in secondary structure element compared to holo, and negative ones a loss. The helical break around residue 74-82 [26–28] is seen to different extents in all Ca2+-binding states of CaM, Fig 6. This break occurs in the linker and yields a compact state with the two lobes in contact. This compact state resembles the binding mode where CaM is wrapped around the target protein and was suggested to be an intermediate conformation during target protein- and Ca2+-binding [69,75]. Furthermore, binding Ca2+ rearranges the helices in the lobes and pins the beta sheets between residues 27/63 and 100/136 instead of 28/62 and 99/137, Fig 6. Apo CaM is more flexible, allowing the beta sheets to shift between 26-28/62-64 and 99-101/135-137, Fig 6. This highlights the rigidity of holo compared to apo CaM. Due to the flexibility of apo, a C-CaM state (Fig 5b, state 2) is found where the end of the fourth binding loop is involved in a beta sheet, deforming the binding loop, in agreement with [31]. We hypothesize that this state may inhibit Ca2+-binding. In this state, residues 129-131 are more prone to form a beta sheet with 99-101 and 135-137, S12 Fig. This extra sheet does not form in apo N-CaM because the coiled loop between helices C and D is slightly shorter than the one in apo C-CaM between helices G and H. The flexibility obtained by removing Ca2+ from C-CaM allows these beta sheet formations to occur. To investigate the possibility for conformational selection during the Ca2+-binding process, we searched for potential overlaps in interhelical angle space between apo and holo. We observe holo-like states in apo N-CaM, Fig 5a. (states 3,5 and 8, plotted with colors lime, cyan and cerise), as well as C-CaM holo-like states that have been reported in previous work [31,76], Fig 5b (states 1 and 5, depicted as red and blue). These states could potentially aid Ca2+-binding though conformational selection. The holo N-CaM states 4 and 6 (shown by cyan and magenta dots) overlap with the apo ensemble, Fig 4a. Their secondary structure frequencies are plotted in S13 Fig, which show similar secondary structure frequencies as apo (Fig 6b), with a shift of beta sheet to residues 28/62. The C-CaM interhelical angles of holo and apo overlap more than the N-CaM interhelical angles do. The overlap in interhelical angles and shift of beta sheet in N-CaM, make these states likely intermediate states in Ca2+-binding, involved in conformational selective binding of Ca2+. To assess whether Ca2+-binding allosterically modulates the conformational ensemble of the other lobe, we used the MD, T-REMD and REST datasets of C-holo (Ca2+ in C-CaM) and N-holo (Ca2+ in N-CaM), Table 1. Solvent exposure for each residue in N-holo and C-holo was compared to the apo (Fig 7a and 7c), and holo (Fig 7b and 7d) ensembles. This shows that the Ca2+-free lobe tends to transition to an apo-like conformational ensemble, while the Ca2+-bound lobe stays in the holo-like ensemble. Although there is a small allosteric influence between the lobes, as seen in Fig 6c and 6d, where N-CaM mostly lacks beta sheet in C-holo, the solvent exposure analysis indicates that the overall conformation of one lobe is mostly independent of the conformation of the other. These results indicate that binding of Ca2+ to one of the lobes does not yield a significant population shift in the other lobe and thus that binding is likely not conformationally cooperative between the two lobes. Cooperative binding between sites and lobes has been investigated but the results obtained in different studies have been largely inconsistent, thus making it difficult to validate our results. This is because many different parameter sets in the models used to fit the experimental data perform equally well [77]. To investigate conformational selection aspects of calmodulin binding to target proteins, the relative solvent exposure of the different CaM residues in different states was calculated and mapped to the contacts in a set of CaM complexes, Tables 4 and 5 and Fig 3. The underlying idea is that in a conformational selection mechanism, water molecules around an exposed residue will be replaced by the target peptide. Therefore, the residues that are in contact with the target peptides should be found exposed to solvent in a state involved in a binding mechanism dominated by conformational selection. To then characterize different modes of binding, the distributions of total solvent exposure per state were clustered and divided into classes, as described in the Materials and Methods section. We used molecular dynamics and temperature enhanced MD with an agnostic spectral clustering scheme to search for conformational selection of CaM Ca2+ and target protein binding. We found that the Ca2+-state of one lobe does not significantly influence the conformation of the other lobe, but binding Ca2+ may occur through conformational selection facilitated by a transition state that is visited by both apo and holo CaM. It is observed here by overlapping conformations in interhelical angle space. In N-CaM, the transition between apo and holo likely involves a beta sheet residue shift. The target protein binding modes of apo CaM are more shallow than those of holo CaM, and the charged interactions dominate due to the induced burial of the hydrophobic binding interface from Ca2+-depletion. The residues involved in binding are equally exposed in apo CaM states and thus all observed states are equally likely to initiate binding to the target protein. The notion that not only the linker, but also the binding interface of C-CaM shows configurational flexibility has previously been proposed both through extensive MD simulations [31], but also in structural studies [2,24,25]. The conformation of C-CaM interface is suggested to vary more than the N-CaM interface [25]. Here, holo C-CaM shows distinct states exposing hydrophobic residues that are otherwise screened from water. Such states are absent in N-CaM, which also exhibits less binding heterogeneity than C-CaM in the CaM-complex structures. This is seen as the shallow class of C-CaM show more shallow binding than the shallow class of N-CaM, while the deep binding class of C-CaM shows deeper binding than the deep binding class of N-CaM. The distinct states of holo C-CaM with a clear mapping to bound structures indicate a tendency for C-CaM binding to be selective and dominated by conformational selection, while N-CaM, which lacks this clear mapping, likely binds through a more flexible mechanism involving intermediate states. For general protein-ligand binding, strong and long-ranged interactions have been thought to favor binding dominated by induced fit, while short-ranged interactions would favor conformational selection. [84] Our results support and extend this idea to CaM by showing that weak hydrophobic interactions dominate deep binding modes which indeed tend to be dominated by conformational selection in C-CaM. In the flexible binding of N-CaM, on the other hand, hydrophobic interactions occur less frequently. We hypothesize that the long-ranged electrostatic interactions of the N-CaM charged residues may initiate fast binding while the hydrophobic pocket in C-CaM may account for selectivity. Previously published studies using NMR [23] proposed that C-CaM is selective while N-CaM binds afterwards through induced fit mechanism, or a coupled conformational selection mechanism initiated by C-CaM [60]. This is in line with the argument given here on C-CaM selectivity and N-CaM flexibility. This study opens the way towards understanding the process and mechanisms behind calmodulin Ca2+-sensing. In a next step, the second aspect of binding, induced fit, may be studied. The full mechanism of specific target peptide binding could also be simulated starting from the conformational selection state of CaM identified here.
Calmodulin is a protein involved in the regulation of a variety of cell signaling pathways. It acts by making usually calcium-insensitive proteins sensitive to changes in the calcium concentration inside the cell. Its two lobes bind calcium and allow the energetically unfavorable exposure of hydrophobic residues to the aqueous environment which can then bind target proteins. The mechanisms behind the simultaneous specificity and variation of target protein binding is yet unknown but will aid understanding of the calcium-signaling and regulation that occur in many of our cellular processes. Here, we used molecular dynamics simulations and data analysis techniques to investigate what effect calcium has on the binding modes of calmodulin. The simulations and analyses allow us to observe and differentiate specific states. One domain of calmodulin is shown to be selective with binding involving short-distance interactions between hydrophobic residues, while the other binds target proteins through a more flexible mechanism involving long-distance electrostatic interactions.
chemical characterization
medicine and health sciences
molecular dynamics
protein interactions
electrophysiology
neuroscience
ion channels
molecular motors
actin motors
protein structure
thermodynamics
motor proteins
research and analysis methods
contractile proteins
proteins
chemistry
binding analysis
biophysics
molecular biology
free energy
physics
biochemistry
biochemical simulations
cytoskeletal proteins
cell biology
neurophysiology
physiology
myosins
biology and life sciences
physical sciences
computational chemistry
computational biology
macromolecular structure analysis
Effect of Ca2+ on the promiscuous target-protein binding of calmodulin