Structures and energies of the transition between two conformations of the alternate frame folding calbindin-D_9k protein: a theoretical study

Abstract

We carried out conventional molecular dynamics simulations, targeted molecular dynamics simulations and energy calculations for the two states (AFF-D_9k-N′ and AFF-D_9k-N) of the alternate frame folding (AFF) calbindin-D_9k protein and their conformational transition in the Ca²⁺-free form to address the dynamical transition mechanism from AFF-D_9k-N′ to AFF-D_9k-N states. We found that the structural characteristics of the two stable AFF-D_9k-N′ and AFF-D_9k-N states present the coupled conformations of mutually exclusive folding. The transition from AFF-D_9k-N′ to AFF-D_9k-N states may occur via two transition states and an intermediate with the first rate-controlling barrier of 4.7 kcal mol⁻¹ and the second barrier of 1.7 kcal mol⁻¹. These results showed that the conformational transition pathway from AFF-D_9k-N′ to AFF-D_9k-N is energetically feasible due to the low rate-controlling barrier. Moreover, the stabilities of the AFF-D_9k-N′ and AFF-D_9k-N conformations calculated by the MM_PBSA method are consistent with the available experimental result. The transition mechanism involves the relative movements of the two mutually exclusive folding regions, EF2 and EF2′, and their coupled folding-unfolding. The crucial mediation of the deformation of EF1 hand located at the middle of the AFF calbindin-D_9k protein plays a key role in this transition process. The correlation analysis revealed the allosteric communications between the EF2′ and EF2 regions via the mediation of the EF1 hand during the transition of AFF-D_9k-N′ to AFF-D_9k-N.

---

1. Introduction

It is well known that allostery of a protein is used to control protein function and plays a crucial role in understanding fundamental biological processes of cell signaling, the regulation of metabolism, etc.^1–4 Some proteins can achieve the biological switchable functions through self-structural transformations and the variation of the free energies.^5–7 Switchable proteins with new functionalities can be created through protein engineering after understanding their switching mechanisms.^8,9 As switchable proteins, the calmodulin family and its homologues have been extensively studied because of their ubiquity in signal transductions and regulations of many cellular processes.^10–15 A wild type (WT) calbindin-D_9k, being a small globular protein of the calmodulin superfamily and possessing a pair of helix-loop-helix structures called EF-hand (Lys1-Gln75), plays a key role in calcium binding and calcium buffer regulations as a typical functional unit.^16–19 The X-ray crystal structure of the WT calbindin-D_9k protein was reported by Szebenyi and co-workers in 1986.²⁰ Its high-resolution three-dimensional solution structure in the apo state has also been reported by Chazin and co-workers in 1995.²¹ Unlike other members of this superfamily, such as calmodulin and troponin C, perform their regulatory function through substantial conformational rearrangements, the WT calbindin-D_9k protein is unable to undergo a large conformational change upon to practise regulatory function.^22–25 To achieve switchable function for some proteins, Stewart N. Loh and co-workers developed an alternate frame folding (AFF) strategy from the wild-type fold (N) to a circularly-permuted fold (N′) to build the switchable protein based on conformational changes.^26–31 The protein conformational switches applied by AFF strategy are in response to a signaling event such as ligand binding etc.⁹ Moreover, the AFF approach is the mutually exclusive folding reaction, i.e., unfolding of one segment of the protein is coupled to folding of another segment, which is widely used in the allosteric processes of the switchable proteins.^26,27

The AFF calbindin-D_9k protein is the modified switchable protein from a wild type (WT) calbindin-D_9k protein (including the sequence of Lys1-Gln75), and was generated by duplicating the parent sequence (designated by prime superscripts of Ser44′-Gln75′) at N-terminal of WT calbindin-D_9k protein with a six-amino acid linker inserted between the residue 75′ and this N-terminal (see Fig. 1). The AFF calbindin-D_9k protein could adopt two distinctly different conformations in Ca²⁺-free form, i.e., N′ and N conformations, with similar stabilities in native condition.²⁷ This duplicate parent sequence equates to the original 44-75 sequence which needs not to correspond to a functional domain for a typical protein.²⁶ The folding component of N conformation of the AFF calbindin-D_9k protein is comprised of the WT amino acid sequence, and adopts the WT protein structure; while the folding component of N′ conformation of the AFF calbindin-D_9k protein with the circularly permuted structure consists of the original sequence of residues 1-43, the duplicate parent sequence and the linker to bridge the original N- and C-termini of the WT calbindin-D_9k protein.^26–28 The N conformation of the AFF calbindin-D_9k protein possesses two EF-hands (assigned as EF1 hand of Lys1-Met43 and EF2 hand of Ser44-Gln75), the structure of which is similar to that of the WT calbindin-D_9k protein, an unstructured EF region of Ser44′-Gln75′ (assigned as EF2′ region) and the linker (see Fig. 1). Its mutually exclusive folding N′ conformation possesses the circularly permuted folding structure of EF1 and EF2′ hands, and an unstructured EF2 region (see Fig. 1).^26–28 Especially, a surface loop of Phe36-Met43, connecting the EF1 and EF2 hands, can be as the only possible site for circular permutation in the AFF strategy.^27,32 Considerably, the AFF strategy for the calbindin-D_9k protein describes the coupled folding–unfolding equilibrium between two modified protein regions that fold in a mutually exclusive fashion, which is similar to the HT model developed by Hilser and Thompson.^33,34 Moreover, the barnase protein and ribose binding-protein have been successfully converted to an artificial zymogen and a fluorescent sensor for ribose, respectively, by applying the AFF modification.^9,30,31,35 As a guideline, the applications of the AFF modification to other biomedical and diagnostic proteins have also been attracting significant attention and interest. Because the AFF calbindin-D_9k protein is the modified switchable protein, the conformational transition of their N′ and N conformations plays a key role in practising its switch function.^26–28 It is difficult to shed light on the microscopic process of the conformational transition experimentally. The interconversion of N′ and N conformations of the AFF calbindin-D_9k protein through the kinetic experiment has been reported by Stewart N. Loh and co-workers.²⁷ And it was also suggested that this conversion process does not require a whole-molecule unfolding.²⁷ However, experimental studies on detail conformational variations and required energies for the transitions of N′ and N conformations of the AFF calbindin-D_9k protein at atomic level are very limited so far. Theoretically, the cooperative interactions between the two binding sites of calcium ions in a WT calbindin-D_9k protein have been investigated using molecular dynamics simulations based on detailed atomic models.^24,36 However, theoretical investigations on the stable structures and conformational transitions of the AFF calbindin-D_9k protein have not yet been detailed so far. Therefore, it remains unclear what the dynamical transition mechanism of the two conformations is and what the exact allosteric communication in this transition is.

Fig. 1 The amino acid sequences of AFF calbindin-D_9k protein in Ca²⁺-free. Amino acids are numbered according to the crystal calbindin-D_9k protein.²⁴ The subsequences of EF2′, EF1, EF2 and linker in AFF calbindin-D_9k protein are colored in magenta, cyan, orange and black, respectively.

To explore the conformational transition process in more detail at the atomic level, we carried out the conventional and targeted molecular dynamics simulations combined with the potential of mean force (PMF) simulations to investigate the structures and energies in the transition process of two stable N′ and N conformations of AFF calbindin-D_9k protein in Ca²⁺-free. The results obtained from the simulations provide valuable insights into understanding the detailed conformational transition process concerning how one state of the protein gradually transforms to another one and what the corresponding free energy requirements are.

2. Models and methods

2.1. Initial structures

Based on the previous experimental studies, the two N′ and N conformations of the alternate frame folding (AFF) calbindin-D_9k protein in Ca²⁺-free form (assigned as AFF-D_9k-N′ and AFF-D_9k-N models) were designed from the crystal structure of the wild type (WT) calbindin-D_9k protein (PDB entry 3ICB). The AFF-D_9k-N′ and AFF-D_9k-N models used in the MD simulations correspond to the mutually exclusive folding structures of folded N-terminal vs. unfolded C-terminal and unfolded N-terminal vs. folded C-terminal, respectively. The AFF-D_9k-N′ and AFF-D_9k-N models consist of a conserved EF1 hand (Lys1-Met43), two mutually exclusive folding regions, EF2 (Ser44-Gln75) and EF2′ (Ser44′-Gln75′), and a short linker (Gly76-Met81) between the EF2′ region and EF1 hand with the connections of Lys1 and Met81, Met43 and Ser44, Gln75′ and Gly76 (see Fig. 1). The sequences of the EF2′ and EF2 regions are nearly identical except for the Glu65′ mutated by Gln65′ in the EF2′ region. In order to choose initial conformations of two AFF-D_9k-N′ and AFF-D_9k-N models for MD simulations, the coordinates of sequences of Lys1-Gln75 for the AFF-D_9k-N model and Ser44′-Gln75′-Gly76-Met81-Lys1-Met43 (which is similar to the sequence of the WT calbindin-D_9k protein in Ca²⁺-free) for the AFF-D_9k-N′ model were taken from the crystal structure of the WT calbindin-D_9k protein (PDB entry 3ICB). The initial coordinates of other alternate remaining sequences for the two models were built by using the loop search method in the Swiss-PDBViewer (also known as DeepView, http://spdbv.vital-it.ch/). Total coordinates with the crystal structure of sequence of the WT calbindin-D_9k protein and the unstructured region of the additional alternate sequences form the initial structures of two AFF-D_9k-N′ and AFF-D_9k-N models as AFF calbindin-D_9k protein in Ca²⁺-free for the MD simulations. 12 Na⁺ ions were used to neutralize each of the AFF-D_9k-N′ and AFF-D_9k-N models, and an ionic strength of 50 mM was generated by adding 4 Na⁺ and 4 Cl⁻ ions for these two models.^37,38 Similar counterion processes were applied to the other models. Each of the systems was explicitly solvated by using the TIP3P water potential inside an orthorhombic box of water molecules with a minimum solute-wall distance of 8 Å.

2.2. Conventional molecular dynamics simulation

All conventional molecular dynamics (CMD) simulations for the two models were carried out using the AMBER 9 package³⁹ and ff03 all atom force field parameters.^40–42 The computational details of the CMD procedure are given in the ESI.†

2.3. Targeted molecular dynamics simulation

Generally, the structure transitions are beyond the reach of conventional simulations without biased potential constraint due to transition simulations on the scale of microseconds and longer. Targeted molecular dynamics (TMD) simulation is a method to observe large-scale conformational transition between two known end-point conformations of a molecule at ordinary temperature by applying a time-dependent with the help of a constraint. The transition is enforced independently of the height of energy barriers relative to biased potential constraint.⁴³ A restraint energy term was added to the energy function proportional to the square of the difference which may be characterized as the mass-weighted root-mean-square deviation (RMSD) of the current structure to the target structure in terms of atomic positions.^43,44 The functional form of the restraint energy can be written aswhere, k₁ is the harmonic force constant per atom, N is the number of the restrained atoms, RMSD(t) is the root-mean-square deviation of the simulated structure at time t relative to the target structure, and RMSD₀(t) is the prescribed target RMSD value at time t that decreases to zero linearly with time to drive the system from an initial structure to the target structure. To determine an appropriate value of the harmonic force constant k₁, we tested the use of various k₁ values for the TMD simulations (the data were depicted in Fig. S1 of the ESI†), and chose k₁ = 5.0 kcal (mol⁻¹ Ų) as the lowest harmonic force constant to apply onto all the backbone atoms of two models to bias the trajectories toward the target structure. The computational details are provided in the ESI.†

2.4. Potential of mean force (PMF) simulations

In order to explore the free energy profile for the structural transformation from AFF-D_9k-N′ to AFF-D_9k-N states, the potential of mean force (PMF) was calculated by using umbrella-sampling MD simulations and the weighed histogram analysis method (WHAM).^45–49 Based on the mutually exclusive folding reaction of AFF-D_9k-N′ to AFF-D_9k-N states, a convenient reaction coordinate was chosen as the mass center distance between the EF1 hand and EF2 region that has been tested to satisfy a simple linear variation feature. In umbrella-sampling calculations, the biasing harmonic potentials were introduced to constrain each conformation to a narrow range of reaction coordinates. The functional form of the restraining potential for reaction coordinate r for the ith window in the current umbrella sampling waswhere k₂ is a force constant and r_i0 is a center of the window. A total of 67 windows from the reaction coordinate of the mass center distances of 23.1 to 10.3 Å with the increments of 0.2 Å were generated for the structural transformation from AFF-D_9k-N′ to AFF-D_9k-N states. The analysis of the all trajectories and calculations of the geometric parameters of protein were performed with the PTRAJ module of the AMBER 9 program. After all the umbrella-sampling MD simulations were finished, the data collected from separate simulation windows were combined along the reaction coordinate. These data were then used to calculate the PMF for the whole structural transformation process with the weighed histogram analysis method (WHAM) using the code developed by Alan Grossfield (http://membrane.urmc.rochester.edu/Software/WHAM/WHAM.html).^47,50 The test of PMF calculation convergence based on the different simulation times was performed and is shown in Fig. S2 of the ESI,† indicating that 2 ns simulation for each window of the PMF calculations was sufficiently long for obtaining the converged PMF results. The absence of noteworthy differences between the curves generated from 2 ns and 3 ns runs suggests that a 2 ns simulation may be sufficient to obtain a reasonable estimate of the free energy along the ordering reaction coordinate.

2.5. Calculations of interhelical angle, B-factor and correlation of atomic motions

To analyze conformational changes in the relative orientations of any two helices, the program INTERHLX (written by Kyoko Yap, available at http://structbio.vanderbilt.edu/chazin/wisdom/interhel.html) was used to calculate the angles between two helices in the EF1 and EF2′ hands in the AFF-D_9k-N′ and AFF-D_9k-N models.^25,51 The B-factor (the temperature factor) of the Cα atom of each residue was estimated from the atomic fluctuations for the EF1 and EF2′/EF2 hands as follows:^52,53where B_i is a B-factor for a Cα atom and (RMSF_i)² is the mean square atomic fluctuation of the Cα atom from its equilibrium (average) position. The dynamic feature of a protein and the extent of correlation of motions in the different regions of a protein were assessed via the calculations of cross-correlation coefficients, C(i,j), given as follows:

C(i,j) = 〈Δr_i × Δr_j〉/(〈Δr_i²〉〈Δr_j²〉)^1/2 (4)

where Δr_i and Δr_j are the displacement vectors for Cα atoms of residues i and j, respectively, and the angle brackets denote the ensemble averages.^39,54 The computational details of the INTERHLX procedure, the B-factors and the motion correlations are given in the ESI.†

3. Results

The root-mean-square deviation (RMSD) values of all the backbone atoms relative to the corresponding starting structures over two trajectories for the AFF-D_9k-N′ and AFF-D_9k-N models were examined to determine the system equilibrium. Plots of RMSDs for the two equilibrium systems, AFF-D_9k-N′ and AFF-D_9k-N models, over simulation times are shown in Fig. 2. As illustrated in Fig. 2, the AFF-D_9k-N′ and AFF-D_9k-N systems reached equilibrium after 35 ns, with stable energies during the remainder of each simulation. Therefore, the trajectory analysis of the two systems yielded the equilibrated conformations between 40 ns and 50 ns simulation times, recording 2500 snapshots at every 4 ps time-interval of each trajectory. Fig. 3 illustrates the plots for the RMSDs of all backbone atoms in the TMD-simulated structures relative to the final structure of transition of AFF-D_9k-N′ to AFF-D_9k-N. The plots depicted in Fig. 3 reveal that the backbone atoms in the TMD simulations all reached the target conformations within 8 ns. The detailed structural transition results obtained from the CMD and TMD simulations are discussed below in detail.

Fig. 2 Root mean square deviation (RMSD) values of all backbone atoms of the AFF-D_9k-N′ and AFF-D_9k-N models with respect to the corresponding starting structures for their simulations.

Fig. 3 Root mean square deviation (RMSD) values of all backbone atoms and all atoms of the whole protein from the targeted molecular dynamics simulation are shown for the AFF-D_9k-N′ to AFF-D_9k-N transition.

3.1. The stable conformations of the AFF-D_9k-N′ and AFF-D_9k-N models

Based on the CMD simulations, we obtained two dynamically stable states, AFF-D_9k-N′ and AFF-D_9k-N models that respectively correspond to N-terminal-folded and N-terminal-unfolded conformations of the AFF calbindin-D_9k protein. The corresponding average structures of the two AFF-D_9k-N′ and AFF-D_9k-N models extracted from the CMD trajectories of their simulations were analyzed and are shown in Fig. 4(a) and (b), respectively. The superpositions for these two average structures and the WT calbindin-D_9k protein in Ca²⁺-free, shown in Fig. 4(c) and (d), represent the structural characteristics of the experimental design of the AFF calbindin-D_9k protein.²⁶ The features from Fig. 4 are that both structures with an identical sequence mainly present one conserved EF1 hand, two mutually exclusive folding regions (EF2 and EF2′) and a linker between the EF1 hand and EF2′ region, i.e., for the AFF-D_9k-N′ model, the EF1 hand (including A helix (Lys1-Ala15), loop_A–B region (Lys16-Ser24), B helix (Lys25-Glu35)), loop_B–C region (Phe36-Met43), the folded EF2′ hand (including C′ helix (Ser44′-Leu53′), loop_C′–D′ region (Asp54′-Ser62′), D′ helix (Phe63′-Gln75′)), the unfolded EF2 region (including the same sequence as in the folded EF2′ hand excepted for the Glu65Gln mutation, C segment (Ser44-Leu53), loop_C–D region (Asp54-Ser62), and D segment (Phe63-Gln75)) and the linker (Gly76-Met81); for the AFF-D_9k-N model, the couple structure with the same sequence as in the AFF-D_9k-N′ model including the conserved EF1 hand located at the middle of the AFF-D_9k-N model, the unfolded EF2′ region, the folded EF2 hand and the linker. The distances between the Cα atoms of fluorescence position residues Ser44′ and Met43 in the AFF-D_9k-N′ and AFF-D_9k-N models are respectively 5.0 Å and 41.5 Å, which supports the experimental results of the fluorescence data.^26,27 Moreover, the mass center distances between EF2 and EF2′ regions are 21.2 Å and 26.3 Å in the AFF-D_9k-N′ and AFF-D_9k-N models, respectively, which predicts the similarity of the couple mutually exclusive folding structures.

Fig. 4 CMD-simulated average structures for (a) the AFF-D_9k-N′ state and (b) the AFF-D_9k-N state; (c) and (d) show the superpositions for these two average structures and the WT calbindin-D_9k protein (white, PDB-ID 3ICB) in Ca²⁺-free, respectively. The four identical components of EF2′, EF1, EF2 and linker in AFF-D_9k-N′ and AFF-D_9k-N states are colored in magenta, cyan, orange and black, respectively.

The free energy difference between the AFF-D_9k-N′ and AFF-D_9k-N states using the MM_PBSA methodology from the MD simulations only is 0.98 kcal mol⁻¹, which indicates that the unfolding and folding structures of EF2/EF2′ region in the two models insignificantly influences the stability of two states due to the coupling of folding and unfolding segments. This result also supports the experimental observation of the approximately isoenergetic existence for the two N′ and N conformations of the AFF calbindin-D_9k protein.²⁷ Details and equations of free energy calculations are shown in the ESI and Fig. S3.†

3.2. Free energy landscape and pathway of the conformational transition between the AFF-D_9k-N′ and AFF-D_9k-N models

Based on the structural features of the AFF-D_9k-N′ and AFF-D_9k-N models, we performed the TMD simulation to explore the conformational transition process between the two models. We first computationally examined the conformational transition from AFF-D_9k-N′ to AFF-D_9k-N starting from the equilibrated AFF-D_9k-N′ structure. The transition process from AFF-D_9k-N′ to AFF-D_9k-N by TMD simulation involves the shifts of the EF2 region and EF2′ hand firstly, then their corresponding folding and unfolding. As noted above, the TMD-simulated protein structures were used as the initial structures for further umbrella-sampling MD simulations in order to determine the free energy profile. The forward PMF calculations by the WHAM program after the umbrella-sampling MD simulations for the transformation from AFF-D_9k-N′ to AFF-D_9k-N were determined and are depicted in Fig. 5(a). The distance between the mass centers of the EF1 hand and EF2 region was used as the reaction coordinate along the transition process for the PMF calculations. Because both conformations of the protein are comprised of a well-folded domain and a full-unfolded region with the identical sequence, the AFF-D_9k-N′/AFF-D_9k-N conformational change belongs to a mutually exclusive folding reaction, i.e. the folding of one segment of the protein is coupled to unfolding of another segment. The reverse transition from the AFF-D_9k-N to AFF-D_9k-N′ still involves the shifts of the EF2 hand and EF2′ region firstly, then their corresponding unfolding and folding, which is shown in Fig. S4 of the ESI.† The reverse process is the reaction cycle as similar to the forward process. Therefore, the PMF calculation for the reverse process of AFF-D_9k-N to AFF-D_9k-N′ was expected as similar to that for the forward one.

Fig. 5 (a) Calculated free energy profile for the transition process of AFF-D_9k-N′ to AFF-D_9k-N. The reaction coordinate was defined as the distances between the mass centers of the EF1 and EF2 hands. The roseate letters denote the different conformations. (b) The three-dimensional average structures for the corresponding conformations of AFF-D_9k-N′, E, TS1, Im, TS2 and AFF-D_9k-N presented in (a).

According to the free energy profile depicted in Fig. 5(a), the conformational transition from AFF-D_9k-N′ to AFF-D_9k-N mainly involves two allosteric transition steps via two free energy barriers and a local minimum separating them. Fig. 5(b) shows the four important conformations E, TS1, Im and TS2 in the free-energy landscape corresponding to the disordered reactant of AFF-D_9k-N′, the conformation at the free energy top in the first barrier, the intermediate in the local minimum and the another conformation at the free energy top in the second barrier, respectively. In the first transition step from AFF-D_9k-N′ to Im via E and TS1 conformations, the unstructured tail of EF2 region in AFF-D_9k-N′ reactant gradually move up; then the EF2′ hand moves from the right side of EF1 hand to its left side through the back of the EF1 hand to produce the conformation E, which favors the movement of EF2 region up to the right side of EF1 hand (see Fig. 5(b)). Moreover, the C′ helix in the EF2′ hand occurs the slight disorder during this mutually movement. As shown in Fig. 5(a), the free energy from AFF-D_9k-N′ to E only increases 1.0 kcal mol⁻¹ for this position shift. Then, the C′ helix of EF2′ hand in the E conformation extends to a loop structure; simultaneously, the EF2 region coils partially itself, which causes the deformation of EF1 hand with the decrease of the hydrogen bonds in the EF2′ region and EF1 hand, to produce the first transition state TS1. As shown in Fig. 5(a), the free energy increases by 4.7 kcal mol⁻¹ from the reactant AFF-D_9k-N′ to TS1 to form the first barrier. Then, the EF2 region continually coils itself to form partial helix with the formations of some new hydrogen bonds, simultaneously the slight recovery of the deformation of the middle EF1 hand, to produce the intermediate (Im) associated with a local minimum on the free energy profile. The free energy from TS1 to Im lowers 3.2 kcal mol⁻¹ due to the partial helix formation of the EF2 region and the recovery of the middle EF1 hand structure. In the second transition step from Im to the product AFF-D_9k-N via TS2, the partial-coiled EF2′ region unfolds continually; simultaneously, the A helix in EF1 hand unfolds partially in Im to form the second transition state (TS2) with a free energy increase of 1.7 kcal mol⁻¹. Then, the partial-coiled EF2 region in TS2 coils continually to form the EF2 hand to produce the product AFF-D_9k-N with the formations of more hydrogen bonds and the decrease of free energy of 2.4 kcal mol⁻¹ (see Fig. 5(a)). Notably, the free energy difference of 0.61 kcal mol⁻¹ calculated by the PMF method between the AFF-D_9k-N′ and AFF-D_9k-N states is consistent with that of 0.98 kcal mol⁻¹ by the MM_PBSA calculations. This result predicted that the current calculation protocols are sufficiently accurate for describing the kind of systems. Summarily, the transition process from AFF-D_9k-N′ to AFF-D_9k-N involves: the shifts of EF2 region and EF2′ hand firstly with the slight increase of free energy; then the extension of the C′ helix in EF2′ hand and the deformation of EF1 hand causing the formation of the rate-controlling barrier of 4.7 kcal mol⁻¹ that is comparable with the energy requirement involving the interactions of calcium with C-lobe of the EF-hand in calmodulin protein;⁵⁵ then the coiling of the EF2 region to form an Im; the unfolding of EF2′ region and the A helix in EF1 hand resulting in the formation of the second barrier of 1.7 kcal mol⁻¹ relative to Im. These results predict that the extension of the EF2′ hand in TS1 causing the considerable deformation of middle EF1 hand, the partial helix formation in the EF2 region of Im with the slight recovery of EF1 hand and the continual unfolding of the partial-coiled EF2′ region in TS2 resulting in the deformation of A helix of EF1 hand again indicate the mediating function of the deformations of the EF1 hand as a mediating agent during the transition from AFF-D_9k-N′ to AFF-D_9k-N. That is, the transition from AFF-D_9k-N′ to AFF-D_9k-N is a mutually exclusive folding reaction via a coupled mediated mechanism. Considerably, the deformation of the middle EF1 hand during this conversion instead of a whole-molecule unfolding further supports the previous experimental kinetic data.²⁷

3.3. Dynamical fluctuation and correlation analyses of the transition from AFF-D_9k-N′ to AFF-D_9k-N

To further interpret the conformational transition characteristics of the AFF calbindin-D_9k protein in Ca²⁺-free, the dynamics of every residue was investigated by residue fluctuation and correlation methods. To test the consistency of AFF calbindin-D_9k protein with the WT calbindin-D_9k protein in Ca²⁺-free, the B-factor (the temperature factor) values for the EF1 and EF2′/EF2 hands (the components of the WT calbindin-D_9k protein) computed from the root-mean-square fluctuations (RMSFs) and the simulation trajectories for the two AFF-D_9k-N and AFF-D_9k-N′ models were analyzed and are shown in Fig. 6(a) and (b), respectively. The B-factor values for the EF1 and EF2 hands in the AFF-D_9k-N model in Fig. 6(a) are qualitatively similar to those from the experimental crystal structure of the WT calbindin-D_9k protein due to the complete identity of the sequences involved in the EF1 and EF2 hands in the AFF-D_9k-N model and the WT calbindin-D_9k protein. For example, the peak between Phe36 and Met43 residues predicts the loop fluctuation between two EF hands; and another small peak between Asp54 and Ser62 residues shows the slight fluctuation of the loop between C and D helices in the EF2 hand. These results are consistent with those from the experimental crystal structure. For the AFF-D_9k-N′ model, the unidentities of the sequences of the EF1 and EF2′ hands in the AFF-D_9k-N′ model and the WT calbindin-D_9k protein include the non-connection of two residues of Met43 and Ser44′ between EF1 and EF2′ hands, N-terminal located at the C′ helix and the connection of Gln75′ residue (that is the C-terminal in the WT calbindin-D_9k protein) to Gly76 of the additional linker. Therefore, based on ignoring these insignificant residues (Phe36-Met43) nearing the connecting loop of EF1 hand and EF2 region, and the additional linker, the B-factor values for the EF1 and EF2′ hands in the AFF-D_9k-N′ model in Fig. 6(b) are similar to those from the experimental crystal structure of the WT calbindin-D_9k protein. Especially, the peak between Asp54 and Ser62 indicating the considerable fluctuation of loop at the middle of C′ and D′ helices in EF2′ hand, is higher than that in the WT calbindin-D_9k protein due to the loop region of Asp54-Ser62 nearby the N-terminal of AFF-D_9k-N′ model. Moreover, the B-factors nearing Gln75′ residue in the AFF-D_9k-N′ model are greatly smaller than those in the WT calbindin-D_9k protein due to the change of C-terminal in the AFF-D_9k-N′ model. Other residues have low B-factor values and are very stable during the simulation. No significant differences in the B-factors are detectable among the EF1 and EF2′/EF2 hands for these two models.

Fig. 6 The calculated residue B-factors of (a) AFF-D_9k-N (red) and (b) AFF-D_9k-N′ (blue) in comparison with the experimental B-factors of the crystal calbindin-D_9k protein (black), respectively. The peaks with the number of residues at the loop regions are marked in grey.

To explore the allosteric communications of mutually exclusive folding process in the EF1 hand and EF2′/EF2 region for the transition from the AFF-D_9k-N′ to AFF-D_9k-N conformations, we constructed and analyzed the motion correlations for all Cα atoms of AFF-D_9k-N′ to AFF-D_9k-N from the simulation trajectories. They are displayed in Fig. 7(a), (b), (c) and (d) for the stable AFF-D_9k-N′ model, the stable AFF-D_9k-N model, TS1 and Im, respectively. These maps show the high motion correlations between the residues. For the AFF-D_9k-N′ model, the motions of C′ and D′ helices in the EF2′ hand significantly correlate with the motions of A and B helices in the EF1 hand represented by the black squares in Fig. 7(a). The uncorrelations of the EF2 region with other regions present its unfolding conformation and its position far from the two EF hands in the AFF-D_9k-N′ model. Similarly, the correlations and uncorrelations of the AFF-D_9k-N model occur between the C/D helix in the EF2 hand and A/B helix in the EF1 hand represented by the black squares in Fig. 7(b), and between the EF2 region and other regions, respectively. However, the correlated motions of TS1 conformation in Fig. 7(c) are insignificant due to the loose structure of all EF regions in TS1. Furthermore, because from TS1 transiting to Im the unfolding conformation of TS1 become more coiled structure in Im, the corresponding correlations of Im in Fig. 7(d) become more considerable. For example, the motions of the D region in EF2 significantly correlate with the motions of the B helix in the EF1 hand represented by the black square in Fig. 7(d) due to the partial coil of D region in Im. As illustrated above, all correlations in these states involve the correlated motions of the EF1 hand with other mutually exclusive folding regions, which supports that the EF1 hand plays a crucial role in mediating the transition process of AFF-D_9k-N′ to AFF-D_9k-N.

Fig. 7 Dynamical cross-correlation maps for (a) the AFF-D_9k-N′ model, (b) the AFF-D_9k-N model, (c) the TS1 conformation and (d) the Im conformation with the key subregions squared in black.

4. Discussion

4.1. Conformational characteristics for the AFF-D_9k-N′ and AFF-D_9k-N models

Base on the average structure and the trajectory analysis, the two stable conformations of AFF-D_9k-N′ and AFF-D_9k-N models for the AFF calbindin-D_9k protein in Ca²⁺-free form have mutually exclusive folding structures. To further investigate the detail conformational characteristics at the conserved middle region of the EF1 hand of the two models compared with that of the WT calbindin-D_9k protein in Ca²⁺-free, the interhelical angles and the mass center distances between A and B helices in the EF1 hand in the AFF-D_9k-N′ and AFF-D_9k-N models have been quantitatively examined by using the program INTERHLX and are shown in Fig. 8. The calculated data indicate that the deviation values of the average interhelical angles and the mass center distances between A and B helices for the AFF-D_9k-N′, AFF-D_9k-N models compared with that of the WT calbindin-D_9k protein are 2°, 3° and 0.2 Å, 0.4 Å, respectively, which suggested that the EF1 hand in the two models present the conserved characteristics that reproduces the experimental observation. Simultaneously, the results suggest the two designed conformations of AFF-D_9k-N′ and AFF-D_9k-N using the AFF approach and MD simulations are reasonable.

Fig. 8 (a) Scheme of the interhelical angle and distance between A and B helices in EF1 hand; the interhelical angles (b) and distances (c) between A and B helices in EF1 hand along with respective integrated distributions for the AFF-D_9k-N′ (black) model, AFF-D_9k-N (red) model, and the experimental values (blue lines) of the WT calbindin-D_9k protein.

4.2. Allosteric process analysis from AFF-D_9k-N′ to AFF-D_9k-N

To understand fully the allosteric communication of mutually exclusive folding reaction from the AFF-D_9k-N′ to AFF-D_9k-N models, we explored the relationship between the motion correlations and their structural changes for this allosteric process mediated by the EF1 hand. It can be seen from the motion correlation analysis in Fig. 7 that the C′/D′ helix of the EF2′ hand correlate to the A/B helix of the EF1 hand with large cross-correlation coefficients at the beginning of the proceeding of mutually exclusive folding (see Fig. 7(a)); then, at the end of this proceeding the strong correlations change to the C/D helix of the EF2 hand with the A/B helix of the EF1 hand (see Fig. 7(b)). As expected, with the proceeding of this allosteric process from the reactant to TS1 then to Im finally to the product, the changes of correlations gradually occur from the correlations of the C′/D′ helix in the EF2′ hand and the A/B helix of the EF1 hand changing to that of the C/D helix in the EF2 hand and the A/B helix of the EF1 hand (see Fig. 7(a)–(d)). Such correlation transition reasonably represents the allosteric communication of structures from the C′ and D′ helices of the EF2′ hand to the C and D helices of the EF2 hand through the mediation of A/B helix of the EF1 hand.

To explore the details of the structural changes along the transition process from AFF-D_9k-N′ to AFF-D_9k-N, we used the middle EF1 hand with the conservation of mass center as the reference point to calculate the variations in mass center distances between the EF2 region and the EF1 hand, between the EF2′ and the EF1 hands, and between the EF2′ hand and EF2 region. The corresponding data are shown in Fig. 9. The occurrences of hydrogen bonds of the EF2′, EF1 and EF2 regions in the AFF-D_9k-N′, TS1, Im and TS2 conformations were also analyzed and are shown in Table 1 by calculating the percentages of times on the sampling CMD simulations of the corresponding conformations. It can be seen during the allosteric process that the distances between the EF2′ and EF1 hands gradually increase from 10.5 Å to 19.2 Å; and that the distances between the EF2 region and the EF1 hand gradually decrease from 23.1 Å to 10.3 Å (see the black and blue lines in Fig. 9). Such mobility of the EF2 region and EF2′ hand was also observed from the average structures of AFF-D_9k-N′, E, TS1, Im, TS2 and AFF-D_9k-N states during the conformational transition in Fig. 5(b). Moreover, the distances between the EF2′ hand and EF2 region change from 21.2 Å at the beginning time of the transition simulation to the top value of 30.8 Å in the middle of E and TS1 conformations to 25.9 Å at the end of transition process (see the red line in Fig. 9), which represents the allosteric characteristics of the rotating shifts of EF2′ and EF2 along the ellipse axis of the EF1 center as discussed below.

Fig. 9 Mass center distances between the EF2 region and the EF1 hand (blue), between the EF2′ and the EF1 hands (black), and between the EF2′ hand and EF2 region (red) along the transition pathway from the simulation of AFF-D_9k-N′ to AFF-D_9k-N; the roseate letters included in the grey areas denote the AFF-D_9k-N′, E, TS1, Im, TS2 and AFF-D_9k-N conformations.

Table 1 Occupancies (%) of hydrogen bonds for the AFF-D_9k-N′, TS1, Im and TS2 conformations

Hydrogen bond	AFF-D9k-N′	TS1	Im	TS2
(Leu53′)N–H⋯O(Phe50′)	93.74	—	—	—
(Leu53′)N–H⋯O(Leu49′)	—	—	99.97	81.27
(Lys71′)N–H⋯O(Gln67′)	99.33	99.98	—	—
(Leu69′)N–H⋯O(Phe66′)	90.41	—	—	—
(Leu69′)N–H⋯O(Gln65′)	—	—	99.67	—
(Gln33)N–H⋯O(Lys29)	99.87	—	—	—
(Leu31)N–H⋯O(Glu27)	99.87	—	—	99.91
(Glu11)N–H⋯O(Lys7)	98.67	—	—	—
(Phe10)N–H⋯O(Lys7)	—	—	99.90	—
(Ile9)N–H⋯O(Glu5)	98.67	—	—	—
(Phe36)N–H⋯O(Leu32)	96.94	—	—	—
(Lys7)N–H⋯O(Pro3)	95.74	—	—	100
(Glu35)N–H⋯O(Leu31)	100	—	—	—
(Thr34)N–H⋯O(Leu30)	91.08	—	—	—
(Gly8)N–H⋯O(Glu4)	86.15	—	99.87	—
(Thr34)N–H⋯O(Leu31)	69.24	—	99.99	99.98
(Lys7)N–H⋯O(Glu4)	—	—	99.99	—
(Leu31)N–H⋯O(Leu28)	—	—	98.68	—
(Leu30)N–H⋯O(Glu26)	96.94	—	—	—
(Glu52)N–H⋯O(Phe50)	—	—	—	94.45
(Asn56)N–H⋯O(Leu53)	—	—	71.68	99.99
(Ser74)N–H⋯O(Val70)	—	—	57.55	96.67

---

In detail, for the first conformational transition step from AFF-D_9k-N′ to Im via E and TS1 conformations, the distances between the EF2′ and EF1 hands first keep at 10.5 Å during the simulation beginning time; oppositely, the distances between the EF2 region and EF1 hand decrease from 23.1 Å in AFF-D_9k-N′ to 17.0 Å in E during simulation times of 0∼3 ns due to the unmoving of stable EF2′ hand and the up-moving of the flexible EF2 loop structure at the transition beginning. Then, the distances between the EF2′ and EF1 hands increase to 14.6 Å at E spending 1.5 ns simulation time, which presents the fast moving of EF2′ hand from the right side of the EF1 hand to its left side through the back of the EF1 hand. Otherwise, the distances between the EF2′ hand and EF2 region fast increase from 19.8 Å to 26.8 Å at E spending 1 ns simulation time also supports the fast moving of EF2′ hand (see red line in Fig. 9). Moreover, the average interhelical angles of the C′ and D′ helices in the EF2′ hand at the two average conformations of AFF-D_9k-N′ and E were measured with the INTERHLX program and are shown in Fig. 10. The difference value of this average interhelical angle is 59° for the two average conformations of AFF-D_9k-N′ and E, which suggested the rotating shift of EF2′ hand from the right side of the EF1 hand to its left side. Then in the transition of E to TS1, although the increase value of distances between the EF2′ and EF1 hands and the decrease value of the distances between the EF2 region and EF1 hand are similar, the distances between the EF2′ hand and EF2 region are still changed from 26.8 Å at E to the top of 30.8 Å and to 29.0 Å in TS1, which presents the rotating shifts of EF2′ and EF2 along the ellipse axis of the EF1 center due to the deformation mediation of EF1 hand, and the loose structure of TS1. Simultaneously, the deformation mediation of EF1 hand and unfolding of EF2′ hand result in the decrease of hydrogen bonds by 93% total occupancy and the increase of free energy from AFF-D_9k-N′ to TS1. For example, the hydrogen bond between the N–H group of Leu53′ and the O atom of Phe50′ in C′ helix of the EF2 hand is maintained with the occupancy of 93.74% of simulation times in the AFF-D_9k-N′ model, while it disappeared in TS1 (see Table 1). From TS1 to Im, the formation of the new partial-coiled helix at the EF2 region and the slight recovery of the deformational A and B helices at the middle EF1 hand cause the increase of hydrogen bonds by 87% total occupancy. For example, the new hydrogen bonds between the N–H group of Ser74 and the O atom of Val70, and between N–H group of Asn56 and the O atom of Leu53 in EF2 region are formed with the occupancy of 57.55% and 71.68% of simulation times in Im, respectively (see Table 1). For the second conformational transition step from Im to the product AFF-D_9k-N via TS2 conformation, the continually unfolding of EF2′ hand from Im to TS2 causing insignificantly the changes of the distances between the EF2′ and EF1 hands indicates the unfolding balance of EF2′ hand; while the fluctuation of distances between the EF2′ hand and EF2 region presents the unfolding characteristics of the A helix in the EF1 hand as a mediated factor and the partial coiling of EF2 region. These structural changes cause the decrease of hydrogen bonds by 20% total occupancy and the increase of free energy from Im to TS2. For example, the hydrogen bond between the N–H group of Gly8 and the O atom of Glu4 in the A helix of the EF1 hand is maintained with the occupancy of 99.87% of simulation times in Im; while it disappeared in TS2 (see Table 1). From TS2 to the product AFF-D_9k-N, the decrease of distances between EF2 region and the EF1 hand from 13.0 Å in TS2 to 10.3 Å in AFF-D_9k-N further supports the considerable coiling of EF2 region to form EF2 hand to produce the product AFF-D_9k-N with the increase of hydrogen bonds and the decrease of free energy. Moreover, the calculated hydrogen bonds in the reactant AFF-D_9k-N′ involving the residues Phe50′, Ser62′ and Phe66′ in EF2′ hand reproduce the experimental results.²⁸

Fig. 10 The difference of the interhelical angles between the C′ and D′ helices of EF2′ hand for the AFF-D_9k-N′ (colored by white) and E (colored by magenta) states in the transition process.

5. Conclusions

Conventional molecular dynamics simulations, targeted molecular dynamics simulations and energy calculations from the potential of mean force have been performed to address the stable structure characteristics for the AFF-D_9k-N′ and AFF-D_9k-N conformations of the AFF calbindin-D_9k protein in Ca²⁺-free form, and the dynamic transition mechanism of the two conformations. The structural characteristics of the two stable AFF-D_9k-N′ and AFF-D_9k-N states represent the coupled conformations of mutually exclusive folding. The investigated results for the transition mechanism demonstrated that the structural transition from AFF-D_9k-N′ to AFF-D_9k-N could occur via two transition states and an intermediate with the first rate-controlling barrier of 4.7 kcal mol⁻¹ and the second barrier of 1.7 kcal mol⁻¹. Therefore, this mutually exclusive folding reaction is energetically feasible due to the low rate-controlling barrier. The calculated free energies for the two stable AFF-D_9k-N′ and AFF-D_9k-N states by the MM_PBSA method are in agreement with the experimental results of the approximately isoenergetic existence. The energy fluctuations are influenced by the EF2′ hand and EF2 region movements, the mutually exclusive folding of EF2′ hand and EF2 region, and the deformation of the EF1 hand. Analysis of the hydrogen bonds with the total occurrences of simulation times decreased by 93% in TS1 relative to the reactant AFF-D_9k-N′ supports the free energy increase of 4.7 kcal mol⁻¹ as the rate-controlling barrier. On the other hand, the total occurrence of hydrogen bonds in TS2 decreased by 20% relative to Im corresponds to forming the second barrier of 1.7 kcal mol⁻¹. The structural change characteristics along the transition pathway indicate the relative movements of the two mutually exclusive folding regions, EF2 and EF2′, at the transition beginning, then the unfolding/folding procedures of EF2′ hand/EF2 region mediated by the deformation of EF1 hand. Moreover, the variations of the calculated mass center distances between the EF hands and the interhelical angles in EF2′ hand also supported the conformational changes during the conformational transition. Especially, the mass center distances between the EF2′ hand and EF2 region changed from 26.8 Å at E to the top of 30.8 Å, and to 29.0 Å in TS1 suggested the rotating shifts of EF2′ and EF2 along the ellipse axis of the EF1 center due to the deformation mediation of EF1 hand. The correlation analysis reveals that the allosteric communications between the C′/D′ helix of the EF2′ hand and A/B helix of the EF1 hand for the reactant AFF-D_9k-N′, between the D helix of the EF2 hand and A/B helix of the EF1 hand at Im, and between the C/D helix of the EF2 hand and A/B helix of the EF1 hand at the product AFF-D_9k-N further verified the mediating role of the EF1 hand in this transition process. The present investigations provide useful insights into understanding the dynamics mechanism of the mutually exclusive folding reaction from AFF-D_9k-N′ to AFF-D_9k-N states for the AFF calbindin-D_9k protein in Ca²⁺-free.

Acknowledgements

The authors acknowledge research support from the National Science Foundation of China (Nos 21271029, 21131003, 21073015, and 20973024), the Major State Basic Research Development Programs (Grant No. 2011CB808500).

References

1. X. Luo, Z. Tang, G. Xia, K. Wassmann, T. Matsumoto, J. Rizo and H. Yu, Nat. Struct. Mol. Biol., 2004, 11, 338–345.
2. M. Hilbert, F. Kebbel, A. Gubaev and D. Klostermeier, Nucleic Acids Res., 2011, 39, 2260–2270.
3. T. Kenakin and L. J. Miller, Pharmacol. Rev., 2010, 62, 265–304.
4. A. W. Fenton, Trends Biochem. Sci., 2008, 33, 420–425.
5. A. J. Björkman and S. L. Mowbray, J. Mol. Biol., 1998, 279, 651–664.
6. K. P. Ravindranathan, E. Gallicchio and R. M. Levy, J. Mol. Biol., 2005, 353, 196–210.
7. M. Zhang, T. Tanaka and M. Ikura, Nat. Struct. Biol., 1995, 2, 758–767.
8. M. M. Stratton and S. N. Loh, Protein Sci., 2011, 20, 19–29.
9. J.-H. Ha and S. N. Loh, Chem.–Eur. J., 2012, 18, 7984–7999.
10. M. A. Wilson and A. T. Brunger, J. Mol. Biol., 2000, 301, 1237–1256.
11. K. P. Hoeflich and M. Ikura, Cell, 2002, 108, 739–742.
12. C. K. Johnson, Biochemistry, 2006, 45, 14233–14246.
13. N. Valeyev, D. Bates, P. Heslop-Harrison, I. Postlethwaite and N. Kotov, BMC Syst. Biol., 2008, 2, 48.
14. M. Schumacher, M. Crum and M. Miller, Structure, 2004, 12, 849–860.
15. Z. Grabarek, J. Mol. Biol., 2005, 346, 1351–1366.
16. L. Maler, J. Blankenship, M. Rance and W. J. Chazin, Nat. Struct. Mol. Biol., 2000, 7, 245–250.
17. K. Julenius, J. Robblee, E. Thulin, B. E. Finn, R. Fairman and S. Linse, Proteins, 2002, 47, 323–333.
18. L. A. Svensson, E. Thulin and S. Forsén, J. Mol. Biol., 1992, 223, 601–606.
19. S. Christakos, C. Gabrielides and W. B. Rhoten, Endocr. Rev., 1989, 10, 3–26.
20. D. M. E. Szebenyi and K. Moffat, J. Biol. Chem., 1986, 261, 8761–8777.
21. N. J. Skelton, J. Kördel and W. J. Chazin, J. Mol. Biol., 1995, 249, 441–462.
22. G. M. C. Mitsuhiko Ikura, A. M. Gronenborn, G. Zhu, C. B. Klee and A. Bax, Science, 1992, 256, 632–638.
23. G. B. M. Ikura, C. B. Klee and A. Bax, Cell Calcium, 1992, 13, 391–400.
24. S. Marchand and B. Roux, Proteins, 1998, 33, 265–284.
25. K. L. Yap, J. B. Ames, M. B. Swindells and M. Ikura, Proteins, 1999, 37, 499–507.
26. D. M. M. Margaret, M. Stratton and S. N. Loh, ACS Chem. Biol., 2008, 3, 723–732.
27. M. M. Stratton and S. N. Loh, Proteins, 2010, 78, 3260–3269.
28. M. M. Stratton, S. McClendon, D. Eliezer and S. N. Loh, Biochemistry, 2011, 50, 5583–5589.
29. H. Zheng, J. Bi, M. Krendel and S. N. Loh, Biochemistry, 2014, 53, 5505–5514.
30. J.-H. Ha, S. A. Shinsky and S. N. Loh, Biochemistry, 2013, 52, 600–612.
31. D. M. Mitrea, L. S. Parsons and S. N. Loh, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 2824–2829.
32. T. A. Cutler, B. M. Mills, D. J. Lubin, L. T. Chong and S. N. Loh, J. Mol. Biol., 2009, 386, 854–868.
33. V. J. Hilser and E. B. Thompson, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 8311–8315.
34. T. A. Cutler and S. N. Loh, J. Mol. Biol., 2007, 371, 308–316.
35. J. S. Butler, D. M. Mitrea, G. Mitrousis, G. Cingolani and S. N. Loh, Biochemistry, 2009, 48, 3497–3507.
36. Y. Kiyota and M. Takeda-Shitaka, J. Phys. Chem. B, 2014, 118, 11496–11503.
37. N. F. Barley, S. R. Prathalingam, P. Zhi, S. Legon, A. Howard and J. R. F. Walters, Biochem. J., 1999, 341, 491–500.
38. D. Bucher, B. J. Grant, P. R. Markwick and J. A. McCammon, PLoS Comput. Biol., 2011, 7, e1002034.
39. T. A. D. D. A. Case, T. E. Cheatham III, C. L. Simmerling, J. Wang, R. E. Duke, R. K. M. M. Luo, D. A. Pearlman, M. Crowley, R. C. Walker, W. Zhang, B. Wang, S. A. R. Hayik, G. Seabra, K. F. Wong, F. Paesani, X. Wu, S. Brozell, V. Tsui, H. L. Y. Gohlke, C. Tan, J. Mongan, V. Hornak, G. Cui, P. Beroza, D. H. Mathews, C. W. S. R. Schafmeister and P. A. Kollman, AMBER 9, University of California, San Francisco, CA, 2006.
40. M. C. Lee and Y. Duan, Proteins, 2004, 55, 620–634.
41. Y. Duan, C. Wu, S. Chowdhury, M. C. Lee, G. Xiong, W. Zhang, R. Yang, P. Cieplak, R. Luo, T. Lee, J. Caldwell, J. Wang and P. Kollman, J. Comput. Chem., 2003, 24, 1999–2012.
42. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, J. Comput. Chem., 2004, 25, 1157–1174.
43. J. Schlitter, M. Engels and P. Krüger, J. Mol. Graphics Modell., 1994, 12, 84–89.
44. J. Schlitter, M. Engels, P. Krüger, E. Jacoby and A. Wollmer, Mol. Simul., 1993, 10, 291–308.
45. G. N. Patey and J. P. Valleau, J. Chem. Phys., 1975, 63, 2334–2339.
46. C. Chipot, P. A. Kollman and D. A. Pearlman, J. Comput. Chem., 1996, 17, 1112–1131.
47. B. Roux, Comput. Phys. Commun., 1995, 91, 275–282.
48. G. M. Torrie and J. P. Valleau, J. Comput. Phys., 1977, 23, 187–199.
49. C. Bartels and M. Karplus, J. Comput. Chem., 1997, 18, 1450–1462.
50. S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen and P. A. Kollman, J. Comput. Chem., 1992, 13, 1011–1021.
51. J. B. A. Kyoko, L. Yap, M. B. Swindells and M. Ikura, Methods Mol. Biol., 2002, 173, 317–324.
52. G. A. Petsko and D. Ringe, Annu. Rev. Biophys. Bioeng., 1984, 13, 331–371.
53. U. D. Priyakumar and A. D. MacKerell Jr, J. Mol. Biol., 2010, 396, 1422–1438.
54. S. K. Sadiq and G. De Fabritiis, Proteins, 2010, 78, 2873–2885.
55. M. Kohagen, M. Lepšík and P. Jungwirth, J. Phys. Chem. Lett., 2014, 5, 3964–3969.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5ra11234f

---