1. Biochemistry and Chemical Biology
  2. Structural Biology and Molecular Biophysics
Download icon

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

  1. Yong Wang
  2. Elena Papaleo
  3. Kresten Lindorff-Larsen  Is a corresponding author
  1. University of Copenhagen, Denmark
Research Article
Cite this article as: eLife 2016;5:e17505 doi: 10.7554/eLife.17505
12 figures, 6 videos and 7 tables


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

Figure 2 with 2 supplements
Free energy landscape of the L99A variant of T4L.

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

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.

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(Spath,Zpath) of L99A sampled by a 667 ns PathMetaD simulation. (B) The two-dimensional free energy surface F(Spath,Zpath) of the triple mutant sampled by a 961 ns PathMetaD simulation. (C) The free energy profiles as a function of Spath of both L99A (black) and the triple mutant (blue).

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 Spath=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−1 for L99A and −1.9 kcal mol−1 for the triple mutant) are shown as yellow lines.

Figure 4 with 5 supplements
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 Spath 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 Spath variable (cyan).

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 (Hα, HN, N, C and Cα) 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.

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 EROSETTA as the reference structure, we calculated the chemical shifts using different algorithms and compared the correlation coefficients and RMSD between them.

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.

Figure 4—figure supplement 4
Effect of changing the force constant and number of replicas in CS-restrained simulation of L99A.

(A) N = 4, ϵCS=24 KJ · mol−1. (B) N = 2, ϵCS=24 KJ · mol−1. (C) N = 2, ϵCS=12 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 EROSETTA structure.

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.

Figure 5 with 3 supplements
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 I0.36 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.

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 RHF-HI and RF114-HD. is defined as the Cα distance between E108 and R137, while RF114-HD is defined as the distance between the Cδ4 atom of F114 and the Cα atom of Y88.

Figure 5—figure supplement 3
Solvent accessible surface area (SASA) calculation of the side chain of F114.

The figure suggests in the direct GE transitions (Trj1 and first half of Trj3) F114 can rotate its side chain inside the protein core. In contrast, in the GI0.36E 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.

Figure 6 with 4 supplements
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 (HF) and helix I (HI). 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).

Figure 6—figure supplement 1
A transiently formed tunnel from the solvent to the cavity forms in the I0.36 state.

(A) Typical structures from the I0.36 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 I0.36, resulting in the internal cavity to be open. The tunnel#1 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#1. Lines in different colours represent different structures. Grey dotted line represents the average tunnel radius.

Figure 6—figure supplement 2
Representative structures of the cavity region in the I0.36 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.

Figure 6—figure supplement 3
Tunnel clustering analysis on I0.36 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#1) populates 27% of the I0.36 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#1. The colour map represents the radius of the tunnel#1 along the tunnel length. (B) Average tunnel radius and minimal tunnel radius of individual structures belonging to tunnel#1 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#1) has the highest population and is highlighted in yellow. (D) A typical structure of I0.36 with an open tunnel#1. HE, HF and HG are coloured in blue, F114 is coloured in red, and tunnel#1 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.

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.

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 Spath > 0.55 and Zpath < 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.

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.

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.

Appendix 1—figure 4
Characteristic transition times between G and E states of the triple mutant.

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

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.

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 GE (A) and EG (B) transition times as a function of the size of a subsample of transition times randomly extracted from the main complete sample.



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.

Video 2
Trajectory of the G-to-E conformational transition observed in Trj2, corresponding to the blue 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.

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

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.

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 Spath and Zpath sampled in a 667 ns PathMetaD simulation of L99A.

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.



Table 1

Free energy differences and rates of conformational exchange.

τGE (ms)τEG (ms)ΔG (kcal mol−1)
InMetaD175 ± 561.4 ± 0.62.9 ± 0.5
InMetaD2.0 ± 1.714.3 ± 8.3-1.2 ± 1.1
Appendix 1—table 1

Simulation details.

Methodlabelsysteminitlengthforce fieldCVsparametersT (K)software version
PathMetaDRUN1L99AG667 nsCHARMM22*Spath,Zpathγ=20, τD=1 ps,τG=1 ps§, ω0=1.5298PLUMED2.1, GMX4.6
PathMetaDRUN2L99A, G113A, R119PG961 nsCHARMM22*Spath,Zpathγ=20, τD=1 ps, τG=1ps, ω0=1.5298PLUMED2.1, GMX5.0
CS-restrained MDRUN3L99AE252 nsCHARMM22*N#=4, ϵCS=24**298PLUMED2.1, ALMOST2.1, GMX5
CS-restrained MDRUN4L99AE233 nsCHARMM22*N=2, ϵCS=24298PLUMED2.1, ALMOST2.1, GMX5
CS-restrained MDRUN5L99AE221 nsCHARMM22*N=2, ϵCS=12298PLUMED2.1, ALMOST2.1, GMX5
CS-restrained MDRUN6L99AE125 nsAMBER99sb*ILDNN=4, ϵCS=24298PLUMED2.1, ALMOST2.1, GMX5
CS-restrained MDRUN7L99AE190 nsAMBER99sb*ILDNN=2, ϵCS=12298PLUMED2.1, ALMOST2.1, GMX5
CS-restrained MDRUN8L99A, G113A, R119PG204 nsCHARMM22*N=4, ϵCS=24298PLUMED2.1, ALMOST2.1, GMX5
CS-restrained MDRUN9L99A, G113A, R119PG200 nsCHARMM22*N=2, ϵCS=24298PLUMED2.1, ALMOST2.1, GMX5
Reconnaissance MetaDRUN10L99AG120 nsCHARMM22*2,3,4,5ΔT=10000††, Δt=250298PLUMED1.3, GMX4.5
Reconnaissance MetaDRUN11L99AG85 nsCHARMM22*2,3,4,510000, 250298PLUMED1.3, GMX4.5
Reconnaissance MetaDRUN12L99AG41 nsCHARMM22*2,3,4,5,7,8100000, 500298PLUMED1.3, GMX4.5
PT-WT-MetaDRUN13L99AG961 nsCHARMM22*1,2,3,7,8γ=20, ω0=0.5297 298‡‡ 303 308 316 325 341 350 362PLUMED 1.3, GMX4.5
PT-WT-MetaDRUN14L99AG404 nsCHARMM22*1,2,3,6,7γ=20, ω0=0.5297 298 303 308 316 325 333 341 350PLUMED 1.3, GMX4.5
PT-WT-MetaDRUN15L99AG201 nsCHARMM22*1,2,3,6γ=20, ω0=0.5297 298 303 308 316 325 333 341 350PLUMED 1.3, GMX4.5
PT-WT-MetaDRUN16L99AG511 nsCHARMM22*1,2,3,4,5γ=5, ω0=0.5295 298 310 325 341 358 376PLUMED 1.3, GMX4.5
PT-WT-MetaDRUN17L99AG383 nsCHARMM22*1,2,3,4,5γ=20, ω0=0.5298 305 313 322 332 337 343PLUMED 1.3, GMX4.5
PT-WT-MetaDRUN18L99AE365 nsCHARMM22*1γ=20, ω0=0.5298 302 306 311 317 323 330 338PLUMED 1.3, GMX4.5
plain MDRUN19L99AG400 nsCHARMM22*298GMX5.0
plain MDRUN20L99AE400 nsCHARMM22*298GMX5.0
plain MDRUN21L99AG400 nsAMBER99sb*ILDN298GMX5.0
plain MDRUN22L99AE400 nsAMBER99sb*ILDN298GMX5.0
plain MDRUN23L99A,G113A,R119PG400 nsCHARMM22*298GMX5.0
plain MDRUN24L99A,G113A,R119PE400 nsCHARMM22*298GMX5.0
plain MDRUN25L99A,G113A,R119PG400 nsAMBER99sb*ILDN298GMX5.0
plain MDRUN26L99A,G113A,R119PE400 nsAMBER99sb*ILDN298GMX5.0
InMetaDRUN27-68§§L99AGCHARMM22*Spath,Zpathγ=20, τG=80 ps, ω0=1.0, σS=0.016¶¶, σZ=0.002##298PLUMED2.1, GMX5
InMetaDRUN69-104***L99AECHARMM22*Spath,Zpathγ=20, τG=80 ps, ω0=1.0, σS=0.016, σZ=0.002298PLUMED2.1, GMX5
InMetaDRUN105-119L99A,G113A,R119PGCHARMM22*Spath,Zpathγ=15, τG=100ps, ω0=0.8, σS=0.016¶¶, σZ=0.002##298PLUMED2.2, GMX5.1.2
InMetaDRUN120-134L99A,G113A,R119PECHARMM22*Spath,Zpathγ=15, τG=100 ps, ω0=0.8, σS=0.016¶¶, σZ=0.002##298PLUMED2.2, GMX5.1.2
ABMDRUN135-154L99A-BenzeneboundCHARMM22*RMSDbenzeneK=50 KJ · mol−1 · nm−2, Starget=4.0 nm298PLUMED2.2, GMX5.1.2
ABMDRUN155-174L99A-BenzeneboundCHARMM22*RMSDbenzeneK=20 KJ · mol−1 · nm−2, Starget=4.0 nm298PLUMED2.2, GMX5.1.2
  1. γ is the bias factor.

  2. τD is the characteristic decay time used for the dynamically-adapted Gaussian potential.

  3. §τG is the deposition frequency of Gaussian potential.

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

  5. # N is the number of replicas.

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

  7. ††ΔT means RUN_FREQ and Δt means STORE_FREQ. Other parameters: BASIN_TOL=0.2, EXPAND_PARAM=0.3, INITIAL_SIZE=3.0.

  8. ‡‡ Replica at 298 K is the neutral replica without energy bias.

  9. §§42 trajectories collected for G-to-E transitions.

  10. ¶¶σS is the Gaussian width of Spath.

  11. ##σZ is the Gaussian width of Zpath (in nm2).

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

Appendix 1—table 2

Definition of collective variables.

1total energybin=500enhance energy fluctuations
2dihedral angle of Cα atoms of consecutive residues F104-Q105-M106-G107σ=0.1
3dihedral angle of Cα atoms of consecutive residues G113-F114-T115-N116σ=0.1
4QG, distance in contact map space to the GXray structureσ=0.5
5QE, distance in contact map space to the EROSETTA structureσ=0.5
6distance between QG and QEσ=0.5
7number of backbone hydrogen bonds formed between M102 and G107σ=0.1
8dihedral correlation between the Cα dihedral angles of consecutive residues in segment N101-G107σ=0.1
9global RMSD to the whole proteinwall potentialavoid sampling unfolding space
Appendix 1—table 3

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

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

Minimum probe radius0.9 Å
Shell depth4
Shell radius3
Clustering threshold3.5
Starting point optimization
Maximum distance3 Å
Desired radius5 Å
Appendix 1—table 5

Clustering analysis of tunnels (top three listed).

IndexPopulationMaximal bottleneck radius (Å)Average bottleneck radius (Å)
Appendix 1—table 6

Unbinding Pathways Explored by ABMD (RMSDBNZ as CV).

k=20 kJ/(mol · nm−2)k=50 kJ/(mol · nm−2)
RUN156 nsP127 nsP2
RUN236 nsP278 nsP1
RUN343 nsP16 nsP1
RUN443 nsP135 nsP1
RUN577 nsP210 nsP1
RUN6176 nsP144 nsP1
RUN741 nsP118 nsP1
RUN8106 nsP115 nsP1
RUN972 nsP17 nsP1
RUN10107 nsP12 nsP1
RUN1161 nsP120 nsP2
RUN1258 nsP226 nsP1
RUN1364 nsP131 nsP2
RUN14173 nsP220 nsP1
RUN15172 nsP134 nsP1
RUN1674 nsP222 nsP1
RUN1720 nsP117 nsP1
RUN1834 nsP135 nsP2
RUN1991 nsP121 nsP2
RUN2061 nsP118 nsP1
Cost1.6 μs0.5 μs
P175% (15/20)75% (15/20)
P225% (5/20)25% (5/20)

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)

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

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