Non-equilibrium effects of molecular motors on polymers

Abstract

We present a generic coarse-grained model to describe molecular motors acting on polymer substrates, mimicking, for example, RNA polymerase on DNA or kinesin on microtubules. The polymer is modeled as a connected chain of beads; motors are represented as freely diffusing beads which, upon encountering the substrate, bind to it through a short-ranged attractive potential. When bound, motors and polymer beads experience an equal and opposite active force, directed tangential to the polymer; this leads to motion of the motors along the polymer contour. The inclusion of explicit motors differentiates our model from other recent active polymer models. We study, by means of Langevin dynamics simulations, the effect of the motor activity on both the conformational and dynamical properties of the substrate. We find that activity leads, in addition to the expected enhancement of polymer diffusion, to an effective reduction of its persistence length. We discover that this effective “softening” is a consequence of the emergence of double-folded branches, or hairpins, and that it can be tuned by changing the number of motors or the force they generate. Finally, we investigate the effect of the motors on the probability of knot formation. Counter-intuitively our simulations reveal that, even though at equilibrium a more flexible substrate would show an increased knotting probability, motor activity leads to a marked decrease in the occurrence of knotted conformations with respect to equilibrium.

---

1 Introduction

Molecular motors that work out-of-equilibrium are key elements for the viability of cells. Typical examples include myosin and kinesin motors,¹ RNA and DNA polymerases,² helicases³ and condensin.^4,5 Although many of these protein complexes are biochemically well characterised, we are far from having a full understanding of their collective action in vivo. In some cases, the behaviour may be tuned by inter-motor interactions that are mediated by the substrate itself, a germane example being the change in polymerase–DNA binding affinity due to substrate supercoiling.^6,7

Most previous theoretical and computational work on non-equilibrium forces within live cells has focused on modelling the propagation of stresses on biopolymer networks^8,9 or active fields which mimic cytoskeleton dynamics.¹⁰ While molecular motors can now be readily investigated using single-molecule techniques,^11–13 the emergent properties of systems where a collection of motors interact on, and through, their target polymer substrate remain poorly studied.^14,15

Here we propose a coarse-grained computational model to describe the action of generic molecular motors on polymers with the aim of improving our understanding of the effects of ATP-consuming translocating machineries on the conformational and dynamic properties of their target substrates. Within our model, motors undergo free diffusion in 3-D until they encounter the polymer, at which point they become “bound” to it via an attractive interaction. While bound, a motor experiences an active force propelling it along the polymer; at the same time the polymer is subjected to an equal but opposite active force. The result is that the motor displays relative motion with respect to the substrate. The motor can unbind from the polymer, for example due to thermal motion or because it reaches the polymer's end. While our model is generic, and our aim is to understand its underlying physics, it could be viewed as a model for DNA or RNA polymerases moving along DNA or kinesin molecules stepping on microtubules.²

There has been a great deal of recent simulation work exploring active polymers, and our model differs from these in a number of ways. To our knowledge our work is the first to explicitly treat discrete motors that track along a polymer in a fully 3-D system. A number of recent studies^14–20 have considered a chain of beads, where each bead experiences an active force, either directed along the chain bonds^17,18,20 or in a direction independent of the chain.^14,19 These works have mainly focused on regimes of different behaviour at different levels of activity (e.g., transitions from active translation to rotational motion¹⁷ or to spiral configurations¹⁸ have been observed – for a review see ref. 21). In the present work the polymer only experiences an active force at the points where the explicit motors are bound (and these vary in time), and we focus on the effect on the conformational and dynamical properties of the polymer. Another class of models have more specifically focused on motor proteins found in the cytoskeleton, like myosin and kinesin which have provided a model experimental system for active materials.^22,23 These can be modelled at a field theoretic level, but there have been a number of simulation studies which resolved the motors explicitly.^24–27 Unlike the current work, those studies focused on dense suspensions of motors and their substrates, and/or 2-D systems where the motors are fixed to a surface, and have often considered joined pairs of motors which can crosslink different substrate molecules to drive them into relative motion. In contrast, here we consider a single, semi-flexible polymer which is free to move in 3-D, and is subject to the action of multiple motors.

Our simulations reveal that motor activity affects the steady state conformations of the polymer substrate; we find that these are more crumpled with respect to their passive counterparts. We attribute this effective softening of the substrate to the emergence of double-folded branches, or hairpins; the effect can be tuned by changing, e.g., the strength of motor activity. We also find that, for large enough motor activity, the centre of mass of the polymer substrate can display super-diffusive behaviour at short times, with a return to diffusive behaviour at long times. Again this effect can be controlled by varying the activity or number of motors. We conclude our work by investigating the effect of motors on the propensity of the polymer to form knots. This is relevant for many biological processes involving DNA, where knotting would be detrimental (e.g., knots and tangles might hinder transcription and replication, repair of double strand breaks, or segregation of chromosomes). We find that the motor activity effectively reduces the steady state probability of finding knotted conformations. This is surprising since at equilibrium, in a good solvent, and in the absence of confinement, softer substrates show an enhanced probability of knotting.²⁸ Thus, our findings point to an intriguing non-equilibrium effect of molecular motors that could be important in many biological systems.

2 The model

We perform Langevin dynamics simulations of a semi-flexible bead-and-spring polymer consisting of L beads in solution with N motors, as described schematically in Fig. 1. We use a standard polymer model which includes finitely-extensible non-linear elastic (FENE) springs connecting consecutive beads, a Kratky–Porod potential for triplets of beads providing bending rigidity, and a Weeks–Chandler–Anderson (WCA) potential to regulate steric interactions (see Appendix A for full details). The intrinsic stiffness of the polymer is controlled by the bending energy K_BEND which appears as a parameter in the Kratky–Porod potential. Each polymer bead also displays an orientation, which is tracked by a vector û_i as shown in Fig. 1, and we include a potential to orient this vector along the local direction of the polymer backbone (if the position of bead i is r_i and the orientation is û_i, then the potential acts to minimise the angle ϕ_i between û_i and |r_i+1 − r_i|).

Fig. 1 Schematic of our model for molecular motors acting on a generic polymer substrate. Diffusing motors (blue) can interact with the polymer backbone (grey); when bound, the motors (red) and the polymer experience opposite forces which cause relative motion. On the right hand side of the figure we show the polymer in more detail: the vector û defines an orientation for each bead and ϕ_i is then defined as the angle between û_i and the vector joining beads i and i + 1; an orientation potential acts to minimise ϕ_i (see text and Appendix A for details).

The motors are represented by freely diffusing beads which interact with each other sterically (via a WCA potential) and via a Morse-like attractive potential with the polymer beads. The latter is given by

Here r is the separation of the motor and the polymer bead, and K_MOT is the interaction strength. The potential has a minimum when the motor and polymer beads overlap. We consider a motor to be “bound” to a polymer bead when the centre-to-centre separation of the bead and the motor is less than r_m cut, which we set to be the same as the motor diameter. In addition to this attraction, when a motor is bound to a polymer bead, it experiences an active force of magnitude f in the direction of the orientation vector û_i of the bead; the polymer bead experiences an opposite active force −fû_i (in compliance with Newton's third law). It is possible for a motor to be within a distance r_mcut of more than one polymer bead at the same time: in this case each interacting polymer bead will experience an active force, and the total active force on the motor will be the sum of that due to each interacting polymer bead. (In practice motors spend over 95% of their “bound” time interacting with two consecutive polymer beads.)

The result of the active force is that the motors and the polymer are driven into relative motion; due to the connectivity of the chain, the polymer tends to move less compared to motors, so the latter effectively “track” along the polymer substrate. Motors can become unbound from the polymer due to thermal fluctuations (controlled by the parameter K_MOT), and also active effects (e.g. tracking off the end of the polymer – this is discussed in more detail below). Since each monomer has a preferred orientation, i.e., the polymer is polar, motors will move in the parallel or anti-parallel direction depending on the sign of f. In the present work we consider only the case where all motors are the same and f ≥ 0.

For simplicity motors and polymer beads are taken to have the same diameter, σ, mass, m, and unless otherwise stated they experience the same friction ξ due to the implicit solvent (though in reality a motor such as RNA polymerase would be much bigger). We use the LAMMPs molecular dynamics software²⁹ to perform the simulations: a Langevin equation for each polymer bead and motor is solved using a velocity-Verlet algorithm. A separate Langevin equation is solved for each monomer orientation û_i. Full details are given in Appendix A. Throughout we quote energies and times in units of k_BT and respectively; how these map to physical units will depend on the system of interest, as detailed in Appendix A.

3 Results

As discussed above, although a motor and the polymer bead it is bound to experience forces of equal magnitude (but acting in opposite directions), the chain connectivity means that the motors can move more freely, and tend to track along the polymer substrate. The three control parameters which we vary in our simulations are the number of motors N, the interaction strength K_MOT, and the force acting on the motors f. However, in what follows we also consider 〈n〉, corresponding to the mean number of motors bound to the polymer at any one time. We note that these parameters are not independent: for large forces, motors have a shorter residency time due to the finite size of the polymer substrate and if both N and f are large, collisions between motors can also lead to a decrease in residency time. That is to say, changing f at fixed N leads to a change in 〈n〉, and this relationship depends on N itself – this is examined in more detail in Appendix B.

As expected, the tracking speed of the motors in general increases with f, though the relationship is complex (as noted, changing f leads to a change in 〈n〉 which can affect motor speed due to collisions). Within the range of parameters investigated, typically a motor will move from one polymer bead to the next in a few τ_LJ (e.g., for f = 20 and N = 200 we measure a speed ∼0.4 beads/τ_LJ). For an L = 500 bead polymer, motors typically spend on the order 50–100 τ_LJ bound to it (increasing with K_MOT and decreasing with increasing f and 〈n〉).

The motor action also changes the dynamics of the polymer. At a global scale this can be quantified by measuring the decorrelation time for the polymer end-to-end vector. For an L = 500 bead polymer in the absence of motor action (f = 0) a typical decorrelation time is of the order 10⁵τ_LJ. Again, there is a complicated dependence on f and N, but the motor effect leads to a decrease in this decorrelation time by around an order of magnitude. Nevertheless, the polymer relaxation time remains much longer than the motor time-scales detailed above, i.e., the positions of motors on the polymer change quickly with respect to the polymer dynamics. This implies that, after some transient time, the properties of the polymer will settle into a non-equilibrium steady-state (NESS) which is independent of initial conditions (and the system is ergodic). In the following sections we examine these steady state polymer properties by averaging measurements over time and over an ensemble of 50 independent simulations, unless otherwise stated.

Motor activity leads to polymer softening

In this section we study the role of motors on the steady state conformations of the polymer substrate. We examine the effect of different values of 〈n〉 and f, evolving the system until a NESS is reached. By visual inspection of the simulated trajectories we notice that, for large 〈n〉, the typical polymer conformations become more crumpled and display kinks which do not appear in passive (f = 0) systems (Fig. 2(a and b)).

Fig. 2 (a and b) Comparison of conformations for a polymer of length L = 500 with 〈n〉 = 0 and 〈n〉 = 17 motors attached at any time. In (b) hairpins are indicated by arrows and the purple beads mark the presence of an active motor bound to the polymer backbone. For clarity of visualisation, motors are depicted as having a bigger size with respect to the polymer beads. (c) Persistence length L_p as a function of the average number of motors attached 〈n〉 for different values of f. Inset: The persistence length for different motor forces can be scaled on a master curve L_p(x = 〈n〉f²) ∼ x^α with α = −0.108 ± 0.007. (d) Similar plot showing R_g as a function of 〈n〉. Again the inset shows the collapse of the curves when plotted as a function of 〈n〉f² with coefficient α = −0.100 ± 0.009. Error bars show the standard error in the mean (though these are often smaller than the points), and the lines connecting symbols are always a guide for the eye.

To quantify the extent of polymer crumpling we measure the radius of gyration R_g defined by

where r_i is the position of the i-th polymer bead, R_cm is the centre of mass of all polymer beads, and 〈⋯〉 denotes ensemble and time average (after the system reaches a steady state). We also measure the persistence length L_p as defined from the bond–bond correlation function

〈b_i·b_i+s〉 = e^−s/L_p, (3)

where b_i = (r_i+1 − r_i)/|r_i+1 − r_i| represents the normalised tangent to the polymer at bead i. Here the average is additionally evaluated over bead index i. To calculate L_p we perform an exponential fit of the bond–bond correlation function (even for cases in which this displays negative values – see below for more on this point).

As shown in Fig. 2(c and d), both R_g and L_p exhibit a marked decrease with increasing 〈n〉, in agreement with our initial qualitative observation. Furthermore, the decrease in effective persistence length L_p is found to have a power-law scaling

L_p ∼ (〈n〉f²)^α,

with α ≈ −0.1, which can also be used to collapse the different curves for R_g onto a universal master curve (Fig. 2(d), inset). In other words we find that the decrease in average size of the polymer can be well described by a scaling function of the parameter 〈n〉f².

In order to shed more light onto the effect of the motors on the metric exponent ν relating the polymer size to its length (R_g ∼ L^ν), in Fig. 3 we plot R_g a function of L for different combinations of the number of motors N and force f. We observe that passive systems (N = 0 or f = 0) closely follow the expected self-avoiding statistics (ν ≈ 0.588). More interestingly, larger values of N or f lead to polymer behaviours that are better fitted by ν ≃ 0.43. This value lies between that of an ideal chain ν = 1/2 and that of a fully collapsed conformation ν = 1/3. It is intriguing to note that a simple scaling argument can partially explain this finding using the value of α measured in Fig. 2 (see Appendix C).

Fig. 3 Radius of gyration of the polymer substrate as a function of substrate lengths L, for (a) a fixed motor force (fσ/k_BT = 20) and different values for the total number of motors N; and (b) a fixed number of motors (N = 200) and different values of the motor force f.

To conclude this section we examine the effect of the motors on the global “shape” of the polymer by measuring two quantities: the asphericity Δ and prolateness S, which quantify the degree of spherical symmetry and the oblateness of the polymer respectively (see Appendix D for formal definitions). In Fig. 4, one can readily observe that for increasing 〈n〉 and f both shape descriptors decrease. This indicates that the coil is becoming more spherical, in agreement with our previous observations that the polymer effectively crumples under the action of the motors.

Fig. 4 (a) Polymer prolateness S as a function of 〈n〉, for three different values of the motor force f. (b) Polymer asphericity Δ shows a similar trend. Insets show curve collapse when plotted as a function of 〈n〉f² in a double logarithmic scale.

Effective softening is explained by hairpin formation

Having established a crumpling effect on the global polymer conformation, we now examine the microscopic mechanism driving this process. We measure the polymer bond auto-correlation function

β(s) ≡ 〈b_i·b_i+s〉,

whereas before angle brackets denote the time and ensemble average, and average over i. For an equilibrium polymer β(s) decays exponentially to zero (the bond orientation becomes uncorrelated for contour separations much larger than the persistence length). If β(s) is close to one for a given s this indicates that bonds with that separation are on average oriented parallel to each other; whereas if β(s) → −1, then bonds with that separation are on average anti-parallel. As shown in Fig. 5(a), β(s) drops below zero and displays a minimum – this indicates a non-negligible probability that bonds separated by about 75 beads are anti-parallel. In other words, the polymer conformations display hairpins, or “branches”³⁰ and the distance s for which β(s) has a minimum approximately corresponds to the length of the hairpins that are formed on the polymer. The minimum β_min becomes deeper as both 〈n〉 and f increase (Fig. 5(b)). In other words, higher values of the motor force cause a more significant anti-correlation of the bonds, i.e., more pronounced hairpins. A fit of β_min to the function f(x) = ax^b gives us its scaling as a function of 〈n〉, with coefficients a = −0.047 ± 0.003 and b = 0.27 ± 0.03. An analogous scaling law can be found for β_min as a function of the force f, this time with coefficients a = −0.024 ± 0.006 and b = 0.39 ± 0.08.

Fig. 5 (a) Polymer bond correlation as a function of the distance along the polymer backbone s, at fixed motor force f = 20. The arrow indicates the direction of increasing 〈n〉. (b) Plot showing the minimum value of β(s) as a function of 〈n〉. Inset: Scaling of the minimum value of this functions, β_min, with 〈n〉f². (c) Polymer bond correlation as a function of the distance along the polymer backbone s, at fixed motor force f = 20, 〈n〉 = 17 and different (initial) persistence lengths L_p. Inset: Position of the bond correlation minimum, L_h (defined as the value which satisfies β(s = L_h) = β_min), as a function of the polymer persistence length L_p. (d) Plot of the polymer bond correlation for the case of motor force f = 20 and different values of the friction ξ experienced by the motors. Inset: Position of the bond correlation minimum L_h as a function of the friction ξ.

A possible physical mechanism underlying hairpin formation is that the motors exert localised tangential forces on the semiflexible polymer they move on. As a result the portion of the polymer “behind” the motor (i.e., the polymer region through which the motor has just moved) is under compression, while the portion ahead is under tension – the situation is reminiscent of forced translocation through a pore.³¹ The conformation of the polymer could then be computed by minimising its bending energy in the presence of such a localised force. This problem may be viewed as a generalisation of the Euler buckling problem.^32,33 Although a detailed calculation is outside the scope of this work, it is reasonable to expect the combination of tensile and compressive forces to result in local hairpin formation.

In Fig. 5(c) we analyse the dependence of the bond correlation function on the intrinsic stiffness of the polymer set by K_BEND (recall that in the absence of motor activity, the persistence length L_p ≃ K_BENDσ/k_BT). To this end, we perform three sets of simulations with K_BEND = 10, 20 and 30k_BT, for the case of 〈n〉 = 17 and f = 20. Plotting β(s) we observe that the position of the minimum becomes progressively shifted towards higher values of s as K_BEND increases, i.e., the hairpins become larger (Fig. 5(c), inset). Additionally, we find that the value of β(s) at the minimum, β_min, gets closer to zero for smaller values of K_BEND.

In Fig. 5(d) we study how changing the relative impact of the active force on the motors and polymer beads affects the hairpin formation. This is achieved by changing the friction experienced by the motors (i.e., we change the viscous drag only for the motors – this changes their effective hydrodynamic size, or could reflect a difference in interactions with macromolecular crowders compared to the polymer). Increasing the friction for the motors means that they will move less as a result of the active force (which is kept at a constant magnitude). Fig. 5(d) shows that as the friction for the motors is increased relative to that of the polymer the minimum of β(s) shifts to shorter s, and gets deeper and narrower. This implies that greater motor friction leads to shorter and tighter hairpins. We can understand this by considering the opposite limit of very low motor friction; in that case the motors will move more quickly relative to the polymer. Since a motor spends less time at a given polymer bead, that bead will experience the active force for less time and we can expect the effect of the motors to be diminished. The polymer will, however, still experience the active force, so we still expect a ballistic component to its motion.

Motor activity gives rise to anomalous or enhanced diffusion

We now turn our attention to the role of molecular motors on the overall polymer dynamics. To this end, we measure the mean square displacement of the polymer centre of mass

MSD_CM(t) = 〈[R_cm(t₀ + t) − R_cm(t₀)]²〉,

where the average is over start times t₀, and R_cm(t) is the position of the centre of mass at time t.

Fig. 6 shows that at short times the system is diffusive as the propulsion exerted by the motors needs a characteristic time in order to provide a displacement comparable to a Brownian displacement. Following the Rouse model³⁴ MSD_CM(t) = (6D₀/L)t, where D₀ = k_BT/ξ is the diffusion coefficient of a single monomer. This behaviour is recovered at all times for 〈n〉 → 0.

Fig. 6 (a) Mean squared displacement of the polymer centre of mass as a function of time. Each set of points shows results for a different value of 〈n〉 at fixed f = 20 k_BTσ⁻¹. Higher values of 〈n〉 induce a faster diffusivity of the polymer. (b) Similar plot but now each set of points shows results for a different value of f, with the total number of motors fixed at N = 400.

At intermediate times, the polymer exhibits super-diffusive motion, i.e., MSD_CM(t) ∼ t^α with α ∼ 1.6. Finally, at large times we observed enhanced diffusion, i.e., α ≈ 1 but the diffusion constant increases with motor activity (f or 〈n〉). In Table 1 we report the measured diffusion constant D for the different cases considered. One can readily observe that the exponent α is roughly independent of 〈n〉 and f, whereas the diffusion coefficient increases with 〈n〉. Moreover, by comparing fits at short and long times for each value of 〈n〉, we find that the crossover between the super-diffusive and diffusive behaviour occurs at progressively larger times as 〈n〉 decreases. Together, this implies that the motor action leads to a persistent motion of the centre of mass of the polymer; at longer times the velocity becomes uncorrelated and the motion is again diffusive.

Table 1 Diffusion properties of the centre of mass of the polymer are measured by fitting to the MSD for different numbers of motors. Fits are evaluated using the function f(t) = a·x^α, and the diffusion coefficient (or mobility μ) is then D = a/6 following the Rouse model

Intermediate times: super-diffusive
〈n〉	α	μ
17	1.720 ± 0.002	3.1 × 10^−5 ± 0.5 × 10^−6
9	1.690 ± 0.001	1.7 × 10^−5 ± 0.2 × 10^−6
5	1.542 ± 0.006	2.6 × 10^−5 ± 1.3 × 10^−6

---

Long times: diffusive
〈n〉	α	D
17	1.025 ± 0.001	2.7 × 10^−2 ± 2.6 × 10^−4
9	1.020 ± 0.002	2.0 × 10^−2 ± 0.06 × 10^−4
5	1.090 ± 0.0001	0.4 × 10^−2 ± 0.08 × 10^−4

---

These three regimes are typical of other active systems, e.g., in polymers where an active force directed along the backbone is applied to each polymer bead^14,20 and in Active Brownian Particles (ABPs).³⁵ At intermediate times, the motion of an ABP is dictated by a force which typically maintains its direction for a characteristic time, leading to super-diffusion. In our model, a correlation on the scale of the whole polymer is set by the end-to-end vector.

Untying knots with molecular motors

Finally, we discuss a possible application of our model by studying the knotting probability of a polymer in non-equilibrium conditions. We start from the observation that a freely diffusing linear polymer in a good solvent may display a physical knot – i.e., a knot within an open segment of the polymer contour – with probability P_knot. It has recently been shown^28,36,37 that P_knot displays a non-monotonic dependence on the stiffness of the substrate with a maximum at around L_p ≃ 5–10σ. In light of the results presented above – i.e., motors drive a decrease in the apparent persistence length of the substrate – we reason that the presence of the motors should lead to a shift in the knotting probability as a function of substrate persistence length. Specifically, we expect that, due to the action of the motors, the maximum in P_knot will be located at larger values of the intrinsic persistence length of the substrate, l_p (to be clear, the lower case l_p refers to the persistence length the polymer would have in the absence of motors).

In order to test this hypothesis we perform simulations of a freely diffusing linear polymer 500 beads long, and measure the probability of forming physical knots as a function of l_p by adopting the algorithm described in ref. 38 and 39 and publicly available at http://kymoknot.sissa.it/.⁴⁰ By comparing P_knot(l_p) for cases with and without motors we find that, surprisingly, the location of the maximum of P_knot is unchanged by the action of the motors, i.e., l_p ≃ 5–10σ. The most apparent effect of the motor activity is a global reduction in the knotting probability across the full spectrum of bending rigidities studied (Fig. 7). In light of this finding we argue that the effect of the motors goes beyond a simple softening of the polymer backbone, and that subtler non-equilibrium effects must be at play. We also speculate that the reduced knotting probability due to motor activity may be relevant for the case of DNA under the action of motors such as polymerase. According to our model, the action of these motors would naturally have the beneficial consequence of reducing the occurrence of knots and self-entanglements in vivo.⁴¹

Fig. 7 Molecular motors acting on a polymeric substrate reduce the knotting probability P_knot(l_p) of a free chain in a good solvent. Here l_p is the intrinsic persistence length which would be measured in the absence of motors. The simulations are performed on a 500 bead long chain and the points are averages over 10⁴ independent configurations (error bars show the standard error in the mean).

4 Conclusions

In this work we performed Langevin dynamics simulations of a coarse grained model to study the effect of molecular motors on a polymer substrate. We examined both the static and dynamic properties of polymers in dilute conditions while varying the number of motors and the magnitude of the force which they generate.

Our main result was that the effect of the motors is to reduce the overall size of the polymer coil by increasing the likelihood of double-folded segments, or hairpins. This “softening” of the backbone manifests as a reduction of the measured persistence length of the polymer. We found that the data collapse onto a single curve when plotted as a function of 〈n〉f², revealing a power-law scaling R_g ∼ (〈n〉f²)^−0.1 and L_p ∼ (〈n〉f²)^−0.1 (insets of Fig. 2(c and d)). We also performed simulations of polymers with different lengths and obtained the explicit values of the metric exponent ν (where R_g ∼ L^ν) as a function of 〈n〉 and f. These confirm that at equilibrium the polymers obey self-avoiding statistics, whereas motor activity leads to an exponent between those expected for a crumpled globule and an ideal chain (Fig. 3). A reduction in R_g and L_p has also been observed in models without explicit motors, where an active force directed tangentially to the polymer is added directly to each polymer bead.^14,20 In that case the softening arises because, for an initially straight segment of the polymer (e.g., of length equal to the persistence length), small fluctuations in backbone direction will lead to components of the tangential force which are perpendicular to the end-to-end vector for the segment; this then drives the polymer into large smooth curves.¹⁴ Those studies, however, did not find hairpins (the decay in bond correlations remained exponential²⁰), suggesting the present case involves a different mechanism.

We also found that the motors affect the dynamics of the polymer. At intermediate times the centre of mass displays super-diffusive behaviour while at longer times it follows an enhanced diffusion. Both regimes are controlled by motor activity (number of bound motors and force applied). While a super-diffusive regime was to be expected, our simulation have provided a detailed quantification of how that can be tuned through the action of the motors (Fig. 6 and Table 1). The effect of the motors on the knotting probability of the polymer was subtler. Unexpectedly, P_knot(l_p) did not simply show a shift due to the change in effective persistence length, but rather there was a global reduction in the steady state knotting probability across a range of substrate rigidities (Fig. 7). We can speculate that since the motor action promotes polymer configurations with hairpin bends, but suppresses knotted ones, then these two types of configuration are incompatible. An interesting future study would be to examine knotting probabilities in other active polymer models where hairpins are not observed.

Throughout this paper we have limited ourselves to the simple case in which all of the motors travel in the same direction along the substrate. While this is the case for some systems, e.g., kinesin on microtubules, in others the motors can travel in either direction. For example, RNA polymerase follows the direction of genes, depending on the specifically oriented DNA sequences in gene promoters. This could be readily incorporated into our model in the future by either considering a mixture of two species of motors which move in different directions on the substrate, or by setting some directional binding sites. Also, we have considered the case of motors which are the same size as the polymer beads – in reality RNA polymerase is a large molecule of around 30 nm in size. However, as our results are mostly affected by the magnitude of the force, we do not expect motor size to have a qualitative affect. Nevertheless it would be interesting to examine systems with larger motors in the future.

Another interesting scenario would be to consider a mixture of passive and active polymers. Indeed, an activity-driven phase separation has been observed in another non-equilibrium polymer system – where two species of polymer were held at different temperatures.⁴² It would be interesting to see whether such demixing behaviour can be recovered when the system is driven away from equilibrium by a more biologically relevant mechanism, such as that presented here.

Our model naturally lends itself to the study of the response of polymer substrates to the action of molecular motors. We therefore expect that in the future is could be employed to address more specific biological questions, such as understanding the role of motors in the buckling of actin networks,^8,9 or in the bending and active collective motion of microtubules^43,44 by suitably tuning the rigidity of the substrate and processivity of the motors.

Finally we stress that the model presented here is among the first to explicitly account for the presence of physical motors; as a consequence, it faithfully captures the feedback between the action of the motors and the change in substrate–motor interaction due to, for instance, alterations in local polymer conformation. For this reason, we believe that our model will be pivotal in studying complex active collective phenomena where interactions between motors are mediated through the target substrate.

Conflicts of interest

There are no conflicts to declare.

Appendix A: full model details

We use the LAMMPS software²⁹ to perform Brownian dynamics simulations, where the position of bead i is determined by the stochastic differential equationwhere m_i is the mass of the bead, ξ_i is the friction it experiences due to an implicit solvent, and η_i is a vector representing random uncorrelated noise which obeys the following relations

〈η_α(t)〉 = 0 and 〈η_α(t)η_β(t′)〉 = δ_αβδ(t − t′). (5)

The potential U_i is a sum of interactions between bead i and all other beads, as described below. F_i is the active force on bead i due to motor–polymer interactions. Such force, applied to the polymer by all the motors in equal measure, is directed along the polymer backbone. For simplicity we assume that all polymer and motor beads in the system have the same mass m_i ≡ m = 1, and (except for the case of Fig. 5d) friction ξ_i ≡ ξ = 2. The orientation of each polymer bead is described by three orthogonal unit vectors f̂_i,v̂_i and û_i which form a right-handed set of axes. A similar stochastic differential equation describes the dynamics of the orientationwhere ω_i is the angular velocity of bead i, and ξ_r,i and I_i are the rotational friction and moment of inertia respectively. The latter is set according to the Stokes–Einstein relation assuming a solid sphere. The rotational noise vector η_r,i(t) has the same properties as the translation case, and the torque T_i is due to the potential given below. The equations of motion are solved in LAMMPS using a standard Velocity-Verlet algorithm. Polymer beads are connected via FENE bonds according to the potentialwhere r_i,i+1 = |r_i − r_i+1| is the separation of the beads, and the first term is the Weeks–Chandler–Andersen (WCA) potentialwhich represents a steric interaction preventing adjacent beads from overlapping. In eqn (8)d_ij is the mean of the diameters of beads i and j, which for simplicity is the same for motors and polymer beads, and denoted σ. In eqn (7)R₀ and K_FENE are the maximum extension and strength of the bond respectively, and we use R₀ = 1.6σ and K_FENE = 30k_BT throughout.

The bending rigidity of the polymer is introduced via a Kratky–Porod potential for every three adjacent DNA beads

U_BEND(θ) = K_BEND[1 − cos(θ)], (9)

where θ is the angle between the three beads as given byand K_BEND is the bending energy. In the absence of motors the persistence length in units of σ is given by L_p = K_BEND/k_BT. Unless otherwise stated we set K_BEND= 20k_BT, i.e. we consider a semi-flexible polymer; this has previously been used as a model for DNA – if we take the bead diameter σ = 2.5 nm, then the persistence length (in the absence of motors) is L_p = 50 nm, which is relevant for DNA in typical physiological conditions. Once the simulation length scale has been fixed, and using the energy unit of k_BT with temperature T = 298 K, one can map simulation time to physical units through the Brownian time τ_B = σ²/D, using the Stokes–Einstein relation to find the diffusion constant for a sphere D = k_BT/(3πησ). For a DNA system a suitable value of fluid viscosity results in a Brownian time τ_B = 36 ns; with our choice of m = 1 and ξ = 2, the simulation time unit would then be τ = τ_B/2 = 18 ns.

The torque on the polymer beads is due to a potential which constrains bead orientations to lie along the polymer tangent, given by

U_ORIENT = K_ORIENT[1 − cos(ϕ_i)], (11)

where ϕ_i is the angle between û_i and the direction of the polymer backbone such thatWe set the orientation energy K_ORIENT = 30k_BT. Steric interactions between non-adjacent DNA beads are also given by the WCA potential [eqn (8)]. Note that this is the polymer model described in ref. 45 but without torsional rigidity.

The motors are represented by freely diffusing beads which interact with each other via a WCA potential and via the potential given in eqn (1). A motor is said to be “bound” to the polymer when its distance to a polymer bead is less than r_mcut = σ. In addition to this interaction, when a motor is bound to a polymer bead, it experiences an active force of magnitude f which is directed along the orientation vector of the polymer bead. As a consequence of Newton's third law, the polymer bead experiences an equal and opposite force. The result is that for non-zero f the motors track along the polymer substrate while driving the motion of the polymer. In order to avoid accumulation and jamming we want motors to detach once they have reached the end of the substrate. To achieve this we impose a large force f_end ≫ K_MOT/σ to motors when bound to the last bead of the chain. Though large, this force does not appreciably influence the dynamics nor the conformation of the chain. Since the polymer is polar, motors can be made to move in the parallel or anti-parallel direction depending on the sign of f. In the present work we consider the case where all motors are the same and f ≥ 0.

Appendix B: motor statistics

The average number of motors attached 〈n〉 is a key quantity, as it determines the average, total active force applied on the polymer. It is also connected to the average time a motor spends attached to the polymer, the residence time. In equilibrium (f = 0), the motors do not move along the polymer and the residence time scales as e^K_MOT/k_BT. We report the average number of attached motors 〈n〉 in Table 2, for the values of the total number of motors N, the interaction strength K_MOT and the motor force f considered in our simulations. We observe that 〈n〉 strongly depends on f, as well as N and K_MOT. For fσ/2 < K_MOT, 〈n〉 increases with K_MOT as the residence time is still proportional to the equilibrium one; we call this the thermodynamic unbinding regime. If fσ/2 ≈ K_MOT, then 〈n〉 decreases; we call this the activity induced unbinding regime, and propose three possible mechanisms. First, for larger forces the motors move faster, and will “fall off the end” of the polymer after a time shorter than the equilibrium residence time; second, the local bends induced by the motors may provide a negative feedback, as a motor is more likely to detach as it moves around a sharp bend; and third, at very high values of 〈n〉, motor collisions might also result in a decreased residence time.

Table 2 Average number of motors attached 〈n〉, as measured from the simulations, for the different values of K_MOT, N and f considered

K MOT = 10			K MOT = 20			K MOT = 30
N	f	〈n〉	N	f	〈n〉	N	f	〈n〉
100	6	9	200	6	165	200	6	200
100	40	1	200	40	12	200	40	160
200	6	18	300	6	250	300	6	300
200	40	3	300	40	20	300	40	250
400	6	34	400	6	300	400	6	400
400	40	6	400	40	30	400	40	320

---

In Fig. 8a we show how the linear density of bound motors 〈n〉/L depends on L and f, at constant N = 200. For small motor force (the thermodynamic unbinding regime) we find that 〈n〉/L increases roughly linearly with increasing L, as would be expected at equilibrium; at high f longer polymers are more populated than shorter ones. This suggests that motors moving off the end of the polymer play a more important role in the latter regime. In Fig. 8b we plot the fraction of motors which are bound, 〈n〉/N, as a function of f at fixed L = 500; data for different N collapse onto the same curve, implying that, even in the large f regime, motor crowding does not have a large effect. Interestingly, if we plot the average “used” motor force f〈n〉N (Fig. 8b, inset) we observe non-monotonic behaviour: the total force is maximal at f ∼ 15k_BT/σ.

Fig. 8 (a) Average fraction of occupied sites as a function of the motor force for different substrate lengths L and fixed N = 200. (b) Average fraction of bound motors for different N and fixed L = 500. Inset: Average “used” motor force as function of the “nominal” motor force for different N and fixed L = 500. In all cases, K_MOT = 10k_BT.

Appendix C: scaling argument for R_g

In the main text we have reported that the radius of gyration scales as R_g ∼ (〈n〉f²)^−α with α = 0.1. We can link this observation with the scaling of R_g with the contour length via the following effective blob picture: at steady state, the only length scale that can be extracted from the parameters of the system and that depend on the motor force is ξ_b = k_BT/f. In analogy with traditional scaling arguments, we will take this to be the typical blob size.⁴⁶ Within a single blob the steric repulsions are not screened and the chain behaves like a self-avoiding walk (we expect the blob size to scale with the number of monomers per blob as ∼g^3/5). On a larger scale, the size of the chain will grow as R_g ∼ m^ν where m is the number of independent blobs and ν is the metric exponent. In order to find the number of blobs in a chain we first need to compute the number of monomers per blob, i.e.,As mentioned above, the size of the polymer then scales aswhere L is the total number of monomers in the chain. By combining these equations, we arrive at R_g ∼ f^(5/3)ν−1. Comparing this with the numerical result from the simulations, R_g ∼ f^−0.2 (Fig. 2(b), inset), we obtainor ν = 0.48 ≈ 1/2, i.e., the exponent of an ideal random walk. It is intriguing to notice that this value is not far from the one obtained by performing explicit simulations with different polymer lengths L in Fig. 3.

Appendix D: measurements from simulated trajectories

As detailed in the text, the main quantities we measured from the polymer were the radius of gyration R_g, the bond–bond correlation function β(s) (from which we obtain the effective persistence length L_p), and the mean squared displacement of the centre of mass of the polymer. In each case the system was initiated with the polymer in a random walk conformation; the dynamics were evolved in the absence of motor interactions for 10⁶ time steps to reach an equilibrium polymer configuration. After switching on the motor interactions, the system was allowed to reach a steady state, and then quantities were determined by averaging over 3 × 10⁸ time steps in each of 50 independent simulations.

Additionally we measured two polymer shape parameters, the asphericity and prolateness. These are calculated from the principal moments of the gyration tensor, which is defined as

where rⁱ_n is the nth element of the position vector r_i of the ith monomer, and R^cm_n is the nth element of the centre of mass position vector R_cm. Specifically the asphericity is definedand the prolatenesswhere λ₁ > λ₂ > λ₃ are the principal moments. The radius of gyration can also be defined as R_g² = λ₁² + λ₂² + λ₃².⁴⁷

Acknowledgements

This research was partly funded by the European Training Network COLLDENSE (N2020-MCSA-ITN-2014) grant number 642774. DMa, DMi and CAB acknowledge funding from ERC (CoG 648050 THREEDCELLPHYSICS). DMi, CNL and EL would also like to acknowledge the contribution and networking support by the “European Topology Interdisciplinary Action” (EUTOPIA) CA17139.

References

1. T. Sanchez, D. T. N. Chen, S. J. Decamp, M. Heymann and Z. Dogic, Nature, 2012, 491, 1–5.
2. B. Alberts, A. Johnson, J. Lewis, D. Morgan and M. Raff, Molecular Biology of the Cell, Taylor & Francis, 2014, p. 1464.
3. C. R. Calladine, H. Drew, F. B. Luisi, A. A. Travers and E. Bash, Understanding DNA: the molecule and how it works, Elsevier Academic Press, 1997, vol. 1.
4. T. B. Le, M. V. Imakaev, L. A. Mirny and M. T. Laub, Science, 2013, 342, 731–734.
5. T. Terakawa, S. Bisht, J. M. Eeftens, C. Dekker, C. H. Haering and E. C. Greene, Science, 2017, 358, 672–676.
6. C. A. Brackley, J. Johnson, A. Bentivoglio, S. Corless, N. Gilbert, G. Gonnella and D. Marenduzzo, Phys. Rev. Lett., 2016, 117, 018101.
7. C. Naughton, N. Avlonitis, S. Corless, J. G. Prendergast, I. K. Mati, P. P. Eijk, S. L. Cockroft, M. Bradley, B. Ylstra and N. Gilbert, Nat. Struct. Mol. Biol., 2013, 20, 387–395.
8. M. Lenz, Phys. Rev. X, 2014, 4, 1–9.
9. C. P. Broedersz and F. C. MacKintosh, Rev. Mod. Phys., 2014, 86, 995–1036.
10. E. Tjhung, D. Marenduzzo and M. E. Cates, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 12381–12386.
11. K. M. Herbert, W. J. Greenleaf and S. M. Block, Annu. Rev. Biochem., 2008, 77, 149–176.
12. J. Michaelis and B. Treutlein, Chem. Rev., 2013, 113, 8377–8399.
13. A. M. Ganji, I. A. Shaltiel, S. Bisht, E. Kim, A. Kalichava, C. H. Haering and C. Dekker, Science, 2018, 7831, 1–9.
14. V. Bianco, E. Locatelli and P. Malgaretti, Phys. Rev. Lett., 2018, 121, 217802.
15. D. Saintillan, M. J. Shelley and A. Zidovska, Proc. Natl. Acad. Sci. U. S. A., 2018, 115, 11442–11447.
16. J. Harder, C. Valeriani and A. Cacciuto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 062312.
17. H. Jiang and Z. Hou, Soft Matter, 2014, 10, 1012–1017.
18. R. E. Isele-Holder, J. Elgeti and G. Gompper, Soft Matter, 2015, 11, 7181–7190.
19. T. Eisenstecken, G. Gompper and R. G. Winkler, Polymers, 2016, 8, 304.
20. S. K. Anand and S. P. Singh, Phys. Rev. E, 2018, 98, 042501.
21. R. G. Winkler, J. Elgeti and G. Gompper, J. Phys. Soc. Jpn., 2017, 86, 101014.
22. V. Schaller, C. Weber, C. Semmrich, E. Frey and A. R. Bausch, Nature, 2010, 467, 73.
23. T. Sanchez, D. T. N. Chen, S. J. Decamp, M. Heymann and Z. Dogic, Nature, 2012, 491, 1–5.
24. F. Nedelec and D. Foethke, New J. Phys., 2007, 9, 427.
25. S. L. Freedman, S. Banerjee, G. M. Hocky and A. R. Dinner, Biophys. J., 2017, 113, 448–460.
26. A. Ravichandran, G. A. Vliegenthart, G. Saggiorato, T. Auth and G. Gompper, Biophys. J., 2017, 113, 1121–1132.
27. N. Gupta, A. Chaudhuri and D. Chaudhuri, Phys. Rev. E, 2019, 99, 042405.
28. L. Coronel, E. Orlandini and C. Micheletti, Soft Matter, 2017, 13, 4260–4267.
29. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
30. D. Michieletto, Soft Matter, 2016, 12, 9485–9500.
31. T. Saito and T. Sakaue, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 042606.
32. L. D. Landau, A. Kosevich, L. Pitaevskii and E. Lifshitz, Theory of Elasticity, Butterworth-Heinemann, Oxford, 1986.
33. A. Wynveen and C. N. Likos, Soft Matter, 2010, 6, 163–171.
34. P. E. Rouse, Jr, J. Chem. Phys., 1953, 21, 1272.
35. M. E. Cates, Rep. Prog. Phys., 2012, 75, 042601.
36. R. Matthews, A. a. Louis and C. N. Likos, ACS Macro Lett., 2012, 1352–1356.
37. P. Poier, C. N. Likos and R. Matthews, Macromolecules, 2014, 47, 3394–3400.
38. L. Tubiana, E. Orlandini and C. Micheletti, Prog. Theor. Phys. Suppl., 2011, 191, 192–204.
39. L. Tubiana, E. Orlandini and C. Micheletti, Phys. Rev. Lett., 2011, 107, 188302.
40. L. Tubiana, G. Polles, E. Orlandini and C. Micheletti, Eur. Phys. J. E: Soft Matter Biol. Phys., 2018, 41, 72.
41. A. Bates and A. Maxwell, DNA topology, Oxford University Press, 2005.
42. J. Smrek and K. Kremer, Phys. Rev. Lett., 2017, 118, 098002.
43. D. a. Head, W. Briels and G. Gompper, BMC Biophys., 2011, 4, 18.
44. S. P. Pearce, M. Heil, O. E. Jensen, G. W. Jones and A. Prokop, Bull. Math. Biol., 2018, 80, 3002–3022.
45. C. A. Brackley, A. N. Morozov and D. Marenduzzo, J. Chem. Phys., 2014, 140(13), 135103.
46. P. G. de Gennes, Scaling Concepts in Polymer Physics, Cornell University Press, 1979.
47. A. Narros, A. J. Moreno and C. N. Likos, Macromolecules, 2013, 46, 3654–3668.

---