Identification and Validation of Reaction Coordinates Describing

Subscriber access provided by UNIVERSITY OF THE SUNSHINE COAST

Article

Identification and validation of reaction coordinates describing protein functional motion: Hierarchical dynamics of T4 Lysozyme Matthias Ernst, Steffen Wolf, and Gerhard Stock J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00571 • Publication Date (Web): 15 Sep 2017 Downloaded from http://pubs.acs.org on September 15, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 14

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment



Page 2 of 14

Page 3 of 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tory and numerous targeted MD simulations, we here pursue various strategies of dimensionality reduction to construct reaction coordinates that reveal the molecular mechanism underlying the hinge-bending motion of T4L. To this end we (i) perform various principal component analyses to determine important coordinates (and associated residues and motions), and (ii) employ targeted MD pulling simulations to select a suitable twodimensional reaction coordinate. Apart from the obvious open-closed transition, we thus identify a locking transition where the side chain of Phe4 changes from a solvent-exposed to a hydrophobically-buried state. This mechanism is found to stabilize the open and closed states of T4L and thereby causes their relatively long life time of ∼ 10 µs. In extension of the usual two-state picture, we suggest a four-state model of the functional motion of T4L, which also accounts for the locking transition. Performing a mean first passage time analysis, we find a hierarchical coupling of the fast (∼ 100 ns) opening-closing motion and the slow (∼ 10 µs) locking transition. The significance of the key residues of the mechanism (Phe4, Phe67 and Phe104) is tested by additional MD simulations of in total six mutants, which are found to considerably modify the energy landscape and the overall kinetics of T4L.

2 2.1

Methods MD Simulations

All simulations were performed using the Gromacs package (version 4.6.7), 49 employing the Amber ff99sb*ILDN force field 50–52 and the TIP3P water model. 53 Following earlier work by Hub and de Groot, 24 we considered the M6I mutant of T4L (PDB 150L, 34 chain D; residues 163-164 omitted as they are not resolved in the crystal structure). Employing a triclinic box with a NaCl salt concentration of 150 mmol L−1 resulted in ∼29 400 atoms in total. Use of the LINCS algorithm 54 to constrain bonds including hydrogen atoms allowed for an integration time step of 2 fs. Particle-Mesh Ewald summation 55 (PME) was used for the calculation of electrostatic interactions. All cut-offs (neighbor search, Verlet scheme, Lennard-Jones interactions, and the real space grid of PME) were set to 1.2 A. Following steepest descent energy minimization of T4L in vacuo to remove sterically unfavorable interactions, the solvation box was constructed and energy minimization in solvent was performed. Equilibration included a 100 ps NVT run using position restraints and the Bussi thermostat 56 at 300 K, followed by a 1 ns NPT run using position restraints and the Berendsen barostat, 57 a 5 ns free NPT simulation (of which the last 4 ns were used to calculate the averaged box volume), and a 10 ns free NVT simulation. Finally an NVT production run with a total length of 50 µs at 300 K was performed, saving atom coordinates every 1 ps. Visual inspection of the trajectories was carried out with VMD. 58

Mutations were introduced into the structure of lowest free energy of the closed state during the first 30 µs of the equilibrium run using PyMOL. 59 The simulation protocol was the same as above with a production run length of 10 µs per mutation.

2.2

Principal (PCA)

component

analysis

In a PCA, the correlated internal motion of a system with N degrees of freedom r = (r1 , . . . , rN )T is described by the covariance matrix σmn = h(rm − hrm i) (rn − hrn i)i ,

(2)

where h. . . i represents the average over all sampled conformations. Diagonalization of this covariance matrix results in N eigenvectors {v(i) } and eigenvalues {λi }, which describe the modes of the collective motion and their respective amplitudes. The principal components (PCs) are then given by the projections of the coordinates r onto the eigenvectors xi = v(i) · r.

(3)

Considering the first few PCs with highest eigenvalues, we may construct a reaction coordinate x = (x1 , . . . , xd )T , which accounts for a large part of the system’s fluctuations. 13,14 Rather than performing a PCA based on the covariance [Eq. (2)] which points out coordinates with high variance, it may be advantageous to consider the correlation (i.e., the normalized covariance) which emphasizes correlated motion. To this end, we normalize each coordinate by its standard deviation, r˜n = rn /σrn , calculate the resulting covariance matrix {˜ σmn } (which is equal to the correlation matrix of {rn }) ˜ (i) , and obtain the corresponding and its eigenvectors v (i) ˜ · ˜r. PCs via x ˜i = v While this strategy appears straightforward, several issues need to be considered to obtain useful reaction coordinates via a PCA. First, the analysis depends crucially on the input coordinates used in the dimensionality reduction method. While Cartesian coordinates are convenient to handle, a Cartesian coordinate PCA is known to break down in the case of large-amplitude motion (as occurring, e.g., in a folding process), since structural dynamics of flexible molecules necessarily results in a mixing of overall and internal motion. 60,61 To circumvent this problem, internal coordinates such as (φ, ψ) backbone dihedral angles 62–64 or distances and contacts between atoms and protein residues 65–70 may be used. While backbone dihedral angle PCA has proven quite powerful to model the folding dynamics of peptides, RNA and small proteins, 63,64 distance-based PCAs also take into account the structure and dynamics of side chains and may therefore be advantageous to describe the functional dynamics of larger proteins. In both cases, it has been found important to perform a preselection of the coordinates to be included in

ACS Paragon Plus Environment 3


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the analysis. That is, coordinates reflecting irrelevant motion (e.g., uncorrelated motion of flexible terminal residues) or hardly any motion (such as essentially rigid secondary structures of the protein core) should be excluded in order to minimize the noise in the analysis. 2.2.1

2.2.2

threshold) of all identified contacts for all frames of the trajectory. This data set of distances is subsequently used as input for the contact PCA according to Eq. (2).

2.3

Dihedral angle PCA

To perform a PCA on circular variables such as angles φ, we may change to sine/cosine-transformed coordinates (r1 = cos φ, r2 = sin φ) to obtain a linear coordinate space with the usual Euclidean distance as induced metric. 63 Alternatively, it has recently been suggested 71 to directly calculate the covariance matrix of angular variables by (i) associating the circular distance as the inner arc between two angles on the unit circle and (ii) defining the circular mean hφi over various observations {φn } as the projection onto the unit circle of the average PN of the sine and cosine projections (x = N1 n=1 sin φn , P N y = N1 n=1 cos φn ) as obtained by the atan2 function, hφi = atan2(y, x). When we project the data φ onto the resulting eigenvectors v(i) according to Eq. (3), we again need to account for the periodicity of the data. To this end, we introduce a cut of the periodic space (say, between φmin = −π and φmax = π) and shift all angles (φ → φ + δ) such that populated sections of the angle distribution remain connected. In the case of backbone dihedral angles (φ, ψ), the cut is naturally introduced at the (hardly crossed) barriers at φ ≈ 0° and ψ ≈ 20°. See ref. 71 for details on this recently proposed method termed dPCA+. Contact PCA

Numerous options exist to perform a PCA based on interatomic distances or contacts. 65–69 Following our recent work, 70 here we defined a contact to be formed if the distance Dij between the closest non-hydrogen atoms of two residues i and j is shorter than 4.5 A, Dij = min(|r i,k − r j,l |) ≤ 4.5 A,

(4)

with the indices k and l running over all heavy atoms of the selected residue pair. In addition, we discard contacts between residues less than four residues apart in sequence, effectively omitting short-range contacts as, e.g., in helical structure elements or turns. Employing this definition, first the contacts of T4L need to be determined from a suitable reference structure. As neither the open nor the closed state of T4L contains all relevant contacts, we considered the contacts of both states, using the (energy minimized) crystal structure for the open state and the MD structure with the lowest radius of gyration for the closed state. The resulting contact matrix shown in Fig. 1c reveals in total 402 contacts, including stabilizing hydrogen bonds of the α and β secondary structures and numerous tertiary contacts. Using the MDAnalysis framework, 72 we calculated the contact distance Dij (now without

Page 4 of 14

Targeted (TMD)

Molecular

Dynamics

To test if some chosen reaction coordinate x accounts for the mechanism of the open-closed transition in T4L, we studied the molecule’s response to external pulling along x. To this end, we used TMD simulations 31–33 which constrain the velocity in direction of x and provide an easy means to calculate the free energy profile ∆G(x). Restricting ourselves to a 1D reaction coordinate x, we drive the system from x(0) = x0 to x(tend ) = xend by constraining coordinate x(t) to xc (t) = x0 + vc t with a constant velocity vc via the constraint function 2

Φ(x(t)) = (x(t) − xc (t)) = 0.

(5)

The constraint is realized via the force Fc = λ

dΦ(x(t)) = 2λ (x(t) − xc (t)) , dx

(6)

where λ is a Lagrange parameter. In the present application, the pulling coordinate was chosen as the distance between the centers of mass of two atom groups of T4L. To account for different masses of the individual atoms, p the force acting on atom i was scaled by a factor mi /mav , where mi and mav represent the mass of atom i and the average atom mass of the individual pull group, respectively. The TMD simulations were performed employing the PULL code as implemented in GROMACS 4.6.7, 49 using the “constraint” mode and ns−1 for pulling a constraint velocity of vc = 0.125 A −1 whole domains and vc = 0.5 A ns for pulling selected side chains. To calculate the free energy change ∆GTMD = G(xend ) − G(x0 ) from the TMD simulations, we employed a thermodynamic integration scheme 1 Z xend Z xend dG hFc ix dx dx = ∆GTMD = dx x0 x0 ≈

N X

hFc ixi ∆x,

(7)

i=1

where hFc ix represents an equilibrium time average of the constraining force at point x. The averages were calculated from 200 ns long TMD runs at fixed position xi (i.e., for vc = 0), using only the last 20 ns to evaluate hFc ixi . Typically N = 16 equidistant points were used to approximate the integral.

2.4

Mean first passage time analysis

Considering a system propagating in a discrete state space, the average transition time τ between two states i and j can be estimated via the mean first passage time


Page 5 of 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(MFPT). 73 Unlike to τ calculated by counting transitions within a certain lag time, the MFPT is more general as it avoids any Markovian assumption. In the simplest case of a system with two states (say, 0 and 1) and a trajectory with a single transition event 0 → 1, which spends n01 time steps δt in initial state 0 before it jumps to state 1, the MFPT is defined as 73 mfpt τ0→1 =

n01 1 X n01 + 1 δt. l δt = n01 2

(8)

l=1

That is, the MFPT represents the sum of all first passage times divided by the number of time steps spent in the initial state. For a two-state system, the MFPT is therefore just half of the average time waited in the state. We now generalize to the case of more than two states (i = 1, 2, . . .) and consider a trajectory that exhibits mul(k) tiple (m > 1) transition events i → j. Defining nij as the number of time steps spent in initial state i occurring in event k, the total number of time steps the system spends in state i before it jumps to state j is Pm (k′ ) nij = k′ =1 nij . Hence the MFPT amounts to n

mfpt τi→j

(k)

m ij 1 XX l δt. = nij

(9)

k=1 l=1

The outcome of a MFPT analysis depends crucially on the partitioning of the continuous MD data into discrete states. Using appropriate reaction coordinates (e.g., p and x defined in section 3.6), to this end we first applied density-based geometric clustering 28 which yields well defined microstates separated by local free energy barriers. As the projection on a low-dimensional space may induce spurious transitions in the vicinity of energy barriers, in a second step we identify core regions of the microstates and count transitions only if the core region of the other state is reached. 30,74 Effectively, this procedure generates a state trajectory i(t) with clear-cut state boundaries, on the basis of which a appropriate MFPT analysis can be carried out. We used an elliptical shape for the core regions (Fig. S2), which was chosen to match the shape of the minima of the free energy landscape (Fig. 6a). Varying the size of the core regions resulted in no significant changes of the MFPT, see Table S1.

3 3.1

Results and Discussion 1D order parameters can provide a qualitative picture of protein motion

To give an overview of the functional dynamics of T4L, we first consider the radius of gyration RG as a commonly used one-dimensional (1D) order parameter. Figure 2a displays the time evolution of RG , obtained

from the 50 µs MD trajectory. Despite large fluctuations, RG (t) clearly reveals several prominent openclosed transitions of T4L, e.g., at t ≈ 1, 2, 3 and 7 µs. On average, we find that the open state exists for about twice the time of the closed state. In agreement with experimental investigations, 43 we furthermore find that both states exists on a time scale of 5 to 20 µs. Other choices of 1D order parameters, such as the root mean square distance (RMSD) to the (open or closed) reference structure gave similar results (data not shown). While the radius of gyration provides a qualitative picture of the functional dynamics of T4L (such as the overall amplitude and the time scale), it may not tell much about the molecular mechanism causing the process. This is indicated by the free energy curve along RG (right panel of Fig. 2a), showing two minima associated with the open and closed states, that are separated by an energy barrier ∆GRG . 1kB T . Modeling the transition rate by the standard expression 73 k = 1/τ = k0 e−∆GRG /kB T ,

(10)

it is obvious that a small energy barrier ∆GRG ≈ 1 kB T cannot account for the long (∼ 10 µs) observed transition time τ of T4L. That is, the rate is almost solely caused by the prefactor k0 , which depends critically on the diffusion coefficients of the system. 75 Rather we find that the projection of the many-dimensional molecular motion onto the 1D coordinate causes large fluctuations of RG , which result in a small barrier. Hence the radius of gyration reflects the consequences rather than the origin of the functional motion of T4L. Although it is possible to further analyze 1D order parameters by identifying all atoms involved in its time evolution (e.g., via functional mode analysis 24,25 ), it is therefore not clear to what extent such an analysis sheds light on the underlying mechanism. For a complex system as such as T4L, 1D order parameters are apparently not adequate to uncover microscopic reaction pathways.

3.2

Cartesian PCA may not work for multi-event trajectories

As discussed in the Introduction, dimensionality reduction methods such as principal component analysis 13,14 (PCA) or independent component analysis 15,17,18 may provide an approach to identify suitable reaction coordinates in a systematic manner. With this end in mind, we first consider a standard PCA performed on the Cartesian coordinates of the backbone atoms of T4L. The time evolution of the resulting first PC, x1 (t) shown in Fig. 2b is overall quite similar to the evolution of the radius of gyration (Fig. 2a). This is in line with the observation that similar atoms are prominently involved in the motion along x1 (t) and along RG (t). Exhibiting less fluctuations than the latter, x1 (t) results in a doublewell free energy curve ∆G(x1 ) with a barrier of about 2 kB T (Fig. Fig. 2b). As shown in Fig. S3, however, the free energy curves ∆G(xi ) of all higher PCs are struc-




Page 6 of 14

Page 7 of 14





Page 8 of 14

Page 9 of 14





Page 10 of 14

Page 11 of 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in the locking coordinate, reflecting the many possible conformations of Phe4 in the free state. In intermediate state 2 the mouth is closed, but the side chain of Phe4 still is buried in the hydrophobic cavity. Hence there is strain in the system which drives it back to the open state, if the transition in p cannot be completed. In intermediate state 3, on the other hand, the side chain of Phe4 is found outside of the cavity and therefore in principle is free to move. While this also holds for state 4 and although these states are connected by only a minor barrier, there are nonetheless distinct structural differences between 3 and 4. They manifest themselves in terms of the alignment of the Phe4 side chain, which needs to find a position on the protein surface that marks the entrance into the hydrophobic cavity (state 3), but frequently is sterically blocked by Arg8 (state 4). As displayed in Fig. S17, the phenyl ring of Phe4 tries to reduce its nonpolar solvent accessible surface area in both states, by sticking to a hydrophobic patch at the outside of the protein and bridging a gap between Asn2, Ile3, Phe67, Asn68 and Val71. In state 4, Arg8 is found over the entrance to the hydrophobic binding pocket, forming a partially bidentate salt bridge with Glu64. In state 3, the side chain of Phe4 is positioned over the entrance, while Arg8 is displaced and the salt bridge with Glu64 perturbed, not allowing the highly stable bidentate connection. Besides this local change in connectivity, state 3 also exhibits a lower total number of intraprotein hydrogen bonds than state 4. Table 1: Mean first passage times (in units of µs) of the four-state model of T4L. to\from

1

2

3

4

1 2 3 4

0.7 5.4 6.5

0.07 3.2 4.8

0.9 0.5 0.04

1.2 0.7 0.02 -

The free energy landscape ∆G(x, p) indicates that transitions between the main states 1 and 4 require subsequent closing and unlocking of T4L. To elucidate this process, Fig. S18 shows the time evolution of trajectory pieces during and shortly before/after the 1↔4 transitions. As expected from the appearance of the energy landscape, we find that most transitions proceed along the route 1→ 2→3→4 and back along 4→3→2→1. To estimate the time scales associated with the individual steps of the process, we performed a mean first passage time (MFPT) analysis of the four-state system (see Methods). The resulting MFPT are comprised in Table 1, and were shown to be quite robust against variation of the cores used for assigning the states (Table S1). The next-neighbor MFPTs are also included in the illustrative scheme in Fig. 6b, which provides the basis for the following discussion. Let us first consider the 1→2→ 3→4 pathway. We

find a MFPT of 6.5 µs for the overall 1→4 transition, which corresponds to an average waiting time of 13.0 µs in state 1 (see Methods). As expected from the discussion above, the locking transition 2→3 with a MFPT of 3.2 µs represents the slowest step of the route. This is about one and two orders of magnitude slower than the other two steps 1→2 and 3→4, respectively. For the other direction of the process, 4→3→ 2→1, we find a significantly shorter overall MFPT of 1.2 µs, which is mostly caused by the relatively fast (0.5 µs) unlocking step 3→2. Out of in total 16 transitions between states 1 and 4, eleven follow directly the proposed “L”-shaped route 1↔2↔3↔4, while the remaining five transitions follow this route closely but may leave out states 2 or 3 (cf. Fig. S18). Hence we find that the relatively long MFPT of 6.5 µs for the overall 1→4 transition is caused by a hierarchical process that consists of fast opening and slow locking motion. Starting in the open state 1, the system attempts every ∼ 1 ns to close its mouth, but –due to frequent (0.07 µs) back transitions– requires a MFPT of 0.7 µs to actually accomplish the transition. This relatively fast process is the prerequisite for the subsequent slow locking step 2→3, before the system rapidly relaxes to the final state 4.

3.8

Mutation studies

To test our hypothesis that Phe4 is a key residue for the functional dynamics of T4L, we performed additional MD simulations of various mutations of our reference system M6I (see Methods). Apart from Phe4, we chose to also mutate residues Phe67 and Phe104. This is because the side chain of Phe67 forms a steric barrier that Phe4 has to overcome when moving in or out of the hydrophobic pocket. By mutating one of these two side chains, we intented to directly modify the main barrier in the opening-closing process. Phe104 was chosen because it was found to show correlated motion with Phe4 and Phe67 38 and because it exhibits major contributions to the first contact PCA eigenvector. Moreover, it is part of an inner hydrophobic core (consisting of Leu7, Glu11, Ile29, Asp70, Ala71, Ala74, Ile100, Asn101, Val103 and sometimes Asp10, Phe67 or Gln105) and is solvent accessible at the inside of the mouth, probably with the task of properly aligning the substrate. We chose mutations Phe→Ala to change the large phenyl side chain to a small hydrophobic one, as well as Phe→Asn to change from a hydrophobic to a hydrophilic side chain of similar size. In total we carried out six MD runs with a trajectory length of 10 µs each, referred to as F4A, F4N, F67A, F67N, F104A and F104N. Showing the time evolution of the radius of gyration and the corresponding free energy profiles of all mutations, Fig. 7 reveals that the various mutations may significantly change the energy landscape and the conformational dynamics of the system. As a thermodynamic effect, we find changes of the positions and values of




Page 12 of 14

Page 13 of 14

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


to the first few contact PCA eigenvectors, we have tested the response of T4L to pulling along various other choices of the reaction coordinate. In this way, we have discovered a prominent role of residue Phe4, which is located at the jaw joint of T4L. It occurs in two states (Fig. 5): In the open conformation, the side chain of Phe4 is buried inside a hydrophobic cavity which restricts its motion (locked state), while in the closed conformation, the side chain is outside of the cavity and exposed to the solvent (free state). By pulling apart Phe4 and the residues forming the lock cavity, TMD simulations have been shown to successfully induce the openclosed transition of T4L (Fig. 4). That is, although this transition is not enforced per se, it happens as a consequence of enforcing Phe4 to bypass Phe67. The locking mechanism explains the origin of the main barrier of the process, which is caused by steric hindrance due to the large size of the phenyl rings. Moreover, the mechanism results in an entropic stabilization of both open and closed states of T4L. In excellent agreement with recent experimental data, 43 we obtain a relatively long life time of ∼ 10 µs for both states. To verify the significance of the key residues of the mechanism, we have performed additional MD simulations of in total six mutants, which were found to considerably modify the energy landscape and the overall kinetics of T4L (Fig. 7). In extension of the usual two-state picture of the open-close transition, the above TMD/PCA study suggests a four-state model (Fig. 6) which also accounts for the locking transition. Performing a mean first passage time analysis of the various steps of the process, we have found a hierarchical coupling of the fast openingclosing motion and the slow locking transition. That is, starting in the open state, the system attempts every ∼ 1 ns to close its mouth, but (due to frequent back transitions) requires hundreds of nanoseconds to accomplish this transition. This relatively fast process is the prerequisite for the subsequent slow locking step, which occurs on a microsecond timescale. While the concept of hierarchical protein energy landscapes has been discussed for some time, 21–23 this study has made a first attempt to concretely specify the various subprocesses, timescales and hierarchical couplings for a two-domain protein. Acknowledgement We thank Florian Sittel and Mithun Biswas for numerous instructive and helpful discussions. Furthermore, we gratefully acknowledge computing resources provided by the Black Forest Grid initiative at the University of Freiburg, as well as the bwGRiD, bwUniCluster and bwForCluster initiatives of the State of Baden-W¨ urttemberg. This work has been supported by the Deutsche Forschungsgemeinschaft. Supporting Information Available: Details of various PCAs and MFPT analyses (Figs. S1-S10), additional TMD simulations (Figs. S11-S15), SASA (Fig. S16), structural details of states 3 and 4 (Fig. S17), transition pathways in the four-state landscape (Fig.

S18) and comparison of mutation studies with free energy landscapes of conPCA and four-state model of equilibrium simulation (Fig. S18). This material is available free of charge via the Internet at http: //pubs.acs.org/.

References (1) Berendsen, H. J. C. Simulating the Physical World; Cambridge University Press: Cambridge, 2007. (2) Onuchic, J. N.; Schulten, Z. L.; Wolynes, P. G. Theory of protein folding: The energy landscape perspective. Annu. Rev. Phys. Chem. 1997, 48, 545–600. (3) Dill, K. A.; Chan, H. S. From Levinthal to Pathways to Funnels: The ”New View” of Protein Folding Kinetics. Nat. Struct. Bio. 1997, 4, 10–19. (4) Wales, D. J. Energy Landscapes; Cambridge University Press: Cambridge, 2003. (5) Best, R. B.; Hummer, G. Coordimate-dependent diffusion in protein folding. Proc. Natl. Acad. Sci. USA 2010, 107, 1088 – 1093. (6) Krivov, S. V. On Reaction Coordinate Optimality. J. Chem. Theory Comput. 2013, 9, 135–146. (7) Krivov, S. V.; Karplus, M. Hidden complexity of free energy surfaces for peptide (protein) folding. Proc. Natl. Acad. Sci. USA 2004, 101, 14766–14770. (8) Altis, A.; Otten, M.; Nguyen, P. H.; Hegger, R.; Stock, G. Construction of the free energy landscape of biomolecules via dihedral angle principal component analysis. J. Chem. Phys. 2008, 128, 245102. (9) Maisuradze, G. G.; Liwo, A.; Scheraga, H. A. How Adequate are One- and Two-Dimensional Free Energy Landscapes for Protein Folding Dynamics? Phys. Rev. Lett. 2009, 102, 238102. (10) Rohrdanz, M. A.; Zheng, W.; Clementi, C. Discovering Mountain Passes via Torchlight: Methods for the Definition of Reaction Coordinates and Pathways in Complex Macromolecular Reactions. Annu. Rev. Phys. Chem. 2013, 64, 295–316. (11) Noe, F.; Clementi, C. Collective variables for the study of longtime kinetics from molecular trajectories: theory and methods. Curr. Opin. Struct. Biol. 2017, 43, 141 – 147. (12) Das, P.; Moll, M.; Stamati, H.; Kavraki, L. E.; Clementi, C. Low-dimensional, free-energy landscapes of protein-folding reactions by nonlinear dimensionality reduction. Proc. Natl. Acad. Sci. USA 2006, 103, 9885–9890. (13) Jolliffe, I. T. Principal Component Analysis; Springer: New York, 2002. (14) Amadei, A.; Linssen, A. B. M.; Berendsen, H. J. C. Essential dynamics of proteins. Proteins 1993, 17, 412–425. (15) Hyv¨ arinen, A.; Karhunen, J.; Oja, E. Independent Component Analysis; John Wiley & Sons: New York, 2001. (16) Lange, O. F.; Grubm¨ uller, H. Generalized Correlation for Biomolecular Dynamics. Proteins 2006, 62, 1053–1061. (17) Lange, O. F.; Grubm¨ uller, H. Full correlation analysis of conformational protein dynamics. Proteins 2008, 70, 1294 – 1312. (18) Perez-Hernandez, G.; Paul, F.; Giorgino, T.; De Fabritiis, G.; Noe, F. Identification of slow molecular order parameters for Markov model construction. J. Chem. Phys. 2013, 139, 015102. (19) Schaudinnus, N.; Bastian, B.; Hegger, R.; Stock, G. Multidimensional Langevin modeling of nonoverdamped dynamics. Phys. Rev. Lett. 2015, 115, 050602. (20) Chipot, C.; Pohorille, A. Free Energy Calculations; Springer: Berlin, 2007. (21) Frauenfelder, H.; Sligar, S.; Wolynes, P. The energy landscapes and motions of proteins. Science 1991, 254, 1598–1603. (22) Metzler, R.; Klafter, J.; Jortner, J. Hierarchies and logarithmic oscillations in the temporal relaxation patterns of proteins and other complex systems. Proc. Natl. Acad. Sci. USA 1999, 96, 11085–11089. (23) Buchenberg, S.; Schaudinnus, N.; Stock, G. Hierarchical biomolecular dynamics: Picosecond hydrogen bonding regulates microsecond conformational transitions. J. Chem. Theo. Comp. 2015, 11, 1330–1336. (24) Hub, J. S.; de Groot, B. L. Detection of functional modes in protein dynamics. PLoS Comput. Biol. 2009, 5, e1000480. (25) Krivobokova, T.; Briones, R.; Hub, J. S.; Munk, A.; de Groot, B. L. Partial Least-Squares Functional Mode Analysis: Application to the Membrane Proteins AQP1, Aqy1, and CLC-ec1. Biophys. J. 2012, 103, 786–796.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(26) Keller, B.; Daura, X.; van Gunsteren, W. F. Comparing geometric and kinetic cluster algorithms for molecular simulation data. J. Chem. Phys. 2010, 132, 074110. (27) Rodriguez, A.; Laio, A. Clustering by fast search and find of density peaks. Science 2014, 344, 1492–1496. (28) Sittel, F.; Stock, G. Robust Density-Based Clustering to Identify Metastable Conformational States of Proteins. J. Chem. Theo. Comp. 2016, 12, 2426–2435. (29) Prinz, J.-H.; Wu, H.; Sarich, M.; Keller, B.; Senne, M.; Held, M.; Chodera, J. D.; Sch¨ utte, C.; Noe, F. Markov models of molecular kinetics: generation and validation. J. Chem. Phys. 2011, 134, 174105. (30) Jain, A.; Stock, G. Identifying metastable states of folding proteins. J. Chem. Theo. Comp. 2012, 8, 3810 – 3819. (31) Schlitter, J.; Engels, M.; Kr¨ uger, P.; Jacoby, E.; Wollmer, A. Targeted Molecular Dynamics Simulation of Conformational Change-Application to the T ↔ R Transition in Insulin. Mol. Simul. 1993, 10, 291–308. (32) Schlitter, J.; Engels, M.; Kr¨ uger, P. Targeted Molecular Dynamics - A New Approach for Searching Pathways of Conformational Transitions. J. Mol Graph. 1994, 12, 84–89. (33) Schlitter, J.; Swegat, W.; M¨ ulders, T. Distance-type reaction coordinates for modelling activated processes. J. Mol. Model. 2001, 7, 171–177. (34) Zhang, X.-J.; Matthews, B. Conservation of solvent-binding sites in 10 crystal forms of T4 lysozyme. Protein Sci. 1994, 3, 1031–1039. (35) Krichevsky, O. T4 lysozyme as a Pac-Man: how fast can it chew? Biophys. J. 2012, 103, 1414–1415. (36) Weaver, L.; Matthews, B. Structure of bacteriophage T4 lysozyme refined at 1.7 ˚ A resolution. J. Mol. Biol. 1987, 193, 189 – 199. (37) Remington, S.; Anderson, W.; Owen, J.; Eyck, L.; Grainger, C.; Matthews, B. Structure of the lysozyme from bacteriophage T4: An electron density map at 2.4 A resolution. J. Mol. Biol. 1978, 118, 81 – 98. (38) Dixon, M.; Nicholson, H.; Shewchuk, L.; Baase, W.; Matthews, B. Structure of a hinge-bending bacteriophage T4 lysozyme mutant, Ile3 → Pro. J. Mol. Biol. 1992, 227, 917 – 933. (39) Zhang, X.; Wozniak, J. A.; Matthews, B. W. Protein Flexibility and Adaptability Seen in 25 Crystal Forms of T4 Lysozyme. J. Mol. Biol. 1995, 250, 527 – 552. (40) Xu J,; Baase W A,; Baldwin E,; Matthews B W, The response of T4 lysozyme to large-to-small substitutions within the core and its relation to the hydrophobic effect. Protein Sci. 1998, 7, 158–177. (41) Cellitti, J.; Bernstein, R.; Marqusee, S. Exploring subdomain cooperativity in T4 lysozyme II: Uncovering the C-terminal subdomain as a hidden intermediate in the kinetic folding pathway. Protein Sci. 2007, 16, 852–862. (42) Baase Walter A,; Liu Lijun,; Tronrud Dale E,; Matthews Brian W, Lessons from the lysozyme of phage T4. Protein Sci. 2010, 19, 631–641. (43) Yirdaw Robel B,; Mchaourab Hassane S, Direct Observation of T4 Lysozyme Hinge-Bending Motion by Fluorescence Correlation Spectroscopy. Biophys. J. 2012, 103, 1525–1536. (44) Arnold, G. E.; Ornstein, R. L. Protein hinge bending as seen in molecular dynamics simulations of native and M6I mutant T4 lysozymes. Biopolymers 1997, 41, 533–544. (45) de Groot, B. L.; Hayward, S.; van Aalten, D. M.; Amadei, A.; Berendsen, H. J. Domain motions in bacteriophage T4 lysozyme: a comparison between molecular dynamics and crystallographic data. Proteins 1998, 31, 116–127. (46) Peng, J.; Zhang, Z. Simulating Large-Scale Conformational Changes of Proteins by Accelerating Collective Motions Obtained from Principal Component Analysis. J. Chem. Theory Comput. 2014, 10, 3449–3458. (47) Harada, R.; Kitao, A. Parallel cascade selection molecular dynamics (PaCS-MD) to generate conformational transition pathway. J. Chem. Phys. 2013, 139, 035103. (48) Zheng, W.; Glenn, P. Probing the folded state and mechanical unfolding pathways of T4 lysozyme using all-atom and coarsegrained molecular simulation. J. Chem. Phys. 2015, 142 . (49) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. Journal of Chemical Theory and Computation 2008, 4, 435–447. (50) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins 2006, 65, 712–725. (51) Best, R. B.; Hummer, G. Optimized Molecular Dynamics Force Fields Applied to the Helix-Coil Transition of Polypeptides. J. Phys. Chem. B 2009, 113, 9004–9015.

Page 14 of 14

(52) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins 2010, 78, 1950–1958. (53) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926. (54) Hess, B. P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 116– 122, PMID: 26619985. (55) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577–8593. (56) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 0141011–0141017. (57) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Dinola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684. (58) Humphrey, W.; Dalke, A.; Schulten, K. VMD – Visual Molecular Dynamics. J. Mol. Graph. 1996, 14, 33–38. (59) Schr¨ odinger, LLC, The PyMOL Molecular Graphics System, Version 1.8. 2015, (60) Mu, Y.; Nguyen, P. H.; Stock, G. Energy Landscape of a Small Peptide Revealed by Dihedral Angle Principal Component Analysis. Proteins 2005, 58, 45 – 52. (61) Sittel, F.; Jain, A.; Stock, G. Principal component analysis of molecular dynamics: On the use of Cartesian vs. internal coordinates. J. Chem. Phys. 2014, 141, 014111. (62) van Aalten, D. M. D.; de Groot, B. L.; Finday, J. B. C.; Berendsen, H. J. C.; Amadei, A. A comparison of techniques for calculating protein essential dynamics. J. Comput. Chem. 1997, 18, 169–181. (63) Altis, A.; Nguyen, P. H.; Hegger, R.; Stock, G. Dihedral angle principal component analysis of molecular dynamics simulations. J. Chem. Phys. 2007, 126, 244111. (64) Riccardi, L.; Nguyen, P. H.; Stock, G. Free energy landscape of an RNA hairpin constructed via dihedral angle principal component analysis. J. Phys. Chem. B 2009, 113, 16660 – 16668. (65) Abseher, R.; Nilges, M. Are There Non-trivial Dynamic Crosscorrelations in Proteins? J. Mol. Biol. 1998, 279, 911. (66) L¨ atzer, J.; Shen, T.; Wolynes, P. G. Conformational Switching upon Phosphorylation: A Predictive Framework Based on Energy Landscape Principles. Biochem. 2008, 47, 2110–2122. (67) Hori, N.; Chikenji, G.; Berry, R. S.; Takada, S. Folding energy landscape and network dynamics of small globular proteins. Proc. Natl. Acad. Sci. USA 2009, 106, 73–78. (68) Allen, L. R.; Krivov, S. V.; Paci, E. Analysis of the Free-Energy Surface of Proteins from Reversible Folding Simulations. PLoS Comput. Biol. 2009, 5, e1000428. (69) Kalgin, I. V.; Caflisch, A.; Chekmarev, S. F.; Karplus, M. New Insights into the Folding of a beta-Sheet Miniprotein in a Reduced Space of Collective Hydrogen Bond Variables: Application to a Hydrodynamic Analysis of the Folding Flow. J. Phys. Chem. B 2013, 117, 6092–6105. (70) Ernst, M.; Sittel, F.; Stock, G. Contact- and distance-based principal component analysis of protein dynamics. J. Chem. Phys. 2015, 143, 244114. (71) Sittel, F.; Filk, T.; Stock, G. Principal component analysis on a torus: Theory and application to protein dynamics. submitted 2017, (72) Michaud-Agrawal, N.; Denning, E. J.; Woolf, T. B.; Beckstein, O. MDAnalysis: A toolkit for the analysis of molecular dynamics simulations. J. Comput. Chem. 2011, 32, 2319–2327. (73) Zwanzig, R. Nonequilibrium Statistical Mechanics; Oxford University: Oxford, 2001. (74) Buchete, N.-V.; Hummer, G. Coarse master equations for peptide folding dynamics. J. Phys. Chem. B 2008, 112, 6057– 6069. (75) Hinczewski, M.; von Hansen, Y.; Dzubiella, J.; Netz, R. R. How the diffusivity profile reduces the arbitrariness of protein folding free energies. J. Chem. Phys. 2010, 132 . (76) McLachlan, A. D. Comparison of methods of matching protein structures. Acta. Cryst. A 1972, 28, 656. (77) Jain, A.; Stock, G. Hierarchical folding free energy landscape of HP35 revealed by most probable path clustering. J. Phys. Chem. B 2014, 118, 7750 – 7760. (78) Li, W.; Ma, A. Recent developments in methods for identifying reaction coordinates. Mol. Sim. 2014, 40, 784–793. (79) Peters, B. Reaction Coordinates and Mechanistic Hypothesis Tests. Annu. Rev. Phys. Chem. 2016, 67, 669–690.


Identification and Validation of Reaction Coordinates Describing

Recommend Documents