Role of water-bridged interactions in metal ion coupled protein allostery

Allosteric communication between distant parts of proteins controls many cellular functions, in which metal ions are widely utilized as effectors to trigger the allosteric cascade. Due to the involvement of strong coordination interactions, the energy landscape dictating the metal ion binding is intrinsically rugged. How metal ions achieve fast binding by overcoming the landscape ruggedness and thereby efficiently mediate protein allostery is elusive. By performing molecular dynamics simulations for the Ca2+ binding mediated allostery of the calmodulin (CaM) domains, each containing two Ca2+ binding helix-loop-helix motifs (EF-hands), we revealed the key role of water-bridged interactions in Ca2+ binding and protein allostery. The bridging water molecules between Ca2+ and binding residue reduces the ruggedness of ligand exchange landscape by acting as a lubricant, facilitating the Ca2+ coupled protein allostery. Calcium-induced rotation of the helices in the EF-hands, with the hydrophobic core serving as the pivot, leads to exposure of hydrophobic sites for target binding. Intriguingly, despite being structurally similar, the response of the two symmetrically arranged EF-hands upon Ca2+ binding is asymmetric. Breakage of symmetry is needed for efficient allosteric communication between the EF-hands. The key roles that water molecules play in driving allosteric transitions are likely to be general in other metal ion mediated protein allostery.


Well-tempered Metadynamics algorithm
The metadynamics algorithm was used to enhance conformational sampling in simulations so that Free Energy Surfaces (FES) could be reliably computed. In complex systems, it is necessary to identify physically reasonable collective variables (CVs) to describe large scale conformational transitions. In the canonical metadynamics algorithm, the system is biased through the addition of a sum of repulsive Gaussian terms to the system potential energy function. The biasing potentials are usually applied to the CVs. The total biasing potential applied on a set of CV, , at time t during the simulation can be expressed as [1]: where h is the height of each Gaussian potential, δ y i is the Gaussian width associated with the CV Y i (x), δt determines the time step at which the biasing potentials are added. In order to accelerate the convergence of metadynamics, the well-tempered algorithm was designed [2], in which the height of the Gaussian potentials is systematically attenuated during the course of the simulations using e −Vmeta(Y (x),t)/∆T , where ∆T represents a characteristic energy [2]. The form of the biasing potential used in our simulations is where h 0 is the initial Gaussian height. In practice, the rate of decrease is controlled by the bias factor η = (T + ∆T )/T [3], where T is the system temperature. Combining the replica exchange methodology and metadynamics, the bias-exchange metadynamics was designed to enhance the sampling efficiency [4]. As an inherent feature of the algorithm, the converged metadynamics runs enable the calculation of the canonical probability distribution of the CVs directly. However, the statistics of other degrees of freedom is affected by the bias, and require "reweighting" techniques to reconstruct the unbiased canonical probability distribution [5]. In particular, for well-tempered metadynamics, the reconstruction algorithm is described in Ref. [6], and can be found in the PLUMED package [7], which is used in our work.

Technical details of the simulations
During the metadynamics simulations, we applied restraining potentials to the values of the MSD (mean square deviation) of a conformation with respect to the folded states to limit the accessible region of phase space sampled by the protein. The use of such restraining potentials physically prevents unreasonable unfolding of nCaM. The functional form of the restraining potential is [3]: where the energy constant κ = 0.8kJ/Å 4 , d(X α i , X(t)) is the MSD from the closed (i = 1) or open (i = 2) structure of the α th (α = 1, 2) EF-hand for a given conformation X(t), the limit d α 0 (α = 1, 2) is chosen as d 1 0 = 41Å 2 for EF 1 and d 2 0 = 36Å 2 for EF 2 . The terminal residues with high flexibility were not included in the calculations of the MSD values. For better convergence, we prevent Ca 2+ ions from leaving too far from the binding sites by applying a lower limit to the coordination number, as one of the reaction coordinates.

Convergence of metadynamics simulations
We sampled the solution conformations of nCaM and the associated Ca 2+ binding using welltempered bias-exchange metadynamics with 5 replicas (in four of which biasing potentials were applied). The simulations were initiated from random conformations, and were continued for 800ns. We assessed the convergence of the metadynamics simulations by ensuring that the freeenergy differences between all the major populated regions for each CV reach a steady state.
For the four biased replicas, this condition was satisfied after about 400ns of sampling (see Fig J in S1 Text). In addition, as an embedded property of well-tempered metadynamics, the automatically rescaled Gaussian height guarantees convergence, and avoids fluctuations around the correct free energy value [2,8]. Therefore, we checked the height of the biasing Gaussian potentials during the simulations. As can be seen in Fig

Coordination probability
Binding pathways are quantified using the probability, P (i, N α ), that the i th residue is coordinated to Ca 2+ when the total ligand number of EF-hand α is N α (N α = 1, 2, . . . , 6). We define where N α and k α i are defined in equation (1) and equation (2) (in the main text), respectively, and θ Nα (N α ) is a step function: is a two dimensional probability distribution of the variables k α i and N α , which cannot be directly calculated as the average over time from the metadynamics. We use the reweighting techniques described in Ref. [6] to reconstruct the unbiased distribution. We calculated the two-dimensional probability distributions on a 60 × 60 grid. Then we constructed bins on the     N 2 (B). The relevant residues in the hydrophobic core are Phe16, Phe19, Ile27, Leu32 for EF 1 , and Ile52, Val55, Ile63, Phe68 for EF 2 . (C) The hydrophobic cluster in the nCaM shows that it is densely packed.  Potential of mean force for the binding/unbinding process of Glu31 to Ca 2+ as a function of the minimum distance between side-chain oxygen atoms of Glu31 and Ca 2+ using modified parameter set of LJ potential function described by Kahlen and coworkers [9]. The black line illustrates the results with the water molecules being realistically treated. The red line shows the results with electrostatic interactions involving the bridging water molecules being turned off. The error bars represent the difference between the results from two groups of simulations, each group containing 10 independent metadynamics simulations. Compared to Fig 6A shown in the main text, the native coordination states corresponding to the free energy minimum around the distance of ∼0.24 nm are less stable. However, after switching off the electrostatic interactions between the bridging water and the Glu31 side-chain oxygen atoms in the control simulation, the free energy barrier for the ligand exchange becomes much higher, suggesting that the bridging water tends to reduce the ruggedness of the ligand exchange free energy landscape, which is consistent with the results shown in Fig 6A based on intact OPLS-AA force field.  Supplementary Equation (S2)) in the metadynamics as functions of simulation time. In our well-tempered metadynamics simulations, the initial value of the Gaussian height was set to 0.3kJ/mol, while the bias factor η = 200 was used to control the rate of potential decrease. After 200 ns, all the biasing Gaussians are smaller than 0.05kJ/mol, which further decrease to values less than 0.03kJ/mol after 300ns. These correspond to only about 1/10 of the initial values.