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

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.

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, , 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 Å3) 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 (; 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 ), but perturbing the structure in regions that the experiments demonstrated to undergo conformational change in a way so that the final model () 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 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 (). 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.

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 (; 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 (; 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 ) rotates and fuses with helix G () into a longer helix, and the side chain of phenylalanine at position 114 () rotates so as to occupy part of the cavity. As the cavity is inaccessible in both the and structures it has been hypothesized that ligand entry occurs via a third ‘cavity open’ state (Merski et al., 2015).

DOI: http://dx.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 ( and ) (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 variable describes the progress of the conformational transition between the and structures with additional ‘interpolation’ using an optimal ‘reference’ path in a simplified model (see details in Appendix and Figure 2—figure supplement 1), while measures the distance to this reference path. In this way, the two-dimensional free energy landscape along and 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.

Appendix 1—table 1.

Simulation details.

DOI: http://dx.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−1 · nm−2, S_target=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−1 · nm−2, S_target=4.0 nm	298	PLUMED2.2, GMX5.1.2

† is the bias factor.

‡ is the characteristic decay time used for the dynamically-adapted Gaussian potential.

§ is the deposition frequency of Gaussian potential.

¶ is the height of the deposited Gaussian potential (in KJ · mol−1).

# N is the number of replicas.

** is the strength of chemical shift restraints (in KJ · mol−1 · ppm−2).

†† means RUN_FREQ and 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.

¶¶ is the Gaussian width of .

## is the Gaussian width of (in nm2).

***36 trajectories collected for E-to-G transitions.

Figure 2—figure supplement 1.

Approximately equidistant frames along the reference path.

The plot reveals a ‘gullwing’ shape of the matrix of pairwise RMSDs of the frames of the reference path, indicating that frames along the reference path are approximately equidistant. We used 31 structures to discretize the path.

DOI: http://dx.doi.org/10.7554/eLife.17505.005

Projecting the sampled free energy landscape along (upper panel of Figure 2) reveals a deep, narrow free energy basin around (labeled by red sphere and corresponding to the G state), and a broader, shallow free energy basin with 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(, ), 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 and correspond to the structures of and , respectively.

Figure 2.

Free energy landscape of the L99A variant of T4L.

In the upper panel, we show the projection of the free energy along , 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 and 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 ranging from 0.55 to 0.83. In the lower panel, we show the three-dimensional negative free energy landscape, -F(, ), 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 and  nm2, which we denote , is labeled by a yellow dot.

DOI: http://dx.doi.org/10.7554/eLife.17505.004

The plot reveals a ‘gullwing’ shape of the matrix of pairwise RMSDs of the frames of the reference path, indicating that frames along the reference path are approximately equidistant. We used 31 structures to discretize the path.

DOI: http://dx.doi.org/10.7554/eLife.17505.005

(A) The two-dimensional free energy surface F(,) of L99A sampled by a 667 ns PathMetaD simulation. (B) The two-dimensional free energy surface F(,) of the triple mutant sampled by a 961 ns PathMetaD simulation. (C) The free energy profiles as a function of of both L99A (black) and the triple mutant (blue).

DOI: http://dx.doi.org/10.7554/eLife.17505.006

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 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 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 and states, we also found a free energy minimum around and  nm2 (denoted as 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 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 and 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 =0.7 and a higher energy conformation around =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.

Figure 2—figure supplement 2.

One and two dimensional free energy landscape of L99A and the triple mutant.

(A) The two-dimensional free energy surface F(,) of L99A sampled by a 667 ns PathMetaD simulation. (B) The two-dimensional free energy surface F(,) of the triple mutant sampled by a 961 ns PathMetaD simulation. (C) The free energy profiles as a function of of both L99A (black) and the triple mutant (blue).

DOI: http://dx.doi.org/10.7554/eLife.17505.006

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, , 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 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 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 =3.5 kcal mol−1 (Figure 3, black line). This value can be compared to the value of 2.1 kcal mol−1 obtained from NMR relaxation dispersion experiments (Mulder et al., 2001), revealing reasonably good, albeit not exact, agreement with the experiments.

Figure 3.

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 (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 () from their populations. The time evolution of 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−1 for L99A and −1.9 kcal mol−1 for the triple mutant) are shown as yellow lines.

DOI: http://dx.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−1 (Figure 3, blue line), in excellent agreement with the experimental measurement (−1.9 kcal mol−1) (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−1 in simulations compared to 4.0 kcal mol−1 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.

Appendix 1—table 2.

Definition of collective variables.

DOI: http://dx.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

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 state or 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 1 kcal mol−1. 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).

Appendix 1—figure 1.

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 > 0.55 and Z < 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.

DOI: http://dx.doi.org/10.7554/eLife.17505.030

Appendix 1—figure 2.

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 Spath < 0.28 and Zpath < −0.01.

DOI: http://dx.doi.org/10.7554/eLife.17505.031

Table 1.

Free energy differences and rates of conformational exchange.

DOI: http://dx.doi.org/10.7554/eLife.17505.008

	τ_G→E (ms)	τ_E→G (ms)	Δ⁢G (kcal mol−1)
L99A
NMR	20	0.7	2.1
InMetaD	175± 56	1.4± 0.6	2.9± 0.5
PathMetaD			3.5
L99A/G113A/R119P
NMR	0.2	4	-1.9
InMetaD	2.0± 1.7	14.3± 8.3	-1.2± 1.1
PathMetaD			-1.6

Appendix 1—figure 3.

Characteristic transition times between G and E states of L99A.

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

DOI: http://dx.doi.org/10.7554/eLife.17505.032

Appendix 1—figure 4.

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

The figure shows the characteristic transition time τ (left panel) and τ (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.

DOI: http://dx.doi.org/10.7554/eLife.17505.033

Appendix 1—figure 5.

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.

DOI: http://dx.doi.org/10.7554/eLife.17505.034

Appendix 1—figure 6.

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.

DOI: http://dx.doi.org/10.7554/eLife.17505.035

The ability to calculate forward and backward rates between and 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.90.5 kcal mol−1 and −1.21.1 kcal mol−1 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−1 and −1.6 kcal mol−1) and experiment (2.1 kcal mol−1 and −1.9 kcal mol−1). 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−1) 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 [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 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 state, is negligible in comparison). Thus, we compared the experimental CS values of the minor state with the values calculated using the 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 (, , , and ) 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.

Figure 4—figure supplement 1.

Equilibrium sampling of conformational regions of the E state of L99A by CS-restrained replica-averaged simulation.

We calculated the RMSD between the experimental CSs and the values back-calculated by CamShift (Kohlhoff et al., 2009) as implemented in ALMOST (Fu et al., 2014). We projected a 250 ns MD trajectory sampled using the CHARMM22* force field (RUN3 in Appendix 1—table 1) onto the RMSDs. The average RMSDs for the five measured nuclei (, , , and ) are 0.23 ppm, 0.38 ppm, 1.97 ppm, 0.83 ppm and 1.06 ppm, respectively (Appendix 1—table 2), which are close to the inherent uncertainty of the chemical shift calculation (Figure 4—figure supplement 2). This indicates the simulation yielded an ensemble in good agreement with experiments.

DOI: http://dx.doi.org/10.7554/eLife.17505.010

Figure 4—figure supplement 2.

Estimation of the inherent uncertainty of the chemical shift calculation by different algorithms: CamShift (Kohlhoff et al., 2009), ShiftX (Neal et al., 2003) and Sparta+ (Shen and Bax, 2010).

Using as the reference structure, we calculated the chemical shifts using different algorithms and compared the correlation coefficients and RMSD between them.

DOI: http://dx.doi.org/10.7554/eLife.17505.011

Appendix 1—table 3.

Average root-mean-square deviation ( in units of ppm) between experimental CSs and those from the CS-restrained replica-averaged simulations.

DOI: http://dx.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

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 and 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.

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 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 variable (cyan).

DOI: http://dx.doi.org/10.7554/eLife.17505.009

We calculated the RMSD between the experimental CSs and the values back-calculated by CamShift (Kohlhoff et al., 2009) as implemented in ALMOST (Fu et al., 2014). We projected a 250 ns MD trajectory sampled using the CHARMM22* force field (RUN3 in Appendix 1—table 1) onto the RMSDs. The average RMSDs for the five measured nuclei (, , , and ) are 0.23 ppm, 0.38 ppm, 1.97 ppm, 0.83 ppm and 1.06 ppm, respectively (Appendix 1—table 2), which are close to the inherent uncertainty of the chemical shift calculation (Figure 4—figure supplement 2). This indicates the simulation yielded an ensemble in good agreement with experiments.

DOI: http://dx.doi.org/10.7554/eLife.17505.010

Using as the reference structure, we calculated the chemical shifts using different algorithms and compared the correlation coefficients and RMSD between them.

DOI: http://dx.doi.org/10.7554/eLife.17505.011

The chemical shifts of the E state of L99A (BMRB 17604) were used. (A) The simulation with CHARMM22* force field. (B) The simulation with Amber ff99SB*-ILDN force field.

DOI: http://dx.doi.org/10.7554/eLife.17505.012

(A) N = 4, KJ · mol−1. (B) N = 2, KJ · mol−1. (C) N = 2, KJ · mol−1. N refers to the number of replicas that the CS values are averaged over. The CHARMM22* force field was used in these simulations. The results also support the conclusion that the conformational space of the minor (E) state covers a relatively wide set of structures including the structure.

DOI: http://dx.doi.org/10.7554/eLife.17505.013

Chemical shift restraints were from BMRB 17,603 and CHARMM22* force field was used.

DOI: http://dx.doi.org/10.7554/eLife.17505.014

Figure 4—figure supplement 3.

Force field dependency of the replica averaged MD simulations of L99A with chemical shift restraints.

The chemical shifts of the E state of L99A (BMRB 17604) were used. (A) The simulation with CHARMM22* force field. (B) The simulation with Amber ff99SB*-ILDN force field.

DOI: http://dx.doi.org/10.7554/eLife.17505.012

Figure 4—figure supplement 4.

Effect of changing the force constant and number of replicas in CS-restrained simulation of L99A.

(A) N = 4, KJ · mol−1. (B) N = 2, KJ · mol−1. (C) N = 2, KJ · mol−1. N refers to the number of replicas that the CS values are averaged over. The CHARMM22* force field was used in these simulations. The results also support the conclusion that the conformational space of the minor (E) state covers a relatively wide set of structures including the structure.

DOI: http://dx.doi.org/10.7554/eLife.17505.013

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 (Figure 4).

Figure 4—figure supplement 5.

Replica-averaged CS-restrained MD simulation of a T4L triple mutant (L99A/G113A/R119P).

Chemical shift restraints were from BMRB 17,603 and CHARMM22* force field was used.

DOI: http://dx.doi.org/10.7554/eLife.17505.014

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(,) (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 state or indirectly via it.

Figure 5—figure supplement 1.

Complete G-to-E transitions of L99A obtained by reconnaissance metadynamics simulations.

The state-specific fraction of contacts (Wang et al., 2012), and , were employed to monitor the conformational transitions to G and E state, respectively. Trajectories Trj1, Trj2 and Trj3 are from the simulations RUN10, RUN11 and RUN12 (Appendix 1—table 1), respectively.

DOI: http://dx.doi.org/10.7554/eLife.17505.016

Figure 5.

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 state or indirectly via it.

DOI: http://dx.doi.org/10.7554/eLife.17505.015

The state-specific fraction of contacts (Wang et al., 2012), and , were employed to monitor the conformational transitions to G and E state, respectively. Trajectories Trj1, Trj2 and Trj3 are from the simulations RUN10, RUN11 and RUN12 (Appendix 1—table 1), respectively.

DOI: http://dx.doi.org/10.7554/eLife.17505.016

Trajectories Trj1 (magenta), Trj2 (blue) and Trj3 (green and orange) are from the simulations RUN10, RUN11 and RUN12 (Appendix 1—table 2), respectively. The steepest descent path (SDP, black) used to define the initial path in PathMetaD is also shown as a reference. To measure the distance between helix F and helix I, and between F144 and helix D, we employed order parameters and . is defined as the distance between E108 and R137, while is defined as the distance between the atom of F114 and the atom of Y88.

DOI: http://dx.doi.org/10.7554/eLife.17505.017

The figure suggests in the direct transitions (Trj1 and first half of Trj3) F114 can rotate its side chain inside the protein core. In contrast, in the route (Trj2 and second half of Trj3) the side chain of F114, which occupies the cavity in the E state, gets transiently exposed to solvent during the transition.

DOI: http://dx.doi.org/10.7554/eLife.17505.018

In the context of coarse-grained kinetic models the results above would suggest at least two possible mechanisms operate in parallel: or . 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 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 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).

Figure 5—figure supplement 2.

Conformational transitions between the G and E states monitored by other order parameters.

Trajectories Trj1 (magenta), Trj2 (blue) and Trj3 (green and orange) are from the simulations RUN10, RUN11 and RUN12 (Appendix 1—table 2), respectively. The steepest descent path (SDP, black) used to define the initial path in PathMetaD is also shown as a reference. To measure the distance between helix F and helix I, and between F144 and helix D, we employed order parameters and . is defined as the distance between E108 and R137, while is defined as the distance between the atom of F114 and the atom of Y88.

DOI: http://dx.doi.org/10.7554/eLife.17505.017

Video 1.

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.

DOI: http://dx.doi.org/10.7554/eLife.17505.019

Video 4.

Trajectory of the E-to-G conformational transition observed in Trj3, corresponding to the yellow 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.

DOI: http://dx.doi.org/10.7554/eLife.17505.022

Figure 5—figure supplement 3.

Solvent accessible surface area (SASA) calculation of the side chain of F114.

The figure suggests in the direct transitions (Trj1 and first half of Trj3) F114 can rotate its side chain inside the protein core. In contrast, in the route (Trj2 and second half of Trj3) the side chain of F114, which occupies the cavity in the E state, gets transiently exposed to solvent during the transition.

DOI: http://dx.doi.org/10.7554/eLife.17505.018

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

DOI: http://dx.doi.org/10.7554/eLife.17505.020

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

DOI: http://dx.doi.org/10.7554/eLife.17505.021

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 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.

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

The figure shows the time evolution of the free energy surface as a function of and sampled in a 667 ns PathMetaD simulation of L99A.

DOI: http://dx.doi.org/10.7554/eLife.17505.023

Figure 6—figure supplement 1.

A transiently formed tunnel from the solvent to the cavity forms in the state.

(A) Typical structures from the state sampled in the simulation (between 430 ns and 447 ns) are mapped onto the free energy surface, and represented by yellow points. (B) The cavity-related regions (helix E, F and G) are coloured in blue, while F114 is coloured in red. F114 tends to be partially solvent exposed in , resulting in the internal cavity to be open. The tunnel connecting the internal cavity and protein surface is coloured in yellow, and has a peanut-shell like shape. (C) shows the radius along the tunnel of structures belong to the cluster of tunnel. Lines in different colours represent different structures. Grey dotted line represents the average tunnel radius.

DOI: http://dx.doi.org/10.7554/eLife.17505.025

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) having a entrance located at the groove between and . A typical representative structure of is shown in Figure 6A. The radii along the structures in cluster vary, but share an overall shape (Figure 6—figure supplement 1), and we find that the maximal bottleneck radius is 2.5 Å, the average bottleneck radius is 1.3 Å, and the average length 11.2 Å.

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).

DOI: http://dx.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 Å

Figure 6—figure supplement 2.

Representative structures of the cavity region in the state.

The figure shows six representative structures of the cavity region revealing multiple tunnels that connect the cavity with the protein surface. The different colours correspond to different tunnels observed, and a structure can have different tunnels with different widths present at the same time. The colours represent the relative size with yellow, purple and green corresponding to tunnels of decreasing width.

DOI: http://dx.doi.org/10.7554/eLife.17505.026

Figure 6—figure supplement 3.

Tunnel clustering analysis on state.

The clustering of tunnels was performed using the CAVER Analyst software (Kozlikova et al., 2014) and the average-link hierarchical algorithm based on the calculated matrix of pairwise tunnel distances. We found that the most weighted tunnel (denoted as tunnel) populates of the basin. The second and third tunnels populate 20% and 15%, respectively, but their maximal bottleneck radii are 1.4 and 1.3 Å, in contrast to the maximal bottleneck radius of tunnel#1 of 2.5 Å. (A) Heat map visualization of the tunnel profile of tunnel. The colour map represents the radius of the tunnel along the tunnel length. (B) Average tunnel radius and minimal tunnel radius of individual structures belonging to tunnel cluster. Note that the gaps indicate these snapshots do not have tunnels. (C) The tunnel radius as a function of the tunnel index which is ranked by the average radius (R). The second widest tunnel (tunnel) has the highest population and is highlighted in yellow. (D) A typical structure of with an open tunnel. H, and are coloured in blue, F114 is coloured in red, and tunnel is coloured in yellow. (E) The figure shows the location of an alkylbenzene (magenta) in a crystal structure of L99A T4L (PDB ID: 4W59). The figure further shows (in yellow) the tunnel induced in the structure by the alkyl chain, as revealed by CAVER3 when applied to the structure after removing the ligand. Because the tunnel overlaps with the alkyl chain of the ligand, only the phenyl moiety of the ligand is visible.

DOI: http://dx.doi.org/10.7554/eLife.17505.027

Appendix 1—table 5.

Clustering analysis of tunnels (top three listed).

DOI: http://dx.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

Figure 6.

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), that has an entrance located at the groove between helix F () and helix 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).

DOI: http://dx.doi.org/10.7554/eLife.17505.024

(A) Typical structures from the state sampled in the simulation (between 430 ns and 447 ns) are mapped onto the free energy surface, and represented by yellow points. (B) The cavity-related regions (helix E, F and G) are coloured in blue, while F114 is coloured in red. F114 tends to be partially solvent exposed in , resulting in the internal cavity to be open. The tunnel connecting the internal cavity and protein surface is coloured in yellow, and has a peanut-shell like shape. (C) shows the radius along the tunnel of structures belong to the cluster of tunnel. Lines in different colours represent different structures. Grey dotted line represents the average tunnel radius.

DOI: http://dx.doi.org/10.7554/eLife.17505.025

The figure shows six representative structures of the cavity region revealing multiple tunnels that connect the cavity with the protein surface. The different colours correspond to different tunnels observed, and a structure can have different tunnels with different widths present at the same time. The colours represent the relative size with yellow, purple and green corresponding to tunnels of decreasing width.

DOI: http://dx.doi.org/10.7554/eLife.17505.026

The clustering of tunnels was performed using the CAVER Analyst software (Kozlikova et al., 2014) and the average-link hierarchical algorithm based on the calculated matrix of pairwise tunnel distances. We found that the most weighted tunnel (denoted as tunnel) populates of the basin. The second and third tunnels populate 20% and 15%, respectively, but their maximal bottleneck radii are 1.4 and 1.3 Å, in contrast to the maximal bottleneck radius of tunnel#1 of 2.5 Å. (A) Heat map visualization of the tunnel profile of tunnel. The colour map represents the radius of the tunnel along the tunnel length. (B) Average tunnel radius and minimal tunnel radius of individual structures belonging to tunnel cluster. Note that the gaps indicate these snapshots do not have tunnels. (C) The tunnel radius as a function of the tunnel index which is ranked by the average radius (R). The second widest tunnel (tunnel) has the highest population and is highlighted in yellow. (D) A typical structure of with an open tunnel. H, and are coloured in blue, F114 is coloured in red, and tunnel is coloured in yellow. (E) The figure shows the location of an alkylbenzene (magenta) in a crystal structure of L99A T4L (PDB ID: 4W59). The figure further shows (in yellow) the tunnel induced in the structure by the alkyl chain, as revealed by CAVER3 when applied to the structure after removing the ligand. Because the tunnel overlaps with the alkyl chain of the ligand, only the phenyl moiety of the ligand is visible.

DOI: http://dx.doi.org/10.7554/eLife.17505.027

The figure shows how ABMD simulations allow us to observe the ligand benzene escape from the internal binding site. We performed two sets of 20 simulations using two different force constants for the ABMD (upper: 50 KJ · mol−1 · nm−2; lower: 20 KJ · mol−1 · nm−2); note also the different time scales on the two plots. The simulations used the RMSD of the ligand to the bound state as the reaction coordinate, but are here shown projected onto the distance between the benzene molecule and the side chain of F114. The three structures in the bottom panel provide representative structures.

DOI: http://dx.doi.org/10.7554/eLife.17505.028

Video 6.

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

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

DOI: http://dx.doi.org/10.7554/eLife.17505.029

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 structures, the location of the tunnel exit is consistent with the dominant tunnel in (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 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 (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.

Figure 6—figure supplement 4.

Ligand unbinding pathways revealed by ABMD simulations.

The figure shows how ABMD simulations allow us to observe the ligand benzene escape from the internal binding site. We performed two sets of 20 simulations using two different force constants for the ABMD (upper: 50 KJ · mol−1 · nm−2; lower: 20 KJ · mol−1 · nm−2); note also the different time scales on the two plots. The simulations used the RMSD of the ligand to the bound state as the reaction coordinate, but are here shown projected onto the distance between the benzene molecule and the side chain of F114. The three structures in the bottom panel provide representative structures.

DOI: http://dx.doi.org/10.7554/eLife.17505.028

Appendix 1—table 6.

Unbinding Pathways Explored by ABMD ( as CV).

DOI: http://dx.doi.org/10.7554/eLife.17505.041

	k=20 kJ/(mol · nm−2)		k=50 kJ/(mol · nm−2)
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)

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−1. 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 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

Our simulations were initiated in the experimentally determined structures of the ground state of L99A (; PDB ID code 3DMV) or minor, E state (; 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 (). Details can be found in the Appendix.

Initial reaction path

Taking and 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

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

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 and were used as the starting conformations. CS data of and were obtained from the BMRB database (Ulrich et al., 2008) as entries 17604 and 17603, respectively.

Reconnaissance metadynamics simulations

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 using Gromacs4.5.5 with PLUMED1.3 plugin. See parameter details in Appendix 1—table 1.

Calculating kinetics using infrequent metadynamics

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.

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Mapping transiently formed and sparsely populated conformations on a complex energy landscape" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom, Yibing Shan, is a member of our Board of Reviewing Editors, and another is James S Fraser (Reviewer #3), and the evaluation has been overseen by Michael Marletta as the Senior Editor. The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

In this manuscript, the authors employ state-of-art enhanced sampling techniques, including multiple variants of metadynamics, replica-averaged restrained MD and well-tempered replica-exchange MD, to conduct a systematic characterization of conformational dynamics of the L99A variant of T4L. This simulation study, which is shown to be in agreement with experimental measures of transition rates and kinetics, provides insights on the possible mechanism of conformational exchange between the ground state and the excited state. In particular, the excited state is shown to be conformationally broad and heterogeneous, and the transition path between the two states to be diverse. The work also identifies an key intermediate state for the transition, although its confirmation requires further experimental evidence.

Essential revisions:

Appreciating the thoroughness and the technical sophistication of the computation work and its robust connections with experimental data, the reviewers recognize the manuscript's value in benchmarking state-of-art simulation techniques and especially, in benchmarking conformational free energy calculations.

One of the more interesting parts of the of this study concerns the calculation of kinetic rate of the transition for L99A and the results is very encouraging. Because of its potential importance, this work would benefit greatly from additional calculations on mutants besides L99A. Specifically, the transition kinetics of L99A/G113A and L99A/G113A/R119P have been experimentally characterized. Calculating the kinetic rates of these two mutants and comparing the results with the experimental data would much strengthen case for the validity of such calculations.

The manuscript's attempt to elucidate the mechanism of ligand entrance is relatively weak. Further studies, for example, using simulations with ligand placed at the putative entrance to the channel in a partially open conformation would strengthen this part of the manuscript. If unsuccessful, the authors should consider weaken the related claims and de-emphasize the part of the work.

Relatedly, the reviewers think the manuscript's value lies more in benchmarking and validating state-of-art free energy calculation, than in providing significant new insight into the mechanism and function of T4L conformational transition. We would like to see revisions that de-emphasize the T4L specific claims and additional discussions with respect to benchmarking quantitative characterization of free energy landscapes using atomistic simulations.

Essential revisions:

Appreciating the thoroughness and the technical sophistication of the computation work and its robust connections with experimental data, the reviewers recognize the manuscript's value in benchmarking state-of-art simulation techniques and especially, in benchmarking conformational free energy calculations.

One of the more interesting parts of the of this study concerns the calculation of kinetic rate of the transition for L99A and the results is very encouraging. Because of its potential importance, this work would benefit greatly from additional calculations on mutants besides L99A. Specifically, the transition kinetics of L99A/G113A and L99A/G113A/R119P have been experimentally characterized. Calculating the kinetic rates of these two mutants and comparing the results with the experimental data would much strengthen case for the validity of such calculations.

We thank the reviewers for their thoughtful comments on our work. As discussed in more detail below, we agree that the benchmarking aspects of our work are important. We have therefore extended our work and performed additional simulations and analyses. In particular, we have calculated the kinetics of exchange of the L99A/G113A/R119P triple mutant. As in the case of the L99A variant, we found results that were in good agreement with experiment. Specifically, the rates of exchange are within a factor of 10 from the experimental rates, and the free energy difference estimated from the rates (-1.2 kcal/mol) is very similar to the value (-1.6 kcal/mol) obtained from the equilibrium sampling. These results are presented in the revised manuscript.

As suggested by the reviewers we have also performed extensive simulations of the L99A/G113A double mutant. At the end of a 1 microsecond long PathMetaD simulation, the free energy difference obtained is ~ -0.5 kcal/mol, in apparently good agreement with experiment (0.3 kcal/mol). In contrast to the results from L99A and the triple mutant we, however, find that this value appears not to be converged with substantial fluctuations and a (non-converged) free energy landscape that differs in shape from that of L99A and L99A/G113A/R119P. Thus, the calculated free energy difference is (in contrast to L99A and the triple mutant) more sensitive to how the major and minor states are separated. Thus, the good agreement that we observe could very well be fortuitous, and we have therefore opted not to include them in the revised manuscript. Instead we mention that we have performed these simulations, but that they did not converge.

The manuscript's attempt to elucidate the mechanism of ligand entrance is relatively weak. Further studies, for example, using simulations with ligand placed at the putative entrance to the channel in a partially open conformation would strengthen this part of the manuscript. If unsuccessful, the authors should consider weaken the related claims and de-emphasize the part of the work.

We agree that the link between the transiently opened tunnel and the pathway for ligand binding and release was weak (as indeed we stated). While we believe that a more quantitative study lies outside the scope of this work, we have now performed and include the results from a set of simulations aimed to demonstrate a more direct link between the tunnel observed in the apo protein and the pathway for ligand (un)binding. In particular, we have used so-called adiabatic biased molecular dynamics (ABMD) to study the mechanism by which the ligand escapes its internal binding site. In these simulations, a bias is introduced for the ligand to leave its binding site, but leaving the specific pathway to be determined by the molecular simulations. Of the 20 ligand unbinding events we observed in such simulations, 15 take a path located in the same region as the transiently opened tunnel we found in the apo protein. To examine the extent to which this result depends on the amount of bias used, we performed these simulations with two different strengths of the biasing force, and reassuringly found the same results in both cases (i.e. 15/20 in both cases). These new results are included in the manuscript. We now also explicitly show the tunnel induced by binding of an alkyl-benzene to the protein.

Relatedly, the reviewers think the manuscript's value lies more in benchmarking and validating state-of-art free energy calculation, than in providing significant new insight into the mechanism and function of T4L conformational transition. We would like to see revisions that de-emphasize the T4L specific claims and additional discussions with respect to benchmarking quantitative characterization of free energy landscapes using atomistic simulations.

We agree that a key result of our work lies in the testing and benchmarking of the ability of molecular simulations to map conformational free energy landscapes and determine the free energy and kinetics associated with conformational exchange. Prompted by the suggestion of the reviewers we have thus strengthened these aspects in several parts of the manuscript. This also includes the newly described kinetics for the triple mutant, as well as the fact that for both L99A and the triple mutant we find good agreement between the kinetic and equilibrium estimate of the free energy difference. In addition, we have toned down some of the more specific points that mostly relate to T4L. For example, we have now completely removed the section and results describing our studies of domain-domain motions in the minor state of L99A T4L. We have also shortened somewhat the description of the results relating to the mechanism of conformational exchange, and finally we have reorganized the sections to bring together the validation parts more coherently. We believe that these changes together make the focus of the paper stronger.