Mapping transiently formed and sparsely populated conformations on a complex energy landscape

Abstract
eLife digest
Introduction
Results and discussion
Materials and methods
Appendix
References
Article and author information
Metrics

Abstract

Determining the structures, kinetics, thermodynamics and mechanisms that underlie conformational exchange processes in proteins remains extremely difficult. Only in favourable cases is it possible to provide atomic-level descriptions of sparsely populated and transiently formed alternative conformations. Here we benchmark the ability of enhanced-sampling molecular dynamics simulations to determine the free energy landscape of the L99A cavity mutant of T4 lysozyme. We find that the simulations capture key properties previously measured by NMR relaxation dispersion methods including the structure of a minor conformation, the kinetics and thermodynamics of conformational exchange, and the effect of mutations. We discover a new tunnel that involves the transient exposure towards the solvent of an internal cavity, and show it to be relevant for ligand escape. Together, our results provide a comprehensive view of the structural landscape of a protein, and point forward to studies of conformational exchange in systems that are less characterized experimentally.

https://doi.org/10.7554/eLife.17505.001

eLife digest

Proteins are the workhorses of cells, where they perform a wide range of roles. To do so, they adopt specific three-dimensional structures that enable them to interact with other molecules as necessary. Often a protein needs to be able to shift between different states with distinct structures as it goes about its job. To fully understand how a protein works, it is important to be able to characterize these different structures and how the protein changes between them.

Many of the experimental techniques used to study protein structure rely on isolating the individual structural forms of a protein. Since many structures only exist briefly, this can be very difficult. To complement experimental results, computer simulations allow researchers to model how atoms behave within a molecule. However, a number of factors limit how well these models represent what happens experimentally, such as the accuracy of the physical description used for the modeling.

Wang et al. set out to test and benchmark how well computer simulations could model changes in structure for a protein called T4 lysozyme, which has been studied extensively using experimental techniques. T4 lysozyme exists in two different states that have distinct structures. By comparing existing detailed experimental measurements with the results of their simulations, Wang et al. found that the simulations could capture key aspects of how T4 lysozyme changes its shape.

The simulations described the structure of the protein in both states and accurately determined the relative proportion of molecules that are found in each state. They could also determine how long it takes for a molecule to change its shape from one state to the other. The findings allowed Wang et al. to describe in fine detail – down to the level of individual atoms – how the protein changes its shape and how mutations in the protein affect its ability to do so. A key question for future studies is whether these insights can be extended to other proteins that are less well characterized experimentally than T4 lysozyme.

https://doi.org/10.7554/eLife.17505.002

Introduction

Proteins are dynamical entities whose ability to change shape often plays essential roles in their function. From an experimental point of view, intra-basin dynamics is often described via conformational ensembles whereas larger scale (and often slower) motions are characterized as a conformational exchange between distinct conformational states. The latter are often simplified as a two-site exchange process, $G ⇌ E$ , between a highly populated ground (G) state, and a transiently populated minor (or ‘excited’, E) state. While the structure of the ground state may often be determined by conventional structural biology tools, it is very difficult to obtain atomic-level insight into minor conformations due to their transient nature and low populations. As these minor conformations may, however, be critical to protein functions, including protein folding, ligand binding, enzyme catalysis, and signal transduction (Mulder et al., 2001; Tang et al., 2007; Baldwin and Kay, 2009) it is important to be able to characterize them in detail. While it may in certain cases be possible to capture sparsely populated conformations in crystals under perturbed experimental conditions, or to examine their structures by analysis of electron density maps (Fraser et al., 2009), NMR spectroscopy provides unique opportunities to study the dynamical equilibrium between major and minor conformations (Baldwin and Kay, 2009; Sekhar and Kay, 2013) via e.g. chemical-exchange saturation transfer (Vallurupalli et al., 2012), Carr-Purcell-Meiboom-Gill (CPMG) relaxation dispersion (Hansen et al., 2008), or indirectly via paramagnetic relaxation enhancement (Tang et al., 2007) or residual dipolar coupling (Lukin et al., 2003) experiments. In favourable cases such experiments can provide not only thermodynamic and kinetic information (i.e. the population of G and E states and the rate of exchange between them), but also structural information in the form of chemical shifts (CS), that can be used to determine the structure of the transiently populated state (Sekhar and Kay, 2013).

Despite the important developments in NMR described above, it remains very difficult to obtain structural models of minor conformations, and a substantial amount of experiments are required. Further, it is generally not possible to use such experiments to infer the mechanisms of interconversion, and to provide a more global description of the multi-state free energy landscape (Zhuravlev and Papoian, 2010; Wang et al., 2012). In the language of energy landscape theory (Onuchic et al., 1997), free energy basins and their depths control the population and stability of functionally distinct states, while the relative positions of basins and the inter-basin barrier heights determine the kinetics and mechanism of conformational exchange. As a complement to experiments, such functional landscapes can be explored by in silico techniques, such as molecular dynamics (MD) simulations, that may both be used to help interpret experimental data and provide new hypotheses for testing (Karplus and Lavery, 2014; Eaton and Muñoz, 2014). Nevertheless, the general applicability of simulation methods may be limited by both the accuracy of the physical models (i.e. force fields) used to describe the free energy landscape and our ability to sample these efficiently by computation. We therefore set out to benchmark the ability of simulations to determine conformational free energy landscapes.

The L99A variant of lysozyme from the T4 bacteriophage (T4L) has proven an excellent model system to understand protein structure and dynamics. Originally designed a ‘cavity creating’ variant to probe protein stability (Eriksson et al., 1992b) it was also demonstrated that the large (150 Å³) internal cavity can bind hydrophobic ligands such as benzene (Eriksson et al., 1992a; Liu et al., 2009). It was early established that the cavity is inaccessible to solvent in the ground state, but that ligand binding is rapid (Feher et al., 1996), suggesting protein dynamics to play a potential role in the binding process. This posts a long-standing question of how the ligands gain access to the buried cavity (Mulder et al., 2000; López et al., 2013; Merski et al., 2015; Miao et al., 2015).

NMR relaxation dispersion measurements of L99A T4L demonstrated that this variant, but not the wild type protein, displayed conformational exchange on the millisecond timescale between the ground state and a minor state populated at around 3% (at room temperature) (Mulder et al., 2001). Such small populations generally lead only to minimal perturbations of ensemble-averaged experimental quantities making structural studies difficult, and hence it was difficult to probe whether the exchange process indeed allowed for ligand access to the cavity. A series of additional relaxation dispersion experiments, however, made it possible to obtain backbone and side chain CSs of the minor E state of L99A (Mulder et al., 2002; Bouvignies et al., 2011). The backbone CS data were subsequently used as input to a CS-based structure refinement protocol (CS-ROSETTA) to produce a structural model of the E state ( $E_{R O S E T T A}$ ; Figure 1) of the L99A mutant (Bouvignies et al., 2011). This model was based in part on the crystal structure of the ground state of L99A (referred to in what follows as $G_{X r a y}$ ), but perturbing the structure in regions that the experiments demonstrated to undergo conformational change in a way so that the final model ( $E_{R O S E T T A}$ ) agrees with experiments. The structure was further validated by creating and solving the structure of a triple mutant variant that inverts the populations of the G and E states. The $E_{R O S E T T A}$ structure revealed substantial local rearrangements in T4L L99A, in particular near the cavity which gets filled by the side chain of a phenylalanine at position 114 ( $F_{114}$ ). Because the cavity is filled and solvent inaccessible in the E-state, the structure did, however, not reveal how ligands might access the cavity.

Figure 1

Download asset Open asset

Structures of the major G and minor E states of L99A T4L and the hidden state hypothesis.

The X-ray structure of the G state ( $G_{X r a y}$ ; PDB ID code 3DMV) has a large internal cavity within the core of the C-terminal domain that is able to bind hydrophobic ligands. The structure of the E state ( $E_{R O S E T T A}$ ; PDB ID code 2LC9) was previously determined by CS-ROSETTA using chemical shifts. The G and E states are overall similar, apart from the region surrounding the internal cavity. Comparison of the two structures revealed two remarkable conformational changes from G to E: helix F (denoted as $H_{F}$ ) rotates and fuses with helix G ( $H_{G}$ ) into a longer helix, and the side chain of phenylalanine at position 114 ( $F_{114}$ ) rotates so as to occupy part of the cavity. As the cavity is inaccessible in both the $G_{X r a y}$ and $E_{R O S E T T A}$ structures it has been hypothesized that ligand entry occurs via a third ‘cavity open’ state (Merski et al., 2015).

https://doi.org/10.7554/eLife.17505.003

In an attempt to benchmark the ability of simulations to map conformational free energy landscapes, we have here employed a series of in silico experiments designed to probe the structure and dynamics of L99A T4L and have compared the results to NMR measurements. We used enhanced-sampling MD simulations in explicit solvent and with state-of-the-art force fields to map the free-energy landscape including the exchange between the major and minor conformations of the protein. We used a series of recently developed metadynamics methods (Laio and Parrinello, 2002) to sample the conformational exchange process and associated structure and thermodynamics, as well as to determine the kinetics and mechanisms of exchange. We obtained additional insight into the structural dynamics of the E state using simulations that employed the experimental CSs as replica-averaged restraints. Our results provide a coherent picture of the conformational dynamics in L99A and extend the insights obtained from recent simulations of a triple mutant of T4L (Vallurupalli et al., 2016), by providing new information about the mechanisms of exchange and the transient exposure of the internal cavity. Together with previous results for Cyclophilin A (Papaleo et al., 2014) the results described here reiterate how simulation methods have now reached a stage where they can be used to study slow, conformational exchange processes such as those probed by NMR relaxation dispersion even in cases where less information is available from experiments.

Results and discussion

Mapping the free-energy landscape

As the average lifetime of the G and E states are on the order of 20–50 ms and 1 ms, respectively (Mulder et al., 2001, 2002; Bouvignies et al., 2011), direct and reversible sampling of the G-E transition at equilibrium would be extremely demanding computationally. Indeed, a recent set of simulations of a triple mutant of T4L, which has a substantially faster kinetics, was able only to observe spontaneous transitions in one direction (Vallurupalli et al., 2016). We therefore resorted to a set of flexible and efficient enhanced sampling methods, collectively known as ‘metadynamics’ (Laio and Parrinello, 2002), that have previously been used in a wide range of applications. In metadynamics simulations, a time-dependent bias is continuously added to the energy surface along a small number of user-defined collective variables (CVs). In this way, sampling is enhanced to reach new regions of conformational space and at the same time allows one to reconstruct the (Boltzmann) free-energy surface. The success of the approach hinges on the ability to find a set of CVs that together describe the slowly varying degrees of freedom and map the important regions of the conformational landscape.

We first performed a set of metadynamics simulations in the well-tempered ensemble (Barducci et al., 2008) using so-called path CVs ( $S_{p a t h}$ and $Z_{p a t h}$ ) (Branduardi et al., 2007; Saladino et al., 2012) with the aid of recently developed adaptive hills to aid in the convergence of the sampling (Branduardi et al., 2012; Dama et al., 2014) (see details in Appendix and Appendix 1—table 1). In short, the $S_{p a t h}$ variable describes the progress of the conformational transition between the $G_{X r a y}$ and $E_{R O S E T T A}$ structures with additional ‘interpolation’ using an optimal ‘reference’ path in a simplified model (see details in Appendix and Figure 2—figure supplement 1), while $Z_{p a t h}$ measures the distance to this reference path. In this way, the two-dimensional free energy landscape along $S_{p a t h}$ and $Z_{p a t h}$ provides a useful description on conformational exchange between ground and excited states that does not assume that the initial reference path describes perfectly the actual path(s) taken.

Projecting the sampled free energy landscape along $S_{p a t h}$ (upper panel of Figure 2) reveals a deep, narrow free energy basin around $S_{p a t h} = 0.2$ (labeled by red sphere and corresponding to the G state), and a broader, shallow free energy basin with $S_{p a t h}$ ranging from 0.6 to 0.8 (labeled by blue sphere and corresponding to the E state). Additional information is obtained from the two-dimensional landscape (shown as a negative free energy landscape, -F( $S_{p a t h}$ , $Z_{p a t h}$ ), in the lower panel of Figure 2) which reveals a complex and rough landscape with multiple free energy minima (corresponding to mountains in the negative free energy landscape). Subsequently, structural inspection of these minima identified that the conformations in the basins around $S_{p a t h} = 0.2$ and $S_{p a t h} = 0.75$ correspond to the structures of $G_{X r a y}$ and $E_{R O S E T T A}$ , respectively.

Figure 2 with 2 supplements see all

Download asset Open asset

Free energy landscape of the L99A variant of T4L.

In the upper panel, we show the projection of the free energy along $S_{p a t h}$ , representing the Boltzmann distribution of the force field employed along the reference path. Differently colored lines represent the free energy profiles obtained at different stages of the simulation, whose total length was 667ns. As the simulation progressed, we rapidly found two distinct free energy basins, and the free energy profile was essentially constant during the last 100 ns of the simulation. Free energy basins around $S_{p a t h} = 0.2$ and $S_{p a t h} = 0.75$ correspond to the previously determined structures of the G- and E-state, respectively (labelled by red and blue dots, respectively). As discussed further below, the E-state is relatively broad and is here indicated by the thick, dark line with $S_{p a t h}$ ranging from 0.55 to 0.83. In the lower panel, we show the three-dimensional negative free energy landscape, -F( $S_{p a t h}$ , $Z_{p a t h}$ ), that reveals a more complex and rough landscape with multiple free energy minima, corresponding to mountains in the negative free energy landscape. An intermediate-state basin around $S_{p a t h} = 0.36$ and $Z_{p a t h} = 0.05$ nm², which we denote $I_{0.36}$ , is labeled by a yellow dot.

https://doi.org/10.7554/eLife.17505.004

The broad nature of the free energy landscape in the region of the minor state is consistent with the observation that our MD simulations initiated from $E_{R O S E T T A}$ display significant conformational fluctuations (RUN20 and RUN22 in Appendix 1—table 1). Furthermore, our metadynamics simulations revealed multiple local free energy minima adjacent to the $E_{R O S E T T A}$ basin, together composing a wider basin (highlighted by the black curve in Figure 2). Thus, these simulations suggest that the E state displays substantial conformational dynamics, a result corroborated by simulations that have been biased by the experimental data (see section ‘Simulations of the minor state using chemical shift restraints’).

In addition to free-energy minima corresponding to the $G$ and $E$ states, we also found a free energy minimum around $S_{p a t h} = 0.36$ and $Z_{p a t h} = 0.05$ nm² (denoted as $I_{0.36}$ and labeled by a yellow sphere in Figure 2) that is located between the G and E states on the one-dimensional free-energy surface. We note, however, that it is difficult to infer dominant reaction pathways from such free energy surfaces, and so from this data alone, we cannot determine whether $I_{0.36}$ occurs as an intermediate in G-E conformational transitions. Indeed, it appears from the two-dimensional surface that there exist multiple possible pathways between G and E, as illustrated by white lines along the mountain ridges of the negative free energy landscape in the lower panel of Figure 2. (We also explored the mechanism of exchange by reconnaissance metadynamics simulations (Tribello et al., 2011), the results of which are described and discussed further below.)

Effect of mutations on the free energy landscape

Based on the encouraging results above for L99A T4L, we examined whether simulations could also capture the effect of mutations on the free energy landscape. Using Rosetta energy calculations on the $G_{X r a y}$ and $E_{R O S E T T A}$ structures it was previously demonstrated that two additional mutations, G113A and R119P, when introduced into the L99A background, cause an inversion in the populations of the two states (Bouvignies et al., 2011; Vallurupalli et al., 2016). Indeed, NMR data demonstrated that the triple mutant roughly inverts the populations of the two states so that the minor state structure (of L99A) now dominates (with a 96% population) the triple mutant. We repeated the calculations described above for L99A also for the triple mutant. Remarkably, the free energy profile of the triple mutant obtained using metadynamics simulations reveals a free energy landscape with a dominant minimum around $S_{p a t h}$ =0.7 and a higher energy conformation around $S_{p a t h}$ =0.15 (Figure 2—figure supplement 2). Thus, like our previous observations for a ‘state-inverting mutation’ in Cyclophilin A (Papaleo et al., 2014), we find here that the force field and the sampling method are sufficiently accurate to capture the effect of point mutations on the free energy landscape. Further, we note that the barrier height for the conformational exchange in the triple mutant is very similar to the value recently estimated using a completely orthogonal approach (Vallurupalli et al., 2016). Finally, we attempted to determine the free energy landscape of the L99A,G113A double mutant, which has roughly equal populations of the two states (Bouvignies et al., 2011), but this simulation did not converge on the simulation timescales at which the two other variants converged.

Calculating conformational free energies

With a free-energy surface in hand and a method to distinguish G- and E-state conformations, we calculated the free energy difference, $Δ G$ , between the two conformational states, and compared with the experimental values. We divided the global conformational space into two coarse-grained states by defining the separatrix at $S_{p a t h} = 0.46$ which corresponds to a saddle point on the free energy surface, on the basis of the observations above that the E state is relatively broad. Although a stricter definition of how to divide the reaction coordinate certainly helps the precise calculation, here we just used this simple definition to make an approximate estimation of the free energy difference. Further, since the barrier region is sparsely populated, the exact point of division has only a modest effect on the results. By summing the populations on the two sides of the barrier we calculated $Δ G$ as a function of the simulation time (Figure 3). Initially during the simulations the free energy profile varies substantially (Figure 2) and the free energy difference equally fluctuates. As the simulations converge, however, the free energy difference between the two states stabilize to a value at approximately $Δ G$ =3.5 kcal mol⁻¹ (Figure 3, black line). This value can be compared to the value of 2.1 kcal mol⁻¹ obtained from NMR relaxation dispersion experiments (Mulder et al., 2001), revealing reasonably good, albeit not exact, agreement with the experiments.

Figure 3

Download asset Open asset

Estimation of free energy differences and comparison with experimental measurements.

We divided the global conformational space into two coarse-grained states by defining the separatrix at $S_{p a t h} = 0.46$ (0.48 for the triple mutant) in the free energy profile (Figure 2—figure supplement 2) which corresponds to a saddle point of the free energy surface, and then estimated the free energy differences between the two states ( $Δ G$ ) from their populations. The time evolution of $Δ G$ of L99A (upper time axis) and the triple mutant (lower axis) are shown as black and blue curves, respectively. The experimentally determined values (2.1 kcal mol⁻¹ for L99A and −1.9 kcal mol⁻¹ for the triple mutant) are shown as yellow lines.

https://doi.org/10.7554/eLife.17505.007

Similar calculations using the simulations of the triple mutant also converge, in this case to about −1.6 kcal mol⁻¹ (Figure 3, blue line), in excellent agreement with the experimental measurement (−1.9 kcal mol⁻¹) (Bouvignies et al., 2011). Combining these two free energy differences we find that the G113A, R119P mutations cause a shift in the G-E free energy of 5.1 kcal mol⁻¹ in simulations compared to 4.0 kcal mol⁻¹ obtained by experiments. Thus, we find that the simulations with reasonably high accuracy are able to capture the thermodynamics of the conformational exchange between the two states. While the generality of such observations will need to be established by additional studies, we note here that comparably good agreement was obtained when estimating the effect of the S99T mutations in Cyclophilin A (Papaleo et al., 2014).

In our previous work on Cyclophilin A (Papaleo et al., 2014) we sampled the conformational exchange using parallel-tempering metadynamics simulations (Bonomi and Parrinello, 2010) using four CVs that we chose to describe the structural differences between the G and E states in that protein. We note here that we also tried a similar approach here but unfortunately failed to observe a complete G-to-E transition, even in a relatively long trajectory of about 1 μs per replica (CVs summarized in Appendix 1—table 2, parameters shown in Appendix 1—table 1). This negative result is likely due to the CVs chosen did not fully capture the relevant, slowly changing degrees of freedom, thus giving rise to insufficient sampling even with the use of a parallel tempering scheme.

Calculating the rates of conformational exchange

Enhanced-sampling simulations such as those described above provide an effective means of mapping the free-energy landscape and hence the structural and thermodynamic aspects of conformational exchange. While the same free-energy landscape also determines the kinetics and mechanisms of exchange it may be more difficult to extract this information from e.g. path-CV-based metadynamics (PathMetaD) simulations. To examine how well simulations can also be used to determine the rates of the G-to-E transitions, quantities that can also be measured by NMR, we used the recently developed ‘infrequent metadynamics’ method (InMetaD, see details in Appendix) (Tiwary and Parrinello, 2013; Salvalaglio et al., 2014; Tiwary et al., 2015a; 2015b). Briefly described, the approach calculates first-passage times for the conformational change in the presence of a slowly-added bias along a few CVs, here chosen as the path CVs also used to map the landscape. By adding the bias slowly (and with lower amplitude) we aim to avoid biasing the transition-state region and hence to increase the rate only by lowering the barrier height; in this way it is possible to correct the first-passage times for the bias introduced.

Using this approach on L99A T4L we collected 42 and 36 independent trajectories with state-to-state transition starting from either the $G$ state or $E$ state, respectively (Appendix 1—figure 1 and Appendix 1—figure 2 ). The (unbiased) rates that we calculated (Table 1 and Appendix 1—figure 3) are in good agreement with the experimental rates (Mulder et al., 2001; Bouvignies et al., 2011) (within a factor of 10), corresponding to an average error of the barrier height of $\sim$ 1 kcal mol⁻¹. We also performed similar calculations for the ‘population-inverting’ triple mutant, where we collected 30 transitions (15 for each direction) using InMetaD simulations. As for L99A, we also here find similarly good agreement with experimental measurements (Vallurupalli et al., 2016) (Table 1 and Appendix 1—figure 4). We estimated the reliability of this computational approach using a Kolmogorov-Smirnov test to examine whether the first-passage times conform to the expected Poisson distribution (Salvalaglio et al., 2014), and indeed the results of this analysis suggest good agreement (Table 1, Appendix 1—figure 5 and Appendix 1—figure 6).

Table 1

Free energy differences and rates of conformational exchange.

https://doi.org/10.7554/eLife.17505.008

	$τ_{G \to E}$ (ms)	$τ_{E \to G}$ (ms)	$Δ G$ (kcal mol⁻¹)
L99A
NMR	20	0.7	2.1
InMetaD	175 $\pm$ 56	1.4 $\pm$ 0.6	2.9 $\pm$ 0.5
PathMetaD			3.5
L99A/G113A/R119P
NMR	0.2	4	-1.9
InMetaD	2.0 $\pm$ 1.7	14.3 $\pm$ 8.3	-1.2 $\pm$ 1.1
PathMetaD			-1.6

The ability to calculate forward and backward rates between $G$ and $E$ provided us with an alternative and independent means to estimate the free energy difference between the two states (Table 1), and to test the two-state assumption used in the analysis of the experimental NMR data. We therefore calculated the free energy difference from the ratio of the forward and backward reaction rates. The values obtained (2.9 $\pm$ 0.5 kcal mol⁻¹ and −1.2 $\pm$ 1.1 kcal mol⁻¹ for L99A and the triple mutant, respectively) are close both to the values obtained above from the equilibrium free energy landscape (3.5 kcal mol⁻¹ and −1.6 kcal mol⁻¹) and experiment (2.1 kcal mol⁻¹ and −1.9 kcal mol⁻¹). In particular, the relatively close agreement between the two independent computational estimates lends credibility both to the free energy landscape and the approach used to estimate the kinetics. The observation that both values for L99A are slightly larger than the experimental number suggests that this discrepancy (ca. 1 kcal mol⁻¹) can likely be explained by remaining force field deficiencies rather than lack of convergence or the computational approach used.

Simulations of the minor state using chemical shift restraints

While the simulations described above used available structural information of G and E states to guide and enhance conformational sampling, the resulting free energy surfaces represent the Boltzmann distributions of the force field and are not otherwise biased by experimental data. To further refine the structural model of the E state we used the relaxation-dispersion derived CSs that were used to determine of $E_{R O S E T T A}$ [BMRB (Ulrich et al., 2008) entry 17604] as input to restrained MD simulations. In these simulations, we used the experimental data as a system-specific force field correction to derive an ensemble of conformations that is compatible both with the force field and the CSs. Such replica-averaged simulations use the experimental data in a minimally-biased way that is consistent with the Principle of Maximum Entropy (Pitera and Chodera, 2012; Roux and Weare, 2013; Cavalli et al., 2013; Boomsma et al., 2014).

We performed CS-restrained MD simulations of the $E$ state of L99A averaging the CSs over four replicas. Although the number of replicas is a free parameter, which should in principle be chosen as large as possible, it has been demonstrated that four replicas are sufficient to reproduce the structural heterogeneity accurately (Camilloni et al., 2013) without excessive computational requirements. The agreement between calculated and experimental CSs was quantified by the root-mean-square deviation between the two (Figure 4—figure supplement 1). In particular, it is important not to bias the agreement beyond what can be expected based on the inherent accuracy of the CS prediction methods (we assumed that the error in the experimental CS measurement, even for the $E$ state, is negligible in comparison). Thus, we compared the experimental CS values of the minor state with the values calculated using the $E_{R O S E T T A}$ structure as input to CamShift (Kohlhoff et al., 2009), Sparta+ (Shen and Bax, 2010) and ShiftX (Neal et al., 2003) (Figure 4—figure supplement 2). The average RMSDs for five measured nuclei ( $H_{α}$ , $H_{N}$ , $N$ , $C^{'}$ and $C_{α}$ ) are 0.2, 0.4, 2.0, 0.8 and 1.1ppm, respectively (Appendix 1—table 1), which are close to the inherent uncertainty of the CS back-calculation, indicating that the level of agreement enforced is reasonable.

To compare the results of these experimentally-biased simulations with the experimentally-unbiased simulations described above, we projected the CS-restrained MD trajectories onto either one (Figure 4) or both (Figure 4—figure supplement 3) of the $S_{p a t h}$ and $Z_{p a t h}$ variables used in the path-variable-driven simulations (PathMetaD). The distribution of conformations obtained using the E-state CSs as restraints is in good agreement with the broad free energy profile of the E-state obtained in the metadynamics simulations that did not include any experimental restraints. To ensure that this observation is not merely an artifact of both simulations using the same force field (CHARMM22*), we repeated the biased simulations using the Amber ff99SB*-ILDN force field and obtained comparable results. We also verified that the conclusions obtained are reasonably robust to other variables such as the number of replicas and the strength of restraints (Figure 4—figure supplement 4).

Figure 4 with 5 supplements see all

Download asset Open asset

Conformational ensemble of the minor state as determined by CS biased, replica-averaged simulations.

We determined an ensemble of conformations corresponding to the E-state of L99A T4L using replica-averaged CSs as a bias term in our simulations. The distribution of conformations was projected onto the $S_{p a t h}$ variable (orange) and is compared to the free energy profile obtained above from the metadynamics simulations without experimental biases (black line). To ensure that the similar distribution of conformations is not an artifact of using the same force field (CHARMM22*) in both simulations, we repeated the CS-biased simulations using also the Amber ff99SB*-ILDN force field (magenta) and obtained similar results. Finally, we used the ground state CSs of a triple mutant of T4L, which was designed to sample the minor conformation (of L99A) as its major conformation, and also obtained a similar distribution along the $S_{p a t h}$ variable (cyan).

https://doi.org/10.7554/eLife.17505.009

As a final and independent test of the structural ensemble of the minor conformation of L99A we used the ground state CSs of the triple mutant (BMRB entry 17603), which corresponds structurally to the E state of L99A, as restraints in replica-averaged CS-biased simulations (Figure 4—figure supplement 5). Although not fully converged, these simulations also cover roughly the same region of conformational space when projected along $S_{p a t h}$ (Figure 4).

Thus, together our different simulations, which employ different force fields, are either unbiased or biased by experimental data, and use either dispersion-derived (L99A) or directly obtained (triple mutant) CS all provide a consistent view of the minor E-state conformation of L99A. We also note that the CS-derived ensembles of the E-state support the way we divided the G- and E-states when calculating conformational free energy differences between the two states.

Mechanisms of conformational exchange

Having validated that our simulations can provide a relatively accurate description of the structure, thermodynamics and kinetics of conformational exchange we proceeded to explore the molecular mechanism of the G-to-E transitions. We used the recently developed reconnaissance metadynamics approach (Tribello et al., 2010), that was specifically designed to enhance sampling of complicated conformational transitions and has been employed to explore the conformational dynamics of complex systems (Tribello et al., 2011; Söderhjelm et al., 2012).

We performed three independent reconnaissance metadynamics simulations of L99A starting from the G state (summarized in Appendix 1—table 1) using the same geometry-based CVs that we also used in the parallel-tempering simulations described above. We observed complete conformational transitions from the G to E state in the reconnaissance simulations in as little as tens of nanoseconds of simulations (Figure 5—figure supplement 1)— at least 1–2 orders of magnitude faster than standard metadynamics. These G-to-E and E-to-G transitions, although biased by the CVs, provide insight into the potential mechanisms of exchange. To ease comparison with the equilibrium sampling of the free energy landscape we projected these transitions onto the free energy surface F( $S_{p a t h}$ , $Z_{p a t h}$ ) (Figure 5). The results reveal multiple possible routes connecting the G and E states, consistent with the multiple gullies found on the free energy surface (Figure 2). The trajectories also suggested that the G-to-E interconversion can either take place directly without passing the $I_{0.36}$ state or indirectly via it.

Figure 5 with 3 supplements see all

Download asset Open asset

Mechanisms of the G-E conformational exchanges explored by reconnaissance metadynamics.

Trajectories labeled as Trj1 (magenta), Trj2 (blue) and Trj3 (green and orange) are from the simulations RUN10, RUN11 and RUN12 (Appendix 1—table 1), respectively. There are multiple routes connecting the G and E states, whose interconversions can take place directly without passing the $I_{0.36}$ state or indirectly via it.

https://doi.org/10.7554/eLife.17505.015

In the context of coarse-grained kinetic models the results above would suggest at least two possible mechanisms operate in parallel: $G ⇌ E$ or $G ⇌ I_{0.36} ⇌ E$ . Further inspection of the structures along these different kinetics routes (see the trajectories of other order parameters in Figure 5—figure supplement 2 and Videos 1–4) suggested an interesting distinction between the two. In the $G ⇌ I_{0.36} ⇌ E$ route the side chain of F114, which occupies the cavity in the E state, gets transiently exposed to solvent during the transition, whereas in the direct $G ⇌ E$ transitions F114 can rotate its side chain inside the protein core (see also the solvent accessible surface area calculation of F114 in Figure 5—figure supplement 3).

Video 1

Download asset

posterframe for video — Trajectory of the G-to-E conformational transition observed in Trj1, corresponding to the red trajectory in Figure 5.

The backbone of L99A is represented by white ribbons, Helices E, F and G are highlighted in blue, while F114 is represented by red spheres.

https://doi.org/10.7554/eLife.17505.019

Video 2

Download asset

Video 3

Download asset

Video 4

Download asset

A potential pathway for ligand binding and escape

As the internal cavity in L99A T4L remains buried in both the G and E states (and indeed occupied by F114 in the E state) it remains unclear how ligands access this internal cavity and how rapid binding and release is achieved. Visual inspection of our trajectories and solvent-accessible surface area analysis revealed structures with transient exposure of the internal cavity towards the solvent. The structures were mostly found in a region of conformational space that mapped onto the $I_{0.36}$ basin (Figure 2), and the events of that basin mostly took place between 430 ns and 447 ns (see Video 5). Thus, we mapped these structures to the free energy surface (Figure 6—figure supplement 1) and analysed them. Overall, the structure is more similar to the G- than E-state, though is more loosely packed. The similarity to the G-state is compatible with rapid binding and position of F114 in this state.

Video 5

Download asset

We used CAVER3 (Chovancova et al., 2012) (see parameters in Appendix 1—table 4) to analyse the structures and found multiple tunnels connecting the cavity with protein surface (Figure 6—figure supplement 1 and 2). The tunnels are relatively narrow with the typical radius of the bottleneck (defined as the narrowest part of a given tunnel) between ∼1 Å − ∼2 Å. We used CAVER Analyst1.0 (Kozlikova et al., 2014) (see details in Appendix and parameters in Appendix 1—table 4) to separate the tunnels into different clusters (Figure 6—figure supplement 3 and Appendix 1—table 5) with the dominant cluster (denoted tunnel $# 1$ ) having a entrance located at the groove between $H_{F}$ and $H_{I}$ . A typical representative structure of $I_{0.36}$ is shown in Figure 6A. The radii along the structures in cluster $# 1$ vary, but share an overall shape (Figure 6—figure supplement 1), and we find that the maximal bottleneck radius is $\sim$ 2.5 Å, the average bottleneck radius is $\sim$ 1.3 Å, and the average length $\sim$ 11.2 Å.

Figure 6 with 4 supplements see all

Download asset Open asset

A transiently formed tunnel from the solvent to the cavity is a potential ligand binding pathway.

(A) We here highlight the most populated tunnel structure (tunnel $# 1$ ), that has an entrance located at the groove between helix F ( $H_{F}$ ) and helix I ( $H_{I}$ ). Helices E, F and G (blue) and F114 (red) are highlighted. (B) The panel shows a typical path of benzene (magenta) escaping from the cavity of L99A, as seen in ABMD simulations, via a tunnel formed in the same region as tunnel #1 (see also Video 6).

https://doi.org/10.7554/eLife.17505.024

Interestingly, a series of structures of L99A were recently described, in which the internal cavity where filled with eight congeneric ligands of increasing size to eventually open the structure size (Merski et al., 2015). We performed a comparable tunnel analysis on those eight ligand-bound structures (PDB ID codes: 4W52 – 4W59), revealing the maximal bottleneck radius of 1.8 Å (bound with n-hexylbenzene, 4W59). Although the size of the tunnel in these X-ray structures is slightly smaller than that in $I_{0.36}$ structures, the location of the tunnel exit is consistent with the dominant tunnel $# 1$ in $I_{0.36}$ (Figure 6—figure supplement 3). We note, however, that the tunnels observed in our simulation and in the ligand-induced cavity-open X-ray structure (4W59), are too narrow to allow for unhindered passage of e.g. benzene with its a van der Waals’ width of 3.5 Å (Eriksson et al., 1992a). Thus, we speculate that the transient exposure in $I_{0.36}$ might serve as a possible starting point for ligand (un)binding, which would induce (Koshland, 1958; López et al., 2013; Wang et al., 2012) further the opening of the tunnel.

As an initial step towards characterizing the mechanism of ligand binding and escape we used adiabatic biased molecular dynamics (ABMD) simulations (Marchi and Ballone, 1999; Paci and Karplus, 1999) to study the mechanism of how benzene escapes the internal cavity (see Appendix for details). In ABMD the system is perturbed by a ‘ratcheting potential’, which acts to ‘select’ spontaneous fluctuations towards the ligand-free state. In particular, the biasing potential is zero when the reaction coordinate (here chosen to be the RMSD of the ligand to the cavity-bound state) increases, but provides a penalty for fluctuations that brings the ligand closer to the cavity. In this way, we were able to observe multiple unbinding events in simulations despite the long lifetime (1.2 ms) of the ligand in the cavity. Most of trajectories (15 of the 20 events observed) reveal that benzene escapes from the cavity by following tunnel $# 1$ (Figure 6—figure supplement 4 and Appendix 1—table 6). A typical unbinding path is shown in the right panel of Figure 6 (see also Video 6). Because the ABMD introduces a bias to speed up ligand escape, we ensured that the observed pathway was the same at two different values of the biasing force constants (Figure 6—figure supplement 4 and Appendix 1—table 6). Future work will be aimed to perform a more quantitative analysis of the ligand binding and unbinding kinetics.

Video 6

Download asset

Conclusions

The ability to change shape is an essential part of the function of many proteins, but it remains difficult to characterize alternative conformations that are only transiently and sparsely populated. We have studied the L99A variant of T4L as a model system that displays a complicated set of dynamical processes which have been characterized in substantial detail. Our results show that modern simulation methods are able to provide insight into such processes, paving the way for future studies for systems that are more difficult to study experimentally.

Using a novel method for defining an initial reference path between two conformations, we were able to sample the free energy landscape described by an accurate molecular force field. In accordance with experiments, the simulations revealed two distinct free energy basins that correspond to the major and minor states found by NMR. Quantification of the free energy difference between the two states demonstrated that the force field is able to describe conformational free energies to an accuracy of about 1 kcal mol⁻¹. This high accuracy is corroborated by previous studies of a different protein, Cyclophilin A, where we also calculated conformational free energies and compared to relaxation dispersion experiments and found very good agreement. For both proteins we were also able to capture and quantify the effect that point mutations have on the equilibrium between the two states, and also here found good agreement with experiments. We note, however, that comparable simulations of the L99A/G113A mutant did not reach convergence.

Moving a step further, we here also calculated the kinetics of conformational exchange using a recently developed metadynamics method. For both the L99A variant and a population-inverting triple mutant we find that the calculated reaction rates are in remarkably good agreement with experiments. The ability to calculate both forward and backward rates provided us with the opportunity to obtain an independent estimate the calculated free energy difference. The finding that the free energy differences estimated in this way (for both L99A and the triple mutant) are close to those estimated from the free energy landscape provides an important validation of both approaches, and we suggest that, when possible, such calculations could be used to supplement conventional free energy estimates.

The free-energy landscape suggested that the E state is relatively broad and contains a wider range of conformations. To validate this observation, we used the same chemical shift information as was used as input to Rosetta and performed replica-averaged CS-restrained simulations. The resulting ensemble demonstrates that the experiments and force field, when used jointly, indeed are compatible with a broader E state. Thus, we suggest that the $E_{R O S E T T A}$ structure and CS-restrained ensemble jointly describe the structure and dynamics of the E state.

While NMR experiments, in favourable cases, can be used to determine the structure, thermodynamics and kinetics of conformational exchange, a detailed description mechanism of interconversion remains very difficult to probe by experiments. We explored potential mechanisms of conformational exchange between the two states, finding at least two distinct routes. One route involved a direct transition with the central F114 entering the cavity within the protein, whereas a different possible mechanism involves transient partial-loosening of the protein. In both cases, the mechanism differs from the reference path that we used as a guide to map the free energy landscape, demonstrating that high accuracy of the initial guess for a pathway is not absolutely required in the metadynamics simulations, suggesting also the more general applicability of the approach.

Finally, we observed a set of conformations with a transiently opened tunnel that leads from the exterior of the protein to the internal cavity, that is similar to a recently discovered path that is exposed when the cavity is filled by ligands of increasing size. The fact that such a tunnel can be explored even in the absence of ligands suggests that intrinsic protein motions may play an important role in ligand binding, and indeed we observed this path to be dominant in simulations of ligand unbinding.

In total, we present a global view of the many, sometimes coupled, dynamical processes present in a protein. Comparison with a range of experimental observations suggests that the simulations provide a relatively accurate description of the protein, demonstrating how NMR experiments can be used to benchmark quantitatively the ability of simulations to study conformational exchange. We envisage that future studies of this kind, also when less is known about the structure of the alternative states, will help pave the way for using simulations to study the structural dynamics of proteins and how this relates to function.

Materials and methods

System preparation

Request a detailed protocol

Our simulations were initiated in the experimentally determined structures of the ground state of L99A ( $G_{X r a y}$ ; PDB ID code 3DMV) or minor, E state ( $E_{R O S E T T A}$ ; 2LCB). The structure of the ground state of the L99A, G113A, R119P triple mutant, corresponding to the E state of L99A was taken from PDB entry 2LC9 ( $G_{R O S E T T A}^{T r i p l e}$ ). Details can be found in the Appendix.

Initial reaction path

Request a detailed protocol

Taking $G_{X r a y}$ and $E_{R O S E T T A}$ as the models of the initial and final structures, we calculated an initial reaction path between them with the MOIL software (Elber et al., 1995), which has been used to explore the mechanism of conformational change of proteins (Wang et al., 2011). Further details can be found in the Appendix and in refs. (Majek et al., 2008; Wang et al., 2011).

Path CV driven metadynamics simulations with adaptive hills

Request a detailed protocol

The adaptive-hill version of metadynamics updates the Gaussian width on the fly according to the local properties of the underlying free-energy surface on the basis of local diffusivity of the CVs or the local geometrical properties. Here, we used the former strategy. Simulation were performed using Gromacs4.6 (Pronk et al., 2013) with the PLUMED2.1 plugin (Tribello et al., 2014). See parameter details in Appendix 1—table 1.

Replica-averaged CS-restrained simulations

Request a detailed protocol

We performed replica-averaged CS restrained MD simulations by using GPU version of Gromacs5 with the PLUMED2.1 and ALMOST2.1 (Fu et al., 2014) plugins. Equilibrated structures of $E_{R O S E T T A}$ and $G_{R O S E T T A}^{T r i p l e}$ were used as the starting conformations. CS data of $E_{R O S E T T A}$ and $G_{R O S E T T A}^{T r i p l e}$ were obtained from the BMRB database (Ulrich et al., 2008) as entries 17604 and 17603, respectively.

Reconnaissance metadynamics simulations

Request a detailed protocol

Reconnaissance metadynamics (Tribello et al., 2010) uses a combination of a machine learning technique to automatically identify the locations of free energy minima by periodically clustering the trajectory and a dimensional reduction technique that can reduce the landscape complexity. We performed several reconnaissance metadynamics simulations with different combinations of CVs starting from $G_{X r a y}$ using Gromacs4.5.5 with PLUMED1.3 plugin. See parameter details in Appendix 1—table 1.

Calculating kinetics using infrequent metadynamics

Request a detailed protocol

The key idea of infrequent metadynamics is to bias the system with a frequency slower than the barrier crossing time but faster than the slow intra-basin relaxation time, so that the transition state region has a low risk of being substantially biased. As the first transition times should obey Poisson statistics, the reliability of the kinetics estimated from InMetaD can be assessed by a statistical analysis based on the Kolmogorov-Smirnov (KS) test (Salvalaglio et al., 2014). See parameter details on Appendix and Appendix 1—table 1.

Appendix

Molecular modeling and system preparation

We used the crystal structure of T4L L99A (PDB ID code 3DMV) as starting point for simulations of the G state of L99A. For the E state of L99A and the G state of the L99A, G113A, R119P triple-mutant we used the CS-ROSETTA structures with PDB ID code 2LCB and 2LC9, respectively.

Each protein was solvated in a dodecahedral box of TIP3P water molecules with periodic boundary conditions. The protein-solvent box had a distance of 10 Å from the solute to the box boundary in each dimension, which results in approximately 10,000 water molecules and more than 32,000 atoms. Chloride counter-ions were included to neutralize the overall electric charge of the system. We used the CHARMM22* force field (Piana et al., 2011) for most of our simulations, but also used the Amber ff99SB*-ILDN (Hornak et al., 2006; Best and Hummer, 2009; Lindorff-Larsen et al., 2010) for some simulations to examine the dependency of the results on the choice of force field.

The van der Waals interactions were smoothly shifted to zero between 0.8 and 1.0 nm, and the long-range electrostatic interactions were calculated by means of the particle mesh Ewald (PME) algorithm with a 0.12 nm mesh spacing combined with a switch function for the direct space between 0.8 and 1.0 nm. The bonds involving hydrogen atoms were constrained using the LINCS algorithm. We employed the V-rescale thermostat (Bussi et al., 2007) to control the temperatureand simulated the system in the canonical ensemble.

Reference transition path and path collective variables

We used path collective variables both to enhance sampling in path driven metadynamics (see below) as well as to represent the conformational landscape sampled by other means. Path collective variables have previously been shown to be very useful in finding free energy channels connecting two metastable states, and also able to construct the global free energy surfaces even far away from the initial path (Branduardi et al., 2007). A reference path is defined by a set of conformations along the path, and the progress along this path can be described mathematically as:

S_{p a t h} (X) = \frac{\sum_{i = 1}^{N} i e^{- λ M_{i} (X)}}{\sum_{i = 1}^{N} e^{- λ M_{i} (X)}}

Here X are the coordinates of the instantaneous protein conformer sampled by MD simulations, N is the number of frames used to describe the reference path (often dependent on the length scale of the conformational transition process), $M_{i} (X)$ is the mean-square deviation (after optimal alignment) of a subset of the atoms from the reference structure of $i$ th frame, and $λ$ is a smoothing parameter whose value should roughly be proportional to the inverse of the average mean square displacement between two successive frames along the reference path. With this definition, $S_{p a t h}$ quantifies how far the instantaneous conformer, $X$ , is from the reactant state and the product state, thus monitoring the progress of the system along the conformational transition channel.

Using $S_{p a t h}$ as the sole CV would assume that the initial reference path contains a sufficient description of the important degrees of freedom between the two states. It is, however, rarely possible to guess such a path because determining the actual pathway taken is a goal of the simulation. Thus, $S_{p a t h}$ is supplemented by a second CV, $Z_{p a t h}$ , which measures the deviation away from the structures on the reference path. I.e. if $S_{p a t h}$ quantifies the progress along the path, $Z_{p a t h}$ measures the distance away from the reference path:

Z_{p a t h} (X) = - \frac{1}{λ} l n \sum_{i = 1}^{N} e^{- λ M_{i} (X)}

The combination of $S_{p a t h}$ and $Z_{p a t h}$ thus maps the entire conformational landscape to a two-dimensional projection, which can also be thought of as a tube connecting the two end states, and where $S$ measures the progress along the tube and $Z$ the width of the tube. The usefulness of the path CVs is, however, dependent on the quality of reference path, which is determined amongst other things by two factors: (1) the relative accuracy of the reference path and (2) how uniform reference structures are distributed along the path. Because of the explosion in number of possible conformations as one progresses along $Z_{p a t h}$ , simulations are mostly enhanced when $S_{p a t h}$ provides a relatively good description of the pathway taken. Further, if reference structures are only placed sparsely along the path, one looses resolution of the free energy surface and also decreases the ability to enhance sampling.

To obtain a good reference path, without knowing beforehand the mechanism of conversion, we here used a method to construct snapshots along a possible initial path. Taking $G_{X r a y}$ and $E_{R O S E T T A}$ as initial and final structures, we calculated the optimal reaction paths between them with the MOIL software (Elber et al., 1995) which has previously been used to explore the mechanism of conformational change of proteins (Wang et al., 2011) . After minimizing endpoint structures, we employed the minimum-energy-path self-penalty walk (SPW) (Czerminski and Elber, 1990) functional embedded in the CHMIN module to obtain an initial guess for the conformational transition path. This path was subsequently optimized in the SDP (steepest descent path) module by minimizing the target function $T$ consisting of two terms $S$ and $C$ . $S$ is an action function that provides approximate most probable Brownian trajectories connecting the reactant and product states, while $C$ is a restraint function aimed to distribute frames approximately uniformly along the path. They can be expressed by:

\begin{aligned} S = \sum_{i = 1}^{N - 1} \sqrt{H_{s} + (\frac{\partial U}{\partial x_{i}})^{2}} | x_{i + 1} - x_{i} | \\ C = λ \sum_{i} (Δ l_{i, i + 1} - < Δ l >)^{2} \end{aligned}

where $N$ is the number of frames along the reference path, $x_{i}$ is the entire vector of conformational coordinates of frame $i$ , $U$ is a potential energy as a function of the mass-weighted coordinate vector, $H_{s}$ is a constant with an arbitrary positive value, which can be tuned to generate the optimal paths with different thermal energies. $Δ l_{i, i + 1} = M^{1 / 2} | x_{i} - x_{i + 1} |$ is the arc-length between consecutive frames. $λ$ controls the strength of the restraint function $C$ .

The SDP or minimum energy path is the limiting path in which the action $S$ is optimized with $H_{s} \to 0$ . An important advantage of an SDP is that it is capable of giving a good guess of the minimal energy path which can reflect the major mechanism, only with inexpensive computation. Further details can be found in Refs. (Wang et al., 2011; Majek et al., 2008).

The SDP was approximated as 31 discrete conformations. Most regions of the $G_{X r a y}$ and $E_{R O S E T T A}$ structures are very similar, with the exception of the cavity-related atoms whose movement determines the conformational transitions between the G and E states. To minimize computations, we used only a subset of heavy atoms around the cavity from amino acid residues 99 to 126 to define $S_{p a t h}$ and $Z_{p a t h}$ , resulting a ‘light version’ of the reference path which included only 212 atoms. We used the $C_{α}$ atoms of the C-terminal domain to align the molecule. It is important to note that by focusing only on atoms surrounding the cavity in the calculation of the path-variables we only enhance sampling of the conformational changes relating to the cavity. We used $λ = 56.0$ based on the consideration that smooth change of the function of $S_{p a t h}$ can be achieved when $e^{- λ < M_{i} (X) >} \geq 0.1$ (Tribello et al., 2014).

We defined $S_{p a t h}$ and $Z_{p a t h}$ using the SDP as described above. The equidistant requirement of the path is satisfied by the penalty function as is evident from the RMSD matrix for path frames which has a gullwing shape (Figure 2—figure supplement 1, indicating that each frame is closest to its neighbor and more different from all other reference frames.

Metadynamics simulations

Metadynamics discourages the system from sampling already visited conformational regions by continuously adding an external history-dependent repulsive potential at the present value of the reaction coordinates or CVs, which are assumed to include the slowly varying degrees of freedom and thus describe the main features of the dynamics (Laio and Parrinello, 2002). The biasing potential in metadynamics results in an artificial (enhanced) dynamics but makes it possible to reconstruct the free energy surfaces by removing the bias introduced. The bias is typically added as Gaussians at regular time intervals, $τ_{G}$ , and is given by:

V_{G} (S, t) = \sum_{k = 1}^{| t / τ_{G} |} ω_{k} δ (S (t), S (k \cdot τ_{G}))

Here $S$ denotes the CVs, $k$ is the index of the individual Gaussians, $ω_{k}$ is the height of the $k$ ’th Gaussian, and $δ (S (t), S (k \cdot τ_{G}))$ is a short-ranged kernel function of the CVs:

δ (S (t), S (k \cdot τ_{G})) = e^{\sum_{j = 1}^{n} - \frac{{| s_{j} (t) - s_{j} (k \cdot τ_{G}) |}^{2}}{2 σ_{k j}^{2}}}

where $n$ is the number of CV, $j$ is the index of a CV, and $σ_{k j}$ and $s_{j} (k \cdot τ_{G})$ are the width and the position of the Gaussian hills, respectively.

We here used a range of different metadynamics approaches to determine the free energy landscape, and the mechanism and kinetics of conformational exchange (see below).

Well-tempered metadynamics

In the well-tempered version of metadynamics, the height of the individual Gaussians, $ω_{k}$ , is decreased as the total bias accumulates over time, in order to improve the convergence of the free energy:

ω_{k} = ω_{0} e^{- \frac{1}{γ - 1} \frac{V_{G} (S, k \cdot τ_{G})}{k_{B} T}}

Here $ω_{0}$ is the initial height, $γ = (T + Δ T) / T$ is referred as the bias factor, which can be tuned to control the speed of convergence and diminish the time spent in lesser-relevant, high-energy states. Thus, the quantity $T + Δ T$ is often referred as the fictitious CV temperature.

Adaptive-width metadynamics

In contrast to standard metadynamics in which the width of the Gaussians is constant, the adaptive-width version of metadynamics updates the Gaussian width $σ_{k j}$ on the fly according to local properties of the underlying free-energy surface on the basis of local diffusivity of the CVs or the local geometrical properties (Branduardi et al., 2012).

In the region of conformational space near the endpoints of the path CVs many conformations are compressed on similar CV values, leading to high-density but low-fluctuation boundaries. It is apparent that the use of a fixed width might give an inaccurate estimation of the free energy profile in the boundaries where the free energy basins of reactant and product states are located, also makes it more difficult for the simulations to converge. Therefore, the feature of shape adaptive of Gaussian potential is particularly helpful for the case of using path CVs that have significant boundary effects.

Metadynamics with path variables (PathMetaD)

We sampled the free energy landscape along $S_{p a t h}$ and $Z_{p a t h}$ , as defined above, using adaptive-width metadynamics, which resulted in a finer resolution and faster convergence of the free energy landscape, in particular near the path boundaries, than standard metadynamics. The production simulations were performed at 298K in well-tempered ensemble (Appendix 1—table 1: RUN1 and RUN2).

Reconnaissance metadynamics

To explore the mechanism of conformational exchange we used reconnaissance metadynamics (Tribello et al., 2010). This is a ‘self-learning’ approach which combines of a machine learning technique to automatically identify the locations of free energy minima by periodically clustering the trajectory and a dimensional reduction technique that can reduce the complex locations to a locally-defined one-dimensional CV by using information collected during the clustering. It has previously been shown that reconnaissance metadynamics makes it possible to determine a path from a large set of input collective variables (Söderhjelm et al., 2012; Tribello et al., 2011).

Infrequent Metadynamics (InMetaD)

We used the recently described ‘infrequent metadynamics’ (InMetaD) to obtain the rates of the conformational exchange process (Tiwary and Parrinello, 2013). In standard metadynamics simulations it is very difficult to obtain kinetic properties because the biasing potential is added both to the free energy basins as well as the barriers that separate them. While it is potentially possible to determine the rates from the height of the free energy barriers, this requires both that the CVs used represent the entire set of slowly varying degrees of freedom, and also a good estimate of the pre-exponential factor to convert barrier height to a rate.

The key idea in InMetaD to circumvent these problems is to attempt to add the bias to the system more slowly than the barrier crossing time but faster than the slow inter-basin relaxation time, so that the transition state region has a lower risk of being biased, and therefore the transitions are less affected. By filling up a free energy basin by a known amount it is possible to determine how much the barrier has been decreased, and hence remove this bias from the rates determined. Thus, as described in more detail below, the approach works by performing a large number of individual simulations to obtain first passage times between the individual basins, which are then corrected by the known enhancement factors to obtain estimates of the unbiased rates. This method has been successfully used to reproduce the kinetics of conformational change of alanine dipeptide (Tiwary and Parrinello, 2013) unbinding of the inhibitor benzamidine from trypsin (Tiwary et al., 2015a) and slow unbinding of a simple hydrophobic cavity-ligand model system (Tiwary et al., 2015b).

In these simulations we used a deposition frequency of 80 or 100 ps (see parameters in Appendix 1—table 1), a value much lower than the deposition frequency of 1 ps used in the PathMetaD simulations described above. In this way we lower the risk of substantially corrupting the transition state region. In addition, a tight upper wall potential on $Z_{p a t h}$ =0.10 nm² is used to confine the sampling based on our converged free energy surface which shows the conformational change mostly occurs within this region.

With these parameters we collected dozens of trajectories that have a state-to-state transition in the G-to-E and E-to-G directions. The passage times observed in each of these were then corrected for the metadynamics bias as follows.

First, we calculate the acceleration factor $α$ from:

α = τ / τ_{M} =< e^{V (s, t) / k T} >_{M}

where the anguluar brackets denote an average over a metadynamics run before the first transition, and $V (s, t)$ is the metadynamics time-depedent bias. The evolution of the acceleration factor $α (t)$ can be expressed by:

α (t) = (1 / t) \int_{0}^{t} d t^{'} e^{V (s, t^{'}) / k T}

Then the observed passage time, $t$ , is reweighted by:

τ_{t r u e} = α (t) * t = \int_{0}^{t} d t^{'} e^{V (s, t^{'}) / k T}

In principle, the transition time should be a Poisson-distributed random variable, and its mean, $μ$ , standard deviation $σ$ and median $t_{m} / l n 2$ all should be equal to each other. In practice, however, they are somewhat sensitive to insufficient sampling (Salvalaglio et al., 2014). So rather than simply calculating averages of the individual times, we estimated the average rate and transition time $τ$ from a fit of the empirical cumulative distribution function (ECDF) with the theoretical cumulative distribution function (TCDF):

T C D F = 1 - e^{- \frac{t}{τ}}

It has previously been shown that $τ$ estimated in this way converges more quickly than the simple average, $μ$ . This is also consistent with our observation, and we find that 10–15 samples appear sufficient to get a reasonably accurate estimation of the transition time. We used a bootstrap approach to estimate the errors.

To examine whether the observed times indeed follow the expected Poisson distribution we used a Kolmogorov-Smirnov (KS) test to obtain a p-value that quantifies the similarity between the empirical and theoretical distributions. Traditionally, a threshold value typically of 0.05 or 0.01 (the significance level of the test) is used to judge if the theoretical (TCDF) and empirical (ECDF) distributions are in agreement. If the p-value is equal to or larger than the threshold value, it suggests that the estimated transition time is quite reliable. If (a) the transition regions were perturbed significantly with infrequent biasing or (b) there are hidden unidentified timescales at play (e.g. the CVs do not capture the slow degrees of freedom) the KS test for time-homogeneous Poisson statistics would fail.

CS-restrained replica-averaged simulation

The simulation methods described above constitute different ways of exploring the thermodynamics (PathMetaD), kinetics (InMetaD) and mechanism (Reconnaissance metadynamics) of conformational exchange. In all of these simulations, sampling is determined by the molecular energy function (force field), and the experimental information on T4L is used only in the construction of the path variables. When additional experimental information is available, one may introduce an additional energy term so as to bias the simulations to be in agreement with this information (Robustelli et al., 2010). As the experimental values are ensemble averages we apply these restraints only to averages calculated over a number of ‘replicas’ that are simulated in parallel. In this way, the information from the experimental data is incorporated into the simulation as a perturbation following the maximum entropy principle (Pitera and Chodera, 2012; Roux and Weare, 2013; Cavalli et al., 2013; Boomsma et al., 2014).

We used this approach to obtain conformational ensembles that include not only information from the molecular force field, but also experimental NMR chemical shifts (CS). In particular, we used either the chemical shifts of the E state of L99A obtained from the analysis of the CPMG experiments, or the native state chemical shifts of a triple mutant that populates the same state as its ground state (Bouvignies et al., 2011). The CS restraints were imposed by adding a pseudo-energy term ( $E_{C S}$ ) to a standard molecular-mechanics force field ( $E_{F F}$ ).

E_{C S} = ϵ_{C S} \sum_{i = 1}^{N} \sum_{j = 1}^{6} {(δ_{i j}^{E x p} - \frac{1}{M} \sum_{k = 1}^{M} δ_{i j k}^{S i m})}^{2}

Here $ϵ_{C S}$ is strength of the CS restraints, $i$ indicates the residue number (total of $N$ ), $j$ indicates each of the six backbone atoms whose chemical shifts were used ( $C_{α}$ , $C_{β}$ , $C^{'}$ , $H_{α}$ , HN and N), $k$ is an index for the total of $M$ replicas, and $δ^{E x p}$ and $δ^{S i m}$ are the experimental and simulated CSs, respectively. The latter quantity, $δ^{S i m}$ , was calculated by CamShift (Fu et al., 2014) as a plugin of PLUMED. The CS values of Pro, Gly, Asp, Glu, His and terminal residues were not included because the accuracy in their predictions are too low to contain sufficient information in this approach (Fu et al., 2014). We set $ϵ_{C S}$ =24 kJ mol⁻¹ ppm⁻² and used M=4 replicas. In principle, the number of replicas is a free parameter that should be set as large as possible when the experimental data and the method for calculating it is noise free (Boomsma et al., 2014). In practice, one uses a finite set of replicas and it has been shown that M=4 replicas is sufficient to capture the dynamics accurately (Camilloni et al., 2013).

We performed replica-averaged CS restrained MD simulations using GROMACS4.6 and the PLUMED2.1 plugin at 298 K. The equilibrated structures of $E_{R O S E T T A}$ and $G_{R O S E T T A}^{T r i p l e}$ were used as the starting configuration (of each of the four replicas) in the CS-restrained simulations of L99A and the L99A,G113A,R119P triple mutant, respectively. The CS data for $E_{R O S E T T A}$ and $G_{R O S E T T A}^{T r i p l e}$ were obtained from the Biological Magnetic Resonance Bank (BMRB) database (Ulrich et al., 2008) with entries 17604 and 17603, respectively.

PT-WT-MetaD failed to get the converged free energy landscape

In practice, the choice of CVs plays a fundamental role in determining the accuracy, convergence and efficiency of metadynamics simulations. If an important CV is missing, the exploration of the free energy landscape will be difficult due to hysteresis. Finding a minimal set of CVs that include all important degrees of freedom is a highly nontrivial task and one often has to proceed by several rounds of trial simulations.

At the beginning, we followed a strategy which had previously been successfully used in the exploration of protein conformational transitions (Sutto and Gervasio, 2013; Papaleo et al., 2014) to design a set of CVs on the basis of static structural comparison between $G_{X r a y}$ and $E_{R O S E T T A}$ . In particular, by comparing these two structures we defined several CVs that described structural differences by individual dihedral angles and hydrogen bonds, as well as dihedral correlation and coordination number (state-specific contact map) (Bonomi et al., 2009) (summarized in Appendix 1—table 1). We combined these in a multiple-replica, parallel tempering approach in the well-tempered ensemble (PT-WT-metaD) (Bonomi and Parrinello, 2010) to further enhance the sampling. In PT-WT-metaD, the energy fluctuations are enlarged by using energy as a biased CV but the average energy is the same as the canonical ensemble, allowing the use of a larger spacing between temperatures and a much fewer number of replicas than normal PT simulations (Barducci et al., 2015). Coordinate exchange with high temperature replicas can enhance the sampling of all the degrees of freedom, even those not included in the biased CVs, and one may include a ‘neutral’ replica (without energy bias, at 298 K). We performed a series of simulations with different combinations of CVs starting from the G state of L99A (Appendix 1—table 2). However, unfortunately, we only observed partial G-to-E transitions, even in a relatively long trajectory of about 1 μs for each replica. This negative result suggested that these manually chosen CVs did not contain all the necessary slow degrees of freedom.

Tunnel analysis

We used CAVER3 (Chovancova et al., 2012) to analyse the structures and CAVER Analyst1.0 (Kozlikova et al., 2014)(http://www.caver.cz/, see also parameters in Appendix 1—table 4) to separate the tunnels into different clusters. CAVER Analyst is a standalone program based on CAVER3.0 algorithm (Chovancova et al., 2012). The settings for the tunnel calculations can be set through the Tunnel Computation window, while the advanced parameters can be set in the Tunnel Advanced Settings window. We used the center of mass of the cavity-related region (residue 93-124) as the position of the starting point. Average-link hierarchical clustering algorithm is performed to build a tree hierarchy of tunnel axes based on their pairwise distances. The size of the resulting clusters is dependent on the Clustering threshold parameter which specifies the level of detail at which the tree hierarchy of tunnel clusters will be cut. We used the default value so that the tree hierarchy of tunnel clusters is cut at the value of 3.5.

Adiabatic bias molecular dynamics

Adiabatic biased molecular dynamics (ABMD) (Marchi and Ballone, 1999; Paci and Karplus, 1999) is an algorithm developed to accelerate the transition from the reactant state to the productive state, here corresponding to the ligand bound state and ligand-free state, respectively. In ABMD the system is perturbed by a ‘ratcheting potential’, which acts to ‘select’ spontaneous fluctuations towards the ligand-free state. The ratcheting potential is implemented in PLUMED2.2 as

V (ρ (t)) = {\begin{cases} 0.5 K (ρ (t) - ρ_{m} (t))^{2}, & ρ (t) > ρ_{m} (t) \\ 0, & ρ (t) \leq ρ_{m} (t) \end{cases}

where

ρ (t) = {(S (t) - S_{t a r g e t})}^{2}

and

ρ_{m} (t) = m i n_{0 \leq τ \leq t} ρ (τ) + η (t)

$K$ is the force constant, $S (t)$ is the instantaneous CV value, $S_{t a r g e t}$ is the target value of the CV and $η (t)$ is an additional white noise acting on the minimum position of $ρ (t)$ . Here, we used the RMSD of the ligand to the cavity-bound state as the CV, and set $S_{t a r g e t} = 4.0$ nm and K=20 kJ · mol⁻¹ · nm⁻² or 50 kJ · mol⁻¹ · nm⁻² to check the dependency of the force constant chosen. The biasing potential is zero when the CV increases but provides a penalty when the CV decreases. In this way, we were able to observe multiple unbinding events in simulations despite the long lifetime (1.2 ms) of the ligand in the cavity.

Appendix 1—figure 1

Download asset Open asset

Two representative InMetaD trajectories of L99A with G to E transitions.

The time point, t′, for the first transition from G to E is identified when the system evolves into conformational region of S_path > 0.55 and Z_path < 0.01. We then calculate the unbiased passage time by multiplying t′ by the corresponding accelerate factor α(t′). Upper panels show the evolution of the reweighted time as a function of metadynamics time. The kinks usually indicate a possible barrier crossing event. Middle panels show the trajectories starting from the G state and crossing the barrier towards the E state. Lower panels show the biasing landscape reconstructed from the deposited Gaussian potential, which can be used to check the extent to which the transition state regions are affected by deposited bias potential.

https://doi.org/10.7554/eLife.17505.030

Appendix 1—figure 2

Download asset Open asset

Two representive InMetaD trajectories of L99A with E to G transitions of L99A.

First transition time for G to E transition is identified when the system evolves into conformational region of S_path < 0.28 and Z_path < −0.01.

https://doi.org/10.7554/eLife.17505.031

Appendix 1—figure 3

Download asset Open asset

Characteristic transition times between G and E states of L99A.

The error bars represent the standard deviation of τ obtained from a bootstrap analysis.

https://doi.org/10.7554/eLife.17505.032

Appendix 1—figure 4

Download asset Open asset

Characteristic transition times between G and E states of the triple mutant.

The figure shows the characteristic transition time τ_G→E (left panel) and τ_E→G (right panel) of the triple mutant as a function of the size of a subsample of transition times randomly extracted from the main complete sample. The error bars represent the standard deviation of characteristic transition times obtained by a bootstrap analysis. The calculated and experimental values of the transition times are shown in blue and red font, respectively.

https://doi.org/10.7554/eLife.17505.033

Appendix 1—figure 5

Download asset Open asset

Poisson fit analysis for G to E transitions and E to G transitions of L99A.

We show the ECDF (the empirical cumulative distribution function) and TCDF (the theoretical cumulative distribution function) in black and blue lines, respectively. The respective p-values are reasonably, albeit not perfectly, well above the statistical threshold of 0.05 or 0.01, indicating the kinetics is not substantially modified by the deposited bias potential in InMetaD. Error bars are the standard deviation obtained by a bootstrap analysis.

https://doi.org/10.7554/eLife.17505.034

Appendix 1—figure 6

Download asset Open asset

Poisson fit analysis for G to E transitions and E to G transitions of the triple mutant.

The figure shows the p-values of the Poisson fit analysis of G → E (A) and E → G (B) transition times as a function of the size of a subsample of transition times randomly extracted from the main complete sample.

https://doi.org/10.7554/eLife.17505.035

Appendix 1—table 1

Simulation details.

https://doi.org/10.7554/eLife.17505.036

Method	label	system	init	length	force field	CVs	parameters	T (K)	software version
PathMetaD	RUN1	L99A	G	667 ns	CHARMM22*	$S_{p a t h}$ , $Z_{p a t h}$	$γ$ =20^†, $τ_{D}$ =1 ps^‡, $τ_{G}$ =1 ps^§, $ω_{0}$ =1.5^¶	298	PLUMED2.1, GMX4.6
PathMetaD	RUN2	L99A, G113A, R119P	G	961 ns	CHARMM22*	$S_{p a t h}$ , $Z_{p a t h}$	$γ$ =20, $τ_{D}$ =1 ps, $τ_{G}$ =1ps, $ω_{0}$ =1.5	298	PLUMED2.1, GMX5.0
CS-restrained MD	RUN3	L99A	E	252 ns	CHARMM22*		N^#=4, $ϵ_{C S} = 24$ **	298	PLUMED2.1, ALMOST2.1, GMX5
CS-restrained MD	RUN4	L99A	E	233 ns	CHARMM22*		N=2, $ϵ_{C S} = 24$	298	PLUMED2.1, ALMOST2.1, GMX5
CS-restrained MD	RUN5	L99A	E	221 ns	CHARMM22*		N=2, $ϵ_{C S} = 12$	298	PLUMED2.1, ALMOST2.1, GMX5
CS-restrained MD	RUN6	L99A	E	125 ns	AMBER99sb*ILDN		N=4, $ϵ_{C S} = 24$	298	PLUMED2.1, ALMOST2.1, GMX5
CS-restrained MD	RUN7	L99A	E	190 ns	AMBER99sb*ILDN		N=2, $ϵ_{C S} = 12$	298	PLUMED2.1, ALMOST2.1, GMX5
CS-restrained MD	RUN8	L99A, G113A, R119P	G	204 ns	CHARMM22*		N=4, $ϵ_{C S} = 24$	298	PLUMED2.1, ALMOST2.1, GMX5
CS-restrained MD	RUN9	L99A, G113A, R119P	G	200 ns	CHARMM22*		N=2, $ϵ_{C S} = 24$	298	PLUMED2.1, ALMOST2.1, GMX5
Reconnaissance MetaD	RUN10	L99A	G	120 ns	CHARMM22*	2,3,4,5	$Δ T$ =10000^††, $Δ t$ =250	298	PLUMED1.3, GMX4.5
Reconnaissance MetaD	RUN11	L99A	G	85 ns	CHARMM22*	2,3,4,5	10000, 250	298	PLUMED1.3, GMX4.5
Reconnaissance MetaD	RUN12	L99A	G	41 ns	CHARMM22*	2,3,4,5,7,8	100000, 500	298	PLUMED1.3, GMX4.5
PT-WT-MetaD	RUN13	L99A	G	961 ns	CHARMM22*	1,2,3,7,8	$γ$ =20, $ω_{0}$ =0.5	297 298^‡‡ 303 308 316 325 341 350 362	PLUMED 1.3, GMX4.5
PT-WT-MetaD	RUN14	L99A	G	404 ns	CHARMM22*	1,2,3,6,7	$γ$ =20, $ω_{0}$ =0.5	297 298 303 308 316 325 333 341 350	PLUMED 1.3, GMX4.5
PT-WT-MetaD	RUN15	L99A	G	201 ns	CHARMM22*	1,2,3,6	$γ$ =20, $ω_{0}$ =0.5	297 298 303 308 316 325 333 341 350	PLUMED 1.3, GMX4.5
PT-WT-MetaD	RUN16	L99A	G	511 ns	CHARMM22*	1,2,3,4,5	$γ$ =5, $ω_{0}$ =0.5	295 298 310 325 341 358 376	PLUMED 1.3, GMX4.5
PT-WT-MetaD	RUN17	L99A	G	383 ns	CHARMM22*	1,2,3,4,5	$γ$ =20, $ω_{0}$ =0.5	298 305 313 322 332 337 343	PLUMED 1.3, GMX4.5
PT-WT-MetaD	RUN18	L99A	E	365 ns	CHARMM22*	1	$γ$ =20, $ω_{0}$ =0.5	298 302 306 311 317 323 330 338	PLUMED 1.3, GMX4.5
plain MD	RUN19	L99A	G	400 ns	CHARMM22*			298	GMX5.0
plain MD	RUN20	L99A	E	400 ns	CHARMM22*			298	GMX5.0
plain MD	RUN21	L99A	G	400 ns	AMBER99sb*ILDN			298	GMX5.0
plain MD	RUN22	L99A	E	400 ns	AMBER99sb*ILDN			298	GMX5.0
plain MD	RUN23	L99A,G113A,R119P	G	400 ns	CHARMM22*			298	GMX5.0
plain MD	RUN24	L99A,G113A,R119P	E	400 ns	CHARMM22*			298	GMX5.0
plain MD	RUN25	L99A,G113A,R119P	G	400 ns	AMBER99sb*ILDN			298	GMX5.0
plain MD	RUN26	L99A,G113A,R119P	E	400 ns	AMBER99sb*ILDN			298	GMX5.0
InMetaD	RUN27-68^§§	L99A	G		CHARMM22*	$S_{p a t h}$ , $Z_{p a t h}$	$γ$ =20, $τ_{G}$ =80 ps, $ω_{0}$ =1.0, $σ_{S}$ =0.016^¶¶, $σ_{Z}$ =0.002^##	298	PLUMED2.1, GMX5
InMetaD	RUN69-104***	L99A	E		CHARMM22*	$S_{p a t h}$ , $Z_{p a t h}$	$γ$ =20, $τ_{G}$ =80 ps, $ω_{0}$ =1.0, $σ_{S}$ =0.016, $σ_{Z}$ =0.002	298	PLUMED2.1, GMX5
InMetaD	RUN105-119	L99A,G113A,R119P	G		CHARMM22*	$S_{p a t h}$ , $Z_{p a t h}$	$γ$ =15, $τ_{G}$ =100ps, $ω_{0}$ =0.8, $σ_{S}$ =0.016^¶¶, $σ_{Z}$ =0.002^##	298	PLUMED2.2, GMX5.1.2
InMetaD	RUN120-134	L99A,G113A,R119P	E		CHARMM22*	$S_{p a t h}$ , $Z_{p a t h}$	$γ$ =15, $τ_{G}$ =100 ps, $ω_{0}$ =0.8, $σ_{S}$ =0.016^¶¶, $σ_{Z}$ =0.002^##	298	PLUMED2.2, GMX5.1.2
ABMD	RUN135-154	L99A-Benzene	bound		CHARMM22*	$R M S D_{b e n z e n e}$	K=50 KJ · mol⁻¹ · nm⁻², $S_{t a r g e t} = 4.0$ nm	298	PLUMED2.2, GMX5.1.2
ABMD	RUN155-174	L99A-Benzene	bound		CHARMM22*	$R M S D_{b e n z e n e}$	K=20 KJ · mol⁻¹ · nm⁻², $S_{t a r g e t} = 4.0$ nm	298	PLUMED2.2, GMX5.1.2

^† $γ$ is the bias factor.
^‡ $τ_{D}$ is the characteristic decay time used for the dynamically-adapted Gaussian potential.
^§ $τ_{G}$ is the deposition frequency of Gaussian potential.
^¶ $ω_{0}$ is the height of the deposited Gaussian potential (in KJ · mol⁻¹).
^# N is the number of replicas.
** $ϵ_{C S}$ is the strength of chemical shift restraints (in KJ · mol⁻¹ · ppm⁻²).
^†† $Δ T$ means RUN_FREQ and $Δ t$ means STORE_FREQ. Other parameters: BASIN_TOL=0.2, EXPAND_PARAM=0.3, INITIAL_SIZE=3.0.
^‡‡ Replica at 298 K is the neutral replica without energy bias.
^§§42 trajectories collected for G-to-E transitions.
^¶¶ $σ_{S}$ is the Gaussian width of $S_{p a t h}$ .
^## $σ_{Z}$ is the Gaussian width of $Z_{p a t h}$ (in nm²).
***36 trajectories collected for E-to-G transitions.

Appendix 1—table 2

Definition of collective variables.

https://doi.org/10.7554/eLife.17505.037

CV	definitions	parameters	purpose
1	total energy	bin=500	enhance energy fluctuations
2	dihedral angle of $C_{α}$ atoms of consecutive residues F104-Q105-M106-G107	$σ$ =0.1
3	dihedral angle of $C_{α}$ atoms of consecutive residues G113-F114-T115-N116	$σ$ =0.1
4	$Q_{G}$ , distance in contact map space to the $G_{X r a y}$ structure	$σ$ =0.5
5	$Q_{E}$ , distance in contact map space to the $E_{R O S E T T A}$ structure	$σ$ =0.5
6	distance between $Q_{G}$ and $Q_{E}$	$σ$ =0.5
7	number of backbone hydrogen bonds formed between M102 and G107	$σ$ =0.1
8	dihedral correlation between the $C_{α}$ dihedral angles of consecutive residues in segment N101-G107	$σ$ =0.1
9	global RMSD to the whole protein	wall potential	avoid sampling unfolding space

Appendix 1—table 3

Average root-mean-square deviation ( $< R M S D >$ in units of ppm) between experimental CSs and those from the CS-restrained replica-averaged simulations.

https://doi.org/10.7554/eLife.17505.038

Nucleus	RUN3	RUN4	RUN5	RUN6	RUN7	RUN8	RUN9
$C^{'}$	0.833	0.655	0.776	0.854	0.793	0.907	0.727
$C_{α}$	1.055	0.879	0.929	1.065	0.940	1.103	0.894
$N$	1.966	1.707	1.771	1.967	1.780	2.011	1.828
$H_{N}$	0.379	0.275	0.291	0.368	0.284	0.414	0.286
$H_{A}$	0.232	0.183	0.186	0.242	0.182	0.246	0.183

Appendix 1—table 4

Parameter set used in tunnel analysis using CAVER3.0 (Chovancova et al., 2012) and CAVER Analyst1.0 (Kozlikova et al., 2014).

https://doi.org/10.7554/eLife.17505.039

Minimum probe radius	0.9 Å
Shell depth	4
Shell radius	3
Clustering threshold	3.5
Starting point optimization
Maximum distance	3 Å
Desired radius	5 Å

Appendix 1—table 5

Clustering analysis of tunnels (top three listed).

https://doi.org/10.7554/eLife.17505.040

Index	Population	Maximal bottleneck radius (Å)	Average bottleneck radius (Å)
#1	27%	2.5	1.3
#2	20%	1.4	1.0
#3	15%	1.3	1.0

Appendix 1—table 6

Unbinding Pathways Explored by ABMD ( $R M S D_{B N Z}$ as CV).

https://doi.org/10.7554/eLife.17505.041

	k=20 kJ/(mol · nm⁻²)		k=50 kJ/(mol · nm⁻²)
Index	Length	Path	Length	Path
RUN1	56 ns	P1	27 ns	P2
RUN2	36 ns	P2	78 ns	P1
RUN3	43 ns	P1	6 ns	P1
RUN4	43 ns	P1	35 ns	P1
RUN5	77 ns	P2	10 ns	P1
RUN6	176 ns	P1	44 ns	P1
RUN7	41 ns	P1	18 ns	P1
RUN8	106 ns	P1	15 ns	P1
RUN9	72 ns	P1	7 ns	P1
RUN10	107 ns	P1	2 ns	P1
RUN11	61 ns	P1	20 ns	P2
RUN12	58 ns	P2	26 ns	P1
RUN13	64 ns	P1	31 ns	P2
RUN14	173 ns	P2	20 ns	P1
RUN15	172 ns	P1	34 ns	P1
RUN16	74 ns	P2	22 ns	P1
RUN17	20 ns	P1	17 ns	P1
RUN18	34 ns	P1	35 ns	P2
RUN19	91 ns	P1	21 ns	P2
RUN20	61 ns	P1	18 ns	P1
Cost	1.6 μs		0.5 μs
Summary
P1	75% (15/20)		75% (15/20)
P2	25% (5/20)		25% (5/20)

References

1. Baldwin AJ
2. Kay LE
(2009) NMR spectroscopy brings invisible protein states into focus
Nature Chemical Biology 5:808–814.
https://doi.org/10.1038/nchembio.238
- Google Scholar
(2008) Well-tempered metadynamics: a smoothly converging and tunable free-energy method
Physical Review Letters 100:020603.
https://doi.org/10.1103/PhysRevLett.100.020603
- Google Scholar
(2015)
Molecular Modeling of Proteins

151–171, Tackling sampling challenges in biomolecular simulations, Molecular Modeling of Proteins, Springer.
- Google Scholar
1. Best RB
2. Hummer G
(2009) Optimized molecular dynamics force fields applied to the helix−Coil transition of polypeptides
The Journal of Physical Chemistry B 113:9004–9015.
https://doi.org/10.1021/jp901540t
- Google Scholar
1. Bonomi M
2. Branduardi D
3. Bussi G
4. Camilloni C
5. Provasi D
6. Raiteri P
7. Donadio D
8. Marinelli F
9. Pietrucci F
10. Broglia RA
11. Parrinello M
(2009) PLUMED: A portable plugin for free-energy calculations with molecular dynamics
Computer Physics Communications 180:1961–1972.
https://doi.org/10.1016/j.cpc.2009.05.011
- Google Scholar
1. Bonomi M
2. Parrinello M
(2010) Enhanced sampling in the well-tempered ensemble
Physical Review Letters 104:190601.
https://doi.org/10.1103/PhysRevLett.104.190601
- Google Scholar
(2014) Combining experiments and simulations using the maximum entropy principle
PLoS Computational Biology 10:e1003406.
https://doi.org/10.1371/journal.pcbi.1003406
- Google Scholar
1. Bouvignies G
2. Vallurupalli P
3. Hansen DF
4. Correia BE
5. Lange O
6. Bah A
7. Vernon RM
8. Dahlquist FW
9. Baker D
10. Kay LE
(2011) Solution structure of a minor and transiently formed state of a T4 lysozyme mutant
Nature 477:111–114.
https://doi.org/10.1038/nature10349
- Google Scholar
(2012) Metadynamics with adaptive gaussians
Journal of Chemical Theory and Computation 8:2247–2254.
https://doi.org/10.1021/ct3002464
- Google Scholar
(2007) From A to B in free energy space
The Journal of Chemical Physics 126:054103.
https://doi.org/10.1063/1.2432340
- Google Scholar
(2007) Canonical sampling through velocity rescaling
The Journal of Chemical Physics 126:014101.
https://doi.org/10.1063/1.2408420
- Google Scholar
(2013) Replica-Averaged Metadynamics
Journal of Chemical Theory and Computation 9:5610–5617.
https://doi.org/10.1021/ct4006272
- Google Scholar
(2013) Molecular dynamics simulations with replica-averaged structural restraints generate structural ensembles according to the maximum entropy principle
The Journal of Chemical Physics 138:094112.
https://doi.org/10.1063/1.4793625
- Google Scholar
1. Chovancova E
2. Pavelka A
3. Benes P
4. Strnad O
5. Brezovsky J
6. Kozlikova B
7. Gora A
8. Sustr V
9. Klvana M
10. Medek P
11. Biedermannova L
12. Sochor J
13. Damborsky J
(2012) CAVER 3.0: a tool for the analysis of transport pathways in dynamic protein structures
PLoS Computational Biology 8:e1002708.
https://doi.org/10.1371/journal.pcbi.1002708
- Google Scholar
1. Czerminski R
2. Elber R
(1990) Self-avoiding walk between two fixed points as a tool to calculate reaction paths in large molecular systems
International Journal of Quantum Chemistry 38:167–185.
https://doi.org/10.1002/qua.560382419
- Google Scholar
(2014) Well-tempered metadynamics converges asymptotically
Physical Review Letters 112:240602.
https://doi.org/10.1103/PhysRevLett.112.240602
- Google Scholar
1. Eaton WA
2. Muñoz V
(2014)
Impact of atomistic molecular dynamics simulations on understanding how proteins fold: an experimentalist’s perspective

Bioinformatics.
- Google Scholar
1. Elber R
2. Roitberg A
3. Simmerling C
4. Goldstein R
5. Li H
6. Verkhivker G
7. Keasar C
8. Zhang J
9. Ulitsky A
(1995) MOIL: A program for simulations of macromolecules
Computer Physics Communications 91:159–189.
https://doi.org/10.1016/0010-4655(95)00047-J
- Google Scholar
(1992) A cavity-containing mutant of T4 lysozyme is stabilized by buried benzene
Nature 355:371–373.
https://doi.org/10.1038/355371a0
- Google Scholar
1. Eriksson A
2. Baase W
3. Zhang X
4. Heinz D
5. Blaber M
6. Baldwin E
7. Matthews B
(1992) Response of a protein structure to cavity-creating mutations and its relation to the hydrophobic effect
Science 255:178–183.
https://doi.org/10.1126/science.1553543
- Google Scholar
(1996) Access of ligands to cavities within the core of a protein is rapid
Nature Structural Biology 3:516–521.
https://doi.org/10.1038/nsb0696-516
- Google Scholar
1. Fraser JS
2. Clarkson MW
3. Degnan SC
4. Erion R
5. Kern D
6. Alber T
(2009) Hidden alternative structures of proline isomerase essential for catalysis
Nature 462:669–673.
https://doi.org/10.1038/nature08615
- Google Scholar
1. Fu B
2. Sahakyan AB
3. Camilloni C
4. Tartaglia GG
5. Paci E
6. Caflisch A
7. Vendruscolo M
8. Cavalli A
(2014) ALMOST: an all atom molecular simulation toolkit for protein structure determination
Journal of Computational Chemistry 35:1101–1105.
https://doi.org/10.1002/jcc.23588
- Google Scholar
(2008) Using relaxation dispersion NMR spectroscopy to determine structures of excited, invisible protein states
Journal of Biomolecular NMR 41:113–120.
https://doi.org/10.1007/s10858-008-9251-5
- Google Scholar
1. Hornak V
2. Abel R
3. Okur A
4. Strockbine B
5. Roitberg A
6. Simmerling C
(2006) Comparison of multiple Amber force fields and development of improved protein backbone parameters
Proteins: Structure, Function, and Bioinformatics 65:712–725.
https://doi.org/10.1002/prot.21123
- Google Scholar
1. Karplus M
2. Lavery R
(2014) Significance of molecular dynamics simulations for life sciences
Israel Journal of Chemistry 54:1042–1051.
https://doi.org/10.1002/ijch.201400074
- Google Scholar
(2009) Fast and accurate predictions of protein NMR chemical shifts from interatomic distances
Journal of the American Chemical Society 131:13894–13895.
https://doi.org/10.1021/ja903772t
- Google Scholar
1. Koshland DE
(1958) Application of a theory of enzyme specificity to protein synthesis
PNAS 44:98.
https://doi.org/10.1073/pnas.44.2.98
- Google Scholar
1. Kozlikova B
2. Sebestova E
3. Sustr V
4. Brezovsky J
5. Strnad O
6. Daniel L
7. Bednar D
8. Pavelka A
9. Manak M
10. Bezdeka M
11. Benes P
12. Kotry M
13. Gora A
14. Damborsky J
15. Sochor J
16. Sochor J
(2014) CAVER Analyst 1.0: graphic tool for interactive visualization and analysis of tunnels and channels in protein structures
Bioinformatics 30:2684–2685.
https://doi.org/10.1093/bioinformatics/btu364
- Google Scholar
1. Laio A
2. Parrinello M
(2002) Escaping free-energy minima
PNAS 99:12562–12566.
https://doi.org/10.1073/pnas.202427399
- Google Scholar
1. Lindorff-Larsen K
2. Piana S
3. Palmo K
4. Maragakis P
5. Klepeis JL
6. Dror RO
7. Shaw DE
(2010) Improved side-chain torsion potentials for the Amber ff99SB protein force field
Proteins: Structure, Function, and Bioinformatics 78:1950–1958.
https://doi.org/10.1002/prot.22711
- Google Scholar
(2009) Halogenated benzenes bound within a non-polar cavity in T4 lysozyme provide examples of I...S and I...Se halogen-bonding
Journal of Molecular Biology 385:595–605.
https://doi.org/10.1016/j.jmb.2008.10.086
- Google Scholar
1. Lukin JA
2. Kontaxis G
3. Simplaceanu V
4. Yuan Y
5. Bax A
6. Ho C
(2003) Quaternary structure of hemoglobin in solution
PNAS 100:517–520.
https://doi.org/10.1073/pnas.232715799
- Google Scholar
(2013) Conformational selection and adaptation to ligand binding in T4 lysozyme cavity mutants
PNAS 110:E4306–E4315.
https://doi.org/10.1073/pnas.1318754110
- Google Scholar
Book
(2008)
Pathways of conformational transitions in proteins

In: Voth G. A, editors. Coarse-Graining of Condensed Phase and Biomolecular Systems ((Ed).). Boca Raton: FL: CRC Press. pp. 185–203.
- Google Scholar
1. Marchi M
2. Ballone P
(1999) Adiabatic bias molecular dynamics: A method to navigate the conformational space of complex molecular systems
The Journal of Chemical Physics 110:3697–3702.
https://doi.org/10.1063/1.478259
- Google Scholar
1. Merski M
2. Fischer M
3. Balius TE
4. Eidam O
5. Shoichet BK
(2015) Homologous ligands accommodated by discrete conformations of a buried cavity
PNAS 112:5039–5044.
https://doi.org/10.1073/pnas.1500806112
- Google Scholar
(2015) Gaussian accelerated molecular dynamics: Unconstrained enhanced sampling and free energy calculation
Journal of Chemical Theory and Computation 11:3584–3595.
https://doi.org/10.1021/acs.jctc.5b00436
- Google Scholar
1. Mulder FA
2. Hon B
3. Mittermaier A
4. Dahlquist FW
5. Kay LE
(2002) Slow internal dynamics in proteins: application of NMR relaxation dispersion spectroscopy to methyl groups in a cavity mutant of T4 lysozyme
Journal of the American Chemical Society 124:1443–1451.
https://doi.org/10.1021/ja0119806
- Google Scholar
1. Mulder FA
2. Hon B
3. Muhandiram DR
4. Dahlquist FW
5. Kay LE
(2000) Flexibility and ligand exchange in a buried cavity mutant of T4 lysozyme studied by multinuclear NMR
Biochemistry 39:12614–12622.
https://doi.org/10.1021/bi001351t
- Google Scholar
1. Mulder FA
2. Mittermaier A
3. Hon B
4. Dahlquist FW
5. Kay LE
(2001) Studying excited states of proteins by NMR spectroscopy
Nature Structural Biology 8:932–935.
https://doi.org/10.1038/nsb1101-932
- Google Scholar
1. Neal S
2. Nip AM
3. Zhang H
4. Wishart DS
(2003)
Rapid and accurate calculation of protein 1H, 13C and 15N chemical shifts

Journal of Biomolecular NMR 26:215–240.
- Google Scholar
(1997) Theory of protein folding: the energy landscape perspective
Annual Review of Physical Chemistry 48:545–600.
https://doi.org/10.1146/annurev.physchem.48.1.545
- Google Scholar
1. Paci E
2. Karplus M
(1999) Forced unfolding of fibronectin type 3 modules: an analysis by biased molecular dynamics simulations
Journal of Molecular Biology 288:441–459.
https://doi.org/10.1006/jmbi.1999.2670
- Google Scholar
(2014) Conformational changes and free energies in a proline isomerase
Journal of Chemical Theory and Computation 10:4169–4174.
https://doi.org/10.1021/ct500536r
- Google Scholar
(2011) How robust are protein folding simulations with respect to force field parameterization?
Biophysical Journal 100:L47–L49.
https://doi.org/10.1016/j.bpj.2011.03.051
- Google Scholar
1. Pitera JW
2. Chodera JD
(2012) On the use of experimental observations to bias simulated ensembles
Journal of Chemical Theory and Computation 8:3445–3451.
https://doi.org/10.1021/ct300112v
- Google Scholar
1. Pronk S
2. Páll S
3. Schulz R
4. Larsson P
5. Bjelkmar P
6. Apostolov R
7. Shirts MR
8. Smith JC
9. Kasson PM
10. van der Spoel D
11. Hess B
12. Lindahl E
(2013) GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit
Bioinformatics 29:845–854.
https://doi.org/10.1093/bioinformatics/btt055
- Google Scholar
(2010) Using NMR chemical shifts as structural restraints in molecular dynamics simulations of proteins
Structure 18:923–933.
https://doi.org/10.1016/j.str.2010.04.016
- Google Scholar
1. Roux B
2. Weare J
(2013) On the statistical equivalence of restrained-ensemble simulations with the maximum entropy method
The Journal of Chemical Physics 138:084107.
https://doi.org/10.1063/1.4792208
- Google Scholar
(2012) Assessing the Performance of Metadynamics and Path Variables in Predicting the Binding Free Energies of p38 Inhibitors
Journal of Chemical Theory and Computation 8:1165–1170.
https://doi.org/10.1021/ct3001377
- Google Scholar
(2014) Assessing the Reliability of the Dynamics Reconstructed from Metadynamics
Journal of Chemical Theory and Computation 10:1420–1425.
https://doi.org/10.1021/ct500040r
- Google Scholar
1. Sekhar A
2. Kay LE
(2013) NMR paves the way for atomic level descriptions of sparsely populated, transiently formed biomolecular conformers
PNAS 110:12867–12874.
https://doi.org/10.1073/pnas.1305688110
- Google Scholar
1. Shen Y
2. Bax A
(2010) SPARTA+: a modest improvement in empirical NMR chemical shift prediction by means of an artificial neural network
Journal of Biomolecular NMR 48:13–22.
https://doi.org/10.1007/s10858-010-9433-9
- Google Scholar
(2012) Locating binding poses in protein-ligand systems using reconnaissance metadynamics
PNAS 109:5170–5175.
https://doi.org/10.1073/pnas.1201940109
- Google Scholar
1. Sutto L
2. Gervasio FL
(2013) Effects of oncogenic mutations on the conformational free-energy landscape of EGFR kinase
PNAS 110:10616–10621.
https://doi.org/10.1073/pnas.1221953110
- Google Scholar
(2007) Open-to-closed transition in apo maltose-binding protein observed by paramagnetic NMR
Nature 449:1078–1082.
https://doi.org/10.1038/nature06232
- Google Scholar
(2015) Kinetics of protein–ligand unbinding: Predicting pathways, rates, and rate-limiting steps
PNAS 112:E386–E391.
https://doi.org/10.1073/pnas.1424461112
- Google Scholar
1. Tiwary P
2. Mondal J
3. Morrone JA
4. Berne BJ
(2015) Role of water and steric constraints in the kinetics of cavity–ligand unbinding
PNAS 112:12015–12019.
https://doi.org/10.1073/pnas.1516652112
- Google Scholar
1. Tiwary P
2. Parrinello M
(2013) From metadynamics to dynamics
Physical Review Letters 111:230602.
https://doi.org/10.1103/PhysRevLett.111.230602
- Google Scholar
(2014) PLUMED 2: New feathers for an old bird
Computer Physics Communications 185:604–613.
https://doi.org/10.1016/j.cpc.2013.09.018
- Google Scholar
(2010) A self-learning algorithm for biased molecular dynamics
PNAS 107:17509–17514.
https://doi.org/10.1073/pnas.1011511107
- Google Scholar
(2011) Exploring the free energy surfaces of clusters using reconnaissance metadynamics
The Journal of Chemical Physics 135:114109.
https://doi.org/10.1063/1.3628676
- Google Scholar
1. Ulrich EL
2. Akutsu H
3. Doreleijers JF
4. Harano Y
5. Ioannidis YE
6. Lin J
7. Livny M
8. Mading S
9. Maziuk D
10. Miller Z
11. Nakatani E
12. Schulte CF
13. Tolmie DE
14. Kent Wenger R
15. Yao H
16. Markley JL
(2008) BioMagResBank
Nucleic Acids Research 36:D402–D408.
https://doi.org/10.1093/nar/gkm957
- Google Scholar
(2012) Studying "invisible" excited protein states in slow exchange with a major state conformation
Journal of the American Chemical Society 134:8148–8161.
https://doi.org/10.1021/ja3001419
- Google Scholar
(2016) Atomistic picture of conformational exchange in a T4 lysozyme cavity mutant: an experiment-guided molecular dynamics study
Chemical Science 7:3602–3613.
https://doi.org/10.1039/C5SC03886C
- Google Scholar
1. Wang J
2. Wang Y
3. Chu X
4. Hagen SJ
5. Han W
6. Wang E
(2011) Multi-scaled explorations of binding-induced folding of intrinsically disordered protein inhibitor IA3 to its target enzyme
PLoS Computational Biology 7:e1001118.
https://doi.org/10.1371/journal.pcbi.1001118
- Google Scholar
1. Wang Y
2. Tang C
3. Wang E
4. Wang J
(2012) Exploration of multi-state conformational dynamics and underlying global functional landscape of maltose binding protein
PLoS Computational Biology 8:e1002471.
https://doi.org/10.1371/journal.pcbi.1002471
- Google Scholar
1. Zhuravlev PI
2. Papoian GA
(2010) Protein functional landscapes, dynamics, allostery: a tortuous path towards a universal theoretical framework
Quarterly Reviews of Biophysics 43:295–332.
https://doi.org/10.1017/S0033583510000119
- Google Scholar

Article and author information

Author details

Yong Wang

Structural Biology and NMR Laboratory, Linderstrøm-Lang Centre for Protein Science, Department of Biology, University of Copenhagen, Copenhagen, Denmark

Contribution
YW, Conception and design, Acquisition of data, Analysis and interpretation of data, Drafting or revising the article

Competing interests
The authors declare that no competing interests exist.

"This ORCID iD identifies the author of this article:" 0000-0001-9156-0377
Elena Papaleo

Structural Biology and NMR Laboratory, Linderstrøm-Lang Centre for Protein Science, Department of Biology, University of Copenhagen, Copenhagen, Denmark

Contribution
EP, Acquisition of data, Analysis and interpretation of data

Competing interests
The authors declare that no competing interests exist.

"This ORCID iD identifies the author of this article:" 0000-0002-7376-5894
Kresten Lindorff-Larsen

Structural Biology and NMR Laboratory, Linderstrøm-Lang Centre for Protein Science, Department of Biology, University of Copenhagen, Copenhagen, Denmark

Contribution
KL-L, Conception and design, Analysis and interpretation of data, Drafting or revising the article

For correspondence
lindorff@bio.ku.dk

Competing interests
The authors declare that no competing interests exist.

"This ORCID iD identifies the author of this article:" 0000-0002-4750-6039

Funding

Novo Nordisk Foundation

Kresten Lindorff-Larsen

Carlsbergfondet

Kresten Lindorff-Larsen

Danish e-Infrastructure Cooperation

Kresten Lindorff-Larsen

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We would like to thank Pratyush Tiwary and Matteo Salvalaglio for the help on the usage and analysis of InMetaD, Gareth Tribello for the help on the usage of reconnaissance metadynamics, Lewis E. Kay, Pengfei Tian, Shuguang Yuan, Micha Ben Achim Kunze for their helpful discussions and comments. Ludovico Sutto and Francesco Gervasio are thanked for input on metadynamics simulations.

Version history

Received: May 4, 2016
Accepted: August 22, 2016
Accepted Manuscript published: August 23, 2016 (version 1)
Version of Record published: October 4, 2016 (version 2)

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

2,614

views
556

downloads
49

citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Mendeley

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Yong Wang
Elena Papaleo
Kresten Lindorff-Larsen

(2016)

Mapping transiently formed and sparsely populated conformations on a complex energy landscape

eLife 5:e17505.

https://doi.org/10.7554/eLife.17505

Share this article

Cite this article

Structures of the major G and minor E states of L99A T4L and the hidden state hypothesis.

Free energy landscape of the L99A variant of T4L.

Estimation of free energy differences and comparison with experimental measurements.

Conformational ensemble of the minor state as determined by CS biased, replica-averaged simulations.

Mechanisms of the G-E conformational exchanges explored by reconnaissance metadynamics.

Trajectory of the G-to-E conformational transition observed in Trj1, corresponding to the red trajectory in Figure 5.

Trajectory of the G-to-E conformational transition observed in Trj2, corresponding to the blue trajectory in Figure 5.

Trajectory of the G-to-E conformational transition observed in Trj3, corresponding to the green trajectory in Figure 5.

Trajectory of the E-to-G conformational transition observed in Trj3, corresponding to the yellow trajectory in Figure 5.

Movie of the calculated two-dimensional free energy landscape of L99A as a function of simulation time.

A transiently formed tunnel from the solvent to the cavity is a potential ligand binding pathway.

A typical trajectory of the benzene escaping from the buried cavity of L99A via tunnel #1 revealed by ABMD simulations.

Two representative InMetaD trajectories of L99A with G to E transitions.

Two representive InMetaD trajectories of L99A with E to G transitions of L99A.

Characteristic transition times between G and E states of L99A.

Characteristic transition times between G and E states of the triple mutant.

Poisson fit analysis for G to E transitions and E to G transitions of L99A.

Poisson fit analysis for G to E transitions and E to G transitions of the triple mutant.

Author details

Yong Wang

Contribution

Competing interests

Elena Papaleo

Contribution

Competing interests

Kresten Lindorff-Larsen

Contribution

For correspondence

Competing interests

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

Categories and tags