Open Access Article
This Open Access Article is licensed under a Creative Commons Attribution-Non Commercial 3.0 Unported Licence

Non-equilibrium effects of molecular motors on polymers

M. Foglino a, E. Locatelli *b, C. A. Brackley a, D. Michieletto ac, C. N. Likos b and D. Marenduzzo a
aSUPA, School of Physics and Astronomy, University of Edinburgh, Edinburgh EH9 3FD, UK
bFaculty of Physics, University of Vienna, Boltzmanngasse 5, A-1090, Vienna, Austria. E-mail:
cMRC Human Genetics Unit, Institute of Genetics and Molecular Medicine, University of Edinburgh, Edinburgh EH4 2XU, UK

Received 7th February 2019 , Accepted 25th June 2019

First published on 11th July 2019

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,1 RNA and DNA polymerases,2 helicases3 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 networks8,9 or active fields which mimic cytoskeleton dynamics.10 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.2

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 studies14–20 have considered a chain of beads, where each bead experiences an active force, either directed along the chain bonds17,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 motion17 or to spiral configurations18 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.28 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 KBEND 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 ri and the orientation is ûi, then the potential acts to minimise the angle ϕi between ûi and |ri+1ri|).
image file: c9sm00273a-f1.tif
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

image file: c9sm00273a-t1.tif(1)
Here r is the separation of the motor and the polymer bead, and KMOT 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 rm[thin space (1/6-em)]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 rmcut 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 KMOT), 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 software29 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 kBT and image file: c9sm00273a-t2.tif 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 KMOT, 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 KMOT 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 105τ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)).
image file: c9sm00273a-f2.tif
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 Lp 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 Lp(x = 〈nf2) ∼ xα with α = −0.108 ± 0.007. (d) Similar plot showing Rg as a function of 〈n〉. Again the inset shows the collapse of the curves when plotted as a function of 〈nf2 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 Rg defined by

image file: c9sm00273a-t3.tif(2)
where ri is the position of the i-th polymer bead, Rcm 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 Lp as defined from the bond–bond correlation function
bi·bi+s〉 = es/Lp,(3)
where bi = (ri+1ri)/|ri+1ri| represents the normalised tangent to the polymer at bead i. Here the average is additionally evaluated over bead index i. To calculate Lp 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 Rg and Lp exhibit a marked decrease with increasing 〈n〉, in agreement with our initial qualitative observation. Furthermore, the decrease in effective persistence length Lp is found to have a power-law scaling

Lp ∼ (〈nf2)α,
with α ≈ −0.1, which can also be used to collapse the different curves for Rg 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 〈nf2.

In order to shed more light onto the effect of the motors on the metric exponent ν relating the polymer size to its length (RgLν), in Fig. 3 we plot Rg 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).

image file: c9sm00273a-f3.tif
Fig. 3 Radius of gyration of the polymer substrate as a function of substrate lengths L, for (a) a fixed motor force (/kBT = 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.

image file: c9sm00273a-f4.tif
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 〈nf2 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) ≡ 〈bi·bi+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”30 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) = axb 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.

image file: c9sm00273a-f5.tif
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 〈nf2. (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 Lp. Inset: Position of the bond correlation minimum, Lh (defined as the value which satisfies β(s = Lh) = βmin), as a function of the polymer persistence length Lp. (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 Lh 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.31 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 KBEND (recall that in the absence of motor activity, the persistence length LpKBENDσ/kBT). To this end, we perform three sets of simulations with KBEND = 10, 20 and 30kBT, 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 KBEND 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 KBEND.

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
MSDCM(t) = 〈[Rcm(t0 + t) − Rcm(t0)]2〉,
where the average is over start times t0, and Rcm(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 model34 MSDCM(t) = (6D0/L)t, where D0 = kBT/ξ is the diffusion coefficient of a single monomer. This behaviour is recovered at all times for 〈n〉 → 0.

image file: c9sm00273a-f6.tif
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 kB−1. 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., MSDCM(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 bead14,20 and in Active Brownian Particles (ABPs).35 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 Pknot. It has recently been shown28,36,37 that Pknot displays a non-monotonic dependence on the stiffness of the substrate with a maximum at around Lp ≃ 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 Pknot will be located at larger values of the intrinsic persistence length of the substrate, lp (to be clear, the lower case lp 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 lp by adopting the algorithm described in ref. 38 and 39 and publicly available at By comparing Pknot(lp) for cases with and without motors we find that, surprisingly, the location of the maximum of Pknot is unchanged by the action of the motors, i.e., lp ≃ 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.41

image file: c9sm00273a-f7.tif
Fig. 7 Molecular motors acting on a polymeric substrate reduce the knotting probability Pknot(lp) of a free chain in a good solvent. Here lp 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 104 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 〈nf2, revealing a power-law scaling Rg ∼ (〈nf2)−0.1 and Lp ∼ (〈nf2)−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 RgLν) 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 Rg and Lp 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.14 Those studies, however, did not find hairpins (the decay in bond correlations remained exponential20), 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, Pknot(lp) 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.42 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 microtubules43,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 software29 to perform Brownian dynamics simulations, where the position of bead i is determined by the stochastic differential equation
image file: c9sm00273a-t4.tif(4)
where mi 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′)〉 = δαβδ(tt′).(5)
The potential Ui is a sum of interactions between bead i and all other beads, as described below. Fi 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 mim = 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 with combining circumflex]i,[v with combining circumflex]i and ûi which form a right-handed set of axes. A similar stochastic differential equation describes the dynamics of the orientation
image file: c9sm00273a-t5.tif(6)
where ωi is the angular velocity of bead i, and ξr,i and Ii 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 Ti 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 potential
image file: c9sm00273a-t6.tif(7)
where ri,i+1 = |riri+1| is the separation of the beads, and the first term is the Weeks–Chandler–Andersen (WCA) potential
image file: c9sm00273a-t7.tif(8)
which represents a steric interaction preventing adjacent beads from overlapping. In eqn (8)dij 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)R0 and KFENE are the maximum extension and strength of the bond respectively, and we use R0 = 1.6σ and KFENE = 30kBT throughout.

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

UBEND(θ) = KBEND[1 − cos(θ)],(9)
where θ is the angle between the three beads as given by
image file: c9sm00273a-t8.tif(10)
and KBEND is the bending energy. In the absence of motors the persistence length in units of σ is given by Lp = KBEND/kBT. Unless otherwise stated we set KBEND= 20kBT, 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 Lp = 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 kBT with temperature T = 298 K, one can map simulation time to physical units through the Brownian time τB = σ2/D, using the Stokes–Einstein relation to find the diffusion constant for a sphere D = kBT/(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

UORIENT = KORIENT[1 − cos(ϕi)],(11)
where ϕi is the angle between ûi and the direction of the polymer backbone such that
image file: c9sm00273a-t9.tif(12)
We set the orientation energy KORIENT = 30kBT. 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 rmcut = σ. 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 fendKMOT/σ 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 eKMOT/kBT. 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 KMOT and the motor force f considered in our simulations. We observe that 〈n〉 strongly depends on f, as well as N and KMOT. For /2 < KMOT, 〈n〉 increases with KMOT as the residence time is still proportional to the equilibrium one; we call this the thermodynamic unbinding regime. If /2 ≈ KMOT, 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 KMOT, 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 fnN (Fig. 8b, inset) we observe non-monotonic behaviour: the total force is maximal at f ∼ 15kBT/σ.

image file: c9sm00273a-f8.tif
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, KMOT = 10kBT.

Appendix C: scaling argument for Rg

In the main text we have reported that the radius of gyration scales as Rg ∼ (〈nf2)α with α = 0.1. We can link this observation with the scaling of Rg 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 = kBT/f. In analogy with traditional scaling arguments, we will take this to be the typical blob size.46 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 ∼g3/5). On a larger scale, the size of the chain will grow as Rgmν 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.,
image file: c9sm00273a-t10.tif(13)
As mentioned above, the size of the polymer then scales as
image file: c9sm00273a-t11.tif(14)
where L is the total number of monomers in the chain. By combining these equations, we arrive at Rgf(5/3)ν−1. Comparing this with the numerical result from the simulations, Rgf−0.2 (Fig. 2(b), inset), we obtain
image file: c9sm00273a-t12.tif(15)
or ν = 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 Rg, the bond–bond correlation function β(s) (from which we obtain the effective persistence length Lp), 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 106 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 × 108 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

image file: c9sm00273a-t13.tif(16)
where rin is the nth element of the position vector ri of the ith monomer, and Rcmn is the nth element of the centre of mass position vector Rcm. Specifically the asphericity is defined
image file: c9sm00273a-t14.tif(17)
and the prolateness
image file: c9sm00273a-t15.tif(18)
where λ1 > λ2 > λ3 are the principal moments. The radius of gyration can also be defined as Rg2 = λ12 + λ22 + λ32.47


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.


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

This journal is © The Royal Society of Chemistry 2019