A Coupled Equilibrium Shift Mechanism in Calmodulin-Mediated Signal Transduction

Summary We used nuclear magnetic resonance data to determine ensembles of conformations representing the structure and dynamics of calmodulin (CaM) in the calcium-bound state (Ca2+-CaM) and in the state bound to myosin light chain kinase (CaM-MLCK). These ensembles reveal that the Ca2+-CaM state includes a range of structures similar to those present when CaM is bound to MLCK. Detailed analysis of the ensembles demonstrates that correlated motions within the Ca2+-CaM state direct the structural fluctuations toward complex-like substates. This phenomenon enables initial ligation of MLCK at the C-terminal domain of CaM and induces a population shift among the substates accessible to the N-terminal domain, thus giving rise to the cooperativity associated with binding. Based on these results and the combination of modern free energy landscape theory with classical allostery models, we suggest that a coupled equilibrium shift mechanism controls the efficient binding of CaM to a wide range of ligands.

(see Supplemental Experimental Procedures), calculated between all structures of the Ca 2+ -CaM and CaM-MLCK ensembles indicate that the structural differences between the Ca 2+ -CaM and CaM-MLCK states are the lowest for residues 29-54 and 101-130, respectively (Table S2). In addition, the distribution of the RMSDDs shown in Figure S4 clearly shows that the structural differences between Ca 2+ -CaM and CaM-MLCK are much more pronounced in the other parts of the structure of CaM. Therefore, structures were first overlaid on residues 29-54 and 101-130, respectively, and then the properties of the non-aligned regions analysed.

Supplemental Experimental Procedures Detailed Description of the Restrained Ensemble-Averaged Molecular Dynamics Protocol
Restrained molecular dynamics techniques were used in combination with ensemble simulations to obtain structural ensembles that satisfy NOE and S 2 restraints. The use of the S 2 order parameters as restraints enables the structural heterogeneity of the native state ensemble to be sampled (Best and Vendruscolo). We have carried out repeated cycles of simulated annealing (Lindorff-Larsen et al.) (see below) in order to increase the efficiency of the sampling. In addition, while NOE distances and S 2 order parameters were initially enforced across sixteen replicas to get structural ensembles (Lindorff-Larsen et al.), a recently established protocol (MUMO) where S2 order parameters are enforced across sixteen replicas but NOE distances, on the other hand, averaged across only two replicas improved the quality of the calculated ensembles (for more details see (Richter et al.)). We used this modified protocol (MUMO) to determine the ensembles of structures representing the Ca 2+ -CaM and CaM-MLCK states, respectively.
The restraint energy was implemented as (Best and Vendruscolo) where X corresponds either to NOE, or S 2 , α X is the force constant associated with each type of restraint, and ρ 0 (t) is defined as (Paci et al.) ρ 0 (t) = min 0≤τ ≤t ρ(τ )

[2]
In case of the NOE restraints where the sum is taken over the number N NOE of experimental NOE distances, and where N rep is the number of molecules (replicas) used in the ensemble-averaged simulations. The where the sum is taken over all equivalent atoms.
where the sum is taken over the number N S 2 of S 2 restraints and S 2,ens values were calculated as: where r i,k is the i To enable large inter-domain motions, the previously reported method (Best and Vendruscolo) (Lindorff-Larsen et al.) of restraining S 2 order parameters was significantly modified. In the procedure described here, molecular motions on the picosecond to nanosecond timescale are restrained within the molecular frame of the individual domains of CaM without altering the overall motion of the domains. At each integration step, a mass-weighted rigid body "ghost"alignment was performed for each domain (backbone atoms), and the non-aligned coordinates of the entire protein were stored. The restraining forces were then calculated individually for each aligned domain. Finally, the forces were rotated by the transposed alignment matrices and applied to the saved non-aligned coordinates. Hence, restraining forces were only applied within the molecular frames of non-aligned individual domains ,i.e., their overall rotational and translational motion was not affected by the restraint. When the individual domains of the saved snapshots were overlaid at the end of the structure calculations as in the "ghost"-alignment procedure, the S 2 order parameters restraints were satisfied.
The structure calculations were initiated with a preparation stage followed by 10 cycles of simulated annealing. In the preparation stage, the structures were first heated and equilibrated at 300 K (α S 2 =1000, α NOE =1000, t=100ps, units for α S 2 and α NOE are kcal/mol and kcal/(molÅ 4 ), respectively). Over the following 640 ps, the force constant α S 2 was progressively increased to 2.1*10 7 and 1.0*10 8 for Ca 2+ -CaM and CaM-MLCK, respectively, and α NOE to 1.7*10 8 and 1.0*10 9 . Subsequently, cycles of simulated annealing were carried out in order to sample conformational space efficiently. During each annealing cycle (lasting 170 ps), the molecules were heated to 500K and then cooled to 300K. After each cycle, 15 structures were extracted at 30 ps intervals, resulting in total of 2400 conformations for further analysis.

Control Simulations
For the control simulation using NOE-derived distances only as restraints, the structure of Ca 2+ -CaM was first heated and equilibrated at 300 K. Over the following 640 ps, the force constant α NOE was progressively increased to 1.7*10 8 . Then simulated annealing cycles were carried out in order to sample conformational space efficiently. In contrast to the standard protocol, a low α NOE (1.0*10 4 ) was used during the annealing cycles. In addition, only one replica could be simulated because no S 2 order parameters were used. After each of the 20 annealing cycles one structure was extracted.
For the normal mode ensemble, a normal mode calculation was carried out on an energy minimized structure of Ca 2+ -CaM. Minimization was done by using cycles of steepest descent and adopted basis Newton-Raphson methods. Harmonic restraints were applied on the backbone atoms and were progressively reduced at each cycle. Finally the restraints were removed and the structure minimized until the energy gradient was smaller than 10 -10 kcal/Å. Electrostatic and Lennard-Jones interactions were force switched. The CHARMM program and the CHARMM22 (Brooks et al.) parameter set were used for the energy minimization and the normal mode analysis.
The structural ensemble was generated by superposition of the first 100 modes. The ensemble properties (S 2 order parameters, RDC Q-factors) were similar when the normal mode analysis was carried out on Ca 2+ -CaM or the NTD and CTD individually.
The setup used for the restrained simulations, i.e. , Ca 2+ -CaM solvated a 8 Å shell of TIP3 water, was taken for the first unrestrained control simulation. A soft boundary potential was used to prevent water molecules from escaping (Beglov and Roux). The Nose-Hoover temperature coupling scheme was used to keep the T= 300K. In the second unrestrained control simulation, Ca 2+ -CaM was solvated in a 82x65x49 Å 3 orthorhombic box and periodic poundary conditions were applied. This unrestrained simulation was carried out at constant N, T, and P, using the Hoover temperature-pressure coupling scheme. All calculations used an atom-based truncation scheme with a list cutoff of 14 Å, a non-bond cutoff of 12 Å, and the Lennard-Jones smoothing function initiated at 10 Å. Electrostatic and Lennard-Jones interactions were force switched.
Molecular dynamics simulations used a 2 fs integration time step and SHAKE for covalent bonds involving hydrogen atoms (Ryckaert et al.). Simulations were carried out for 2 ns.

Helical Axes Definition
To calculate interhelical angles and the spatial disposition of individual helices, their axes were

Vector Geometry Mapping (VGM)
VGM is a method developed specifically to determine the conformational changes in EF-hand motifs (Yap et al.,1999). It compares the position of the exiting helix of an EF-hand with respect to its entering helix by orienting the two helices in a reference coordinate system. A reference EFhand motif is used to define first the coordinate system. We used the EF-hand motifs of the RDCrefined solution structures of Ca 2+ -CaM as reference conformation. Changes in the relative orientations of helices in the query EF-hand are then reported as angles θ, φ, and ω within the coordinate system of the reference. θ is the angle between the exiting helix of the query EF-hand and the Z-axis, φ is the angle between this helix and the x-axis, and ω is the counterclockwise rotation of the helix about its long axis with respect to the exiting helix of the reference.

Distance Difference (DD) Analysis
The DD (Nelson and Chazin) is defined as the difference between the distance of a given atom pair in one structure (Ca 2+ -CaM) and the distance of the same atom pair in another structure (CaM-MLCK): where d ij is the distance between atom i and j. For this analysis only Cα-atoms were used. To get a single value for a whole structure or segments of it, we defined RMSDD ij as: where N is the number of possible atom pairs in the entire protein or the selected segment.

Analysis of Interdomain Motions
To calculate the S 2 global, the structures of the Ca 2+ -CaM ensemble were first aligned on the entire CTD, and then the normalized vector connecting the residue at the centre of the flexible linker (residue 79) with the centre of mass of the NTD was used in Eq. 7.
We analysed the global orientations of the two CaM domains in the Ca 2+ -CaM state in several ways. First, we monitored the distribution of the position of the axis of helix IV whilst keeping helix V fixed ( Figures 6A and 6B). This analysis was repeated by aligning the entire CTD of the different structures instead of just the helix V regions ( Figure S8). This type of analysis allows monitoring very accurately the getting together of the two domains and was subsequently used to determine how compaction of CaM affects internal motion ( Figure 6E-6H). However, to fully assess the motions of the two domains with respect to each other, Euler angles defining the orientations of the NTD with respect to the CTD were calculated. To do so, an orthogonal axis system was first defined. The axis of helix 4 was defined as γ-axis, the β-axis is perpendicular to this axis in the plane defined by the axes of helix 4 and helix 3, and the α-axis is perpendicular to  Average backbone S 2 order parameters in the two domains of Ca  Correlations between experimental and calculated S 2 order parameters are given in parenthesis.

MLCK Ensemble
Segment (      and CaM-MLCK (blue), respectively. The first and second EF hand in the NTD (helices I-II and III-IV, respectively) are shown in A and B, respectively. The first and second EF hand in the CTD (helices V-VI and VII-VIII, respectively) are shown in C and D, respectively. The θ, φ, and ω angles observed in the RDC-refined solution structure of Ca 2+ -CaM (red) and the X-ray structure CaM-MLCK (blue) are shown in the inlay. Overall we find a significant overlap between the Ca 2+ -CaM and CaM-MLCK states of CaM with respect to the tertiary structure of their EF-hands; this overlap is more pronounced in the CTD.