Breaking time-reversal symmetry for ratchet models of molecular machines

TL;DR

This paper explores how to induce time-reversal symmetry breaking in molecular machine models using deterministic protocols over sawtooth potentials, revealing conditions that maximize asymmetry and contrasting externally driven versus autonomous systems.

Contribution

It demonstrates how deterministic driving protocols can break time-reversal symmetry in ratchet models, highlighting the role of potential symmetry and driving parameters in molecular machines.

Findings

Symmetric sawtooth potentials maximize time asymmetry under deterministic driving.

Higher barriers and intermediate width driving potentials increase time asymmetry.

Externally driven systems differ from autonomous ones in symmetry requirements for directed motion.

Abstract

Biomolecular machines transduce free energy from one form to another to fulfill many important roles inside cells, with dissipation required to achieve directed progress. We investigate how to break time-reversal symmetry at a given dissipation cost by using deterministic protocols to drive systems over sawtooth potentials, which have frequently been used to model molecular machines as ratchets. Time asymmetry increases for sawtooth potentials with higher barriers and for driving potentials of intermediate width. For systems driven over a sawtooth potential according to a protocol, we find that symmetric sawtooths maximize time asymmetry, while earlier work examining ratchet models of molecular machines required asymmetric sawtooth potentials to achieve directed behavior. This distinction arises because deterministically driven machines are externally provided with direction, while autonomous machines must generate directed behavior.

Figures

Equations

$$A$$

$$+\frac{1}{2}{⟨\ln \frac{2}{1+\exp (-\beta W[\tilde{X}|\tilde{\Lambda }]-\beta \Delta F)}⟩}_{\tilde{\Lambda }}.$$

$$h[\Lambda ]\equiv \frac{1}{2}\left(\beta ⟨W[X|\Lambda ]{⟩}_{\Lambda }+\beta ⟨W[\tilde{X}|\tilde{\Lambda }]{⟩}_{\tilde{\Lambda }}\right).$$

$$E_{\text{r}}(x)=\{\begin{array}{ll}\Delta E_{\text{r}}\frac{x\text{ mod }\ell }{{\ell }_1},& \text{if }0<x\text{ mod }\ell \le {\ell }_1\\ \Delta E_{\text{r}}\frac{\ell -x\text{ mod }\ell }{{\ell }_2},& \text{if }{\ell }_1<x\text{ mod }\ell \le \ell \end{array},$$

$$\Delta E^{‡}(x_{\text{trap}})$$

$$=\left[\Delta E_{\text{r}}+\frac{1}{2}k{\left(x_{\text{trap}}-\frac{\ell }{2}\right)}^2\right]-\frac{1}{2}k{x}_{\text{trap}}^2$$

$$=\Delta E_{\text{r}}-\frac{1}{2}k\ell \left(x_{\text{trap}}-\frac{\ell }{4}\right).$$

$$t_{\text{cross}}$$

$$=\frac{k{\ell }^4}{16D}\frac{e^{\Delta E_{\text{r}}+k{\ell }^2/8}}{e^{k{\ell }^2/2}-1}.$$

$$E^{\prime }=\frac{1}{2}{k}^{\prime }(x^{\prime }-x_{\text{trap}}^{\prime }{)}^2+\frac{\Delta E_{\text{r}}^{\prime }{x}^{\prime }}{{\ell }_1^{\prime }},$$

$$\frac{E^{\prime }}{H}=\frac{1}{2}\frac{k^{\prime }{\ell }^2}{H}\frac{(x^{\prime }-x_{\text{trap}}^{\prime }{)}^2}{{\ell }^2}+\frac{\Delta E_{\text{r}}^{\prime }}{H}\frac{x^{\prime }}{\ell }\frac{\ell }{{\ell }_1^{\prime }}.$$

$$\Delta E_{\text{r}}=0,$$

$$⟨W{⟩}_{\text{max}}=\frac{1}{Z}{\int }_{-\infty }^{\infty }\Delta E_{\text{trap}}(x)e^{-E[x,x_{\text{trap}}=0]}\mathrm{d}x,$$

Full text

††preprint: APS/123-QED

I Introduction

Biomolecular machines perform a variety of essential tasks inside living cells Alberts (1998), and inspire the development of synthetic molecular machines Zhang et al. (2018). Molecular machines operate out of equilibrium and stochastically Seifert (2012), buffeted by fluctuations at the nanoscale Brown and Sivak (2017a). As with macroscopic machines, molecular machines can be designed to optimize efficiency den Broeck et al. (2012); Schmiedl and Seifert (2008) and speed Wagoner and Dill (2016); Brown and Sivak (2017b, 2018). However, unlike their macroscopic counterparts, design can meaningfully impact the precision of molecular machine operation Barato and Seifert (2015).

Molecular machines can pay a significant fraction of their free energy budget for directed progress Machta (2015); Ariga et al. (2018). To quantify directional behavior, Feng and Crooks Feng and Crooks (2008) proposed time-reversal asymmetry, expressed as the Jensen-Shannon divergence between trajectory distributions for processes driven in forward and reverse directions.

Earlier investigations (unfolding/refolding of ribonucleic acid hairpins in experiment Feng and Crooks (2008) and alanine decapeptide in molecular dynamics simulations Procacci and Marsili (2010)) found the trade-off between time asymmetry (seen as a potential benefit) and dissipation (a cost) remained near the prediction of linear-response theory Marconi (2008). Linear-response behavior for time asymmetry and dissipation is expected for slow processes that remain near equilibrium. For a range of dissipation (approximately 1-8 $k_{\text{B}}T$), there is a substantial gap between the linear-response prediction for time asymmetry and the maximum possible time asymmetry, while for lower and higher dissipation the difference between linear-response and the maximum possible time asymmetry is small. This gap (at medium dissipation) presents an opportunity for significant improvement on linear-response behavior Brown and Sivak (2016). We aim to determine characteristics that allow time asymmetry to increase for a given dissipation cost, or equivalently decrease dissipation for a given time asymmetry.

We investigate the trade-off between time asymmetry and dissipation in the context of molecular machines, using a sawtooth potential that is frequently used to represent a ratchet model of a molecular machine Astumian and Bier (1994); Julicher et al. (1997); Astumian and Derenyi (1998); Parmeggiani et al. (1999); Lipowsky and Harms (2000); Reimann (2002); Suzuki and Munakata (2003); Gabrys et al. (2004); Linke et al. (2005); Lacoste and Mallick (2009); Levien and Bressloff (2015); Kedem et al. (2017); Saakian and Klumper (2018). ‘Flashing’ ratchet models, which switch between multiple potential energy landscapes, have been used to generate autonomous directed motion Suzuki and Munakata (2003). While we use the sawtooth landscape from flashing ratchets, we apply an externally directed protocol (in contrast to autonomous behavior) to produce directed motion. A protocol typically represents the influence of an experimental apparatus, providing a suitable description of externally driven synthetic molecular machines, with more limited direct application to autonomous biomolecular machines. Understanding how molecular machines can achieve higher time asymmetry for a given dissipation cost may point towards design principles for molecular machines. These principles could inform the design of more efficient synthetic molecular machines.

In our previous work Brown and Sivak (2016), we explored how to manipulate time asymmetry for a step potential representing energy transmission, using a harmonic trap to move the system from one side of the step to the other. Intermediate step heights and protocols remaining near the step led to time asymmetries that improved most on linear response.

In contrast to step potentials, we find that high sawtooth potentials provide greater improvement over linear response. Changing from a symmetric to an asymmetric sawtooth causes time asymmetry to decrease, stemming from the greater overlap in initial equilibrium distributions for more asymmetric sawtooths.

II Methods

Time-reversal asymmetry $A$ (here measured in nats, $1/\ln 2$ bits) quantifies the expected information gain about whether a single observed trajectory was produced by a particular nonequilibrium process or its time reversal, assuming a 50% prior probability of the forward or reverse direction Crooks and Sivak (2011). If the forward and reverse trajectory distributions are identical, then any observation provides no information on the protocol direction, defining a minimum $A=0$. If the forward and reverse trajectory distributions are completely distinct, then any observed trajectory can be definitively assigned to the forward or reverse process, defining a maximum $A=\ln 2$ nats, amounting to one bit of information.

The Crooks fluctuation theorem Crooks (1999) quantitatively describes the relative probabilities of forward and reverse trajectories solely in terms of the work required for each trajectory. This permits expression of the time asymmetry $A$ in terms of the forward and reverse work distributions Feng and Crooks (2008),

[TABLE]

$\Lambda$ labels a protocol, the time schedule of a controllable parameter $\lambda(t)$ over $t\in[0,\Delta t]$ for a protocol of duration $\Delta t$. The reverse protocol $\tilde{\Lambda}$ is the time-reversal of $\Lambda$, i.e. $\lambda(\Delta t-t)$. $W[X|\Lambda]$ is the work done during protocol $\Lambda$ when the system follows trajectory $X$. $\Delta F$ is the free energy change during protocol $\Lambda$, and $\beta\equiv(k_{\text{B}}T)^{-1}$, for Boltzmann’s constant $k_{\text{B}}$ and temperature $T$.

Dissipation $h$ is the average excess work for an equal weighting of forward and reverse protocols, and hence equals the average work because the equilibrium free energy changes during $\Lambda$ and $\tilde{\Lambda}$ cancel Feng and Crooks (2008):

[TABLE]

Time asymmetry is mathematically upper bounded by $A\leq h/4$ and $A\leq\ln[2/(1+e^{-h})]$ Feng and Crooks (2008); Taneja (2005). A system which remains near equilibrium responds linearly to small perturbations Chandler (1987) – for a driven ensemble in this linear-response regime, the work distribution is Gaussian, and the Crooks fluctuation theorem imposes a connection between the mean and variance of the work Feng and Crooks (2008); Speck and Seifert (2004), $\langle W\rangle-\Delta F=\frac{1}{2}\sigma_{W}^{2}$.

We investigate time asymmetry and dissipation with an overdamped system of diffusivity $D$ on a potential landscape $E[x,x_{\text{trap}}(t)]=E_{\text{trap}}[x,x_{\text{trap}}(t)]+E_{\text{r}}(x)$. The first component is a quadratic trap, $E_{\text{trap}}[x,x_{\text{trap}}(t)]=\frac{1}{2}k[x-x_{\text{trap}}(t)]^{2}$, with a time-dependent minimum $\lambda(t)=x_{\text{trap}}(t)$ as the control parameter. The second component $E_{\text{r}}(x)$ is a time-independent periodic sawtooth (ratchet) potential (Fig. 1) that enables directional motion for simple models of molecular machines Astumian and Bier (1994); Julicher et al. (1997); Astumian and Derenyi (1998); Parmeggiani et al. (1999); Lipowsky and Harms (2000); Reimann (2002); Suzuki and Munakata (2003); Gabrys et al. (2004); Linke et al. (2005); Lacoste and Mallick (2009); Levien and Bressloff (2015); Kedem et al. (2017); Saakian and Klumper (2018),

[TABLE]

with $\ell=\ell_{1}+\ell_{2}$ the wavelength over which the sawtooth repeats, and $\Delta E_{\text{r}}$ the barrier height. The periodicity of $E_{\text{r}}(x)$ represents a cyclical molecular machine, with $\Delta F=0$ over a cycle. $E$, $k$, $x$, $\Delta E_{\text{r}}$, and $\Delta t$ are nondimensionalized by energy scale $k_{\text{B}}T$, length scale $\ell$, and a timescale such that system diffusivity $D=1$ (see Appendix A for details).

For each protocol, the trap translates at a constant velocity $u=1/\Delta t$ from $x_{\text{trap}}=0$ to $x_{\text{trap}}=1$ for the forward protocol, and in the opposite direction for the reverse protocol. The system begins in equilibrium and diffuses as the trap translates. Appendix B provides further simulation details.

For a given spring constant $k$, ratchet height $\Delta E_{\text{r}}$, sawtooth lengths $\ell_{1}$ and $\ell_{2}$, and trap speed $u$, work distributions are compiled from many repetitions of forward and reverse protocols, from which a single time asymmetry [Eq. (II)] and dissipation [Eq. (2)] are calculated.

III Results

III.1 Symmetric ratchets

We first study symmetric sawtooth potentials, $\ell_{1}=\ell_{2}=0.5$. Fig. 2a displays curves of time asymmetry vs. dissipation, with each curve generated by parametrically varying the trap speed $u$. As trap speed $u$ increases, both the time asymmetry and the dissipation increase in a manner that depends on the barrier height $\Delta E_{\text{r}}$ and spring constant $k$. For each curve, both the time asymmetry and dissipation are negligible for low trap speed $u$ (bottom left of Fig. 2a panels). A slow-moving trap allows the system to stay near equilibrium, leading to forward and reverse work distributions with similar shapes, widths, and mean work magnitudes (Fig. 2b, far left panel). These near-equilibrium trajectories are well-described by linear response, hence staying near the linear-response time asymmetry-dissipation curve.

As $u$ increases (left to right in Fig. 2b), the forward and reverse work distributions become non-Gaussian and more separated, leading to less overlap and causing both time asymmetry and dissipation to increase. For these higher trap speeds, the system departs from linear response (highest dissipations in each panel of Fig. 2a).

At higher trap speeds, individual system trajectories can be roughly categorized as either crossing from the initial side to the far side of the sawtooth potential barrier, or remaining on the initial side without crossing the barrier. An activated process Grabert (1988) occurs at a typical rate $k_{\text{a}}=\tfrac{2D}{(\Delta x^{\ddagger})^{2}}e^{-\Delta E^{\ddagger}}$ for respective distance $\Delta x^{\ddagger}=\ell/2$ and energy difference $\Delta E^{\ddagger}$ from the minimum to the barrier. For $k\leq 2\Delta E_{\text{r}}/\ell^{2}$ the energy minimum on the initial side of the barrier remains at $x=0$ and the energy barrier remains at $x=\ell/2$ regardless of $x_{\text{trap}}$, but the barrier height varies as the trap moves, giving

[TABLE]

For average crossing rate $\langle k_{\text{cross}}\rangle=\tfrac{1}{\ell}\,\int_{0}^{\ell}k_{\text{a}}(x_{\text{trap}})\,\text{d}x_{\text{trap}}$ this gives the barrier crossing timescale for the protocol as a whole:

[TABLE]

The leftmost panel of Fig. 2b ($u=7.81\times 10^{-2}$) corresponds to $t_{\text{cross}}/\Delta t\approx 0.11$, allowing the system to reliably cross the barrier during the protocol duration, leading to a single component in the work distributions.

Faster driving ($u=0.625$ in Fig. 2b) corresponds to $t_{\text{cross}}/\Delta t\approx 0.86$. At this driving speed the system does not reliably cross the barrier during the protocol duration, and the work distributions contain two distinct components, corresponding to two trajectory categories. The wide ‘crossing’ components in Fig. 2, at low work magnitude, represent systems that cross to the far side of the barrier, thus staying closer to $x_{\text{trap}}$, such that the trap does less work on the system. The narrow ‘non-crossing’ peaks, at higher work magnitude, are for systems that do not have sufficient time to cross to the far side of the barrier, keeping the system further from $x_{\text{trap}}$ and increasing the work done by the trap on the system.

As the trap speed $u$ increases further, the system is more likely to remain on the initial side of the barrier without crossing, shifting probability from the crossing component to the non-crossing component (Fig. 2b, $u=5$ and $u=2\times 10^{4}$). At the highest trap speeds, the protocol duration is insufficient for substantial system diffusion, so the system never crosses to the far side of the barrier.

As trap speed $u$ increases (Fig. 2b, left to right), work distributions depart from Gaussian distributions, gaining the secondary components discussed in the preceding paragraphs. These components, with less overlap than the Gaussian work distributions of linear response, cause time asymmetries above linear response, as shown in Fig. 2a.

For higher barrier $\Delta E_{\text{r}}$, time asymmetry $A$ is higher at a given dissipation (Fig. 2a). This can be considered as either increasing the time asymmetry at a given dissipation, or reaching a given time asymmetry at a lower dissipation. Figure 2c shows work distributions for increasing barrier height $\Delta E_{\text{r}}$. As $\Delta E_{\text{r}}$ increases, the system is less likely to cross from the initial to the far side of the barrier, which narrows the work distributions and decreases their overlap, producing higher time asymmetries.

Maximum dissipation is achieved with the fastest trap speeds, approaching instantaneous trap movement, such that the system is unable to substantially diffuse over the protocol duration. The maximum achievable time asymmetry increases with trap stiffness $k$. For a system effectively frozen in its initial position (drawn from the initial equilibrium ensemble), the maximum dissipation $h_{\text{max}}=\tfrac{1}{2}k$ (shown as a vertical dashed line in Fig. 2a) equals the energy change in changing the trap from $x_{\text{trap}}=0$ to $x_{\text{trap}}=1$, averaged over the equilibrium distribution given $x_{\text{trap}}=0$ (see Appendix C).

As linear-response behavior is achieved by relatively slow processes that remain near equilibrium and in earlier empirical investigations of time asymmetry across a range of speeds and processes Feng and Crooks (2008); Procacci and Marsili (2010), we seek driving characteristics that exceed linear-response time asymmetries for a given dissipation. Figure 3 shows time asymmetry $A$ relative to linear-response time asymmetry $A_{\text{LR}}$ at the same dissipation, $A_{\text{rel}}(h)\equiv A(h)-A_{\text{LR}}(h)$, for two speeds $u$ as spring constant $k$ and barrier height $\Delta E_{\text{r}}$ are varied. $A_{\text{rel}}$ quantifies the improvement of time asymmetry upon the linear response time asymmetry at a given dissipation, distinct from the absolute time asymmetry $A$.

The dependence of $A_{\text{rel}}$ on $k$ is unchanged across a wide range of $u$, with $A_{\text{rel}}\approx 0$ for low $k$ and $A_{\text{rel}}$ sharply increasing to substantial positive values at $k\approx 10$, provided $\Delta E_{\text{r}}$ is not too small. The characteristic width $\delta x$ of system fluctuations in the trap scales as $\delta x=k^{-1/2}$, and at $k\approx 10$ is $\delta x/\ell\approx 0.3$.

$A_{\text{rel}}\approx 0$ for low $\Delta E_{\text{r}}$ and sharply increases at a threshold value of $\Delta E_{\text{r}}$, provided $k\gtrapprox 10$. The threshold $\Delta E_{\text{r}}$ value for the increase in $A_{\text{rel}}$ decreases as $u$ increases. For $u=1$, $A_{\text{rel}}$ increases at $\Delta E_{\text{r}}\approx 4$, corresponding to $t_{\text{cross}}/\Delta t\approx 1.4$. This barrier-crossing timescale is slightly longer than the protocol duration and is similar to relative speeds necessary for substantial non-crossing trajectories in Fig. 2. For $u=100$ (protocol duration of 0.01), the $\Delta E_{\text{r}}$ necessary to increase $A_{\text{rel}}$ drops to nearly zero. Even the barrier-crossing timescale at $\Delta E_{\text{r}}=1$ is much longer than this protocol duration.

In Fig. 3, $A_{\text{rel}}>0$ for traps significantly tighter than the sawtooth wavelength and for barrier crossing timescales much longer than the protocol duration. These conditions combine to limit barrier crossing (which implies more distinction between forward and reverse trajectory distributions) and increase dissipation, allowing $A_{\text{rel}}$ to rise.

For larger $\Delta E_{\text{r}}$, there is an optimum $k\approx 18$ which maximizes $A_{\text{rel}}$ (Fig. 3). At low and high dissipation, linear-response time asymmetry nears the maximum possible time asymmetry, limiting $A_{\text{rel}}$. At intermediate dissipation, there is substantial room between linear response and forbidden time asymmetries, permitting larger $A_{\text{rel}}$. Intermediate $k$ values select these intermediate dissipation values that permit large $A_{\text{rel}}$.

After an initial sharp increase of $A_{\text{rel}}$ at a threshold $\Delta E_{\text{r}}$, $A_{\text{rel}}$ plateaus with further increases of $\Delta E_{\text{r}}$. This plateau underlies the diminishing returns in increasing time asymmetry for increasing $\Delta E_{\text{r}}$ in Fig. 2a, caused once $\Delta E_{\text{r}}$ is sufficiently large to prevent nearly all system trajectories from crossing the barrier.

III.2 Asymmetric ratchets

We now explore asymmetric ratchets with $\ell_{1}\neq\ell_{2}$. Figure 4 shows time asymmetry vs. dissipation for various sawtooth asymmetries. Potentials with $\ell_{1}=0.5\pm\delta$, for a given $\delta\in[0,0.5]$, will yield the same time asymmetry, because time asymmetry does not depend on which work distribution is labeled forward and which is labeled reverse. Accordingly, we increase $\ell_{1}$ from 0.5 towards 1 to explore the effect of an asymmetric potential. We find that time asymmetry $A$ at a given dissipation decreases as the asymmetry of the potential increases.

To understand the time asymmetry decrease as $\ell_{1}$ increases, Fig. 5 compares work distributions for symmetric ($\ell_{1}=\ell_{2}$) and asymmetric ($\ell_{1}\neq\ell_{2}$) potentials. For slow trap speed, the work distributions for symmetric and asymmetric ratchets are similar and appear Gaussian (Fig. 5, far left). As trap speed increases, the work distributions for symmetric and asymmetric ratchets become distinct.

A symmetric sawtooth potential and an intermediate trap speed (Fig. 5a, middle) lead to non-Gaussian forward and reverse work distributions, with a clear peak corresponding to trajectories that do not cross the barrier, and a tail of trajectories that cross the barrier. For a comparable asymmetric ratchet, the forward and reverse work distributions are not symmetric (Fig. 5b, middle). The forward-protocol peak is wider, as $\ell_{1}>0.5$ allows the system during the forward protocol to slowly climb up the gradual potential slope (Fig. 1). The primary reverse-protocol peak is narrower, as the system is unable to make the steep climb; a secondary component of the reverse-protocol work distribution occurs at the same work value as the forward-protocol peak, due to trajectories which stay on the gradually sloped side of the sawtooth. The overlap of the forward-protocol peak and reverse-protocol secondary component decrease the time asymmetry compared to the symmetric work distributions. Compared to symmetric potentials, asymmetric sawtooth potentials increase the probability that the system is on the gradual side of the ratchet potential (independent of whether the trap starts on the steep or gradual side), leading to more similar trajectories overall, and thus more difficulty in distinguishing between the forward and reverse protocols.

Fast driving of a symmetric ratchet produces forward and reverse work distributions composed of single, symmetric, non-Gaussian peaks (Fig. 5a, right). For asymmetric potentials, the work distribution peaks nearly recapitulate the sawtooth potential shape, with one side rising gradually and the other quite steeply (Fig. 5b, right). The work distributions for both symmetric and asymmetric potentials are caused by the initial equilibrium distributions, as there is insufficient time for substantial system diffusion. In addition to the primary peak in each distribution, there is a much smaller secondary component at more positive work values – the reverse-protocol secondary component overlaps with the forward-protocol primary peak. As for intermediate speeds, the greater overlap of the work distributions for the asymmetric sawtooth leads to a lower time asymmetry. The cause of this greater overlap in work distributions is a greater overlap in initial equilibrium positions: a system with a trap which begins the protocol centered on the steep side of the potential (reverse protocol) has a non-vanishing equilibrium probability to be on the gradual side of the potential.

Overall, in an asymmetric ratchet the system is more likely to be in the same initial equilibrium position during forward and reverse protocols, compared to symmetric ratchets. Slower protocols that allow time for diffusion eliminate this difference in overlap of initial positions reflected in the work distributions.

IV Discussion

We have investigated the trade-offs of time asymmetry and dissipation for driven systems diffusing in sawtooth potentials, which are frequently used as simple ratchet models for molecular machines Astumian and Bier (1994); Julicher et al. (1997); Astumian and Derenyi (1998); Parmeggiani et al. (1999); Lipowsky and Harms (2000); Reimann (2002); Suzuki and Munakata (2003); Gabrys et al. (2004); Linke et al. (2005); Lacoste and Mallick (2009); Levien and Bressloff (2015); Kedem et al. (2017); Saakian and Klumper (2018). Accordingly, we explored the effect of sawtooth potential (height and asymmetry) and protocol (speed and driving strength) characteristics on time asymmetry. We have described how these factors lead to a relatively high time asymmetry for a given dissipation, or equivalently how to reduce the dissipation necessary to achieve a given time asymmetry.

For high barriers and intermediate spring constants, high trap velocities lead to time asymmetries well above linear response for a sawtooth potential (Fig. 2a). The barrier is difficult to cross in either direction (as the system must climb the barrier), such that for fast protocols the system rarely crosses the barrier during the entire protocol. This finding suggests molecular machines could gain a higher time asymmetry by proceeding through much of each cycle rapidly, consistent with recent work describing the kinesin cycle with slow transitions following rapid events, such as ATP-induced isomerization of kinesin following ATP binding Chang et al. (2013); Isojima et al. (2016), diffusive search for a forward microtubule binding site following neck linker docking Isojima et al. (2016), and ADP release following binding of the kinesin head to the microtubule Hackney (1988); Nogales and Zhang (2016).

Increasing the barrier height of a sawtooth potential increases time asymmetry for a given dissipation, until the barrier is sufficiently high that the system rarely crosses the barrier (Figs. 2 and 3). The probability of crossing the barrier decreases as the barrier height increases, leading to more distinct forward and reverse trajectory distributions.

The linear-response relationship between time asymmetry and dissipation holds for slow processes that allow the system to stay near equilibrium Speck and Seifert (2004). To assess performance for driving the system over a sawtooth potential, we compared time asymmetries for the sawtooth to time asymmetries for linear-response behavior at the same dissipation (Fig. 3), finding a clear optimum spring constant which maximizes the improvement over linear response. At this optimum spring constant, relatively fast protocols produce dissipation in the range ($1k_{\text{B}}T\lesssim h\lesssim 8$ $k_{\text{B}}T$) at which linear-response time asymmetry is far below forbidden time asymmetries. This optimum spring constant suggests that molecular machines may benefit from states that allow some meaningful position fluctuations without high energetic cost, as long as these fluctuations are not so large as to be comparable to the separation between states. These states — that are not too tightly or loosely constrained — are consistent with the somewhat constrained distributions of kinesin position Isojima et al. (2016) and $\text{F}_{1}$-ATPase angle Yasuda et al. (1998).

When driving over a step potential Brown and Sivak (2016), there also is an optimal spring constant that maximizes time asymmetry relative to linear response. However, driving over steps exhibits increased time asymmetry for intermediate step heights, while driving over ratchets maximizes time asymmetry for large barrier heights. Sawtooth and step potentials are similar in that the crossing timescale for either a sufficiently high barrier or step may be longer than the protocol duration, leading to distinct trajectory distributions for forward and reverse protocols. However, moving the trap from the beginning to the end of the protocol for a step potential changes the equilibrium energy, while for a sawtooth potential the equilibrium energy remains unchanged. For a sufficiently large step, although the trap may be located on the high side of the step, most equilibrium probability will be on the low side of the step, allowing very similar trajectories for forward and reverse protocols; this is not possible for a sawtooth potential. This contrast between large steps and sawtooth barriers—regarding how equilibrium probability is located—causes time asymmetry to be maximal for intermediate step heights and large sawtooth barriers.

Previous work has focused on the effect of barrier shape on transition characteristics and molecular machine behavior Wagoner and Dill (2016); Brown and Sivak (2017b, 2018). Accordingly, we examined how an asymmetric sawtooth potential (with a barrier unevenly spaced between neighboring states) affects time asymmetry compared to a symmetric potential (barrier halfway between states). We found that symmetric ratchets maximize time asymmetry (Fig. 4). Symmetric ratchets maximize the distance between the sawtooth potential barrier maximum and the closest neighboring sawtooth potential minimum, decreasing the equilibrium probability to find the system on the opposite side of the barrier from the trap. This increases the chance of distinct forward and reverse trajectories and leads to higher time asymmetries.

Our finding that time asymmetry (at a given dissipation) is maximized for symmetric sawtooth potentials (Fig. 4) contrasts with previous work on ratchet potentials for molecular machines, primarily with flashing ratchets Astumian and Bier (1994); Julicher et al. (1997); Astumian and Derenyi (1998); Parmeggiani et al. (1999); Lipowsky and Harms (2000); Reimann (2002); Suzuki and Munakata (2003); Linke et al. (2005); Kedem et al. (2017); Saakian and Klumper (2018). In particular, flashing ratchets model autonomous Machta (2015) behavior, and thus require spatially asymmetric sawteeth (or more generally a spatially asymmetric potential) to generate directed motion. In contrast, the deterministic protocols used in this work (i.e. the imposed trap translation schedules) externally provide direction, allowing the sawtooth to be symmetric while maintaining directed behavior. This difference between these two types of driving leads the sawtooth shape required for directed motion in flashing ratchets to be different from the sawtooth shape that optimizes the trade-off between time asymmetry and dissipation for deterministic protocols.

Our results most directly apply to synthetic molecular machines, which are typically externally driven Erbas-Cakmak et al. (2015). Although biomolecular machines are autonomous, some aspects of biomolecular machine operation may approximate an externally driven process: rapid changes, such as ATP binding or neck-linker docking for kinesin Muretta et al. (2015), can occur spontaneously and drive a machine transition (distinct from the initial change) in a manner resembling a control protocol, and may yield similar design principles.

Discrete-state models of molecular machines often describe the effect of a load (e.g. a constant resisting force) by slowing down forward transitions, with rate constant $k^{+}\sim k^{+}_{0}e^{-\delta\omega_{\text{load}}}$, or accelerating reverse transitions, $k^{-}\sim k^{-}_{0}e^{-(1-\delta)\omega_{\text{load}}}$. Here $k^{\pm}$ are the transition rate constants, $k^{\pm}_{0}$ are the rate constants absent a load, and $\omega_{\text{load}}$ is the work done against the load over the transition Brown and Sivak (2018). The splitting factor $\delta$ quantitatively describes how the influence of the load is distributed to the forward and reverse transitions. $\delta$ can be thought of as reflecting the asymmetry of the transition-state position between the two metastable states Brown and Sivak (2017b) ($\delta=0$ indicates the transition state is close to the rear state, and $\delta=1$ indicates the transition state is close to the forward state), analogous to how the distance to the transition state affects unfolding rates Elms et al. (2012). Splitting factors have been estimated for a variety of biomolecular machines, and their values indicate that barriers between molecular machine states range from symmetric to very asymmetric Brown and Sivak (2018), suggesting that evolution could steer biomolecular machines towards symmetric potentials.

Acknowledgements.

This work was supported by a Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant (D.A.S.), funds provided by the Faculty of Science, Simon Fraser University through the President’s Research Start-up Grant (D.A.S.), a Tier-II Canada Research Chair (D.A.S.), an NSERC Undergraduate Student Research Award (A.Z.), and WestGrid (www.westgrid.ca) and Compute Canada Calcul Canada (www.computecanada.ca). The authors thank Steven J. Large and Emma Lathouwers (SFU Physics) for useful discussions and feedback.

Appendix A Nondimensionalization

The system energy $E^{\prime}$ at position $x^{\prime}$ in a quadratic trap with spring constant $k^{\prime}$ centered at $x^{\prime}=x_{\text{trap}}^{\prime}$ with sawtooth barrier height $\Delta E_{\text{r}}^{\prime}$ is

[TABLE]

for $x^{\prime}<\ell_{1}^{\prime}$. We nondimensionalize all quantities using energy scale $H$ and length scale $\ell$:

[TABLE]

We set $E=E^{\prime}/H$, $k=k^{\prime}\ell^{2}/H$, $x=x^{\prime}/\ell$, $\Delta E_{\text{r}}=\Delta E_{\text{r}}^{\prime}/H$, $\ell_{1}=\ell_{1}^{\prime}/\ell$, and $\ell_{2}=\ell_{2}^{\prime}/\ell$. We choose $H=k_{\text{B}}T$, so that all energies are in units of $k_{\text{B}}T$. $\ell$ is chosen so that $\ell_{1}+\ell_{2}=1$, i.e. the nondimensionalized wavelength of the sawtooth is one. The same nondimensionalization applies for positions $x^{\prime}>\ell_{1}^{\prime}$. We also nondimensionalize all times by timescale $\tau$ chosen so that dimensionless diffusivity $D$ is unity.

Appendix B Simulation details

The system moves on a one-dimensional spatial lattice with spacing $\Delta x_{\text{step}}=0.01$, with each sawtooth repeat $\ell=1$ in length. Each timestep $\Delta t_{\text{step}}$, the system randomly attempts to move either left or right, with moves accepted according to the Metropolis criterion Metropolis and Ulam (1949). $D=1$ is maintained by choosing the timestep, $\Delta t_{\text{step}}=(\Delta x_{\text{step}})^{2}/2D$. After each attempted system move, the trap takes a step of size $\Delta x_{\text{trap}}=u\,\Delta t_{\text{step}}$, where $u$ is the trap speed.

Appendix C Instantaneous trap translation

The system is on an infinitely repeating static sawtooth potential (sawtooth wavelength is one) and initially the trap is centered at $x_{\text{trap}}=0$. The energy at position $x$ is $E[x,x_{\text{trap}}(t)]=E_{\text{trap}}[x,x_{\text{trap}}(t)]+E_{\text{r}}(x)$. If the trap moves instantaneously, the system does not move. When $x_{\text{trap}}$ is increased by one (the sawtooth wavelength),

[TABLE]

and $\Delta E_{\text{trap}}(x)=k(\tfrac{1}{2}-x)$. The mean work done by this instantaneous trap movement, the maximum that can be achieved for a monotonic protocol connecting the control-parameter endpoints, is

[TABLE]

where $Z\equiv\int_{-\infty}^{\infty}e^{-E[x,x_{\text{trap}}=0]}{\mathrm{d}}x$. For a symmetric sawtooth potential, the ${\mathcal{O}}(x)$ term in Eq. (9) does not contribute to the integral, because the integral bounds are symmetric about $x=0$ and $kx\,e^{-E[x,x_{\text{trap}}=0]}$ is overall an odd function of $x$. Therefore $\langle W\rangle_{\text{max}}=\tfrac{1}{2}k$. The same result is found for instantaneous trap translation from $x=1$ to $x=0$ (reversing the protocol). Thus, using Eq. 2, the maximum dissipation is $h_{\text{max}}=\langle W\rangle_{\text{max}}=\tfrac{1}{2}k$.