⇦Chinese Physics Letters, 2017, Vol. 34, No. 6, Article code 068701⇨ Allosteric Mechanism of Calmodulin Revealed by Targeted Molecular Dynamics Simulation * Qian-Yun Liang(梁倩云)†, Chun-Li Pang(庞春丽)†, Jun-Wei Li(李军委), Su-Hua Zhang(张素花), Hui Liu(柳辉), Yong Zhan(展永)**, Hai-Long An(安海龙)** Affiliations Key Laboratory of Molecular Biophysics of Hebei Province, Institute of Biophysics, School of Sciences, Hebei University of Technology, Tianjin 300401 Received 22 November 2016 ^*Supported by the Natural Science Fund for Distinguished Young Scholars of Hebei Province under Grant Nos C2015202340 and C2013202244, the Fund for Outstanding Talents of Hebei Province under Grant No C201400305, the National Natural Science Fund of China under Grant Nos 11247010, 11175055, 11475053, 11347017, 31600594, 31400711 and 11647121, the Fund for the Science and Technology Program of Higher Education Institutions of Hebei Province under Grant No QN2016113,　the Scientific Innovation Grant for Excellent Young Scientists of Hebei University of Technology under Grant No 2015010, and the Natural Science Foundation of Hebei Province under Grant No C2017202208.
^†These authors contributed to this work equally.
^**Corresponding author. Email: hailong_an@hebut.edu.cn; zhany@hebut.edu.cn Citation Text: Liang Q Y, Pang C L, Li J W, Zhang S H and Liu H et al 2017 Chin. Phys. Lett. 34 068701 Abstract Calmodulin (CaM) is involved in the regulation of a variety of cellular signaling pathways. To accomplish its physiological functions, CaM binds with Ca$^{2+}$ at its EF-hand Ca$^{2+}$ binding sites which induce the conformational switching of CaM. However, the molecular mechanism by which Ca$^{2+}$ binds with CaM and induces conformational switching is still obscure. Here we combine molecular dynamics with targeted molecular dynamics simulation and achieve the state-transition pathway of CaM. Our data show that Ca$^{2+}$ binding speeds up the conformational transition of CaM by weakening the interactions which stabilize the closed state. It spends about 6.5 ns and 5.25 ns for transition from closed state to open state for apo and holo CaM, respectively. Regarding the contribution of two EF-hands, our data indicate that the first EF-hand triggers the conformational transition and is followed by the second one. We determine that there are two interaction networks which contribute to stabilize the closed and open states, respectively. DOI:10.1088/0256-307X/34/6/068701 PACS:87.14.E-, 87.15.ap, 87.15.kp, 87.15.hp © 2017 Chinese Physics Society Article Text Ca$^{2+}$ is the most common second messenger and plays vital physiological roles in variant physiological processes, such as muscle contraction, blood vessel contraction and expansion, secretion of hormones and modulation of neuron activity,[1-4] in both prokaryotes and eukaryotes. Ca$^{2+}$ accomplishes its physiological roles by binding with Ca$^{2+}$ binding proteins (CaBPs) and triggering conformational switching of CaBPs.[5-7] There are dozens of CaBPs involved in numerous-physiological processes.[8-10] Generally speaking, there is more than one Ca$^{2+}$ binding site in CaBPs, which are composed of 6–8 oxygen atoms derived from 5–7 residues. Calmodulin (CaM) which has been discovered for more than four decades is a representative CaBPs. It is involved in the regulation of a variety of cellular signaling pathways.[11,12] Over the past decades, researchers have made great efforts to understand the structure–function relationship of CaM.[13-17] CaM is a single 148-amino-acid polypeptide chain and is highly conserved in the process of evolution. It has been identified that there are four helix-loop-helix EF-hand motifs which form four distinct Ca$^{2+}$ binding sites. A pair of EF-hand motifs forms a globular domain. Two structurally independent globular domains (lobes) connected by a flexible helical linker form the dumbbell-shaped CaM. It is critical for CaM to transit from closed state to opened state induced by Ca$^{2+}$ binding. Although the x-ray crystal structures reveal the functional states of many CaBPs, they neither show the detailed allosteric mechanism of them nor determine the distinct Ca$^{2+}$ binding sites contributing to the state transition. Pang et al. have studied the interactions between Ca$^{2+}$ and Ca$^{2+}$-binding sites in many CaBPs by using molecular dynamics (MD) simulation.[17,18] Li et al. have reported that the targeted molecular dynamics (TMD) simulation is conductive to investigate the conformational change of macromolecular as an ion channel.[19] Here the CaM was used as a model of CaBPs to investigate the Ca$^{2+}$ inducing conformational change and the multisite collaboration by the simulation method. The N-terminal lobe of CaM was selected to be the simulation fragment in this work. A series of MD and TMD were performed on the N-terminal of CaM which includes two distinct EF-hand Ca$^{2+}$ binding sites. Conformational transition from closed state to opened state was achieved in both apo and holo situation. Our data show that the Ca$^{2+}$ binding speeds up the state transition by weakening the interactions between residues in the EF-hand loops. The interactions in the first EF-hand are weaker than those in the second EF-hand. The different patterns of the interactions in the two Ca$^{2+}$ binding sites make unequal contribution to the transition from closed state to opened state of CaM. This work provides an even better understanding of the allosteric dynamics of CaM and a practical simulation method to probe into the molecular mechanism of state transition of CaBPs regulated by Ca$^{2+}$.[20] Four structures of the N-terminal lobe of CaM, closed apo and holo, opened apo and holo, were constructed based on the crystal structures 1CFD[21] (closed and without Ca$^{2+}$) and 1CLL[22] (opened and with four Ca$^{2+}$), respectively. Each structure was immigrated into the water box and NaCl was added for electrical neutralization. Meadionize, a plugin of VMD, was used to add the Ca$^{2+}$ and Cl$^{-}$ into the minima of the electrostatic potential map of the solvated system by replacing water molecule. The electrostatic potential was generated by the 'potential' utility of the macroscopic electrostatics with atomic detail (MEAD) program package.[23,24] The equilibrated simulation box dimensions are 73 Å $\times$ 73 Å $\times$ 63 Å. The total number of atoms in the simulation system is 30862, including the fragment of CaM containing 1172 atoms, 9893 water molecules, and ions. For structural optimization, the MD simulations are performed by using NAMD (version2.9 (http://www.ks.uiuc.edu/Research/namd/))[25] with force field CHARMM27[26] and are finished on a cluster of computers. The integration time step is 2 fs. The simulation time was set to 10 ns. The electrostatic interactions are computed with a cut-off distance 12 Å and the particle mesh Ewald (PME) method. Periodic boundary conditions are introduced to stabilize the simulation space. The temperature coupling method is used to keep the temperature at 310 K. At the beginning of simulations, the system is minimized for 10000 steps to remove bad contacts. The Ca atoms of 40 to 44 amino acids are fixed to prevent the protein from extending out of the water box.[25,27]

Fig. 1. MD simulations of the N-terminal of CaM. (a) N-terminal of CaM in its closed-apo-state. (b) N-terminal of CaM in its closed-holo-state. Ca$^{2+}$ are highlighted as green ball. The left and right sides of the figure are EF-hand 1 and EF-hand 2, respectively. The yellow and red helixes of each EF-hand are E helixes and F helixes, respectively. The blue loop represents the Ca$^{2+}$ binding loop. (c, d) RMSD variations for closed apo and holo structures, respectively.

To achieve conformational change, TMD simulations were carried out,[28] a subset of atoms in the simulation is guided towards a final 'target' structure by means of steering forces. The targeted process is used to reach the final structure by stretching the protein with a harmonic strength applied to selected atoms.[29] The equilibrium closed state and the experimentally solved opened state of the CaM represent the starting and ending coordinates, respectively, and all the Ca atoms of the N-terminal lobe have been chosen as targeted atoms. The number of atoms used to calculate the root mean square deviation (RMSD) from the target structure was set to be the same as the number of restrained atoms. The simulation time of TMD was set to 10 ns. The time evolution of RMSD with respect to the target structure was observed. The targeted RMSD value was decreased monotonically from the initial RMSD to the target structure until it reached near 0 Å at the end of the TMD simulation. Visualisation and data analysis were carried out by VMD. To study the interplay between Ca$^{2+}$ binding and allosteric motion of CaM, we integrated the computational models which are in closed-apo-state (Fig. 1(a)). The Ca$^{2+}$ is placed in the energy minimum point flowing by a 10 ns MD. As shown in Fig. 1(b), there are two ions of Ca$^{2+}$ binding in the two EF-hand like Ca$^{2+}$ binding sites of the N terminal lobe of CaM, respectively. For each site, there are two helices connected by a loop. The RMSD values based on all the Ca atoms were calculated in Figs. 1(c) and 1(d). The RMSD values of the two models are less than 2.5 Å, which indicates that the models achieve the equilibrium states.

Fig. 2. Comparison between the presence and absence of Ca$^{2+}$ in binding site during TMD and MD simulation. The RMSD of the TMD simulation of (a) closed-apo-state and (b) closed-holo-state. (c) RMSD of the N-terminal lobe in closed apo and holo state during MD simulation. The H-bond number of loop 1 and loop 2 after 10 ns simulation of (d) MD and (e) TMD, respectively.

It is known that binding with Ca$^{2+}$ will facilitate the state transition of CaM. To understand the molecular mechanism by which Ca$^{2+}$ induces the conformational transition, we set the structure of open states as target structure and performed TMD on closed-apo-state and closed-holo-state to reproduce the state transition with or without Ca$^{2+}$. We measured the RMSD of the two TMD processes, as shown in Figs. 2(a) and 2(b). Both RMSDs show a pattern with a linear rise stage and a plateau stage. The rise stage is corresponding to the global conformational transition. The plateau stage means the local variation of the conformation. The difference between the two processes is the time spent to achieve the open state. It takes about 6.5 ns and 5.25 ns for the apo- and holo-closed states to reach the opened states, respectively, which indicates that the Ca$^{2+}$ promotes the opening process of CaM. H bonds play critical roles in stabilizing the structure and conformation transitions of proteins. There are two types of H-bonds formed between main chains and side chains of residues, respectively. The H-bonds between the main chains are thought to stabilize the secondary structure of protein molecules. The H-bonds between the side chains contribute to stabilize the tertiary structure and are involved in conformational change. In this work, the H-bonds formed between side chains of residues are analyzed to evaluate the stereo conformational change. Comparing the H-bonds of the apo- and holo-states, we find that binding with Ca$^{2+}$ weakens the interactions in both loop 1 and loop 2 (loop 1: Asp20 to Glu31 in EF-hand 1 and loop 2: Asp 56 to Glu67 in EF-hand 2, Figs. 2(d) and 2(e)), which indicates that the loops become less constrained after binding with Ca$^{2+}$ (Fig. 2(c)). Before Ca$^{2+}$ binding, the CaM is in a stable equilibrium state and the loops in both EF-hands are relatively rigid. After Ca$^{2+}$ is added in the site, electrostatic attraction will occur between Ca$^{2+}$ and the binding residues in EF-hand loops, and the interactions in the loops itself will be weakened. The flexibility of loops increases to cater to the conformational change of CaM. The data indicate that Ca$^{2+}$ binds at its site and triggers the state transition of CaM by weakening the interactions and increasing the flexibility of the loops. Next, we evaluate the detailed contributions to allosteric switching of the two EF-hand Ca$^{2+}$ binding sites. As shown in Fig. 3, the TMD simulations lead to a large-scale conformational change and accomplish an 'open' conformation (see Fig. 3(a)), in which the hydrophobic binding pocket of each EF-hand loop is exposed as shown in Fig. 3(b). The significant difference between the closed and the opened states of CaM is the angles formed by the two helixes of each EF-hand, which are 50$^{\circ}$ and 100$^{\circ}$ for the closed and opened states, respectively. It takes about 5.25 ns and 6.5 ns for the first EF-hand and second EF-hand to achieve their opened states (Figs. 3(c) and 3(d)), which indicates that the first EF-hand triggers the conformational transition and the second EF-hand follows.

Fig. 3. Conformational switching of N-terminal of CaM achieved by TMD. (a) The initial (closed) and final (open) structure of the N-terminal of CaM achieved by TMD, in which the lower image is their overlap. (b) The initial (closed) and final (open) structure of the two loops of CaM achieved by TMD, as well as their overlap, respectively. (c, d) The evolutions of angles formed by E helix and F helix in EF-hand 1 and EF-hand 2, respectively. The black and red lines are corresponding to the MD and TMD simulations. (e) The average value of the transition time.

Fig. 4. The RMSD of the relay-like TMD. (a) The time course of the RMSD. (b) The crystal structure of the CaM. Backbones of EF1 and EF2 are highlighted in blue and red, respectively.

As shown in Fig. 3, the conformational switching of the second EF-hand comes later than that of the first one, which indicates that the two EF hands may contribute unequally to the state transition of CaM. For validation, relay-like TMD is performed in the N-terminal lobe of CaM. Firstly, a 10 ns TMD simulation is carried on by setting the Ca atoms in EF-hand 1 (Fig. 4(b)) as target atoms. Then, another 10 ns TMD follows the last state of the previous TMD by setting the Ca atoms in EF-hand 2 (Fig. 4(b)) as target atoms. Likewise, the second 20 ns TMD is performed in the reverse order of the two EF hands. The time courses of RMSD of the two TMDs show that the CaM achieves the targeted conformation more easily by the first TMD than the second one as shown in Fig. 4(a), which indicates that the state transition of CaM should start at the first EF-hand, and the two EF-hands contribute unequally to the conformational switching. Conversely, since the first EF-hand is more easily opened, the process of its folding is slower. Also we think that the folding of the first EF-hand is the rate-limiting step. To understand the molecular mechanism which determines the contribution of the two EF-hands on state transition, we measure the H-bonds in different cases of the EF-hands in both MD and TMD simulations. As shown in Fig. 5 and Table 1, the patterns of H-bond in different cases are significantly different. Figure 5(a) shows that there are very weak interactions between E and F helixes of the first EF-hand. However, the interactions between E and F helixes of the second EF-hand are much stronger than that of the first one (Fig. 5(b)). The linear fitting data show that the interactions decrease gradually during the state transition, which means that these interactions may stabilize the closed state. Weaker interactions between E and F helixes of the first EF-hand means that the protein more easily achieves conformational switching. The interactions in the loops increase during the state transition (Figs. 5(c) and 5(d)), which means that the interactions in the loops may serve as stabilizing the open conformation.

Fig. 5. The evolution of H-bond number in different cases in TMD simulations. (a, b) The evolutions of numbers of H-bonds formed between residues in E and F helixes of each EF-hand, respectively. (c, d) The evolutions of numbers of H-bonds formed between residues in loops of each EF-hand, respectively. The red lines are the linear fitting data.

**Table 1.** The mean value of H-bond Number in different cases of CaM.
	Protein	EF 1	EF 2	Loop 1	Loop 2
ApoMD	66.80	0.24	4.44	5.80	5.11
ApoTMD	59.14	1.10	1.30	6.80	6.92
HoloMD	60.79	0.20	3.68	4.67	4.45
HoloTMD	59.92	0.32	2.76	6.26	5.20

The variation of H-bond represents that of conformation. Since the conformations of apo-closed state and holo-closed state are similar after the MD simulations, the number of H-bonds changes just slightly from apo-closed state to holo-closed state. However, comparing the apo-closed state with the holo-closed state in both MD and TMD, the same variation tendencies are found in Table 1. Combining MD with TMD, we achieve the dynamical state transition and find the conformational pathway of CaM. To understand the Ca$^{2+}$-dependent state transition of CaM, researchers have applied static structures in diverse state or combined distinct receptors and CaM, which accumulate numerous experimental and computational data,[30,31] then identified four EF-hands Ca$^{2+}$ binding sites in the CaM. To understand the dynamical process of the CaM state transition, we perform the TMD on the closed and opened states of CaM induced by Ca$^{2+}$ binding. It was reported by Wang et al. that Ca$^{2+}$ can bind to the two EF-hands sites differently. One is conformational selected binding and the other is Ca$^{2+}$ inducing conformational change binding,[32] which indicate that the two EF-hands may play different roles in the state transition of CaM. However, they did not discuss the structural basis and the mechanism of the two EF-hands contributing unequally to state transition of CaM. Here by performing the TMD on closed CaM, we identify the allosteric pathway and reveal that Ca$^{2+}$ facilitates the state transition of CaM by weakening the interactions and increasing the flexibility of CaM. The difference of interactions between the E and F helixes of the two EF hands is the basic reason for unequal contribution of the EF-hands. Moreover, we find that the first EF-hand triggers the conformational switching followed by the second one. There is cooperation between the two EF-hand Ca$^{2+}$ binding sites of CaM. Generally, CaBPs have multiple Ca$^{2+}$ binding sites. Four EF-hand Ca$^{2+}$ binding sites are recognized in CaM and other EF-hand proteins. At least two Ca$^{2+}$ binding sites are identified in large-conductance Ca$^{2+}$-activated K$^{+}$ channel[33] and Ca$^{2+}$-activated Cl$^{-}$ channel.[34] To determine whether the distinct Ca$^{2+}$ binding sties integrate to promote the conformational switching of CaBPs, or independently, we perform relay-like TMDs and find that there is a time sequence between the two EF-hands. Our data show that transition is more easily accomplished when the conformational switching is provoked by the first EF-hand than by the second one. It is indicated that the first EF-hand triggers the conformational switching and induces the second one to transit from closed state to opened state. Our work could shed light on understanding of the molecular mechanism of conformational switching of CaBPs modulated by multiple Ca$^{2+}$ binding sites. In conclusion, we have performed a series of MD and TMD simulations on a model protein of CaM to address the Ca$^{2+}$ activating conformational change of CaBPs. The results indicate that the Ca$^{2+}$ binding speeds up the state transition process and the different binding sites contribute to the allosteric change unequally. This work facilitates the understanding of the allosteric dynamics of CaM and the molecular mechanism of other allosteric proteins regulated by Ca$^{2+}$. References Glial calcium and diseases of the nervous systemThe central role of calcium in the effects of cytokines on beta-cell function: Implications for type 1 and type 2 diabetesCalcium - a life and death signalThe versatility and universality of calcium signallingMolecular simulation assisted identification of Ca2+ binding residues in TMEM16ACa2+ Regulation in the Near-Membrane Microenvironment in Smooth Muscle CellsCa2+-Dependent Interaction between FKBP12 and Calcineurin Regulates Activity of the Ca2+ Release Channel in Skeletal MuscleTMEM16A/B Associated CaCC: Structural and Functional InsightsCalmodulin: a prototypical calcium sensorStructure of a Trapped Intermediate of Calmodulin: Calcium Regulation of EF-hand Proteins from a New PerspectiveCalmodulin in ActionGenetic polymorphism and protein conformational plasticity in the calmodulin superfamily: Two ways to promote multifunctionalityCombining fragment homology modeling with molecular dynamics aims at prediction of Ca2+ binding sites in CaBPsMolecular dynamics simulation studies on Ca2+-induced conformational changes of annexin IIdentification of the Conformational transition pathway in PIP2 Opening Kir ChannelsThe 1.0 Å crystal structure of Ca2+-bound calmodulin: an analysis of disorder and implications for functionally relevant plasticitySolution structure of calcium-free calmodulinCalmodulin structure refined at 1.7 Å resolutionElectrostatic calculations of the pKa values of ionizable groups in bacteriorhodopsinScalable molecular dynamics with NAMDAll-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins ^†NAMD2: Greater Scalability for Parallel Molecular DynamicsTargeted molecular dynamics: A new approach for searching pathways of conformational transitionsIdentification of Three Interactions to Determine the Conformation Change and to Maintain the Function of Kir2.1 Channel ProteinCalmodulin as a direct detector of Ca2+ signalsEF-hand calcium-binding proteinsEnergy landscape views for interplays among folding, binding, and allostery of calmodulin domainsStructure of the gating domain of a Ca2|[plus]|-activated K|[plus]| channel complexed with Ca2|[plus]|/calmodulinIncreased Expression of the Human Ca ^{2 +} -activated Cl ⁻ Channel 1 (CaCC1) Gene in the Asthmatic Airway

[1]	Batistič O and Kudla J 2012 BBA-BioMemb. 1820 1283
[2]	Jiao H and Liu S 2011 J. Cell Sci. 93 2029
[3]	Nedergaard M, Rodríguez J J and Verkhratsky A 2010 Cell Calcium 47 140
[4]	Ramadan J W, Steiner S R, O'neill C M and Nunemaker C S 2011 Cell Calcium 50 481
[5]	Berridge M J, Bootman M D and Lipp P 1998 Nature 395 645
[6]	Berridge M J, Lipp P and Bootman M D 2000 Nat. Rev. Mol. Cell. Bio. 1 11
[7]	Pang C L, Yuan H B, Cao T G, Su J G, Chen Y F, Liu H, Yu H, Zhang H L, Zhan Y, An H L and Han Y B 2015 J. Comput. Aid. Mol. Des. 29 1035
[8]	Bazzazi H, Kargacin M E and Kargacin G J 2003 Biophys. J. 85 1754
[9]	Shin D W, Pan Z, Bandyopadhyay A, Bhat M B, Kim D H and Ma J J 2002 Biophys. J. 83 2539
[10]	Pang C L, Yuan H B, Ren S X, Chen Y F, An H L and Zhan Y 2014 Protein Peptide Lett. 21 94
[11]	Chin D and Means A R 2000 Trends Cell Biol. 10 322
[12]	Grabarek Z 2005 J. Mol. Biol. 346 1351
[13]	Carafoli E and Klee C B 1999 Calcium as a Cellular Regulator (New York: Oxford University Press)
[14]	Hoeflich K P and Mitsuhiko I 2002 Cell 108 739
[15]	Mitsuhiko I and Ames J B 2006 Proc. Natl. Acad. Sci. USA 103 1159
[16]	Linda J and Eldik V 1998 Calmodulin and Signal Transduction (New York: Academic Press)
[17]	Pang C L, Cao T G, Li J W, Jia M W, Zhang S H, Ren S X, An H L and Zhan Y 2013 J. Comput. Aid. Mol. Des. 27 697
[18]	Shesham R D, Bartolotti L J and Li Y 2008 Protein Eng. Ees. Sel. 21 115
[19]	Li J W, Lu S Q, Liu Y Z, Pang C L, Chen Y F, Zhang S H, Yu H, Long M, Zhang H L, Logothetis D E, Zhan Y and An H L 2015 Sci. Rep. 5 11289
[20]	Wilson M A and Brunger A T 2000 J. Mol. Biol. 301 1237
[21]	Kuboniwa H, Tjandra N, Grzesiek S, Ren H, Klee C B and Bax A 1995 Nat. Struct. Biol. 2 768
[22]	Chattopadhyaya R, Meador, W E, Means A R and Quiocho F A 1992 J. Mol. Biol. 228 1177
[23]	Bashford D 2006 Lect. Notes. Comput. Sci. 1343 233
[24]	Bashford D and Gerwert K 1992 J. Mol. Biol. 224 473
[25]	Phillips J C, Braun R, Wang W, Gumbart J, Tajkhorshid E, Villa E, Chipot C, Skeel R D, Kalé L and Schulten K 2005 J. Comput. Chem. 26 1781
[26]	Mackerell A D, Bashford D, Bellott M, Dunbrack R L, Evanseck J D, Field M J, Fischer S, Gao J, Guo H, Ha S et al 1998 J. Phys. Chem. 102 3586
[27]	Kalé L, Skeel R, Bhandarkar M, Brunner R, Gursoy A, Krawetz N, Phillips J, Shinozaki A, Varadarajan K and Schulten K 1999 J. Comput. Phys. 151 283
[28]	Schlitter J, Engels M and Krüger P 1994 J. Mol. Grap. 12 84
[29]	Li J W, Xiao S Y, Xie X X, Yu H, Zhang H L, Zhan Y and An H L 2015 Chin. Phys. Lett. 32 028702
[30]	Faas G C, Sridhar R, Lisman J E and Istvan M 2011 Nat. Neurosci. 14 301
[31]	Lewitbentley A and Rety S 2000 Curr. Opin. Struct. Biol. 10 637
[32]	Li W, Wang W and Takada S 2014 Proc. Natl. Acad. Sci. USA 111 10550
[33]	Schumacher M A, Rivard A F, Bächinger H P and Adelman J P 2001 Nature 410 1120
[34]	Hoshino M, Morita S, Iwashita H, Sagiya Y, Nagi T, Nakanishi A, Ashida Y, Nishimura O, Fujisawa Y and Fujino M 2002 Am. J. Resp. Crit. Care 165 1132