Simulating a chemically fueled molecular motor with nonequilibrium molecular dynamics

Most computer simulations of molecular dynamics take place under equilibrium conditions—in a closed, isolated system, or perhaps one held at constant temperature or pressure. Sometimes, extra tensions, shears, or temperature gradients are introduced to those simulations to probe one type of nonequilibrium response to external forces. Catalysts and molecular motors, however, function based on the nonequilibrium dynamics induced by a chemical reaction’s thermodynamic driving force. In this scenario, simulations require chemostats capable of preserving the chemical concentrations of the nonequilibrium steady state. We develop such a dynamic scheme and use it to observe cycles of a particle-based classical model of a catenane-like molecular motor. Molecular motors are frequently modeled with detailed-balance-breaking Markov models, and we explicitly construct such a picture by coarse graining the microscopic dynamics of our simulations in order to extract rates. This work identifies inter-particle interactions that tune those rates to create a functional motor, thereby yielding a computational playground to investigate the interplay between directional bias, current generation, and coupling strength in molecular information ratchets.

Three SI Movies depict representative simulations. All movies show dynamics at moderate fuel concentration with µ FTC = 0.5, µ ETC = −10, and µ C = −10. The file Motor_Operation.gif shows the dynamics of the full simulation cell with frames taken every 50,000 time steps. The movie shows a clear preference for clockwise motion over counterclockwise motion. The Cleave_Close.gif and Attach_Close.gif movie files show representative pathways for the k cleave,close and k attach,close pathways, respectively. In these movies each frame represents 50 individual time steps. The cleave pathway shows that the transitions happens with a simple unbinding of a red C particle from the close catalytic site. The attach pathway shows a transition that occurs as an FTC cluster decomposes into ETC and C at the catalytic site, leaving behind a C particle as a blocking group. These two pathways are not time-reverses of each other.
In all movie files we translated and rotated the system configuration in each frame to give a clear view of the shuttling ring's relative motion in the frame of a stationary large ring. We first translated the motor's center of mass to the center of the frame and then performed rotations so that the (orange) binding sites were at, approximately, the 12 o'clock and 6 o'clock positions and the large ring was in a plane approximately perpendicular to the observer.

Appendix B: Grand Canonical Monte Carlo Details
To target the grand Canonical distribution of Eq. (15), GCMC trial moves are attempted throughout a volume V 0 . Using the Metropolis criterion, the probability of inserting an additional i molecule and accepting this trial move is P acc i addition (r, p → r , p ) = min 1, where P gen is the probability of generating such a trial move. Here we adopt a notation where r, p denotes the positions and momenta, respectively, of all particles in the initial system with N i (r) molecules of species i (where i can be FTC, ETC, or C). The trial move inserts an additional molecule of i and the positions and momenta of all particles in this configuration is given by r , p . The probabilities P (r, p) and P (r , p ) are given by the grand canonical distribution, Eq. (15). The probability of generating the insertion move is the probability of uniformly selecting that the next move would be an insertion of species i times the probability that the inserted species will be given a particular label in the list of i species times the Boltzmann probability of the inserted molecular configuration: Here r 0 and p 0 are the positions and momenta, respectively, of only the inserted molecule's particles. The trial configuration r , p is then r, p ∪ r 0 , p 0 The total energy H(r 0 , p 0 ) is the internal energy of the inserted molecule if it were in isolation, with a normalizing partition function Z 0 i for the single species i in a volume equal to the volume available for GCMC moves, V 0 : Here h is Planck's constant and n i is the number of particles in species i. Practically, to achieve configurations to insert that come from this Boltzmann ensemble we first pre-sample the FTC and ETC species with standard Markov chain Monte Carlo to build a library of configurations [1]. We then uniformly choose a configuration from this library, rotate that configuration into a uniformly sampled orientation, then place it inside the volume V 0 at a uniformly sampled location. Random velocities for all the inserted particles are sampled from the Boltzmann ensemble independently of the positions. The library generation and random rotation are not necessary for the C species since it is a single spherically symmetrical particle with no internal potential energy. The probability of generating the reversal of the insertion move is simply the probability of choosing to remove a species i particle times the probability of selecting the particular one of the N i (r ) particles to delete: Substituting Eq. (15), Eq. (B2), and Eq. (B4) into Eq. (B1) yields P acc i addition (r, p → r , p ) = min 1, Recognizing that kinetic energy is purely additive across particles, K(p ) = K(p) + K(p 0 ), and that A 0 i = − log Z 0 i /β is the Helmholtz free energy of a single molecule of species i in volume V 0 , we can simplify even further, yielding In a similar way, the acceptance probability of a trial move that removes a molecule of species i is given as In both cases we introduce the definition µ i ≡ µ i − A 0 i as a modified chemical potential that is the true external chemical potential less the free energy of the corresponding species at a reference state with a single molecule in volume V 0 and inverse temperature β. We need not calculate this free energy and can instead use µ i as the "knob" to turn to adjust the concentration of species in the simulation box. We would only need to calculate A 0 i if we wanted to know the true absolute external chemical potential µ i . As µ i becomes larger the concentration of species in the simulation box will increase and as µ i decreases the concentration will decrease. By setting a high chemical potential for FTC and low chemical potential for ETC and C, we can set up a NESS as the FTC decomposes to form ETC and C.
To summarize, a GCMC step begins by uniformly selecting one of six trial moves-insertion or deletion of either FTC, ETC, or C. If an insertion move for FTC or ETC is selected then a configuration of the species is selected from a pre-sampled equilibrium Boltzmann library [1]. This configuration is then randomly rotated. In our case we generated libraries of 10,000 samples for FTC and ETC using standard Markov chain Monte Carlo. The inserted molecule is then uniformly placed in the insertion volume V 0 and velocities are assigned according the Boltzmann distribution. If a deletion move is selected, a uniformly chosen molecule of the species is removed from the simulation box (if there is at least one molecule of the species present). In either case, the trial move is conditionally accepted, with rejected moves returning the system to its original state. Following the GCMC move, the system undergoes Langevin dynamics for a certain number of time steps until the GCMC procedure is repeated. With Langevin time step ∆t = 0.005 we find that 100 time steps between GCMC trial moves is sufficient to maintain equilibrated chemostats.
The relationship between the external chemical potential µ i and the internal composition is given by With concentration C i = N i /V and reference concentration C 0 i = 1/V , the relationship between chemical potential and composition simplifies to The non-ideality of the solution is packaged into the activity coefficient γ i . At low concentrations, the solution is ideal (γ i ≈ 1) because intermolecular interactions are small. As the concentration becomes larger these interactions start to dominate and the relationship between chemical potential and composition becomes non-ideal. With our motor system we only perform GCMC moves in an outer volume, away from the motor (see Fig. 1). This is done to keep the natural dynamics of the motor unperturbed by abrupt insertions and deletions. If the volume where GCMC moves take place is V 0 , which is less than the total simulation volume accessible to the molecules V , and we assume that the diffusion of the molecules is fast relative their reaction (the diffusive Damköhler number is low) then the composition is given as From Fig. C2b we see that as µ FTC is cycled between high and low values the composition of FTC reaches a new steady state very quickly, justifying this assumption. The simulations in Fig. C2 switched between µ FTC = 0.5 (shaded) and µ FTC = −10 (unshaded). The total volume V = L 3 outer with L outer = 34 with the GCMC volume V 0 = L 3 outer − L 3 inner and L outer = 30. The inverse temperature β is 2. Assuming a dilute solution γ i ≈ 1 then from Eq. (B10) we would expect N FTC ≈ 8.7 for the shaded areas of Fig. C2a and N FTC ≈ 0 for the unshaded areas, matching the simulation data well. Thus, at these simulation parameters, the assumptions of ideal solution and low Damköhler number are justified.
As further proof that we can assume our simulations are ideal solutions and that the GCMC chemostats are working as expected, we performed tests where the entire simulation volume was available for GCMC moves and no motor was present. Those simulations were carried out with the full potential and also with non-interacting simulations where intermolecular interactions between different molecules in the simulation were not included. The ideal, non-interacting system allows for direct comparison with Eq. (B9) when γ i = 1. As shown in Fig. B1, non-ideality begins to emerge around µ i = 1.

Chemostat Concentrations
The results of Figs. 3a and b demonstrate clear net clockwise cycling of the motor. In Fig. 3a for moderate numbers of FTC in the simulation box ( N FTC ≈ 2−8) both Motor I and Motor II show significantly increased accuracy above the equilibrium random cycling of 50%. Both motors also show significantly increased current in this range of FTC concentration compared to an equilibrium expectation value of 0. To demonstrate that this net directional cycling is coupled to the reaction of FTC to ETC and C we carried out a series of control simulations where the external chemical potential of FTC was set to -10, low enough that no FTC is expected to appear in the simulation. The concentrations of other species were then varied to determine their effect on the motor. We varied the concentrations of ETC, C, and a variant of FTC where the free central particle was no longer free, but constrained with a FENE spring to the center of the tetrahedron. For all intents and purposes, this constrained FTC (FTC(c)) is almost identical to FTC, but cannot react to form ETC and C. The form of the constraint for FTC(c) was where r c is the vector between the central particle and the center of mass of the four tetrahedron particles, k constraint = 120, and the maximum extension, r max,c is given as the maximum distance that defines the FTC molecule as being in a reactant state, 0.25 in our case. The results of these control simulations, presented in Fig. C1, show that for all three species ETC, C, and FTC(c) there is no deviation from an accuracy of 50% and a current of 0. In other words, these simulations are all equilibrium simulations as the decomposition of FTC and its replenishment through the external chemical potential is required to establish a nonequilibrium steady state.

Chemostat Timescales
To realize current, it was necessary that the chemostats behave sufficiently quickly that NESS concentrations of FTC, ETC, and C could be achieved significantly faster than the timescale for shuttling ring diffusion. We check that this criterion is satisfied by periodically switching from µ FTC = 0.5 to µ FTC = −10, while holding fixed µ ETC = µ C = −10. The responses to this time-periodic perturbation, plotted in Fig. C2, show that the FTC and ETC concentrations relax to their NESS values very rapidly, while the concentration of C relaxes more gradually due to its attractive interactions with the catalytic sites of the motor. Even with this slower relaxation, that central particle equilibrates much faster than the typical timescale to complete cycles. Consequently, the periodically driven motor alternates between regimes of fueled directed motion (shaded) and unfueled equilibrium diffusion (unshaded). When FTC is present the motor has a clear bias toward the CW direction and when FTC is not present the motor reverts to unbiased diffusion where it can cycle CW or CCW with equal probability, resulting in no net current.

Appendix D: Forcefield Parameters
The model consists of 11 different particle types: the two binding sites on the large ring (BIND), 22 inert sites on the large ring (INERT), six catalytic sites (2 each of CAT1, CAT2, and CAT3), 12 sites around the shuttling ring (SHUTTLE), four different types of particles on the vertices of the tetrahedra (TET1, TET2, TET3, and TET4), and the central particle (CENT). The forcefield is fully defined by the masses and radii of these particles (Table I), the pairwise interactions (Table II), and the three-body angular interactions along the backbone of the rings (Table III).      (10)) and bonding interactions. Vertices of the tetrahedra are linked with harmonic bonds while neighboring particles around the two rings are bonded together with FENE (finitely extensible nonlinear elastic) bonds in the order: BIND, CAT2, CAT1, CAT3, 11 INERT, BIND, CAT2, CAT1, CAT3, 11 INERT. For harmonic bonds kij is that of Eq. (11) while for FENE it is the kF,ij of Eq. (12). The two bold entries report the value for Motor I | the value for Motor II. An asterisk in the particle name, e.g. TET*, indicates that the value is the same for TET1, TET2, TET3, and TET4.   (13). Note that the angular interactions are symmetric about the center particle such that k A,ijk would be equivalent to k A,kji .