Motor-Like Properties of Non-Motor Enzymes

Molecular motors are thought to generate force and directional motion via nonequilibrium switching between energy surfaces. Because all enzymes can undergo such switching, we hypothesized that the ability to generate rotary motion and torque is not unique to highly adapted biological motor proteins, but is instead a common feature of enzymes. We used molecular dynamics simulations to compute energy surfaces for hundreds of torsions in three enzymes, adenosine kinase, protein kinase A, and HIV-1 protease, and used these energy surfaces within a kinetic model that accounts for intersurface switching and intrasurface probability flows. When substrate is out of equilibrium with product, we find computed torsion rotation rates up ~140 cycle s−1, with stall torques up to ~2 kcal mol−1 cycle−1, and power outputs up to ~50 kcal mol−1 s−1. We argue that these enzymes are instances of a general phenomenon of directional probability flows on asymmetric energy surfaces for systems out of equilibrium. Thus cyclic probability fluxes, corresponding to rotations of torsions and higher-order collective variables, are expected in any chiral molecule driven between states in a non-equilibrium manner. The same considerations apply to synthetic chiral molecules switched in a nonequilibrium manner between energy surfaces by light, redox chemistry, or catalysis.


Text Introduction
A biological molecular motor is an enzyme that uses the free energy of an out-ofequilibrium chemical reaction to drive mechanical motion. This motion must have a specific direction to fulfill the motor's functional role; for example, a helical flagellum must rotate in the appropriate sense to propel the organism. The ability to generate a directional cycle of conformational changes appears to be a highly specialized protein property. However, it has been argued that a key mechanistic aspect of biomolecular motors is that they are switched in a nonequilibrium manner between distinct energy surfaces (1)(2)(3)(4)(5)(6)(7)(8), and non-motor enzymes also switch between at least two conformational free energy surfaces as they bind and release substrate. It is thus of interest to ask whether enzymes not normally thought of as motor proteins have motor-like properties when operating away from chemical equilibrium.

Fig. 1.
Illustrative apo and substrate-bound free energy surfaces of an enzyme (solid green and blue lines, respectively) and the flow of probability along coordinate θ, after the effective potential changes due to substrate binding (green to blue probability distribution curves).
Consider a periodic molecular coordinate, θ, such as a bond torsion or a higherorder collective variable, in an enzyme undergoing a catalytic cycle of substrate-binding, catalysis, and product release. When binding of substrate switches the effective energy surface from that of the apo state, μ0(θ), to that of the bound state, μ1(θ), probability flows along θ towards a new equilibrium distribution corresponding to μ1(θ) (Figure 1).
Such a flow of probability can be modeled with the Fokker-Planck equation, for example (4). Because enzymes are chiral, neither energy landscape, μ0(θ) and μ1(θ), can possess mirror symmetry; i.e., neither can be an exactly even function around any value of θ. Thus, there will be more probability flow in one direction than the other, due to the asymmetry of the starting probability distribution and of the patterns of barriers and troughs along the bound free energy surface. After the substrate is catalytically removed, this process repeats itself, with a probability distribution that has approached equilibrium on the substrate-bound surface diffusing asymmetrically under the potential of the apo free energy surface, μ0(θ). Unless there is an infinitely high energy barrier at the same location on both surfaces, there is no basis to expect that the net flows will cancel identically for a chiral molecule, and the resulting asymmetric net flow of probability during this two-step process constitutes directional motion. Note that any such motion will be superimposed on a background of rapid, non-directional stochastic fluctuations in θ.
As this reasoning is general in nature, we hypothesized that all periodic coordinates of any enzyme catalyzing an out-of-equilibrium reaction must undergo some degree of directional rotations. In order to uncover these rotations and begin to assess their magnitudes, we used molecular simulations and kinetic analysis to analyze the rotational dynamics of bond torsions induced by nonequilibrium, stochastic, stateswitching in three, diverse, non-motor enzymes: adenosine kinase (ADK), protein kinase A (PKA), and HIV-1 protease (HIVP).

Materials and Methods
The one-dimensional free energy surfaces of protein main-and side-chain torsions, discretized into bins, were obtained from equilibrium molecular dynamics (MD) simulations of the enzymes in their apo and substrate-bound states (Supplementary Methods). These data, coupled with literature values for the enzyme kinetic parameters (Table S1), enabled us to define first order rate constants for transitions along and between the free energy surfaces (Fig. S1). The resulting set of rate equations was solved for the non-equilibrium steady state probability distribution, for a given concentration of substrate, and this, in turn, was used to compute the probability flux on each surface. The net flux, , in units of torsional rotational cycles per second, is an indication of directional rotation; e.g., a positive value implies clockwise rotation of the torsion.

Results
For all three enzymes, we find that multiple torsion angles undergo directional rotations, as indicated by nonzero probability flux, when excess substrate is present. Thus, at high substrate concentration, about 40 torsions in ADK and PKA rotate faster than 10 cycle s -1 , and about 140 rotate faster than 1 cycle s -1 (Fig. 2a). Directional rotation is also observed for HIVP, though the rates are lower (Fig. 2b, red). The lower rates largely reflects the lower kcat value of HIVP, relative to PKA and ADK. Thus, artificially assigning cat = 200 s −1 to HIVP leads to substantial increases in the number of torsions with fluxes of at least 10 s -1 and at least 1 s -1 (Fig. 2b, orange and Fig. S3). The tendency toward lower fluxes in HIVP may also reflect the smaller scale of its conformational changes (Fig. S2): a small conformational change may lead to more similar energy surfaces in the two states, and hence less opportunity to generate rotational flux by the mechanisms discussed below.
Although the maximum rotation rates are different for ADK, PKA and HIVP (180,    The mechanism by which high rates of directional flux are generated is illustrated by the 2 torsion of ADK Thr175 (Fig. 4a). This angle has a two-peaked probability distribution in both the bound and apo states, but the peak near + 2 is favored in the apo state, while that near − 2 is favored in the bound state (Fig. 4b,c). In the presence of substrate, the bound-state energy minimum near − 2 is highly occupied (Fig. 4b,c).
Catalytic breakdown of substrate pumps the system to the secondary energy minimum of the apo state at − 2 ( This process parallels fluctuating potential mechanisms previously invoked to explain highly evolved molecular motors (3,4,7,8,(20)(21)(22)(23). Intuitively, high flux can be generated when each energy surface (apo and bound) has a main energy barrier and a main energy well, and the two surfaces are offset, so that probability flow in a given direction on each surface bypasses the barrier on the other surface, as schematized in There are also many torsions whose net directional flux is small (< 1 cycle s -1 ), and whose dynamics are instead dominated by reciprocating motion, in which probability flows in one direction in the bound state and in the other direction in the apo state, within some angular range. These reciprocating probability fluxes correspond to driven, oar-like motion of the atoms controlled by the torsion. The enzymes ADK and PKA, respectively, are predicted to have ~1250 and ~750 torsions with reciprocating motions at rates of at least 1 cycle s -1 (Fig. 2e,f), and minimal directional flux. The maximal reciprocating fluxes are greater than the maximal directional fluxes, and, for ADK and PKA, approach the catalytic rates (Fig. S4). The mechanism of reciprocating motion is illustrated by 2 of Asn 138 in ADK (Fig. 4e), which has a net flux that is nonzero but well below 1 cycle s -1 , yet has reciprocating fluxes of up to 130 cycles s -1 ( Fig. 4e-h). Reciprocating flux occurs when the main energy wells on the two surfaces are in different locations, but little probability flow on each surface bypasses the main barrier on the other surface, as schematized in Fig. S5 c-d. Thus, if the main barriers in the two states are coincident, movement will be dominated by reciprocating motion.

Discussion
We have provided reasoning and evidence that directional rotation that can work against a load is not limited to highly adapted biological motor proteins. Instead, any enzyme operating out of equilibrium should generate directional conformational fluxes.
The magnitudes of these motions will depend on a number of factors, including the levels and locations of the barriers and energy wells on the two surfaces, and the chemical driving forces and kinetic constants. We also observe high reciprocating flux in torsions with energy surfaces that are not well tuned to generate net directional motion.
The present modeling approach uses atomistic modeling to obtain the effective potentials, and then embeds these in a kinetic model. This efficient technique is more broadly applicable, and could be used to study not only proteins but also synthetic molecular motors. This approach can also be extended to account for greater complexity, such as a position-dependent diffusion constant and additional energy surfaces. Although we have focused on torsional motions, the same reasoning applies to motions through higher-order conformational subspaces, including larger-scale protein domain motions, because they, too, will have asymmetric potentials, due to the chirality of these molecules. Such motions might help explain the observation that the translational diffusion coefficients of some enzymes rise with substrate concentration and hence with enzyme velocity (24,25). The present results also indicate that even enzymes from earliest evolutionary time would have had the ability to generate directional motion, and thus would have had a good starting point from which to embark on an evolutionary path to today's highly adapted motor proteins. Using fewer bins would effectively smooth the surfaces, reducing the differences between apo and bound and thus decreasing the computed fluxes and torques.

Acknowledgments
The kinetic model allows transitions between neighboring bins in the same state, with forward and reverse rate constants , → +1 and , +1→ respectively, where = 0 for the apo state and = 1 for the bound state; and transitions between states in the same bin, with rate constants 0→1, and 1→0, for binding and dissociation, respectively.
These intersurface rate constants are set to where on is the macroscopic on-rate for binding of substrate to the enzyme and is the substrate concentration; and where Δ , the binding free energy associated with bin , is given by where , is the free energy associated with bin on surface : We now show that this set of rate constants leads to the correct dissociation constant, Kd, for the full system. In particular, because d = 1 0 ⁄ , where C1 and C0 are the concentrations of bound and free protein, respectively, we want Using the prior equations, we obtain The forward and reverse rate constants are equal to if the energy difference between the two bins is zero, and any energy difference increases and decreases the forward and reverse rate constants, respectively, by the same factor, as expected on physical grounds. The value of is discussed in Assignment of Numerical Parameters, below.
Solving the first order kinetic system described above would yield an equilibrium Boltzmann probability distribution across the bins and states, , eq ∝ exp (− , ) with the following normalization: So far, we have included only one mechanism for transitioning from the bound state to the apo state: dissociative release of substrate. However, there is a second mechanism, which is the one that can put the system out of equilibrium, namely catalytic degradation of the substrate followed by product release. Referring to Eq (2), we account for this by adding another set of first-order rate processes going from the bound state bins to the apo state bins, with rate constants of cat : As noted above, for simplicity, we assume that the concentration of product is negligibly small. When and cat are nonzero, the system is no longer at equilibrium, but instead goes to a non-equilibrium, steady-state probability distribution, , ss . This distribution is obtained by the numerical method described in the following subsection.
Steady-state solution of the kinetic model The net probability flux is the sum across both surfaces, At steady state, is uniform across bins, and we write ≡ in the main text as "directional flux". We also report "reciprocating flux" for each torsion as the peak magnitude of directional flux across either surface, = max(| 0, |, | 1, |).

Applying a load (torque) to a torsion angle
One test of these motors is their ability to do work against a mechanical load. We apply load by tilting the free energy surfaces with a constant slope, corresponding to a torque, τ, that opposes the direction of flux. This is done by supplementing the chemical potential differences in equation (4) on both energy surfaces. That is, the energy difference between adjacent bins used to compute the transition rates (equation (9)) are modified by a factor ±Δ , depending on direction, where Δ > 0 gives a negative torque that is used to oppose positive probability flux, and vice versa. So that the load is correctly imposed across the periodic boundary, we replace ,1 − , with ,1 − , + Δ and , − ,1 with , − ,1 − Δ . The stall torque is the load that brings the next flux to zero; increasing the torque beyond this limit will cause the motor to slip backwards. We identify the stall torque by solving the kinetic system iteratively for a systematic scan of applied torques. The power produced for a given torque is given by = , where is the flux computed for a given applied torque τ. Empirically, the maximum power is found to occur at half the stall torque.

Assignment of numerical parameters
The present model requires numerical values of , on , d = 1/ a , and cat . An initial estimate of was made by simulating a molecule of butane, with the force field torsion and nonbonded terms set to zero to create a flat energy surface, and determining the rotational diffusion constant for this barrierless torsion. The mean value over one hundred 1 ns Langevin dynamics simulations is 3×10 15 degree 2 s -1 . However, this large value led to numerical instability during diagonalization of the transition matrix, so we tested how the computed fluxes vary with D, and found that they are insensitive to the value of once it is above ∼ 10 9 degree 2 s -1 . We used = 3 × 10 12 degree 2 s -1 for the present calculations, as this is well into the regime where the results are independent of , but not so large as to cause numerical problems (Fig. S7). The enzyme kinetic parameters and their basis in the experimental literature are provided in Table S1.

Torsion potentials of mean force from molecular dynamics simulations
Adenylate Kinase PDB accession 4AKE (28) was used as a starting structure for apo ADK; crystallographically ordered water molecules were retained. Hydrogens were added with pdb2pqr (29) at pH 7.0, bringing the net charge of the protein to -4. The substratebound protein was similarly modeled using PDB accession 3HPQ (30), which includes the ligand AP5A, a transition state analog, bound to the active site. AP5A carries a charge of -5 on the five phosphate groups, and the protein again carries a charge of -4 in the bound state. Partial atomic charges for AP5A were determined using the AM1-BCC method (31) in the antechamber program, and the remaining force field parameters were assigned from GAFF (32,33). The apo and substrate-bound simulation systems were neutralized with 4 and 9 sodium ions, respectively. Both protein structures were solvated in a truncated octahedron with 12 Å padding.
Each system was energy-minimized for 20,000 steps, thermalized to 300 K over 1 ns, and equilibrated for 100 ns. Production simulations were then carried out for 1.0 μs using PME electrostatics with a 9 Å cutoff, and hydrogen mass repartitioning (34) with a 4 fs time step, using pmemd.cuda.MPI module of Amber 16 (35). Histograms (60 bins) of the torsion types listed in Table S2 were computed over the entire production simulation, using the cpptraj module (36) of Amber 16.
HIV Protease PDB accession 1HHP (37) was used to model the apo structure of HIV-1 protease; crystallographically ordered waters were retained. To be consistent with Uniprot P03367, we made the following computational mutations: K14R, S37N, R41K, L63P, and I64V. Hydrogens were added with pdb2pqr at pH 7.0, bringing the net charge of the protein to +4. The substrate-bound protein was modeled using PDB accession 1KJF (38) with 10 residues of the Gag protein co-crystallized (RPGNFLQSRP; residues 443-452; fragment p1-p6 or SP2-p6). The following mutations were made to make the apo and bound structures consistent: K7Q and N25D. The charge of the bound protein and peptide was +6. Both models were solvated in a truncated octahedron with 12 Å padding and 4 or 6 chloride ions, respectively, added to neutralize the charge. The simulation procedures were identical to those for ADK.

Protein Kinase A
Simulations were carried out as described in reference (14).

Numerical convergence
Convergence of the results was tested by using increasing amounts of simulation time to derive apo and bound state population histograms and then predicting the directional flux based on those histograms. For a sample angle shown in Fig. S8, after 500 ns, the calculated directional flux is very close to its value at 1000 ns, and the uncertainty in the flux value drops from about ±60 cycle s -1 to about ±1 cycle s -1 going from 500 ns to 1000 ns.  conformations of ADK, rendered as a surface, with substrate AP5A (red) as spheres. In the right two panels, the absolute magnitude of directional flux is mapped onto the apo structure, using a color gradient with thresholds from < 1 cycle s -1 (blue) to > 10 cycle s -1 (red). Middle row: the apo (green; PDB accession 1HHP) and substrate-bound (blue; PDB accession 1KJF) conformations of HIVP, rendered as a surface, with substrate Gag peptide (red) as spheres. In the right two panels, the absolute magnitude of directional flux is mapped onto the apo structure, using the same color gradient, with thresholds < 0.1 cycle s -1 (blue) to > 2 cycle s -1 (red), to account for the lower level of directional flux in HIVP, even at kcat = 200 s -1 . Bottom row: the apo (green; PDB accession 1CMK) and substrate-bound (blue; PDB accession 3FJQ), with ATP (red) as spheres. In the right two panels, the absolute magnitude of directional flux is mapped onto the apo structure, using the same thresholds as for ADK.     Fig. 3c). The probability flux drawn separately for each surface and as a sum (grey points), showing exactly zero directional flux due the symmetry of the two energy surfaces (cf. Fig. 3d).    (39).1xA9k [c] The kinetics of PKA are complex and depend on, among other factors, the presence of divalent ion species and occupancy in the active site. The observed catalytic rate includes two or more conformational changes and lies between 50 s -1 and a fast, burstphase of 500 s -1 . We used the rate-limiting step (ADP · Pi release) from Figure 11 of (13) in this manuscript.
[d] An average of Kd values for ATP found in (12) and (40). The calculated offset is an average of the offset for each Kd. The Kd values for AMP are roughly four times as large. Our model implicitly assumes a single substrate and single on rate.
[g] A literature value of Kd for the specific Gag fragment used in the simulations was not found in the literature. However, reference (41) reports KM = 5.3 × 10 -4 M for the p1/p6 substrate sequence PGNFLQS (the simulated bound peptide sequence was RPGNFLQSRP) and for sufficiently small catalytic rate, d~M . Several other references list values for kcat / KM.
[h] No value could be found in the literature, so the same order of magnitude as ADK is used.
[i] Catalytic rates as low as 0.3 s -1 are reported in (41) for the p1/p6 substrate sequence, but other Gag sequences have catalytic rates on the order of 10 s -1 , see for example (15)(16)(17).  N  PRO  chi2  CD  CG  CB  CA  PRO  chi3  N  CD  CG  CB  SER  phi  C  N  CA  C  SER  psi  N  CA  C  N  SER  chi1  OG  CB  CA  N  THR  phi  C  N  CA  C  THR  psi  N  CA  C  N  THR  chi1  CG2  CB  CA  N  THR  chi2  HG1  OG1  CB  CA  TRP  phi  C  N  CA  C  TRP  psi  N  CA  C  N  TRP  chi1  CG  CB  CA  N  TRP  chi2  CD1  CG  CB  CA  TYR  phi  C  N  CA  C  TYR  psi  N  CA  C  N  TYR  chi1  CG  CB  CA  N  TYR  chi2  CD1  CG  CB  CA  TYR  chi3  HH  OH  CZ  CE1  VAL  phi  C  N  CA  C  VAL  psi  N  CA  C  N  VAL  chi1  CG1  CB  CA  N