Abstract

The microtubule (MT)-associated motor protein kinesin utilizes its conserved ATPase head to achieve diverse motility characteristics. Despite considerable knowledge about how its ATPase activity and MT binding are coupled to the motility cycle, the atomic mechanism of the core events remain to be found. To obtain insights into the mechanism, we performed 38.5 microseconds of all-atom molecular dynamics simulations of kinesin-MT complexes in different nucleotide states. Local subdomain dynamics were found to be essential for nucleotide processing. Catalytic water molecules are dynamically organized by the switch domains of the nucleotide binding pocket while ATP is torsionally strained. Hydrolysis products are 'pulled' by switch-I, and a new ATP is 'captured' by a concerted motion of the α0/L5/switch-I trio. The dynamic and wet kinesin-MT interface is tuned for rapid interactions while maintaining specificity. The proposed mechanism provides the flexibility necessary for walking in the crowded cellular environment.

eLife digest

Motor proteins called kinesins perform a number of different roles inside cells, including transporting cargo and organizing filaments called microtubules to generate the force needed for a cell to divide. Kinesins move along the microtubules, with different kinesins moving in different ways: some ‘walk’, some jump, and some destroy the microtubule as they travel along it. All kinesins power their movements using the same molecule as fuel – adenosine triphosphate, known as ATP for short.

Energy stored in ATP is released by a chemical reaction known as hydrolysis, which uses water to break off specific parts of the ATP molecule. The site to which ATP binds in a kinesin has a similar structure to the ATP binding site of many other proteins that use ATP. However, little was known about the way in which kinesin uses ATP as a fuel, including how ATP binds to kinesin and is hydrolyzed, and how the products of hydrolysis are released. These events are used to power the motor protein.

Hwang et al. have used powerful computer simulation methods to examine in detail how ATP interacts with kinesin whilst moving across a microtubule. The simulations suggest that regions (or 'domains') of kinesin near the ATP binding site move around to help in processing ATP. These kinesin domains trap a nearby ATP molecule from the environment and help to deliver water molecules to ATP for hydrolysis. Hwang et al. also found that the domain motion subsequently helps in the release of the hydrolysis products by kinesin.

The domains around the ATP pocket vary among the kinesins and these differences may enable kinesins to fine-tune how they use ATP to move. Further investigations will help us understand why different kinesin families behave differently. They will also contribute to exploring how kinesin inhibitors might be used as anti-cancer drugs.

To understand the mechanisms of the different kinesins, it is important to search for and elucidate the conserved features of the motor head that are involved in the nucleotide processing events of the motility cycle; that is, ATP binding, hydrolysis, and product release (ADP and Pi, inorganic phosphate) (Figure 1A). One of the most intensely studied family member is Kinesin-1 (Kin-1; hereafter we refer Kin-1 as kinesin). It forms a dimer to walk toward the MT plus-end using one ATP per step (Figure 1B) (Svoboda et al., 1993; Block, 2007). It unbinds from the MT in the ADP state, and after making a step, it releases ADP and enters the nucleotide-free APO state with high MT-affinity (Cross, 2016; Hancock, 2016). Binding of an ATP triggers forward force generation (the ‘power stroke’; Figure 1A) by driving the cover strand (CS) and the neck linker (NL), which are located respectively on the N- and C-terminal ends of the motor head, to fold into a β-sheet named the cover-neck bundle (CNB; Figure 1C) (Rice et al., 1999; Hwang et al., 2008; Khalil et al., 2008). ATP hydrolysis completes a step (Milic et al., 2014; Andreasson et al., 2015).

Overview of kinesin structure and motility cycle.

(A) Diagram of the ATPase cycle of a motor head. Binding of an ATP triggers the conformational change to the post-stroke state. (B) Model of a kinesin dimer bound to the MT. The rear and front heads are in the post- and pre-stroke states, respectively. The neck linker (NL) connects the C-terminal end of the motor head to the α-helical neck stalk. (C,D) Comparison between (C) post- (Gigant et al., 2013) and (D) pre-stroke (Cao et al., 2014) states, defined based on the orientation of α⁢6 relative to α⁢4. CS: cover strand. αH12/βH12: C-terminal helices of α-/β-tubulins that form major contacts with kinesin, mainly with L11, α⁢4, L12, and β⁢5⁢a/b (also called L8). In the pre-stroke state, α⁢6 shortens and its C-terminal end connecting to the NL is positioned behind α⁢4 (red star). This is coupled with the rightward tilting and clockwise rotation about the vertical axis of the motor head (wide arrows). (E) The ATP binding pocket. MT is not shown. Kinesin structures are compared in Supplementary file 1. A complete list of kinesin domain names are in Figure 1—figure supplement 1.

The pre- and post-stroke states differ in the motor head orientation relative to the MT. In the pre-stroke state, the head tilts rightward relative to the MT plus-end direction, and rotates clockwise when viewed from above (wide arrows in Figure 1D). The head rotates in the opposite direction in the post-stroke state (Sindelar and Downing, 2010).

The nucleotide pocket on the rear-left side of the motor consists of the phosphate loop (P-loop), switch-I (sw-I), and switch-II (sw-II) (Figure 1E). These elements are conserved among different proteins including myosin and G-protein (Vale and Milligan, 2000). Compared to an isolated kinesin, a MT-bound kinesin has an at least 10-fold higher ATP hydrolysis rate (Vale, 1996; Ma and Taylor, 1997), which suggests that the nucleotide pocket is allosterically controlled by the interface with the MT. Among kinesin’s MT-facing domains (Figure 1C), L11 and α⁢4 undergo large conformational changes upon binding to the MT. L11 is located after sw-II (Figure 1E), followed by α⁢4 that N-terminally extends by a few turns when the motor head binds to the MT (Supplementary file 1) (Sindelar and Downing, 2010; Atherton et al., 2014; Shang et al., 2014). However, substantial conformational variations are present in these conserved domains, notably in sw-I and L11 (see Supplementary file 1). Also, the extent of the kinesin-MT interface varies depending on experimental conditions (Morikawa et al., 2015).

Thus, it is necessary to identify core features of the motor head that are essential for nucleotide processing. Such information about a single head is a prerequisite for the atomic-level understanding of the motility of a dimer. We characterize these features via multi-microsecond molecular dynamics simulations on the Anton supercomputer (Shaw et al., 2009; Shaw et al., 2014) of a motor head complexed with a tubulin dimer. Compared to previous all-atom simulations that used biasing potentials and were limited in time (Li and Zheng, 2012; Shang et al., 2014; Chakraborty and Zheng, 2015), the unbiased simulations described here reveal the conformational changes of kinesin-MT complexes on a more realistic time scale. We find that the nucleotide binding pocket is conformationally the most dynamic part of the motor head, whose internal motions actively drive the nucleotide processing events. In particular, we show how ATP hydrolysis occurs in a fluctuating environment, and demonstrate the role of the kinesin-MT interface for this process.

The dynamic nature of the kinesin mechanism elucidates how it robustly carries out its motility cycle despite significant conformational perturbations to the motor head and the MT in the crowded cellular environment (Leduc et al., 2012). The source of the chemical energy and the mechanism involved in ‘walking on tracks’ are applicable to other translocating motors such as myosin on actin (Vale and Milligan, 2000; Hwang and Lang, 2009), making the present mechanistic results of general interest.

Results

Simulation overview

We studied Kin-1 in different nucleotide or structural states. The names of the simulated systems are given below in italics. Simulation times and conformational states are in parentheses.

We also carried out simulations of Eg5 (Kin-5 family). However, currently only Kin-1 has atomic-resolution x-ray structures of the motor head complexed with the MT in both the pre- and post-stroke states (Gigant et al., 2013; Cao et al., 2014). Consequently, we focus our analysis on Kin-1, and use Eg5 for comparison.

Conformational dynamics of the motor head.

Related data for systems not shown here are in Figure 2—figure supplement 1. (A,B) Average displacement of Cα atoms from their initial positions (top row) and RMSD (bottom row) during the first/last 400 ns (black/red lines). The central β-sheet of kinesin was used for alignment. Domain names are in Figure 1—figure supplement 1. (A) ATP and (B) APO. (C) Rolling average (100-ns window) of displacements of Cα atoms excluding the CS and NL. Only residues with displacements greater than 1 Å were considered. (D–F) Conformation of domains around the nucleotide pocket at the end of simulation (cf., Figure 1E). (D) ATP. Sw-I lost its pseudo β-hairpin conformation. (E) Kin-only. N-terminal end of α⁢4 (star) and L11 also unfolded. (F) APO. Direction of motion of α⁢0 is marked by a dashed arrow. (G,H) Unfolding of the front end of the motor head. (G) ATP. The initial structure is shown in transparent gray, for comparison. (H) APO. (I) Curvature of the central β-sheet in ATP (yellow) and APO (red). Surfaces take the average curvatures in respective simulations. A higher saddle-point (Gaussian) curvature in APO can be seen by the greater bending of the corners 1, 3 and 2, 4 in opposite directions (arrows). (J) Combined energy landscape parameterized by the mean (M2) and Gaussian (G) curvatures of the central β-sheet.

Functionally important subdomains are mobile

We quantified the conformational motion by measuring the average displacement and root-mean-square deviation (RMSD) of Cα atoms relative to the first frame, during the first and the last 400 ns (Figure 2A,B and Figure 2—figure supplement 1A–D). Displacements represent deformation from the initial structure, and RMSD shows the degree of conformational fluctuation. To find conformational changes over time, we plotted rolling averages of mean Cα displacements that are greater than 1 Å (Figure 2C). To focus on changes in the core domains, the CS and NL, located at the termini of the motor head, were excluded from this calculation. Displacements saturate after 1–3 μs. Kin-only exhibits the greatest displacement, reflecting larger changes without the MT. ADP𝑝𝑟𝑒 had a large displacement between 2–3 μs, which is due to the release of ADP described below. Average displacements of other MT-bound structures were in the 2.5–3 Å range.

For each coordinate frame, we measured the mean curvature M2 (concaveness) and the Gaussian curvature G (saddle-point curvature) of the central β-sheet, and calculated the curvature free energy (potential of mean force; PMF) for each simulation (described in Materials and methods). These two curvatures quantify the bending and twisting of the β-sheet, respectively (Sun et al., 2003). Pre-stroke states had generally higher curvature, especially in G (Figure 2I and Figure 2—figure supplement 1K). Since ATP and ADP+Pi have very similar curvature, neither ATP hydrolysis nor Pi release (see below) is driven by the tendency of the β-sheet to adopt a higher curvature. Similarly, ATP and Kin-only had nearly the same curvature, indicating that binding to the MT does not impose any strain on the central β-sheet (Figure 2—figure supplement 1K).

We obtain information concerning the effect of curvature changes between pre- and post-stroke states by superposing the PMFs for ATP and APO (Figure 2J). APO has a local free energy minimum that is 0.84 kB⁢T (kB⁢T: thermal energy at 300 K) higher than that of ATP. There is also a ∼1.7 kB⁢T energy barrier from ATP towards APO. The pre- and post-stroke states respectively have similar PMFs regardless of the details of individual simulations (Figure 2—figure supplement 1K). Further, the PMF in Figure 2J does not directly represent the properties of the central β-sheet itself, but it implicitly reflects the energetics of the whole system, including domains surrounding the β-sheet, nucleotide and the MT, in controlling the curvature. These free energies are well below the 10-kB⁢T free energy (8 nm step×5 pN stall force) used by kinesin, which can also be seen by the large overlap in individual curvature distributions between the pre- and post-stroke states (Figure 2—figure supplement 1K). By comparison, the rotary motor F1-ATPase has about 5-kB⁢T curvature energy changes (Sun et al., 2003). Therefore, curvature changes in kinesin are not substantial enough to drive ATP hydrolysis nor the transitions between pre- and post-stroke states.

Mobility of kinesin on the MT.

(A,B) Summary of translational and rotational motion. See Figure 3—figure supplement 1A–G for detailed analysis. Direction and magnitude of (A) translation and (B) rotation of major domains between pre- and post-stroke states. (C) Hydration of the kinesin-MT interface. Blue blobs: regions where water density is higher than three times the bulk value. The hydration shell at cutoff equal to the bulk density is in Figure 3—figure supplement 1H. (D) Trajectory of buried area between kinesin and the MT (100-ns rolling average). Yellow: raw data for ADP+Pi, revealing large fluctuation. Raw data for other simulations fluctuate with comparable magnitude.

Kinesin-MT interface is dynamic and hydrated

Next we studied the motion of the motor head relative to the MT. For positional and orientational reference, we used the central β-sheet, α⁢6, α⁢4, and β⁢5⁢a/b (Figure 3A and Figure 3—figure supplement 1A). The central β-sheet, with its low RMSD, represents the overall position and orientation of the motor head. α⁢6 changes its orientation between pre- and post-stroke states (Figure 1C,D). α⁢4 and β⁢5⁢a/b are MT-binding domains. For each domain, translations in longitudinal, transverse, and normal (perpendicular to the MT surface) directions, and rotations about these three directions were measured.

We also calculated the water density map for the kinesin-MT interface during the last 500 ns. The map was visualized with two different density cutoffs. When a cutoff equal to the bulk density (0.0333 Å-1) is used, globular hydration shells surround the interface (Figure 3—figure supplement 1H). With a cutoff equal to three times the bulk density, a collection of blobs appear, which correspond to regions where water oxygens are found with high probability during the simulation (Figure 3C). They are located within the kinesin-MT interface and crevices, for all simulations. Lack of any correlation between the extent of interfacial hydration and the conformational state can also be seen by the buried area within the kinesin-MT interface. It undulates with 80–720-ns correlation times and with instantaneous fluctuations of a few hundred Å2 (e.g., yellow trace in Figure 3D).

We measured the binding energy between kinesin and the MT during the last 500 ns (Figure 3—figure supplement 1I–L). Pre-stroke states interact less with the β-tubulin (Figure 3—figure supplement 1L, squares), which is consistent with its front side (β⁢5⁢a/b) lifting from the MT (Figure 3A). Interaction with α-tubulin differs more across simulations. Overall, ATP and APO have the strongest binding energy (Figure 3—figure supplement 1L, circles), which are in line with experiments where the ATP and APO states have high MT affinity compared to the ADP state (Woehlke et al., 1997). However, our binding energies do not include water-mediated interactions and entropic contributions, which are expected to be comparable to the binding energy in magnitude (Zoete et al., 2005), so that the net binding free energy is much smaller than those in Figure 3—figure supplement 1L. Thus, the calculated binding energies, although they reflect the interaction between kinesin and the MT in different nucleotide states, do not correspond quantitatively to the experimental binding affinities. In any case, the mobility of the motor head and extensive hydration of the interface observed in all simulations suggest that the kinesin-MT interface is highly dynamic. This point is further explored in the analysis of kinesin-MT contacts below.

Nucleotide pocket experiences large changes in intra-kinesin contacts

To understand the conformational behavior of the system at the individual amino acid level, we traced all intra-kinesin and kinesin-MT contacts. Hydrogen bonds (H-bonds; including salt bridges) and nonpolar contacts were considered, majority of which form and break with less than 100% occupancy during the simulation (example occupancy trajectories are in Figure 4—figure supplement 1A). Among intra-kinesin contacts, there were fewer H-bonds (1100–1500) than nonpolar contacts (1800–2400). The occupancy distribution is U-shaped in logarithmic scale, majority of which have lower than 20% occupancy (Figure 4A). The number of contacts with greater than 80% occupancy were 137–153 (H-bond) and 302–339 (nonpolar). Among contacts showing irreversible transitions, we monitored those whose occupancy before breakage or after formation is greater than 80% (Figure 4B; Supplementary file 2). Post-stroke states had more contacts break than form, which occurred mainly within the first 3 μs (Figure 4—figure supplement 1B). Locations at which changes occurred are clustered around the nucleotide pocket and the front side (Figure 4B), that also had high RMSD (Figure 2A). With a bound nucleotide, contacts involving sw-I undergo extensive changes, which is responsible for the greater number of changes in the post- than in the pre-stroke state (Figure 4—figure supplement 1B). This is mainly because in the post-stroke state, contacts between the N- and C-terminal sides of sw-I forming the hairpin break. But contacts between sw-I with ATP and sw-II, necessary for the hydrolysis of ATP, remain intact (Supplementary file 2A–C; see below).

Contact analysis.

(A) Occupancy distribution of intra-kinesin contacts. Filled symbols: nonpolar contacts. Open symbols with solid lines: H-bonds. Figure 4—figure supplement 1A shows examples of occupancy trajectories. (B) Locations of intra-kinesin contacts that broke or formed during simulation (colored red). Bottom view. See Supplementary file 2 for the list of contacts and their transition times. Cumulative numbers of contacts formed or broken over time are in Figure 4—figure supplement 1B. Post-stroke states involve more contacts formed or broken, mainly around the ATP pocket and the front end of the motor head. The first frame of each simulation was used for visualization. In APO, the frontal part of the motor head was unstable from the beginning of the simulation, thus no clear contact changes were identified. (C) Residues forming kinesin-MT contacts with higher than 80% occupancy. Top view. For kinesin, only MT-facing domains are shown. Red: kinesin residues, yellow: MT residues. Stick: residues forming H-bonds, sphere: residues forming nonpolar contacts. In ATP, the three contact positions (rear, middle, front) are marked. Post-stroke states have more contacts with the front part. Supplementary file 3 lists individual kinesin-MT contacts.

Kinesin-MT contacts are plastic

Fewer contacts formed between kinesin and the MT, 170–240 H-bonds and 250–390 nonpolar contacts, of which only 4–10 (H-bond) and 5–16 (nonpolar) had greater than 80% occupancy (Supplementary file 3). Majority of nonpolar contacts are by charged or polar residues so that a hydrated interface is maintained (Figure 3C). In contrast to intra-kinesin contacts, very few contacts formed or broke irreversibly during the simulation (Supplementary file 3). Kinesin-MT contacts can be grouped into the rear (mainly L11 and α⁢4), middle (L12 and α⁢5), and front (β⁢5⁢a-L8b, herein called L8/β⁢5) (Figure 4C). The first two interact respectively with α and β-tubulins, and they are present in all nucleotide states. The front contacts are less robust in the pre-stroke states, lacking high-occupancy nonpolar contacts. This is consistent with the increase of its normal position (Figure 3A), higher (weaker) binding energy with β-tubulin (Figure 3—figure supplement 1L), and also with variations in its MT-binding mode in available structures (Sindelar and Downing, 2007; Morikawa et al., 2015; Shang et al., 2014; Atherton et al., 2014). Moreover, our analysis agrees with previous alanine-scanning experiments (Woehlke et al., 1997). Mutations in L8/β⁢5 (H156, E157, R161) caused marginal changes in the MT binding affinity, while mutating K252, Y274, and R278 in the rear and middle parts of the interface affected the MT affinity more strongly, which form high-occupancy contacts in our simulations (Supplementary file 3).

Variations in contact occupancy suggest that the kinesin-MT interface is maintained by an ensemble of contacts that do not need to be present simultaneously at any given time. In this way, kinesin may be able to bind to the MT quickly without needing to establish a precise combination of contacts. It also allows a certain degree of mobility of the motor head relative to the MT (Figure 3A). Nevertheless, kinesin binds selectively to the cleft on the MT lattice with β-tubulin in front, but not with α-tubulin in front (Figure 1B). Although the two tubulins have similar sequence and structure (Löwe et al., 2001), we found that some of the residues making contacts with kinesin diverge. Especially, H12 of β-tubulin has several contact residues that are non-homologous to those of α-tubulin (Figure 4—figure supplement 1C). Thus, the kinesin-MT interface is tuned so that it permits flexibility in binding, yet it is specific enough to recognize the MT binding site.

Sw-I hairpin unfolds

Our RMSD and contact analyses show that sw-I is among the most mobile kinesin subdomains. In the ATP-state, although its pseudo-hairpin structure has been suggested to be hydrolysis-competent (Kull and Endow, 2002), in all simulations of the post-stroke state, it unfolded (Figure 2D,E and Figure 2—figure supplement 1E; Video 1). The unfolding occurred well after simulation began, at 1.83 μs (ATP), 1.38 μs (Kin-only), and 0.98 μs (ADP+Pi). Even for a 1.73-μs simulation of an isolated Eg5 in the ATP state (Parke et al., 2010), sw-I unfolded at 522 ns (Video 1). Prior to full unfolding, contacts within the hairpin partially broke (Figure 5A,B), and other contacts within the surrounding domains or with sw-I broke or formed even earlier (Supplementary file 2). Thus, unfolding of the sw-I hairpin is a result of gradual changes that accumulate over time, rather than being an isolated event.

Mobility of sw-I in nucleotide processing.

(A) Unfolding of the sw-I hairpin in ATP (Video 1 and Figure 5—figure supplement 1A). Residues forming backbone H-bonds in the initial hairpin, and the rotation of α⁢3 are marked. (B) Occupancy trajectories of contacts for the hairpin in ATP. (C) Orientational angle of α⁢3 measured relative to the first frame of ATP. Red arrow: approximate time at which the sw-I hairpin unfolds. (D) Capturing ATP in the APO state by the ‘α⁢0/L5/sw-I trio’ (a 2.04-μs simulation; Video 2). 960 ns: Adenosine ring is close to its position in the bound state, but the phosphate moiety points outward (star). 1608 ns: Spontaneous formation of an α-helical turn in sw-I is visible. 2040 ns: ATP is positioned behind L5. Major domains that made contacts with the moving ATP are labeled (Figure 5—figure supplement 2A). (E) Pi release in ADP+Pi (Video 3; Figure 5—figure supplement 2B shows Pi release in Eg5). Box: Magnified view of Pi in contact with ADP. At 739.92 ns, sw-I pulls Pi out (arrow), after which it snaps back (arrow in 740.88 ns). (F) ADP release in ADP𝑝𝑟𝑒 (Video 4). Sw-I in an α-helical conformation turns and contacts ADP (200.88 ns). Outward rotation of sw-I moves ADP out of P-loop (668.16 ns). Later, sw-I loses its α-helical conformation (star in 2910 ns). A magnified view is in Figure 5—figure supplement 2E.

Process of ADP release by sw-I.

After unfolding, α⁢3 at the N-terminal side of sw-I rotated outward, increasing the distance between the two ends of sw-I (Figure 5A and Figure 5—figure supplement 1A). In the pre-stroke states, α⁢3 generally points outward (larger θ3 in Figure 5C), suggesting a tendency to move outward in the absence of ATP that holds sw-I. To refold into a hairpin, inward rotation of α⁢3 is necessary. Such rotation requires a broader conformational motion of the motor head that may occur over a time scale longer than that of our simulation. The apparent instability of the sw-I hairpin is at odds with its presence in several crystal structures (Supplementary file 1). In fact, the hairpin is stabilized by crystal contacts in these structures (Figure 5—figure supplement 1B–E). In comparison, the sw-I hairpin in myosin forms extensive contacts with the upper 50 kDa domain (Figure 5—figure supplement 1F). However, unfolding of sw-I in our simulation is only partial, where its N-terminal (outer) side separates from the C-terminal side, while the latter maintains contact with ATP. This was also the case for Kin-only (Figure 2E). In other x-ray structures of kinesin in the ATP (analog) states where sw-I does not adopt a clear hairpin conformation (Supplementary file 1), the C-terminal side in contact with the nucleotide is visible, which agrees with the partial unfolding in our simulation (e.g., Figure 5—figure supplement 1G).

To further examine the partial unfolding of the sw-I hairpin, we aligned the initial structure of ATP and the structure after the hairpin unfolding to high-resolution cryo-EM maps in the ATP states (Figure 6). The central β-sheet and α⁢4 were used as alignment references since their conformation varied little during the simulation (Figure 2A). To highlight differences between these and the cryo-EM structures, we rigidly docked the structures instead of performing flexible fitting. In both structures with the sw-I hairpin folded and unfolded, the N-terminal side deviates more from cryo-EM maps compared to the C-terminal side. Also, the outward rotation of α⁢3 (∼10∘; Figure 5C) is not enough to show any significant deviation from cryo-EM maps. In ATP, the unfolded N-terminal side deviates less from its position in the hairpin state compared to Kin-only (Figure 2D vs. E). Thus, at cryogenic temperatures, the N-terminal side is likely to settle to a hairpin-like state with the C-terminal side as a template, instead of landing in different configurations that lead to low electron density. These findings suggest the mobility of the N-terminal side of sw-I does not contradict existing cryo-EM data.

Sw-I conformation in cryo-EM structures of the ATP-state kinesin-MT complexes.

Yellow/magenta ribbons: Initial structure of ATP (PDB 4HNA) and structure at 3.79-μs that has an unfolded sw-I. Compared to the C-terminal side, the N-terminal side of sw-I aligns less well with cryo-EM maps (arrows), consistent with its mobility in our simulation. References: (A) (Atherton et al., 2014), (B) (Shang et al., 2014), (C,D) (Zhang et al., 2015). Resolutions of the maps are shown in each panel.

Binding of ATP is mediated by the α0/L5/sw-I trio

What would be the functional role of sw-I’s mobility? We first consider binding of an ATP molecule to kinesin in the APO state. For the APO system, we added a free Mg-ATP and performed another 2.04-μs simulation. To prevent ATP from diffusing away, we imposed a 32 Å radius spherical boundary on ATP around the center of mass of kinesin. During the simulation, ATP formed and broke contacts with various parts of kinesin (Figure 5D; Video 2). Nonpolar contacts were dominant, with the adenosine ring of ATP pointing toward kinesin and the charged phosphate moiety pointing away (Figures 5D, 960 ns, and Figure 5—figure supplement 2A). Direct binding of an ATP with the phosphate moiety pointing inward will be unfavorable due to the desolvation penalty for the phosphate moiety and the hydrophobic attraction for the adenosine ring.

ATP binding is more likely a multi-step process orchestrated by the surrounding domains. Among domains whose contact occupancy with ATP was high (labeled in Figures 5D, 2040 ns; Figure 5—figure supplement 2A), sw-I, α⁢0 (including L1a; Figure 1—figure supplement 1), and L5 take the shape of a funnel with the P-loop in the middle. During the simulation, the three domains transiently made contacts with ATP or even held it for a while (Figure 5D; Video 2). This α0/L5/sw-I ‘trio’ act like an antenna that captures nearby ATP and delivers it to the P-loop. Sw-I, the most mobile member of the trio (higher RMSD than the other two; Figure 2B), may be particularly important. When it moves away from the P-loop, it may form contacts with the adenosine ring, so that the phosphate moiety of ATP points towards the P-loop. A closing motion of sw-I will then bring the phosphate moiety in contact with the P-loop. Sw-I’s opening and closing motions have been observed in other simulations described below (cf., Figure 5E,F). Since ATP is amphiphilic, and since the trio domains closely surround the P-loop, their dynamic role should hold even though a complete binding event was not observed in our simulation – alternative scenarios such as ATP approaching kinesin with the phosphate moiety pointing towards the P-loop, or the trio domains not interacting with the incoming ATP despite their proximity, are physically unlikely.

ATP hydrolysis mechanism.

(A) Average water density map near Pγ during the last 500 ns of ATP. Blobs in semi-transparent blue have water density greater than three times the bulk density (cf., Figure 3C). Densities for two catalytic water molecules bridging between Pγ and E236 are in dark blue. A coordinate frame where the two-water bridge is present is displayed, with the catalytic water molecules marked by stars. (B) Water density map for Kin-only. The density close to E236 is broader. (C) Alignment of the structures in (A) (cyan) and (B) (yellow). R203/E236 are in blue (ATP) and red (Kin-only). In Kin-only, downward movement of L11 leads to lowering of R203-E236 (arrow). (D) Number of 2-water bridges between Pγ and E236 during respective simulations. (E) Conformations of Mg-ATP at the end of ATP and ATP-only. Mg2+ forms bidentate (ATP) and tridentate (ATP-only) contacts with O atoms (dotted lines). O atoms of Pβ and Pγ are eclipsed in ATP (blue stars), whereas they are staggered in ATP-only. (F) Angles defined in (E) (avg±std). The dihedral angle ϕγ reflects the eclipsed vs. staggered states.

Catalytic water molecules are dynamically coordinated

A critical question regarding the mobility of sw-I is whether it can support ATP hydrolysis. As noted above, even when the sw-I hairpin is unfolded, its inner side maintains contact with ATP. In its conserved SSR motif (S201-S202-R203) (Vale and Milligan, 2000; Kull and Endow, 2002), S201 and S202 contact Mg-ATP with higher than 99% occupancy in both ATP and Kin-only. Furthermore, R203 contacts E236 of sw-II (Figure 7A and Figure 7—figure supplement 1). We investigated whether this organization is sufficient to coordinate the catalytic water molecules. A previous ab initio calculation on Eg5 suggested a two-water mechanism: The ‘lytic’ water next to Pγ donates an OH group necessary for hydrolysis, and the released H atom travels through the second ‘transfer water,’ arriving at E236 (McGrath et al., 2013). Formation of a two-water bridge between Pγ and E236 is thus necessary for hydrolysis.

For ATP and Kin-only, we calculated the water density map around the phosphate moiety during the last 500 ns. In ATP, high-density blobs corresponding to the two catalytic water molecules were found, whereas in Kin-only, the density for the transfer water broadened (Figure 7A,B). This is because in Kin-only, the absence of the support by MT leads to a downward shift of L11 and R203-E236, so that the channel leading to Pγ widens (Figure 7C). Furthermore, during the simulation period, two-water bridges formed with higher frequency in ATP than in Kin-only (Figure 7D). These suggest that ATP hydrolysis is carried out in a dynamically fluctuating environment where contacts between ATP and sw-I (S201/S202), and between sw-I and sw-II (R203-E236) create a narrow channel that facilitates formation of the two-water bridge between Pγ and E236. Our result also explains the higher ATP hydrolysis rate when kinesin is bound to the MT: Without support from α-tubulin that keeps L11 ordered, R203-E236 move downward and broaden the channel (Figure 7C), thereby reducing the two-water bridge formation. This agrees with the lower but non-zero ATP hydrolysis rate of kinesin in the absence of the MT (Vale, 1996).

Kinesin-bound ATP is torsionally strained

In general, binding of a substrate to an enzyme is believed to induce mechanical strain, thereby lowering the activation energy for cleavage (Williams, 1993; Bustamante et al., 2004). To check for strain on the kinesin-bound ATP, we measured the internal coordinates of the phosphate moiety in ATP and Kin-only. For comparison, we performed a 4-ns simulation of an isolated Mg-ATP in water (named ATP-only). Between isolated and kinesin-bound ATP, the length of the cleaved Oβ–Pγ bond increased from 1.576±0.032 Å (ATP-only; avg±std) to 1.586±0.032 Å (ATP) and 1.586±0.033 Å (Kin-only). In the case of the Ras protein, increase in bond length by 0.01 Å has been shown to affect catalysis of GTP (Klähn et al., 2005). However, for myosin, although the Oβ–Pγ bond elongates in the active site, an energetic analysis reveals no significant destabilization of ATP (Yang and Cui, 2009). Furthermore, in the CHARMM param36 force field (Pavelites et al., 1997), the equilibrium length of the Oβ–Pγ bond is longer, 1.68 Å. Without an ab initio calculation of the energetics, it is difficult to assess the impact of stretching the bond on hydrolysis.

Internal angles of ATP had greater changes (Figure 7E,F). In particular, the dihedral angle ϕγ increased by 21∘–29∘ when ATP is bound to kinesin. This places the O atoms of Pγ in an ‘eclipsed’ (cis) position compared to ATP-only, where they are in a more relaxed, ‘staggered’ (trans) position (Figure 7E). The different torsional states of Pγ also affects contact with Mg2+. In ATP and Kin-only, Mg2+ forms bidentate contacts with two O atoms each from Pβ and Pγ. In ATP-only, it forms tridentate contacts with two O atoms of Pγ and one from Pβ (Figure 7E, dotted lines). Since the increase in the dihedral energy is not substantial (1.2 kcal/mol), torsional angles of the phosphate moiety should readily change when ATP binds to kinesin (calculation using a modified force field that worked well for certain ATP-bound protein structures (Komuro et al., 2014), also yielded only marginal changes in the dihedral energy). But holding only one Oγ atom by Mg2+ will result in a greater electron withdrawal effect compared to the case when the contact is shared between two oxygens. This permits the lytic water to more easily attack Pγ momentarily on the opposite side. Similar dihedral transition and charge redistribution in the phosphate group of GTP upon binding to Ras/GTPase activating protein have been observed (Rudack et al., 2012). For myosin, the eclipsing is present in both pre-powerstroke and post-rigor states. The latter state is incapable of hydrolyzing ATP since critical residues in the switch domains are displaced (Lu et al., 2017). Although a torsion-based ATPase mechanism may hold across nucleotide triphosphatases, there are likely multiple hydrolysis pathways, whose relative energetics may be determined collectively by the ATP conformation, catalytic water coordination, and conformations of residues immediately surrounding ATP as well as remote domains of the motor (Lu et al., 2017).

Release of hydrolysis products is mediated by the mobile sw-I

ADP+Pi models the state after ATP hydrolysis, and Pi was pulled out by sw-I as it moved away from the P-loop (Figure 5E; Video 3). The gap created between sw-I and ADP was sufficient for Pi to exit above (Figure 5E, 739.92 ns). In reality, there may be multiple Pi release paths due to the mobile sw-I. In a similar 2.04-μs simulation of the Eg5-MT complex, Pi released at 701 ns in a rearward direction (Figure 5—figure supplement 2B; Video 3).

Recent experiments suggest that the duration of the ADP+Pi state affects the processivity of a kinesin dimer (Milic et al., 2014; Andreasson et al., 2015; Mickolajczyk et al., 2015; Hancock, 2016). In the above simulations, Pi was monovalent (H2PO-4). In two simulations (3.7 μs and 3.8 μs each) of the Eg5-MT complex with a divalent phosphate (HPO2-4; P2-i), P2-i formed an extensive network of contacts with Mg-ADP and sw-I, and did not release (Figure 5—figure supplement 2C,D). P2-i is a high-energy transition state where the proton released after hydrolysis is added to convert it to Pi (McGrath et al., 2013). Since the proton can instead release into bulk water, the time of conversion from divalent to monovalent phosphate may depend on the time scale of proton transfer and other factors such as conformational fluctuation of kinesin. The phosphate release time also depends on the orientation of the phosphate in the nucleotide pocket. In another 2-μs simulation of Kin-1 with a monovalent Pi (as in ADP+Pi), its lone oxygen atom formed a contact with Mg2+, and release did not happen until the end of the simulation, analogous to the situation in Figure 5—figure supplement 2D. While various factors affect the time scale of Pi release, since sw-I firmly contacts Pi and is mobile, its outward motion is expected to be involved in driving the release.

Sw-I also facilitates ADP release, which was observed in ADP𝑝𝑟𝑒 (Figure 5F). At 201 ns, it swung toward ADP and its conserved N198 formed nonpolar contact with the adenosine ring (Figure 5F). This contact was transient and broke again. After a number of attempts, ADP was gradually pulled out (Figure 5—figure supplement 2E; Video 4). Sw-I then lost its α-helical conformation (Figures 5F, 2910 ns). It is unclear whether the α-helical state of sw-I is required for ADP release. A possible advantage is that the helix is more rigid than a disordered state, and it can exert a lever action for pulling ADP out of the P-loop. In the simulation of ATP binding, sw-I spontaneously formed an α-helical turn (Figures 5D, 2040 ns), which indicates that it can transition between disordered and helical states unless the SSR motif is stabilized by a bound ATP.

Discussion

The present results and previous findings lead us to propose a detailed model of kinesin dimer motility in which subdomain dynamics plays an essential role (Figure 8 and Video 5). It begins with the hydrolysis of the bound ATP in the rear head, which involves the sw-I–II connection, dynamic water coordination, and torsional strain in catalyzing ATP hydrolysis (Figure 8A). After ATP hydrolysis, the rear head changes its position and orientation slightly, which allows the front head to release ADP and fully bind to the MT (Figure 8B). Until the rear head releases Pi and detaches, ATP binding to the front head is prevented (‘gated’; Figure 8C). Once the rear head unbinds, possibly coupled with unfolding of α⁢4 (see below), ATP binding to the front head occurs, assisted by the α⁢0/L5/sw-I trio domains (Figure 8D). The resulting formation of the CNB in the front head generates the power stroke in which the rear head is thrown forward and begins a diffusive search for the next binding site. The E-hook of the MT helps with capturing what is now the front head (Figure 8E).

Motility of a Kin-1 dimer driven by subdomain dynamics.

Video 5 shows an overview of the process. In each panel, only relevant domains are shown. Tubulin dimers are labeled Tub-1–Tub-3. (A) ATP hydrolysis is driven by sw-I motion, sw-I–II connection, dynamic water coordination, and torsional strain (wavy line) on Pγ. The front head may not be fully bound. (B) ATP hydrolysis in the rear head allows full binding of the front head to Tub-2 (Milic et al., 2014), which releases ADP, mediated by sw-I. (C) Until Pi releases from the rear head (assisted by sw-I), ATP binding to the front head is gated (Andreasson et al., 2015; Hancock, 2016), potentially by α⁢0 that is close to the rearward-pointing NL. (D) Rear head (ADP state) unbinds, possibly coupled with unfolding of α⁢4. This allows ATP binding to the front head, assisted by the trio domains. (E) Power stroke is generated by the CNB formation (Hwang et al., 2008). The new front head contacts Tub-3 via interaction with the E-hook (Sirajuddin et al., 2014).

Illustrative model of the motility of a kinesin dimer.

The present work focuses on the properties of a single kinesin head that is likely to form the basis for diverse motility characteristics of different kinesins. In this regard, while the above model of dimer motility describes how it may be achieved by subdomain dynamics, the present model cannot address aspects pertaining specifically to the dimer motility, which would require knowledge of the communication between the two heads. Nevertheless, several general conclusions can be made. Our model highlights the active nature of nucleotide processing, where binding of ATP, hydrolysis, and release of hydrolysis products are mediated by concerted motions of mobile subdomains. The inherent mobility of sw-I is consistent with an experiment based on fluorescence resonance energy transfer (Muretta et al., 2015). It was shown in the paper that sw-I in both Kin-1 and Kin-5 bound to the MT stayed in an ‘open’ state more than 50% (mole fraction), even in the ATP state, as opposed to the ‘closed’ state that was assumed to take the hairpin conformation. The higher mobility or deformation of sw-I for an isolated kinesin compared to the MT-bound kinesin with ATP (Figure 2D,E) is also consistent with a previous electron paramagnetic resonance experiment (Naber et al., 2003). In the ATP state, the mobility of sw-I is limited such that its C-terminal side maintains contacts with ATP and sw-II that are necessary to support the hydrolysis. This was the case even for an isolated kinesin, whose lower ATP hydrolysis rate is likely due to the widening of the nucleotide pocket that reduces coordination of the catalytic water molecules (Figure 7).

The conformational fluctuation of the outer (N-terminal) side of sw-I may allosterically affect the catalytic water coordination and the precise orientation of residues of the SSR motif on the C-terminal side. In the case of myosin, these factors have been shown to affect the energetics of ATP hydrolysis (Lu et al., 2017). The hairpin conformation of sw-I could be more advantageous for hydrolysis than the unfolded state. However, since the hydrolysis reaction itself proceeds over a picosecond time scale (McGrath et al., 2013), the hairpin does not need to stay stably folded. Thus, variations in the sequence of the distal part of sw-I may provide an allosteric mechanism to fine-tune ATP hydrolysis and other kinetic rates by controlling its flexibility. This idea is supported by the behavior of the R190A/D231A mutant (Cao et al., 2014), which is expected to destabilize or prevent the hairpin state (Figure 5—figure supplement 1A). Rather than abolishing hydrolysis, its catalytic rate is 22% of the wild-type value. This demonstrates that the hairpin state may promote hydrolysis, but it is not required. Another illuminating feature of the dynamic role of sw-I is that the ATP hydrolysis rate depends on whether kinesin is bound to the MT filament or to unpolymerized tubulin (Alonso et al., 2007; Gigant et al., 2013). Since sw-I is the most labile element within the nucleotide pocket, its conformational motion may be affected the most when kinesin is bound to an unpolymerized tubulin, so as to influence the catalytic rate.

We propose that α⁢0 is the structural element responsible for ATP gating in the front head (Figure 8C). For the gating, the rearward orientation of the NL, rather than its tension, is essential (Clancy et al., 2011; Andreasson et al., 2015). In x-ray structures of kinesins with a rearward-docked NL, it interacts with the 3-stranded β⁢1 domain, which is linked to the C-terminal end of α⁢0 (Figure 1—figure supplement 1) (Sablin and Fletterick, 2004; Guan et al., 2017). This raises the possibility that the positional fluctuation of α⁢0 (Figure 2) is controlled by the interaction between the NL and β⁢1, thereby affecting the ATP binding. Further tests are needed to validate this proposal.

The relatively low specificity of the hydrated kinesin-MT interface (Figure 3C, Figure 4C, Supplementary file 3) is suited for rapid interaction with MTs in vivo (Leduc et al., 2012). In ADP+Pi, we did not observe unbinding of the motor head after Pi release, as is expected for the ADP-state kinesin. It is significant that ADP+Pi had more high-occupancy contacts with the MT than the other states (Supplementary file 3; this was the case even when only the last 500 ns, well after Pi release in ADP+Pi, were considered), and its MT-binding energy was among the strongest (Figure 3—figure supplement 1L). It appears to us unlikely that extension of the simulation time will lead to detachment of the motor head. To disrupt the kinesin-MT interface, we suggest that a conformational transition must occur. A likely subdomain for this transition is the N-terminus of α⁢4. It unfolded transiently in ADP+Pi and more extensively in Kin-only (Figure 2E). Conformational fluctuation of the unfolded α⁢4 can disrupt the rear part of the kinesin-MT contact (Figure 4C), leading to detachment. This picture agrees with the recent finding that Pi release is required for unbinding of kinesin from the MT (Milic et al., 2014). The presence of Pi in the nucleotide pocket suppresses the unfolding of α⁢4, thereby keeping the motor head bound to the MT. This is an interesting subject for studies beyond the present simulations.

Compared to a static mechanism that requires a specific structure, the dynamic mechanism for nucleotide processing and MT binding found here, provides greater flexibility in fine-tuning time scales and affinities. For example, the conserved residues that contact ATP cannot account for differences in catalytic rates among various kinesins. Altering residues that do not directly contact ATP can affect the conformational dynamics, which could in turn influence the frequency of 2-water bridge formation and charge fluctuation around Pγ, thereby controlling the catalytic rate.

Of interest are experimental means to test the dynamic roles of subdomains by introducing mutations. Since the N-terminal side of sw-I does not contact ATP directly, it may be possible to mutate it to alter the dynamical properties of sw-I without severely impairing motility. For example, a more flexible sw-I may enhance rates for all phases of nucleotide processing. However, caution must be exercised here, since if sw-I is made too flexible, it may disrupt contacts with the nucleotide. Similarly, elongating α⁢0 by lengthening L1a and L1b (Figure 1—figure supplement 1E) may increase the ATP binding rate, but it may also affect the gating behavior (Figure 8C). The design of mutants and their predicted behavior will require careful analysis and simulations.

We also showed that the elastic energy of the curvature of the central β-sheet or deformation of the MT, are unlikely to drive motility (Figure 2 and Figure 2—figure supplement 1). The aspect ratio of the motor head is too small to store any significant deformational energy. Even for myosin V that has a larger aspect ratio, there is little evidence that the twist of its β-sheet has any strong energetic role (Cecchini et al., 2008). Furthermore, the hydrated and dynamic kinesin-MT interface is unlikely to induce substantial strain in the motor head. Instead, small amount of strain may play a role in fine-tuning kinetic rates.

The free energy change upon NL docking in the absence of load (0.7–2.9 kcal/mol) (Rice et al., 2003; Muretta et al., 2015) is much smaller than the maximum work done near the stall force (5.8–8.1 kcal/mol). Additional energy is likely to originate from the CNB formation (Hwang et al., 2008; Khalil et al., 2008) and binding of the front head to the MT (Figure 8E). Since the hydrolysis energy of ATP thermalizes rapidly (on the picosecond time scale), it is unlikely to have any direct role. However, the differential binding energy of ATP and its hydrolysis products are likely to be important in triggering the large transitions of the motility cycle. Thus, while the net free energy change after a motility cycle may be that of hydrolyzing an ATP, there is a large free energy flow between kinesin and the environment during each phase of the cycle. In this regard, kinesin’s mobile subdomains are ‘free energy transducers.’

The present results provide the basis for understanding the role of local subdomain dynamics in the kinesin motility cycle. We emphasize that the extension of the simulations to multiple microseconds for each state, made possible by the use of Anton, played an important role in obtaining converged results. Suggestions are made concerning elements that will have to be tested by future experiments and simulations. This is important because the variation in structural rigidity, hydration, and protein-protein interaction found in the simulations provide a dynamic description of how kinesin works that is significantly different from conclusions based solely on static crystal and cryo-EM structures. Given the commonality among translocating motor proteins (Hwang and Lang, 2009), it is likely that local subdomain dynamics plays active roles for driving conformational changes and reactions, more generally.

Materials and methods

PDB structures used

ADP+Pi: The coordinate frame of ATP at 1043 ns, with ATP converted to ADP and Pi.

ADP𝑝𝑟𝑒: PDB 2P4N (9 Å resolution), a cryo-EM structure of nucleotide-free kinesin-MT complex (Sindelar and Downing, 2007). The fitted kinesin x-ray structure was based on PDB 1BG2 (1.8 Å resolution), which has a bound ADP and sw-I in α-helical conformation (Kull et al., 1996). The missing L11 and the N-terminal part of α⁢4 in the MT-facing domain were modeled after PDB 4HNA.

APOα: Same as ADP𝑝𝑟𝑒, with ADP removed. Sw-I was left in the α-helical conformation.

Among the above, PDB 4HNA and 4LNU are x-ray structures of kinesin-MT complexes respectively in ATP (ATP analogue) and APO states. Although the tubulin dimers in these structures are slightly curved, it has been shown not to affect the kinesin-MT interface (Gigant et al., 2013; Cao et al., 2014). We thus did not straighten the tubulin structure.

System preparation

We constructed the kinesin structure up to the NL (M1–A337), excluding the α-helical stalk. The C-terminal end of a tubulin has 13-aa glutamate-rich E-hook that are invisible in x-ray structures due to its flexibility. For α and β tubulins, we omitted the last 9 (E443–Y451) and 4 (E452–A455) residues of E-hooks, respectively. These truncations render the system size to fit within Anton. The E-hook of α tubulin is located on the minus end side of a tubulin dimer (at the left end of αH12 in Figure 1C) and is away from the kinesin motor head. The E-hook of β tubulin locates on the right side of the motor head (at the end of βH12 in Figure 1C). Being negatively charged and flexible, E-hooks are known not to affect kinesin in the MT-bound state, and it is more important for making non-specific electrostatic contacts with an unbound head (Figure 8E) (Lakämper and Meyhöfer, 2005; Sirajuddin et al., 2014). Truncations of E-hooks are thus unlikely to affect our result for MT-bound kinesins.

For each system, the protein structure was placed in a cubic TIP3P water box of linear size 113–119 Å (for kinesin-MT complex; 88 Å for Kin-only) and it was made electrically neutral by adding ions to about 50 mM concentration. The number of atoms in our systems were in the 150,000–170,000 range (65,000 for Kin-only). A periodic boundary condition was applied to the box.

Preparatory simulation

In preparation for simulations on Anton, the solvated system was simulated using CHARMM (Brooks et al., 2009) on a conventional computer cluster. Initially, a series of energy minimization procedure was done with harmonic constraints applied to proteins and nucleotides, which were gradually reduced to zero in successive 200-step minimization cycles. Next, the system underwent heating (100 ps) and equilibration (200 ps) runs under 1-atm pressure. During heating to 300 K, harmonic constraints were applied to backbone heavy atoms of proteins except for the 4-aa N-terminal end of kinesin’s CS and the C-terminal E-hook domains of MT that are flexible. The spring constant of the harmonic constraint was 1 kcal/mol⋅Å2 during heating, and 0.5 kcal/mol⋅Å2 during equilibration. It was further reduced to 0.25 kcal/mol⋅Å2, with only Cα atoms restrained (excluding those of the flexible domains noted above), and simulation continued for 2 ns using the constant temperature (300 K) and pressure (1 atm) (CPT) dynamics method implemented in CHARMM. The final phase of the preparatory run lasted 2 ns with 0.5-kcal/mol⋅Å2 harmonic constraints applied to Cα atoms of the loops of tubulins that are near the interface with neighboring MT protofilaments (aa 57–61, 83–88, and 279–286, for both tubulins). They are located on the bottom in Figure 1B, and restraining them mimics the effect of the tubulin dimer embedded within a polymerized MT. For Kin-only that lacks the MT, we harmonically restrained the Cα atoms of L229–D231 of β⁢7 (Figure 1—figure supplement 1) with a 0.1-kcal/mol⋅Å2 spring constant. These atoms are located approximately at the center of mass of the motor head, and the weak restraint suppresses translational diffusion of kinesin.

For simulation, the CHARMM param36 force field was used. For ADP+Pi, the monovalent form of Pi (H2PO-4) was constructed based on phosphate parameters in the param22 force field. The SHAKE algorithm was applied to fix the length between hydrogen and its base heavy atom. The integration time step was 2 fs.

Simulation on Anton

We wrote a Python script to convert the CHARMM restart file at the end of the preparatory run to the the Desmond Maestro format file, which was further processed using the Anton software. The CHARMM param36 force field was used through the Viparr utility of Anton. Harmonic restraints of spring constant 0.25 kcal/mol⋅Å2 were applied to Cα atoms of the same residues of the MT loops as in the last phase of the preparatory simulation. SHAKE was applied to hydrogen atoms, with a 2-fs integration time step. The multigrator integration method of Anton was used under a CPT (300 K, 1 atm) condition. Coordinates were saved every 0.24 ns. After simulation, coordinate trajectories were converted to CHARMM DCD format files by using VMD (Humphrey et al., 1996), for analysis using CHARMM. All simulations were carried out on Anton (Shaw et al., 2009) except for APO, which was on the newer Anton-2 machine (Shaw et al., 2014).

Curvature of the central β-sheet of kinesin

We considered seven strands within the central β-sheet of kinesin: β⁢1 (V11–F15), β⁢3 (T80–G85), β⁢4 (I130–Y138), β⁢5 (I142–D144), β⁢6 (S206–K213), β⁢7 (K226–L232), and β⁢8 (T296–C302) (Figure 1—figure supplement 1). This choice excludes regions of β⁢4, β⁢6, and β⁢7 at the front end of the motor head that deformed in some simulations (Figure 2G,H and Figure 2—figure supplement 1G–J). To calculate curvature, the central β-sheet in each coordinate frame was oriented to a reference kinesin structure whose least-square-fit plane was oriented to the x⁢y-plane of the Cartesian coordinate system and the center of mass positioned at the coordinate origin. In this configuration, z-coordinates of Cα atoms within the central β-sheet were parameterized by their x and y coordinates and were fit using the quadratic expansion (Sun et al., 2003)

(1)z⁢(x,y)=a0+a1⁢x+a2⁢x2+a3⁢y+a4⁢x⁢y+a5⁢y2,

where {ai} (i=0to5) are fitting parameters that vary among coordinate frames. Fitting was done using the SciPy package of Python. Examples of fitting surfaces are in Figure 2I. For a given frame, the mean and Gaussian curvatures are given by M=2⁢(a2+a5) and G=a42-4⁢a2⁢a5 (Sun et al., 2003).

To calculate the potential of mean force (PMF) versus the curvature, we calculated a 2-dimensional histogram of M2 and G normalized by the maximum count, ρ⁢(M2,G). PMF is given by -kB⁢T⁢ln⁡ρ (Figure 2—figure supplement 1K). To align PMFs for ATP and APO (Figure 2J), we identified bins of the histogram that have nonzero counts in both simulations. We determined the constant free energy shift Δ that needs to be added to the PMF for APO so that the mean-square difference of the two PMFs in the overlapping bins is minimized. The mean-square difference was calculated weighted by the histogram counts of respective PMFs. Denote the histogram values of ATP and APO in the k-th bin within the overlap region by ρk𝐴𝑇𝑃 and ρk𝐴𝑃𝑂, respectively, and similarly denote their PMFs by Ek𝐴𝑇𝑃 and Ek𝐴𝑃𝑂. Minimizing ∑kρk𝐴𝑇𝑃⁢ρk𝐴𝑃𝑂⁢(Ek𝐴𝑇𝑃-Ek𝐴𝑃𝑂-Δ)2 yields

(2)Δ=∑kρk𝐴𝑇𝑃⁢ρk𝐴𝑃𝑂⁢(Ek𝐴𝑇𝑃-Ek𝐴𝑃𝑂)∑kρk𝐴𝑇𝑃⁢ρk𝐴𝑃𝑂

where the sum is for bins in the overlap region. We added Δ to the PMF for APO, and the merged PMF Ek for the k-th bin in the overlap region was set to

Each coordinate frame was aligned to the first frame of ATP, with the Cα atoms of H12 helices (aa417–432) of α- and β-tubulins as reference for alignment. Let Rα−tub, and Rβ−tub be the centers of masses of H12 in respective tubulins. An orthonormal triad {uL,uN,uT} was constructed in the following way (Figure 3—figure supplement 1A):

Longitudinal direction: uL∝(Rβ−tub−Rα−tub).

Normal direction: Let uα4 be the axis vector of α⁢4 in the reference structure (first frame of ATP) pointing to the right, from the N- to C-termini of α⁢4. The unit vector in the normal direction was set as uN∝(uα4×uL).

Transverse direction: uT=uL×uN.

Let rβ be the center of mass of the central β-sheet. We projected (rβ−Rα-tub) onto the three directions of the triad. Differences of these projections from those of the reference structure were defined respectively as the longitudinal (ΔL), normal (ΔN), and transverse (ΔT) displacements. Displacements of α⁢4, Q320, and β⁢5⁢a/b were measured similarly.

To calculate the orientation of the motor head, we used an approximately rectangular section of β⁢4 (I130–E136), β⁢6 (S206–V212), and β⁢7 (K226–D231). Using their Cα atoms, we calculated the major and normal axes of the least-square-fit plane, vβ1 and vβ2, respectively. We projected vβ1 onto the plane spanned by uL and uN, and measured the forward tilt angle θβ as the angle between this projection and uL. The azimuthal angle ϕβ was measured as the angle between the projection of vβ1 onto the plane spanned by uL and uT, with uL. The transverse tilt angle ωβ was between the projection of vβ2 onto the plane spanned by uN and uT, with uN. Increase in ϕβ is associated with clockwise rotation of the motor head when viewed from top, and for ωβ, it is counterclockwise rotation when viewed from the MT plus end. Since the central β-sheet is tilted to the left, ωβ is typically negative, and in Figure 3—figure supplement 1F, -ωβ was plotted. A larger (less negative) ωβ indicates that the motor head tilts more to the right.

Azimuthal (ϕα) and transverse tilt (ωα⁢6) angles of α⁢6 were similarly measured using the projection of the axis of α⁢6 on respective planes. The orientation angle ψα⁢4 of α⁢4 was measured between its axis and uT.

Hydration analysis

To calculate the water density map for the kinesin-MT interface (Figure 3C), we adopted a method that we developed previously (Ravikumar and Hwang, 2011). Coordinate frames were aligned to the first frame of ATP with Cα atoms of the reference domains consisting of α⁢4 (E250–E270), and parts of H11–H12 of α-tubulin (F395–E420) and H12 of β-tubulin (M425–Y435). A search box was set whose boundary is at least 15 Å away from any atom in the above domains. The box was divided into a cubic grid of linear size 0.7 Å. For each cell in the grid, the fraction of frames where a water oxygen is found was calculated and divided by the volume of the cell (0.73 Å3). The map was saved into an MRC electron density map format file and visualized using UCSF Chimera (Pettersen et al., 2004). The water density map around the phosphate moiety of ATP (Figure 7A,B) was calculated similarly, with Cα atoms of the P-loop (G85–G90) and Pγ as positional reference for aligning coordinate frames.

To calculate the buried area in the kinesin-MT interface (Figure 3D), for each coordinate frame, we used a 1.4 Å probe radius to calculate the solvent accessible surface area of kinesin and the MT, together and separately, and measured their difference.

Binding energy calculation

Calculation of the binding energy was done based on a previously developed method (Zoete et al., 2005). Briefly, for each coordinate frame during the last 500 ns, the following energies were measured: van der Waals (Lennard-Jones), electrostatic, generalized Born solvation free energy, and the nonpolar energy. Calculation of energy terms was done using the Generalized Born with a simple SWitcthing (GBSW) module of CHARMM (Im et al., 2003). To find the binding energy, we used the kinesin motor head, α-tubulin, and β-tubulin, individually or in combination. For example, between the motor head and the α-tubulin, we measured EK for kinesin, Eα for the α-tubulin, and EK⁢α for them together (E denotes an energy term). Their binding energy was defined as EK⁢α-EK-Eα. Similar calculations were done for the binding energy between kinesin and the β-tubulin, and kinesin and the whole tubulin dimer.

Contact analysis

For each coordinate frame in a simulation, hydrogen bonds (H-bonds) were identified with the donor-acceptor distance cutoff of 2.4 Å. A residue pair was considered to form a nonpolar contact if the pair has neutral atoms (partial charge less than 0.3e; e=1.6×10-19 C) that are closer than 3.0 Å. The occupancy of a bond is the fraction of frames over which the bond is formed during the simulation, and its occupancy trajectory is the rolling (running) average with a 96-ns (400 frames) window.

To identify formation and breakage of a bond during simulation, we calculated its occupancy for the first and last 200 frames (48 ns), respectively. If the initial occupancy is less than 0.05, and the final one is greater than 0.5, the bond was regarded to have formed during the simulation, and vice versa for identifying bonds that broke. We traced the trajectory forward (bond formation) or backward (bond breakage) in time and located the time point where the local occupancy became greater than 0.5, as the corresponding transition time (Figure 4—figure supplement 1A).

41–53, Anton 2: raising the bar for performance and programmability in a special-purpose molecular dynamics supercomputer, Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis, 10.1109/SC.2014.9.

Decision letter

Antoine M van Oijen

Reviewing Editor; University of Wollongong, Australia

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.

[Editors’ note: this article was originally rejected after discussions between the reviewers, but the authors were invited to resubmit after an appeal against the decision.]

Thank you for submitting your work entitled "Kinesin Motility Driven by Subdomain Dynamics" for consideration by eLife. Your article has been favorably evaluated by a Senior Editor and three reviewers, one of whom is a member of our Board of Reviewing Editors. The following individuals involved in review of your submission have agreed to reveal their identity: Charles Vaughn Sindelar (Reviewer #2); Qiang Cui (Reviewer #3).

Our decision has been reached after consultation between the reviewers. Based on these discussions and the individual reviews below, we regret to inform you that your work will not be considered further for publication in eLife.

The main concern raised by the reviewers is that the observation of a mobile and unstructured Switch I element in the kinesin motor goes against a number of experimental observations previously reported. Combined with additional concerns surrounding potential artifacts in the simulations, we feel that this study does not provide sufficient evidence to support the claims.

In the extensive online consultation session between the reviewers, a number of additional points were raised that may help you in strengthening the work for submission elsewhere. One of these points was that there may be additional ordered water molecules not present in the simulations that could play a critical role in maintaining the structured nature of Switch I. Another aspect that was discussed is the influence that different force fields may have on the stability of the Switch I element.

Reviewer #1:

This manuscript describes the use of molecular dynamics simulations to investigate the dynamics of the switch I loop in kinesin and its conformational transitions throughout the ATPase and translocation cycles. The molecular dynamics interestingly suggest that elements of the kinesin motor domain are highly mobile and that the that plasticity plays a key role in nucleotide hydrolysis and coupling of the liberated energy to mechanical motion. While the technical aspects of the study are impressive with the extent and duration of the simulations, I have significant concerns about the lack of experimental validation of the proposed models. Especially with the observations described of a highly mobile switch I loop which seem to contradict published structural work – and with one of the co-authors being a specialist of exactly the right type of experimental approaches – one would expect strong experimental evidence be presented in the paper.

Reviewer #2:

In this manuscript, the authors use molecular dynamics simulations to look at mobile elements in the kinesin motor domain that are involved in nucleotide processing. Based on these simulations, they argue that the switch I loop and other elements close to the nucleotide pocket are much more mobile (particularly in the ATP-bound state of the motor) than was evident in prior structural studies performed with cryo-EM or X-ray crystallography. They further argue that this putative high mobility is relevant to the ATP binding, hydrolysis and product release functions of the motor. They also describe attempts to observe nucleotide binding and unbinding directly via simulation; while neither of these processes completed within the timescale of the reported molecular dynamics runs, some additional claims are made regarding structural elements that may be involved in nucleotide binding/unbinding. At the end of the manuscript a model for the kinesin cycle is presented that closely follows the longstanding 'consensus' model for kinesin motility, although the model does not incorporate (nor is reference made to) recent experimental results indicating that a forward step does not come until after ATP is hydrolyzed.

The main messages in this manuscript are centered around the claim that the switch I nucleotide binding loop is highly mobile, even in the microtubule-bound, ATP-bound state where prior structural work has strongly indicated to the contrary. If true, this would certainly change our understanding of the kinesin mechanism, but it is not clear from the manuscript what the mechanistic purpose would be, for example, of having the nucleotide pocket be so unstable in the presumptive catalytically active state. Such a bold and unexpected claim would seem to require independent validation, ideally from experimental methods. However, the only cited support for the switch I mobility claim in the current manuscript is highly indirect (FRET distance measurements from a recent study by (Muretta, Jun et al. 2015)). The evidence to the contrary, on the other hand, is abundant. For example, there are numerous reported cryo-EM structures of kinesin bound to ATP analogs on microtubules where the switch I loop shows little sign of disorder, going all the way back to Kikkawa's sub-nanometer structure in 2006 (Kikkawa and Hirokawa 2006) and continuing up to recent higher-resolution structures (Atherton, Farabella et al. 2014, Shang, Zhou et al. 2014). Perhaps the most striking cryo-EM structure, although not reported as such, is a near-atomic structure from the Nogales lab of microtubule-attached kinesin with a non-hydrolyzing switch II (E236A) mutation, bound with intact ATP ((Zhang, Alushin et al. 2015); structure can be viewed at EMDB ID 6348); again, there is little to no evidence of switch I mobility in this latter structure.

Moreover, the structure of kinesin's switch I loop has been studied by direct labeling with EPR probes (Naber, Minehardt et al. 2003) and found to be highly immobilized in the microtubule-bound, ATP analog bound states (this work is not cited in the current manuscript, but should be). The X-ray co-crystal structure of a kinesin-tubulin complex bound to ATP analogs (Gigant, Wang et al. 2013) also indicate that the switch I loop is highly ordered, without any indication from B-factors of switch I mobility. The current manuscript tries to argue that switch I in this latter structure has been perturbed by a crystal structure contact, but the associated figure (Figure 5—figure supplement 1B) reveals that the contact consists of a single hydrogen bond coming from lysine side chain that extends out from a nearby symmetry mate – I do not find this to be a compelling argument. The switch I loop is already known to be unstable in ATP analog-bound kinesins when tubulin is absent (viz. the previously cited EPR study as well as numerous X-ray studies showing disorder in the presence of ATP analogs; see for example (Nitta, Kikkawa et al. 2004), so it is not surprising to see more extensive crystal contacts stabilizing switch I in the cases where switch I folds to its 'catalytically active' conformation in the absence of tubulin (cf. Figure 5—figure supplement 1C-E). In summary, multiple existing lines of experimental evidence run counter to the authors' claim that switch I is mobile in its catalytically active phase on microtubules.

A major concern here is that the nucleotide pocket in the reported molecular dynamics simulations could be artifactually destabilized by, for example, inaccuracies in the initial conditions or energy function of the simulation. Indeed, the reported simulations seem to indicate that switch I is not only 'mobile' in the ATP-bound, microtubule-bound state, but markedly unstable: it is stated that the switch I loop dissociates from its ATP-coordinating position in multiple simulations, but apparently never re-folds. This feature of the simulation results amplifies the above concerns, because a truly unstable switch I would not show up in experimental density maps (X-ray or EM)- contrary to existing evidence.

The paper makes a number of other potentially interesting claims that are, however, not well substantiated either by the simulations or by prior findings. For example, it is suggested that (1) the initial stages of ATP recruitment might be facilitated by a trio of structural elements of L5, switch I and α 0; or (2) that ADP and Pi products might be 'carried out' of the nucleotide pocket by the switch I loop as a part of nucleotide release. Unfortunately, due to limitations in existing computing power, the timescales of the reported simulations (while impressive) is not sufficient to investigate such claims in detail. The validity of these latter two proposals is also called into question due to the fact that they both involve the switch I loop; the likelihood that the current simulations have problems that may destabilize the structure of switch I (see above) is therefore a significant concern. Moreover, recent findings that the tethered partner head steps only steps forward after hydrolysis, while ADP•Pi is bound in the active site of attached kinesin domain ((Milic, Andreasson et al. 2014, Andreasson, Milic et al. 2015, Mickolajczyk, Deffenbaugh et al. 2015); also reviewed in (Hancock 2016)) seems to indicate that ADP•Pi should persist in kinesin's active site long enough for a successful forward step. In contrast, in the currently reported simulations the active site falls apart (and phosphate dissociates) in less than a microsecond. This seeming discrepancy with the recent experimental evidence is noteworthy, but is never mentioned in the manuscript despite many of the relevant articles having been cited.

I was interested in the observation and discussion of strain in ATP and in kinesin's β sheet. However, these are not strongly emphasized in the manuscript and I am left with questions/concerns, particularly with regards to the β sheet strain. The model for strain is a second-order Gaussian surface, which is characterized by a single curvature parameter. It is not clear to me whether this can accurately capture the 'wrinkle' that appears in kinesin's β sheet when the two halves rotate with respect to each other. The free energy of the strain could therefore be significantly underestimated. For the ATP strain, it is not clear whether the magnitude of the observed strain (1.2kcal/mol) falls in the expected range, or how this would change our view of kinesin's mechanochemistry.

In summary, I do not find the major claims of this paper, regarding the nucleotide pocket and the switch I loop, to be justified by the reported experiments – and indeed I think they likely arise from artifacts related to the simulation methods. The Discussion does not connect the simulation results to the most recent findings regarding kinesin's ADP•Pi state, even where there are seeming contradictions. Some of the other observations in the paper, particularly regarding strain in ATP and/or the β sheet, seem to have more merit. However, on the whole I think this work should be substantially revised and belongs in a more specialized journal.

Reviewer #3:

In this molecular dynamics (MD) simulation study of kinesin, the authors conducted a systematic analysis of conformational dynamics of a single kinesin motor domain, either bound to a microtubule (MT) dimer unit or in solution; several nucleotide binding states (ATP, ADP/Pi, ADP and Apo) have also been examined. Compared to previous MD studies, the current work substantially extended the simulation time scale to multiple microseconds (for each state) by taking advantage of resources available on Anton. This is significant since as the authors showed that in most cases, the conformational properties of interest reach a plateau only after about one or a few microsecond(s). The extensive sampling is also essential to a meaningful estimate of the bending energy of the central β-sheet and an evaluation of its relevance to nucleotide processing.

The MD results clearly highlighted that the motor domain is highly dynamical and features rich variations in structural rigidity, hydration and protein-protein as well as protein-nucleotide contacts. These observations provide a rather different description of the kinesin motor domain from static crystal structures and explain the variation of several structural motifs (e.g., SwI) in different crystal structures. The comparison of different nucleotide states also made it possible to reveal specific structural motifs that are likely important to ATP binding, activation of hydrolysis and facilitation of the release of Pi and ADP. The plastic nature of kinesin/MT interface has not been emphasized in previous studies and is likely important to kinesin/MT recognition in vivo. Overall, the results provided new insights into the motor domain properties that are likely important to the mechanochemical coupling in kinesin, and the study has laid the ground work for better understanding communication between two kinesin domains at a molecular level.

[Editors’ note: what now follows is the decision letter after the authors submitted for further consideration.]

Thank you for resubmitting your work entitled "Kinesin Motility Driven by Subdomain Dynamics" for further consideration at eLife. Your revised article has been favorably evaluated by Arup Chakraborty (Senior Editor) and three reviewers, one of whom is a member of our Board of Reviewing Editors and another who was not involved with the evaluation of the original manuscript.

The manuscript has been improved considerably and reviewers are satisfied with the way in which you have addressed the points raised during the previous round of review. All reviewers, including a new reviewer, agree that the work described in this paper is of high quality. There remain some technical concerns that will need to be addressed before we can reach a firm decision to accept the paper for publication in eLife. These technical points are outlined below:

1) It is stated that kinesin-bound ATP is torsionally strained. There is a comparison with ATP in solution to support this conclusion (For comparison, we performed a 4-ns simulation of an isolated Mg-ATP in water). The authors may wish to have a look at an improved parameterization of the dihedrals of ATP in the CHARMM force field (Komuro et al. JCTC vol 10, 4133-4142, 2014) to verify that this conclusion is robust or force field dependent.

2) The significance of the PMF shown in Figure 2J is not entirely clear. It is stated that "To align PMFs for ATP and APO, we identified bins of the histogram that have nonzero counts in both simulations. We determined the constant free energy shift that needs to be added to the PMF for APO so that the mean-square difference of the two PMFs in the overlapping bins is minimized." The energy barriers on this combined free energy surface are also discussed ("There is also a Ì1.7 kB T energy barrier from ATP towards APO"). In practice, two separate PMFs were extracted from the marginal distributions (-kt*ln(rho)) obtained from two separate simulations that correspond to two very different molecular contexts (one with ATP bound and one without ATP bound). These two PMFs were combined by matching the histograms that have nonzero counts. But that is presuming that the two surfaces exist with respect to the same part of configuration space. Since one of them has a ligand bound and the other does not, why should this be so? Unless we misunderstand, the configurational sampling of the two state do not overlap, so how can the two PMFs be combined? In principle, the relative shift between two free energy surfaces, one ligand-bound and one APO, should depend on the concentration of the ligand. If the ligand concentration is zero, then the PMF of the ligand-bound state shifts to infinity, and likewise, if the ligand concentration is infinite then the PMF of the APO state shifts to infinity. This point needs to be addressed and clarified.

3) In Figure 3E, the binding energies are of the order of -100 to -250 kcal/mol. These values appear to be very large, and are probably not standard binding free energies. The Materials and methods section states that "Calculation of the binding energy was done based on a previously developed method (Zoete et al., 2005)." This method is MM/GB-SA (molecular mechanics-generalized Born surface area) which is an end-point approximation to the binding free energy. There is scepticism that the binding energy estimated by the Generalized Born method is sufficiently accurate to draw meaningful conclusions for the system being considered. How sensitive are the conclusions to this analysis?

Author response

Editors’ note: the author responses to the first round of peer review follow.]

Reviewer #1:

This manuscript describes the use of molecular dynamics simulations to investigate the dynamics of the switch I loop in kinesin and its conformational transitions throughout the ATPase and translocation cycles. The molecular dynamics interestingly suggest that elements of the kinesin motor domain are highly mobile and that the that plasticity plays a key role in nucleotide hydrolysis and coupling of the liberated energy to mechanical motion. While the technical aspects of the study are impressive with the extent and duration of the simulations, I have significant concerns about the lack of experimental validation of the proposed models. Especially with the observations described of a highly mobile switch I loop which seem to contradict published structural work – and with one of the co-authors being a specialist of exactly the right type of experimental approaches – one would expect strong experimental evidence be presented in the paper.

We thank reviewer 1 for recognizing the technical strength and novelty of our work. The consistency of sw-I’s (partial) mobility with existing data is detailed in our response to comments by reviewer 2 below, who had more specific points. Our simulation agrees with the available experimental data, as explained above and throughout the manuscript. Given that the proposed mechanism is novel, we believe it is important to publish it, as a prediction that we hope will stimulate further experiments.

Change made:

Subsection “Concluding Discussion”, last paragraph: Need for future experiments is explained, in light of the conceptual advances made in this work.

Reviewer #2:

In this manuscript, the authors use molecular dynamics simulations to look at mobile elements in the kinesin motor domain that are involved in nucleotide processing. Based on these simulations, they argue that the switch I loop and other elements close to the nucleotide pocket are much more mobile (particularly in the ATP-bound state of the motor) than was evident in prior structural studies performed with cryo-EM or X-ray crystallography. They further argue that this putative high mobility is relevant to the ATP binding, hydrolysis and product release functions of the motor. They also describe attempts to observe nucleotide binding and unbinding directly via simulation; while neither of these processes completed within the timescale of the reported molecular dynamics runs, some additional claims are made regarding structural elements that may be involved in nucleotide binding/unbinding. At the end of the manuscript a model for the kinesin cycle is presented that closely follows the longstanding 'consensus' model for kinesin motility, although the model does not incorporate (nor is reference made to) recent experimental results indicating that a forward step does not come until after ATP is hydrolyzed.

We reply to the comments about the mobility of subdomains in other replies below. In referring to a kinesin dimer, in the caption of Figure 8, we do make reference to the paper that reviewer 2 mentioned regarding “recent experimental results indicating that a forward step does not come until after ATP is hydrolyzed” (Milic et al., 2014). But we mistakenly stated that Pirelease (instead of ATP hydrolysis) completes a step. This is corrected in the revision. However, this does not affect our main findings, which concern the nucleotide processing mechanisms of a single motor head. Figure 8 demonstrates how the dynamic mechanism may be at work in the dimer motility. In reviewer 3’s second point, reviewer 3 recognizes this: “the results provided new insights into the motor domain properties that are likely important to the mechanochemical coupling in kinesin, and the study has laid the ground work for better understanding communication between two kinesin domains at a molecular level.” We clarified these points in the revised manuscript.

While the motility of a kinesin dimer is not the main focus of the present work, the model of dimer motility based on subdomain dynamics provides a possible explanation for how Pirelease triggers detachment of the rear head, which directly aligns with the finding by Milic et al.(2014): With Piin place, conformational transition of α4 that leads to detachment of the motor head (Figure 8D) is suppressed, so that Pirelease is required for unbinding of kinesin from the MT.

Changes made:

Introduction, second paragraph, Figure 8B, C, and its caption are corrected, to be consistent with the model proposed in Milic et al.(2014). Video 5 is also updated.

Subsection “Concluding Discussion”, second paragraph: Scope of the present work is explained (please also see the second paragraph of the Introduction).

The main messages in this manuscript are centered around the claim that the switch I nucleotide binding loop is highly mobile, even in the microtubule-bound, ATP-bound state where prior structural work has strongly indicated to the contrary. If true, this would certainly change our understanding of the kinesin mechanism, but it is not clear from the manuscript what the mechanistic purpose would be, for example, of having the nucleotide pocket be so unstable in the presumptive catalytically active state. Such a bold and unexpected claim would seem to require independent validation, ideally from experimental methods. However, the only cited support for the switch I mobility claim in the current manuscript is highly indirect (FRET distance measurements from a recent study by (Muretta, Jun et al. 2015)). The evidence to the contrary, on the other hand, is abundant. For example, there are numerous reported cryo-EM structures of kinesin bound to ATP analogs on microtubules where the switch I loop shows little sign of disorder, going all the way back to Kikkawa's sub-nanometer structure in 2006 (Kikkawa and Hirokawa 2006) and continuing up to recent higher-resolution structures (Atherton, Farabella et al. 2014, Shang, Zhou et al. 2014). Perhaps the most striking cryo-EM structure, although not reported as such, is a near-atomic structure from the Nogales lab of microtubule-attached kinesin with a non-hydrolyzing switch II (E236A) mutation, bound with intact ATP ((Zhang, Alushin et al. 2015); structure can be viewed at EMDB ID 6348); again, there is little to no evidence of switch I mobility in this latter structure.

As mentioned at the beginning of this reply, unfolding of the sw-I hairpin in the ATP state is only partial: Its N-terminal side breaks backbone hydrogen bonds with the C-terminal side of sw-I (R190–T195 in Figure 5A, B). However, the C-terminal side containing the SSR motif maintains catalytically important contacts with ATP and the sw-II domain throughout the simulation (subsection “Catalytic water molecules are dynamically coordinated”, first paragraph). Thus, reviewer 2’s remark, “the nucleotide pocket be so unstable in the presumptive catalytically active state,” is inapplicable.

To compare the observed mobility of sw-I with cryo-EM structures that reviewer 2 mentioned, we docked the structures before and after unfolding of the sw-I hairpin. We only considered recent high resolution structures of Kinesin-1, and did not consider Kikkawa’s structure for Kinesin-3 (KIF1A). For an unbiased comparison, instead of performing flexible fitting, we rigidly docked atomic structures to the cryo-EM map using the stably folded central β-sheet and α4 as reference. We added this as the new Figure 6. Even for the structure with intact hairpin (yellow in Figure 6), the N-terminal side of sw-I does not match well with the EM density compared to the C-terminal side which firmly contacts ATP (Figure 6, arrows). Given the anchoring of the C-terminal side, the motion of the N-terminal side of sw-I is limited. As Figure 2D vs. E shows, the unfolded state of sw-I is closer to the hairpin conformation compared to the case of an isolated kinesin in the ATP state. At cryogenic temperatures, the N-terminal side is thus likely to settle to a hairpin-like state with the C-terminal side as a template, instead of landing in different configurations that lead to low electron density.

In cryo-EM experiments, in addition to the low temperature, MTs are densely covered with kinesins, so that kinesins would be conformationally more restricted than when individual motors walk on the MT at room or body temperature. Furthermore, kinesins can hydrolyze ATP when they are bound to unpolymerized tubulin dimers, with catalytic rates ranging from lower than to comparable to those bound to MT filaments (Alonso, et al., 2007). Even an isolated kinesin can hydrolyze ATP, though with a lower rate (Introduction, third paragraph). We also demonstrate that even when the N-terminal side of sw-I is detached from the C-terminal side, catalytic water molecules and contacts between sw-I and ATP/sw-II are maintained, which are sufficient for hydrolyzing an ATP. Thus, the “mechanistic purpose” of the flexibility of sw-I would be to fine tune kinetic rates. Since the residues of sw-I contacting ATP and sw-II are highly conserved, they by themselves cannot determine family-specific catalytic rates. By modifying flexibility of the distal domain, the condition under which the hydrolysis reaction occurs can be controlled. We also note that, the hydrolysis reaction occurs over a picosecond time scale, as ab initiocalculations show (McGrath et al., 2013). Even though a fully folded hairpin state of sw-I may support hydrolysis more efficiently, it does not need to stay stably folded throughout the ATP state. Cao et al.(2014) shows that the R190A/D231A mutant has the hydrolysis rate reduced to 22% of the wild-type value. As R190 is at the N-terminus of sw-I, breaking this contact will destabilize the hairpin, but it does not abolish ATP hydrolysis (see reply to reviewer 2’s next pointbelow). This further supports that the propensity for forming hairpin can affect the hydrolysis rate, but it is not required. Thus, the x-ray and cryo-EM data, which are static in nature, provide the foundation for developing the dynamic picture in our study, rather than contradicting it.

As to the FRET paper by Muretta et al. (2015), reviewer 2 states that the measurement is “highly indirect.” However, a main focus of the paper is the mobility of sw-I, as its title suggests: “The structural kinetics of switch-1 and the neck linker explain the functions of kinesin-1 and Eg5.” As explained in our manuscript (subsection “Concluding Discussion”, second paragraph), their finding that sw-I stays ‘open’ in significant mole fraction, even in the ATP state, is consistent with our picture. But we also agree with reviewer 2 that experiments probing the flexibility of sw-I are scant, perhaps because the dynamic nature of sw-I has not been recognized. Since our work is the first systematic study, we hope that our work will lead to further experimental studies probing the dynamic aspect.

Changes made:

Added Figure 6 that compares cryo-EM structures and structures before and after unfolding of the sw-I hairpin, together with discussions in the last paragraph of the subsection “Sw-I hairpin unfolds”. For Figure 6C, D, Nogales’ paper (Zhang et al., 2015) is cited – we thank the reviewer for pointing out this impressive work.

Subsections “Nucleotide pocket experiences large changes in intra-kinesin contacts”; “Sw-I hairpin unfolds”, second paragraph; “Catalytic water molecules are dynamically coordinated”, first paragraph: Explanations about the partial nature of the sw-I unfolding are given.

Subsection “Concluding Discussion”, second and third paragraphs: Expanded discussion about the flexibility of sw-I, including the partial nature of unfolding, mechanistic purpose of fine tuning, works by Alonso, et al. (2007), and the behavior of the R190A/D231A mutant.

Moreover, the structure of kinesin's switch I loop has been studied by direct labeling with EPR probes (Naber, Minehardt et al. 2003) and found to be highly immobilized in the microtubule-bound, ATP analog bound states (this work is not cited in the current manuscript, but should be). The X-ray co-crystal structure of a kinesin-tubulin complex bound to ATP analogs (Gigant, Wang et al. 2013) also indicate that the switch I loop is highly ordered, without any indication from B-factors of switch I mobility. The current manuscript tries to argue that switch I in this latter structure has been perturbed by a crystal structure contact, but the associated figure (Figure 5—figure supplement 1B) reveals that the contact consists of a single hydrogen bond coming from lysine side chain that extends out from a nearby symmetry mate – I do not find this to be a compelling argument. The switch I loop is already known to be unstable in ATP analog-bound kinesins when tubulin is absent (viz. the previously cited EPR study as well as numerous X-ray studies showing disorder in the presence of ATP analogs; see for example (Nitta, Kikkawa et al. 2004), so it is not surprising to see more extensive crystal contacts stabilizing switch I in the cases where switch I folds to its 'catalytically active' conformation in the absence of tubulin (cf. Figure 5—figure supplement 1C-E). In summary, multiple existing lines of experimental evidence run counter to the authors' claim that switch I is mobile in its catalytically active phase on microtubules.

The EPR measurement of Naber et al.(2003) shows that sw-I changes from an open to a closed state when kinesin binds to the MT in the presence of ATP, and it does not address the nucleotide-dependent conformational transition of sw-I in a MT-bound kinesin. As shown in Figure 2D, E, as well as in the displacement/RMSD data, sw-I is deformed more from the initial hairpin conformation in Kin-only than in ATP. Instead of contradicting the EPR measurement, our simulation is thus consistent with it.

For the x-ray structure with ATP-analog in Nitta et al. (2004) that reviewer 2 pointed out (PDB 1VFV), we constructed the corresponding crystal and added it as Figure 5—figure supplement 1G. Sw-I in this case is only partially visible. As expected, it does not have any crystal contact. Notably, the C-terminal side of sw-I including the SSR motif forms contact with the nucleotide, whereas the partially visible N-terminal side is detached from the C-terminal side (shown on the right side of Figure 5—figure supplement 1G). This is entirely consistent with our observation in Kin-only simulation. We have verified similar behaviors of partially visible sw-I in other x-ray structures of kinesin in ATP-analog states listed in Supplementary file 1.

Regarding the crystal contact in PDB 4HNA (Figure 5—figure supplement 1B), the lysine side chain of the neighboring complex makes contact with the backbone oxygen of R190, which is located at the N-terminal end of sw-I. Since the sw-I hairpin opens in a zipper-like manner starting from R190 (Figure 5A, B), constraining it there should be sufficient to stabilize the entire hairpin. For other crystal contacts, although reviewer 2 said “The switch I loop is already known to be unstable in ATP analog-bound kinesins when tubulin is absent,” we feel it is important to explicitly show these contacts, especially since our simulation demonstrates the flexibility of sw-I.

Regarding reviewer 2’s summary, we again emphasize that, unfolding is only partial, and our results elucidate dynamic aspects of kinesin as an ATPase without contradicting previous experiments.

A major concern here is that the nucleotide pocket in the reported molecular dynamics simulations could be artifactually destabilized by, for example, inaccuracies in the initial conditions or energy function of the simulation. Indeed, the reported simulations seem to indicate that switch I is not only 'mobile' in the ATP-bound, microtubule-bound state, but markedly unstable: it is stated that the switch I loop dissociates from its ATP-coordinating position in multiple simulations, but apparently never re-folds. This feature of the simulation results amplifies the above concerns, because a truly unstable switch I would not show up in experimental density maps (X-ray or EM)- contrary to existing evidence.

Through many tests that we performed, we are confident that unfolding of the sw-I hairpin is not a simulation artifact. The unfolding also occurs in currently ongoing simulations of a kinesin dimer on the MT protofilament (not included in our manuscript). Our careful simulation preparation and execution are also recognized by reviewer 3: “The work has been rigorously conducted and analyzed.” Neither is it likely that the unfolding of sw-I is due to the inaccuracy of the force field, since many other aspects of our simulation are consistent with previous experimental data. The force field for ATP cannot be an issue either, since contacts between ATP with the P-loop and the C-terminal part of sw-I were maintained throughout simulations, and ATP does not have any direct contact with the N-terminal part of sw-I. Moreover, the CHARMM param36 force field used in this study is presently the most accurate available, and it has been successfully used in many other multi-microsecond simulations on the Anton supercomputer.

While reviewer 2 mentions that sw-I is “markedly unstable,” stability has a relative meaning. Since the sw-I hairpin unfolding occurred only after at least several hundred nanoseconds from the start of our simulation, it may not be observed in shorter simulations. For example, Chakraborty and Zheng (2015) performed a 400-ns all-atom simulation of a kinesin-MT complex in the ATP state, and did not report unfolding of sw-I. Notably, the 3–4-˚A RMSF of sw-I in their simulation (Figure 2 of Chakraborty and Zheng, 2015) is larger than that during the first 400 ns of our ATP simulation (less than 2.1 ˚A; Figure 2A of our manuscript), which further supports that our simulation system was well-prepared. Although not discussed, the RMSF analysis in Chakraborty and Zheng (2015) also indicates that sw-I is among the most fluctuating subdomains of the motor head, which is consistent with our result.

As explained in Figure 5A–C, refolding of sw-I is hampered by the outward rotation of α3, which keeps the two termini of sw-I apart. Its refolding is thus not simply a local event, but it will involve the conformational motion of the entire motor head that allows inward rotation of α3. Such global-level changes would occur on a time scale longer than that of our simulation, which may be why refolding of the sw-I hairpin was not observed. However, the rate of folding/unfolding of the sw-I hairpin, although it might affect the ATP hydrolysis rate, does not affect our main finding, since we have clearly shown that sw-I can support ATP hydrolysis as long as its C-terminal side maintains proper contacts with ATP and sw-II (Figure 5).

Regarding reviewer 2’s remark, “it is stated that the switch I loop dissociates from its ATP-coordinating position,” we did notmake that statement – even after the hairpin unfolds, catalytically important contacts involving the C-terminal part of sw-I and coordination of catalytic water molecules, are maintained. As explained in response to reviewer 2’s second and third pointsabove, we believe that the revised manuscript clarifies the partial nature of unfolding and its consistency with previous experimental data.

Changes made:

Subsection “Functionally important subdomains are mobile”, second paragraph: Consistency of our result with the mobility of sw-I in Chakraborty and Zheng (2015) is explained.

Subsection “Functionally important subdomains are mobile”, second paragraph: Explains the slow time scale associated with the refolding of the sw-I hairpin.

The paper makes a number of other potentially interesting claims that are, however, not well substantiated either by the simulations or by prior findings. For example, it is suggested that (1) the initial stages of ATP recruitment might be facilitated by a trio of structural elements of L5, switch I and α 0; or (2) that ADP and Pi products might be 'carried out' of the nucleotide pocket by the switch I loop as a part of nucleotide release. Unfortunately, due to limitations in existing computing power, the timescales of the reported simulations (while impressive) is not sufficient to investigate such claims in detail. The validity of these latter two proposals is also called into question due to the fact that they both involve the switch I loop; the likelihood that the current simulations have problems that may destabilize the structure of switch I (see above) is therefore a significant concern. Moreover, recent findings that the tethered partner head steps only steps forward after hydrolysis, while ADP•Pi is bound in the active site of attached kinesin domain ((Milic, Andreasson et al. 2014, Andreasson, Milic et al. 2015, Mickolajczyk, Deffenbaugh et al. 2015); also reviewed in (Hancock 2016)) seems to indicate that ADP•Pi should persist in kinesin's active site long enough for a successful forward step. In contrast, in the currently reported simulations the active site falls apart (and phosphate dissociates) in less than a microsecond. This seeming discrepancy with the recent experimental evidence is noteworthy, but is never mentioned in the manuscript despite many of the relevant articles having been cited.

Observing a complete ATP binding event in an unbiased simulation is beyond the present-day computing power. However, physical factors associated with ATP binding can be found by observing how kinesin interacts with a free ATP without needing to simulate the complete binding process. For this, our simulation of ATP binding provides the following information: 1) The trio domains (α0, L5 and sw-I) surrounding the P-loop are mobile (Figure 2B; their mobility was also observed by their large RMSF in Figure 2 of Chakraborty and Zheng, 2015); 2) Since the trio domains are located further on the outer side of the motor head than the P-loop, when an ATP approaches, it will first interact with the trio (Figure 5—figure supplement 2A); 3) Since ATP is amphiphilic, it is more likely to approach kinesin with the adenosine ring pointing inward while the more hydrophilic triphosphate moiety points outward (Figure 5D). It is thus expected that the trio domains will play important roles for binding of an ATP in a proper orientation. Also, reviewer 2’s concern about the mobility of sw-I was for the hairpin conformation in the ATP state, and it does not apply to the APO-state kinesin, whose sw-I is known to be mobile from previous experiments (e.g., it is invisible in PDB 4LNU, used for the APO simulation) as well as our present simulation. Longer simulation trajectory would of course be useful for providing more details (e.g., how ATP orients, or how its phosphate moiety binds to the P-loop), but the present simulation was sufficient to establish the basic physical picture. The role of L5 for ATP binding has previously been considered (Atherton et al.2014). However, to our knowledge, the functional roles of the mobile sw-I and α0, and the concerted role of the trio domains, are new in the present study. We also would like to add that, though our simulation times are too short for observing the full ATP binding event, their lengths are unprecedented.

Reviewer 2 questions the proposed Pirelease mechanism because of the possibility that the sw-I hairpin may have been made unstable in our simulations. As explained above, not only is this possibility unlikely, but Pirelease occurred with an intact hairpin (Figure 5E and Figure 5—figure supplement 2B). Unfolding of the sw-I hairpin is thus irrelevant. Since the hairpin is by itself mobile (high RMSD during the first 400 ns; Figure 2—figure supplement 1B), and since it maintains contacts with Pi, as sw-I moves, it is expected pull Pi.

Another concern by reviewer 2 is that the sub-microsecond release time of Piin our simulation may be too short compared to the longer lifetime of the ADP·Pistate indicated in recent experiments. But the divalent form of phosphate HPO42-;Pi2- did not release in our simulations, and it may persist in reality (subsection “Release of hydrolysis products is mediated by the mobile sw-I”, second paragraph). Even for the monovalent Pi, in our recent simulations (not included in our original submission), we found that the propensity for its release depends on its orientation. Author response image 1A reproduces the inset of Figure 5E, right before the release of Pi. In this case, a hydrogen atom of Pimakes contact with an Oβatom of ADP. Due to the low charge of this hydrogen (0.33e), the contact is readily broken, leading to the release of Pi. However, we recently found that Pican tumble, and its oxygen atom can form a contact with Mg2+. In this case, Pidid not release till the end of a 2-µs simulation (Author response image 1B). This is due to the stronger electrostatic interaction between the phosphate oxygen (−0.63e) and Mg2+. To estimate the release time, an ensemble of different release pathways must be considered, which is not the focus of the present study. Regardless of the rate of release, it is physically likely that Pirelease is assisted by the mobile sw-I.

Influence of Piorientation on its release.

(A) Re-rendering of Figure 5E, 622.80 ns. Pidoes not have a direct contact with Mg2+. (B) The last frame of a new 2-µs simulation. An oxygen atom of Pimaintains contact with Mg2+(marked by a bent line).

For the simulation of ADP release, we used sw-I in an α-helical conformation. Given its mobility, it is more plausible that it contributes to the release of ADP, for otherwise ADP will have to break its contacts with the P-loop spontaneously, which will be energetically less favorable. Overall, the dynamic roles of sw-I and other subdomains surrounding the ATP pocket make physical sense, and converse scenarios where sw-I remains static through these processes, are less probable.

Changes made:

Abstract: The total simulation time is updated, to account for the new 2-µs simulation mentioned in the second paragraph of the subsection “Release of hydrolysis products is mediated by the mobile sw-I”.

Subsection “Binding of ATP is mediated by the α0/L5/sw-I trio”, last paragraph: Physical plausibility of the dynamic role of the trio domains for ATP binding is explained.

Subsection “Release of hydrolysis products is mediated by the mobile sw-I”, second paragraph: The issue of the lifetime of the ADP+Pistate is mentioned with regard to the recent findings in the references that reviewer 2 mentioned. Citation for Mickolajczyk et al.(2015) is added.

Subsection “Release of hydrolysis products is mediated by the mobile sw-I”, second paragraph: The orientation-dependent lifetime of a bound Piis added. However, Author response image 1B above is not included, since the situation is similar to Figure 5—figure supplement 2D.

I was interested in the observation and discussion of strain in ATP and in kinesin's β sheet. However, these are not strongly emphasized in the manuscript and I am left with questions/concerns, particularly with regards to the β sheet strain. The model for strain is a second-order Gaussian surface, which is characterized by a single curvature parameter. It is not clear to me whether this can accurately capture the 'wrinkle' that appears in kinesin's β sheet when the two halves rotate with respect to each other. The free energy of the strain could therefore be significantly underestimated. For the ATP strain, it is not clear whether the magnitude of the observed strain (1.2kcal/mol) falls in the expected range, or how this would change our view of kinesin's mechanochemistry.

We first address reviewer 2’s question, whether the 1.2 kcal/mol torsional energy “falls in the expected range.” Its smallness compared to the ∼10 kcal/mol free energy of ATP hydrolysis indicates that torsional strain readily develops as ATP binds to kinesin. Yet, this primes the scissile bond for attack by the catalytic water (subsection “Kinesin-bound ATP is torsionally strained”, last paragraph). Thus, energetics of ATP hydrolysis is determined as a combined effect among ATP geometry, catalytic water coordination, configurations of the immediately surrounding residues as well as allosteric behaviors of remote parts of the protein. Thanks to reviewer 3’s suggestions, we have now expanded the discussion on both bond elongation and torsional strain in ATP.

Regarding the twist of the central β-sheet, instead of a Gaussian surface that reviewer 2 referred to (eax2+bxy+cy2), we used a quadratic surface (Equation 1). Perhaps confusion arose because Figure 2J was taken to depict the surface geometry. But as its axis labels and caption indicate, it plots distribution (free energy) of two curvature values. From this, it is also evident that, instead of “a single curvature parameter,” two curvature values are used, the mean and Gaussian curvatures. They respectively describe bending and twisting of curved surfaces in differential geometry. The ‘wrinkle’ that reviewer 2 describes is about twisting (“two halves rotate with respect to each other”), hence it is captured by Gaussian curvature. Indeed, Figure 2J and Figure 2—figure supplement 1K show that Gaussian curvature varies more between pre- and post-stroke states compared to mean curvature. To our knowledge, our work is the first to quantitatively measure the curvature energy of the central β-sheet. It turns out to be small between the pre- and post-stroke states. Although further studies are warranted, we believe that our analyses of the strain in ATP and the curvature of the central β-sheet clearly advance our understanding of kinesin’s mechanochemistry.

Changes made:

Subsection “Kinesin’s central β-sheet does not store enough energy to drive nucleotide processing”, second paragraph: Meaning of mean and Gaussian curvatures as descriptors of the bending and twisting of the β-sheet is explained.

We did not add any additional explanation about the curvature of the central β-sheet, since we believe existing explanations (subsection “Kinesin’s central β-sheet does not store enough energy to drive nucleotide processing”, second and last paragraphs; subsection “Concluding Discussion”, eighth paragraph) are sufficient.

Subsection “Kinesin-bound ATP is torsionally strained”, last paragraph: Implication of the smallness of the torsional energy in kinesin-bound ATP is explained.

Subsection “Kinesin-bound ATP is torsionally strained”, last paragraph: Comparison with the torsional strain of ATP in myosin, and additional discussion about its effect on ATP hydrolysis are provided.

In summary, I do not find the major claims of this paper, regarding the nucleotide pocket and the switch I loop, to be justified by the reported experiments – and indeed I think they likely arise from artifacts related to the simulation methods. The Discussion does not connect the simulation results to the most recent findings regarding kinesin's ADP•Pi state, even where there are seeming contradictions. Some of the other observations in the paper, particularly regarding strain in ATP and/or the β sheet, seem to have more merit. However, on the whole I think this work should be substantially revised and belongs in a more specialized journal.

Through our replies above, we believe that we have sufficiently addressed reviewer 2’s concerns about the mobility of sw-I, as well as the reviewer’s impression that experimental data are inconsistent with our results. We appreciate that reviewer 2 recognizes merits in our analyses of the strain in ATP and the central β-sheet. Given that the proposed mechanisms are novel, and that they have broad relevance to other kinesin families and nucleotide triphosphatases, we believe eLife is an ideal place to publish this work.

Reviewer #3:

In this molecular dynamics (MD) simulation study of kinesin, the authors conducted a systematic analysis of conformational dynamics of a single kinesin motor domain, either bound to a microtubule (MT) dimer unit or in solution; several nucleotide binding states (ATP, ADP/Pi, ADP and Apo) have also been examined. Compared to previous MD studies, the current work substantially extended the simulation time scale to multiple microseconds (for each state) by taking advantage of resources available on Anton. This is significant since as the authors showed that in most cases, the conformational properties of interest reach a plateau only after about one or a few microsecond(s). The extensive sampling is also essential to a meaningful estimate of the bending energy of the central β-sheet and an evaluation of its relevance to nucleotide processing.

The manuscript has been improved considerably and reviewers are satisfied with the way in which you have addressed the points raised during the previous round of review. All reviewers, including a new reviewer, agree that the work described in this paper is of high quality. There remain some technical concerns that will need to be addressed before we can reach a firm decision to accept the paper for publication in eLife. These technical points are outlined below:

1) It is stated that kinesin-bound ATP is torsionally strained. There is a comparison with ATP in solution to support this conclusion (For comparison, we performed a 4-ns simulation of an isolated Mg-ATP in water). The authors may wish to have a look at an improved parameterization of the dihedrals of ATP in the CHARMM force field (Komuro et al. JCTC vol 10, 4133-4142, 2014) to verify that this conclusion is robust or force field dependent.

Upon the reviewer’s suggestion, we used the modified ATP force field from Komuro’s paper and re-calculated the dihedral energy of ATP. In the table below, energies are in kcal/mol units (avg ± std):

Force Field

ATP-only (in solution)

ATP (bound to kinesin)

Original (param36)

28.48 ± 2.38

29.63 ± 1.85

Modified (Komuro)

28.11 ± 2.37

27.59 ± 1.86

Between ATP-only (isolated ATP) and ATP (kinesin-bound), the energy decreases slightly with the modified parameters whereas it increases with the original parameters. As explained in Komuro’s paper, the modified parameters render ATP more flexible, so that its conformation is more stable in certain structures with a bound ATP. However, they also state that the original parameters “have worked reasonably well in simulations performed on many ATP-bound proteins.” They also note: “There seems to be intrinsic difficulty in reproducing all the properties with existing force-field parameters, and different applications may need different adjustments.” Regardless of which parameters are more appropriate for our system, the conclusion that the dihedral energy change does not make a significant contribution, is still valid.

Change made:added:

“[…]torsional angles of the phosphate moiety should readily change when ATP binds to kinesin (calculation using a modified force field that worked well for certain ATP-bound protein structures (Komuro et al., 2014), also yielded only marginal changes in the dihedral energy).”

2) The significance of the PMF shown in Figure 2J is not entirely clear. It is stated that "To align PMFs for ATP and APO, we identified bins of the histogram that have nonzero counts in both simulations. We determined the constant free energy shift that needs to be added to the PMF for APO so that the mean-square difference of the two PMFs in the overlapping bins is minimized." The energy barriers on this combined free energy surface are also discussed ("There is also a Ì1.7 kB T energy barrier from ATP towards APO"). In practice, two separate PMFs were extracted from the marginal distributions (-kt*ln(rho)) obtained from two separate simulations that correspond to two very different molecular contexts (one with ATP bound and one without ATP bound). These two PMFs were combined by matching the histograms that have nonzero counts. But that is presuming that the two surfaces exist with respect to the same part of configuration space. Since one of them has a ligand bound and the other does not, why should this be so? Unless we misunderstand, the configurational sampling of the two state do not overlap, so how can the two PMFs be combined? In principle, the relative shift between two free energy surfaces, one ligand-bound and one APO, should depend on the concentration of the ligand. If the ligand concentration is zero, then the PMF of the ligand-bound state shifts to infinity, and likewise, if the ligand concentration is infinite then the PMF of the APO state shifts to infinity. This point needs to be addressed and clarified.

We thank the reviewer for the insightful comment. The reviewer correctly points out that Figure 2J is not the potential of mean force (PMF) from the simulation of an alchemical system where the Hamiltonian contains both ATP and APO states. The two curvature values as reaction coordinates, do not describe other conditions in individual simulations either (such as the presence or absence of ATP). However, as Figure 2—figure supplement 1K shows, PMFs for individual simulations are similar for a given conformational state of the motor head (pre- or poststroke), and they do not depend directly on the bound nucleotide (ATP or ADP+Pi) or whether or not kinesin is bound to the MT. The large overlap in the curvature distributions between the pre- and post-stroke states (Figure 2—figure supplement 1K) indicates that the central β-sheet cannot store sufficient energy for driving the motility cycle. In this regard, the ‘combined’ PMF in Figure 2J is merely a device that provides a semi-quantitative picture of free energy differences. Even if an alchemical simulation were performed (which would be difficult in practice), the free energies are unlikely to be very different, and our conclusion that the curvature energy does not change significantly between the pre- and post-stroke states, should hold.

Changes made:

Subsection “Kinesin’s central β-sheet does not store enough energy to drive nucleotide processing”, second paragraph: The abbreviation ‘PMF’ is introduced.

Subsection “Kinesin’s central β-sheet does not store enough energy to drive nucleotide processing”, last paragraph: The nature of the PMF in Figure 2J and its relation to Figure 2—figure supplement 1K, are explained.

3) In Figure 3E, the binding energies are of the order of -100 to -250 kcal/mol. These values appear to be very large, and are probably not standard binding free energies. The Materials and methods section states that "Calculation of the binding energy was done based on a previously developed method (Zoete et al., 2005)." This method is MM/GB-SA (molecular mechanics-generalized Born surface area) which is an end-point approximation to the binding free energy. There is scepticism that the binding energy estimated by the Generalized Born method is sufficiently accurate to draw meaningful conclusions for the system being considered. How sensitive are the conclusions to this analysis?

As explained in Zoete’s paper, rather than the issue of the accuracy of the Generalized Born method, the binding energy between proteins without considering entropy tends to be larger than the experimental binding free energy. The binding free energy consists of respectively large enthalpy and entropy terms that subtract to yield experimental values. In Zoete’s work, the binding energy for an insulin dimer was −39 kcal/mol, where an insulin monomer is 51 residues in size. As a kinesin motor head is more than 6 times larger, the calculated binding energy at about −240 kcal/mol, as in Figure 3E in our previous submission, is reasonable. However, although the differences in the binding energies are consistent with experimental observations and with our other analyses regarding the kinesin-MT interface (subsection “Kinesin-MT interface is dynamic and hydrated”, last paragraph and subsection “Kinesin-MT interface is dynamic and hydrated”, first paragraph), we agree that only considering the binding energy (enthalpy) is not sufficient for a quantitative comparison with experiment. We thus explained limitations of this calculation and moved the figure to Figure 3—figure supplement 1L. Given the large size of the system, calculation of the full binding free energy would be the subject of another project.

“However, our binding energies do not include water-mediated interactions and entropic contributions, which are expected to be comparable to the binding energy in magnitude (Zoete et al., 2005), so that the net binding free energy is much smaller than those in Figure 3—figure supplement 1L. Thus, the calculated binding energies, although they reflect the interaction between kinesin and the MT in different nucleotide states, do not correspond quantitatively to the experimental binding affinities.”

National Institutes of Health (R01GM087677)

Pittsburgh Supercomputing Center (Anton Supercomputer)

Texas Advanced Computing Center

Texas A&M Supercomputing Facility

Texas A and M University (Open Access to Knowledge Fund)

Wonmuk Hwang

CHARMM Development Project

Martin Karplus

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

Acknowledgements

This work was supported in part by the National Institutes of Health (NIH) grant R01GM087677. The work done at Harvard University was supported by the CHARMM Development Project. Anton/Anton-2 computer time was provided by the Pittsburgh Supercomputing Center (PSC) through the NIH grant R01GM116961. Anton at PSC was generously made available by DE Shaw Research. We also used machines at the Texas A&M Supercomputing Facility and Texas Advanced Computing Center at UT Austin. We thank the reviewers for their helpful comments. The open access publishing fees for this article have been covered by the Texas A&M University Open Access to Knowledge Fund (OAKFund), supported by the University Libraries and the Office of the Vice President for Research.

eLife is a non-profit organisation inspired by research funders and led by scientists. Our mission is to help scientists accelerate discovery by operating a platform for research communication that encourages and recognises the most responsible behaviours in science.eLife Sciences Publications, Ltd is a limited liability non-profit non-stock corporation incorporated in the State of Delaware, USA, with company number 5030732, and is registered in the UK with company number FC030576 and branch number BR015634 at the address:
eLife Sciences Publications, Ltd
1st Floor, 24 Hills Road
Cambridge CB2 1JP
UK