Mathematical model suitable for efficient simulation of thin semi-flexible polymers in complex environments

We present an alternative approach to simulations of semi-flexible polymers. In contrast with the usual bead-rod compromise between bead-spring and rigid rod models, we use deformable cylindrical segments as basic units of the polymer. The length of each segment is not preserved with end points diffusing under constraints keeping the polymer chain nature intact. The model allows the simulation of tension transport and elasticity properties. In particular we describe a new cooperative regime in the relaxation of the polymer from its fully elongated configuration.


Introduction
Ever increasing interest in biological systems from the physical community has been pushing the boundaries of both theoretical and numerical methods. One specific area that is undergoing such a development is the physics of biopolymers with DNA as its most famous example. Other important examples of biopolymers include tubulin and F-actin. A common characteristic of the polymers mentioned above is that they belong to the class of thin long semi-flexible polymers, e.g. the width of F-actin is d = 6 nm, its typical length being L E 1 mm and its persistence length l p = 17.2 mm, for other examples see Andrews. 1 Especially the small width-length ratio poses a remarkable challenge for simulations of systems consisting of such polymers. Moreover, these polymers naturally occur in dense highly connected active networks in the cytoskeleton, an environment that is rich in molecular motors, which exert forces exceeding the thermal forces due to collisions with particles in the solution. This leads to additional requirements on the stability of the used models. In this paper we present a model suitable for simulating polymers under these conditions.
The traditional approach to simulate polymers is based on bead-spring models like the Gaussian chain, the Zimm or Rouse model. [1][2][3][4] The main advantage of such models is their relative simplicity and, at least in case of the Gaussian chain, many results can be obtained analytically. 2,5 Despite many successful applications of such models, e.g. Jendrejack et al. 6 and Mohammadinejad et al., 7 the main disadvantage is the large number of degrees of freedom necessary for the correct description of long thin polymers. In complex environments, the bead size and spacing between the beads has to be comparable with the diameter of the polymer, e.g. for L = 1 mm long actin filaments we need to have approximately 150 beads to represent a single polymer. This fact renders such models infeasible for the simulation of systems consisting of a large number of polymers. An interesting attempt to address this weakness of bead-spring models was presented by Panja et al., 8 who decompose the dynamics of a polymer into the dynamics of its normal modes. However the open question remains how to include large local forces, e.g. caused by molecular motors, in such a framework as they become non-local in the representation by normal modes. Another attempt at improving bead-spring models include modifications of the integration scheme. 9 Another approach used to simulate large polymer networks is to consider the polymers as rigid rods. 10 The main advantage of such an approach is to have faster simulations at the expense of neglecting deformations of single polymers. However as forces produced by molecular motors are of the order of JFJ p pN, 11,12 i.e. they exceed the usual thermal fluctuations, the assumption of rigidity of the polymers based solely on their persistence length may not suffice to correctly describe the behaviour of polymers in more violent environments, like in the cytoskeleton.
The class of bead-rod models is a natural compromise between bead-spring models and rigid rod models. It allows to preserve some details of the dynamics of single polymers, while it reduces the amount of necessary degrees of freedom. The simplest and most commonly used bead-rod model is a bead-spring model without hydrodynamical interactions and where collisions between polymers are treated based on a cylindrical rather than a spherical geometry. Random forces and friction applied to the beads also take into account the cylindrical geometry of the segments. 12,13 Another issue with this approach is that it is necessary to properly renormalise the spring constant with the increasing length of segments. 14 A different approach is taken in Somasi et al. 4 or in Montesi et al. 15 where instead of using bead-spring models, constraints are introduced to preserve the length of each segment.
In this paper we present an extensible cylinder chain model for polymers as an alternative to bead-rod models. The main difference with the bead-rod models is that we associate the degrees of freedom with extensible cylindrical segments rather than with beads, even though we represent the segments by the positions of their end points. The equations of motion for those points correspond to overdamped diffusion of a single cylinder, 16,17 which can be seen by the presence of two independent thermal noise components per cylinder segment corresponding to both translational and rotational-extensional degrees of freedom. This is fundamentally different from bead-rod models. 4,12,13,15 Lastly these segments are connected into chains by constraints posed on the positions of the end points of neighbouring segments. The length of a single segment is not preserved during the time evolution in our model, hence our model is more closely related to models used by Kim et al. 13 or Gordon et al. 12 with a finite tension propagation speed.
The paper is organized in the following manner: In Section 2 we derive from first principles the overdamped diffusion of a single stiff cylinder representing a single segment of the chain representing the whole polymer. In the next section, Section 3, we introduce constraints on the motion of single segments and the bending model thus concluding the presentation of the polymer model. Section 4 is the summary of our model in a form of a simulation scheme. Section 5 contains some of the basic tests of our model, including a comparison to the wormlike-chain (WLC). In Section 6 we apply our model to study the relaxation of polymers towards equilibrium and observe a new cooperative regime previously not discussed in the literature.

Overdamped diffusion of thin stiff cylinder
Our starting point is underdamped diffusion of rigid bodies submerged in an environment with constant temperature and viscosity. Even though there have been numerous studies on this subject, also in the context of rod-like polymers, 2,16,17 here we present an alternative approach suited for performing simulations of very thin polymers, where the diameter of the polymer is almost negligible compared with the segments length, in an environment where large forces will be generated locally and/or high polymer densities are realised.
We start with an overview of diffusion of rigid bodies, namely with the discretisation of a rigid body O into elements of constant volume V i , see Fig. 1. Then we identify the relevant macroscopic degrees of freedom and extend the description to include deformations. In the end we take the continuous limit V i -0 + to obtain macroscopic equations describing the underdamped diffusion of stiff bodies followed by taking the overdamped limit.
Each discrete element i is characterized by its position x i , velocity v i , volume V i , mass m i and, in case it belongs to the surface of the rigid body qO, by its normal vector to the surface n i .
The time evolution of each internal element, i e qO, is governed by the set of Lagrange's equations of the first kind dx i = v i dt, where F tot i ({x j },t) is the total force applied on the i-th element. It contains all the external forces applied to that particular volume element and also all virtual forces associated with constraints representing the rigidity of the body, i.e. 8i, j: Jx i À x j J = d ij , where the distance d ij is time-independent.
We assume that the time evolution of surface elements, i A qO, can be described by the Langevin equation where F tot i ({x j },t) is again the total force applied to the element i, g i 0 ({x j },{n j }) is the friction matrix associated with the portion of the Stokes drag corresponding to the given element. Although the Stokes drag is non-local as it depends on the overall shape and orientation of the body moving through the fluid, we can still associate locally applied pressure to each element with the corresponding portion of the drag. The thermal white noise dW i 0 (x i ,n i ,t) corresponds to fluctuations caused by collisions with molecules of the surrounding liquid. Consequently, the noise applied to the element i always assumes the direction of the inner normal Àn i . We further assume that the noise applied to different discrete elements is uncorrelated, hence We remark here that in the whole paper, we are working within the framework of the Itô formalism. 18 Note that the only difference between rigid and deformable bodies on this level of description is in the exchange of virtual forces representing the constraints by actual forces representing the elastic properties of deformable bodies.
Technical remark. From now on, in order to make the text more clear, we will not explicitly list all the dependencies for derived quantities as in this introduction, but only the relevant ones, e.g.

Underdamped diffusion of a stiff cylinder
We apply the general framework introduced above to a thin, stiff cylinder. First, we consider the cylinder to be rigid and we will later release this constraint. The cylinder has an obvious set of symmetries: it is centrally symmetric as well as rotationally symmetric with respect to the axis t. We assume the cylinder to be homogeneous, hence its centre of mass has to be at the same position as the symmetry centre X, i.e. it has to lie on the axis at one half of the cylinder's length l. For these reasons, the natural macroscopic degrees of freedom consist of those describing the movement of the centre of mass and rotations around it. By considering the cylinder to be thin we mean that any torque applied to the cylinder is negligible compared to the friction opposing it, so we can simplify our model just to include rotations around any axis perpendicular to the axis of the cylinder. In this assumption we deviate from the approach taken by Dhont 17 and we will see later that it has a great effect on the relative simplicity of the eqn (14), describing the overdamped diffusion of a cylinder.
A consequence of this assumption is that we have to discretise the cylinder only along its axis, see Fig. 2. In that case the time evolution of all elements needs to be described by (1), although we have to replace the noise dW i 0 with the one integrated along the element's surface dW i and subsequently the friction matrix g i 0 by its effective counterpart g i . From symmetry then follows where by x i E AElt/2 we denote the first and last discrete element representing the cylinder in the direction of t.
Now we extend our description to include deformations. As we are trying to describe a single segment much shorter than the persistence length of the overall polymer l { l p we can neglect the bending of the cylinder and we do not need to apply the renormalisation procedure for the segment stiffness introduced by Gutjahr et al. 14 as its effect will be negligible. The stiffness of the cylinder, characterized by its high Young's modulus E, is entering our equations twice. Firstly, the deformations in the direction perpendicular to the cylinder's axis are considered negligible as the width of the cylinder is small in comparison with the cylinders length. Secondly, the relaxation of the tension s along the axis of the cylinder to a uniform profile occurs on a shorter time scale than the relaxation of velocities and positions, i.e. t relax (s) { t relax (v i ) { t relax (x i ). Hence we assume that the tension profile is homogeneous at every instant of time, so we need to add just a single degree of freedom in order to describe the longitudinal deformations, namely the relative stretching velocity l defined as the longitudinal component of the relative end points velocity per unit length, The velocity of each element v i can then be expressed as where V is the velocity of the centre of mass X and x is the angular velocity perpendicular to the axis of the cylinder, xÁt = 0. Such decomposition of the velocity automatically obeys all the necessary constraints on the thin stiff cylinder, hence we do not need to add any virtual forces to our description any more. In order to obtain a set of equations for the macroscopic degrees of freedom we insert (2) into (1) and sum the set of equations with weights P dv i , P (x i À X) Â dv i , P (x i À X)Ádv i . These weights naturally decompose the motion of the cylinder into the translational, rotational and deformational degrees of freedom respectively. Let us define the total force F which is responsible for translations of a segment and the total Fig. 2 The discretisation of a homogeneous thin cylinder into cylindrical slices along the cylinder's axis determined by the orientation vector t pointing in the direction from the Àend to the +end of the cylinder. The length of the cylinder is denoted as l and the position of its centre of mass is denoted as X.
Rotations around the axis are neglected hence the angular velocity x is always assumed to be perpendicular to the axis of the cylinder, tÁx = 0. effective ''torque'' T responsible for a segment's rotations and deformations as where r i is the distance from the centre of the mass, r i = (x i À X)Át. We need to stress here that the effective ''torque'' T is not a torque in its standard physical meaning. Its components are parallel with the applied force and, as will be demonstrated later in eqn (4), it fully characterizes both the moment of force and stress in the longitudinal direction. We also define effective friction matrices where g trans is the friction matrix with respect to translations. Due to the symmetry of cylinders it has only two independent components g > trans and g J trans . The deformational-rotational friction matrix g dr has a perpendicular component g > rot which is given by the friction coefficient for rotations around the axis perpendicular to the cylinder's axis and a parallel component g J def that corresponds to the damping of the degrees of freedom associated with longitudinal deformations. The fact that those seemingly uncoupled degrees of freedom, namely rotations and deformations, can be assembled into a single friction matrix g dr stems from the assumption that the width of the cylinder and its associated degrees of freedom are considered to be negligible. For more details see Appendix A, where we have obtained the resulting equations of motion where M is the total mass of the cylinder and W (a) t , W (b) t are independent Wiener processes. The first equation links the velocity associated with translational degrees of freedom to the total force F. The second equation links the perpendicular component of the angular velocity x to the perpendicular component of the effective ''torque'' T. The last equation links the deformations to the parallel component of the effective ''torque'' T. We can identify the second term in square brackets on the right-hand side in eqn (4b) as the Coriolis force contribution. Similarly, the second term in square brackets on the right-hand side in eqn (4c) corresponds to the centrifugal force.
To have a closed set of equations we also need the following equations see Appendix A again.
At this point we can substitute the friction coefficients obtained from the discretisation with the ones based on the actual Stokes drag on the cylinder 19 where Z is the viscosity of the environment, d is the diameter of the cylinder and c J , c > , c r are end-correction terms. 2,[20][21][22] The remaining friction coefficient g J def will later, in the overdamped regime, be connected with the relaxation time of longitudinal perturbations t l .

Stiffness model
It is easy to determine the resulting total force F and effective ''torque'' T from external forces acting on particular points along the cylinder. The situation is a bit more complicated for internal elastic forces between the cylinder's elements. Here we demonstrate a general method of how to incorporate internal forces to the model in the very simple case of quadratic potentials between two neighbouring discrete elements where N is the number of discrete elements, k is the stiffness, l 0 is the equilibrium length of the cylinder and x N À x 0 = lt. The stiffness k is related to the Young's modulus E by where r is the cylinder radius. As we will show further on the scaling of the stiffness k by the total number of elements N is necessary for having a finite non-zero total elastic potential in the continuous limit N -N. First let us write the contribution from the elastic forces F el i to a given element i If we assume that cylinders are homogeneous at each time again, i.e. the corresponding relaxation time for non-homogeneous longitudinal deformations is much shorter than the typical time scale on which we observe our system, we obtain a much simpler formula for the potential = i V iiÀ1 = À= i V ii+1 = k(l À l 0 )t and consequently for the elastic forces Hence elastic forces do not contribute to the total force applied on the cylinder, but they are included in the effective ''torque'' This result is in good agreement with (4), as the elastic potential, without the help of other external force fields, cannot cause a displacement of the centre of mass nor a rotation of the object.
If we take the continuous limit of the total elastic potential we can see that the contribution to the effective ''torque'' can also be obtained from Furthermore, this expression holds for any other potential for which the discretisation is homogeneous along the cylinder, e.g. FENE-type potentials.

Representation by end points
Although eqn (4) and (5) fully describe the underdamped diffusion of a single segment, it is more convenient for simulations to represent the motion of the whole segment by the motion of its end points, here denoted as ''plus'' x + and ''minus'' x À end points, see Fig. 2. It follows from (5) that The dynamical variables are then given in terms of positions and velocities of the segments' end points as where we can see that x by definition obeys xÁt = 0.
Using these expressions along with (5) we can combine all equations of motion of a segment (4) into those describing the motion of its end points where l, g trans (t) and g dr (t) are now function of the positions of the end points (10). We see that there is no longer a contribution from centrifugal and Coriolis forces as we are now describing the segments in an inertial frame of reference. Hence, the fact that the two points belong to the same segment is represented only by the friction terms and by autocorrelations between the velocities V = (v + , v À ), In this section we have derived the underdamped diffusion of stiff segments in terms of velocities (11) and positions (9) of its end points. The biggest advantage of this approach is that the resulting equations fully describe the system without any additional constraints, hence it is very convenient to use it as a basis for numerical simulations.

Overdamped limit
The simplest method to obtain an overdamped limit is by considering the acceleration terms Mdv X to be negligible in (11). By this we obtain balance equations for all forces and moments which form a set of linear equations that can be solved with respect to the velocities. 17 However such a simple method does not yield the correct result as our friction matrices depend on the orientation of the segment. Thus we need to proceed in a fashion similar to that of diffusion in a nonhomogeneous environment, 23 where the inhomogeneity is caused by the shape of the cylinder instead of external factors in the environment. A more refined method, 24 suitable to more complex environments or configurations, starts from the assumption that canonical momenta have much shorter relaxation times than positions, t p { t x . Then by observing the system on the time scale of positions, i.e. the slow degrees of freedom, we can expand the evolution equation in the degree of separation e = t p /t x and obtain an autonomous Markovian dynamics for slow degrees of freedom up to linear order in e, for technical details see Appendix C. Moreover in many cases the estimate for the relaxation time of canonical momenta can be obtained during the process, thus providing an explicit lower bound for the time step in simulations, see eqn (51).
Following the second method and using that g trans and g dr commute (3) and are invertible we obtain the Smoluchowski equation for the joint probability distribution m t (x + ,x À ) where = + is the gradient with respect to the ''plus'' coordinate x + and = À is the gradient with respect to the ''minus'' coordinate x À . The Smoluchowski equations can be equivalently represented as a set of overdamped Langevin equations, see Appendix B, The set of Langevin eqn (14) describing overdamped diffusion of a single stiff thin segment is our first main result. The additional terms in the Langevin eqn (14) with respect to the underdamped case (11) can be interpreted as a thermal effective ''torque'' whose sign depends on the ratio of rotational friction with respect extensional friction. The thermal ''torque'' can be understood as an analogy to the term =ÁD in inhomogeneous diffusion in the Itô integration scheme, [25][26][27][28][29][30] where in our case the main source of inhomogeneity is the dependency of the friction matrices on the orientation of the segment. Moreover this contribution cannot be obtained in traditional approaches to derive overdamped diffusion based on the balance of forces and moments. 17 Similar to the velocities in the underdamped regime (12) we can see that positions X = (x + ,x À ) in the overdamped limit are autocorrelated This is the main difference with bead-rod models presented by Kim et al. 13 and by Gordon et al., 12 where the noises applied to the beads are uncorrelated. Although bead-rod models 4,15 preserve correlations to some degree after the segments' length preserving constraints are applied, they do not agree with (16). Notably, in case the polymer consists of single segments in the bead-rod models 4,15 the '+' ends' and 'À' ends' displacements are correlated only in the direction along the axis of the rod as can be seen in their model's version of the autocorrelation matrix where z is the friction matrix used by Somasi et al. 4 In Montesi et al. 15 the z corresponds to z ii . 2.4.1 Average length. In order to connect the last remaining unknown friction coefficient g J def with the relaxation time of longitudinal perturbations t l we briefly analyse the time evolution of the average length of the segments hli. It follows from the Smoluchowski eqn (13) that the time evolution of the mean value of length l = Jx + À x À J is given by where the first term corresponds to the applied tension and the second term to the relaxation toward equilibrium with the ''relaxation time'' which is negative. Hence without any potential opposing it, the segments always tend to increase their length indefinitely independently of the sign of the ''thermal torque'', (15). Notice also that the mean value is taken with respect of the full probability distribution m(x + ,x À ).
In the case where the only component of the effective ''torque'' parallel with the segment's axis is the contribution from the elastic potential (8), we can conclude that an equilibrium distribution of the length r(l) is given by where Z is the partition function.
Quadratic stiffness potential. When there is a potential that is actually preserving a finite length of the segment l, the estimate for the relaxation time can be unrealistic as the first term on the right-hand side in eqn (18) cannot be neglected. To demonstrate this fact we provide a more realistic estimate for a quadratic potential, as was introduced in Section 2.2. First, we assume that we are in the low temperature regime where the system relaxes quickly to thermal equilibrium with equilibrium length l 0 . Consequently, we can assume that any fluctuation that arises in the system is small and hence the probability distribution is close the equilibrium one. This enables us to decompose the full probability distribution m t into the equilibrium distribution (19) and a small correction Dm t , m t = r + Dm t . Also, we can assume that the correction Dm t is well localized around l 0 . Hence, we can expand eqn (18) up to linear order in the length difference from which a new estimate for the relaxation time immediately follows This also demonstrates that the relation between the relaxation time of longitudinal perturbations and the friction coefficient is not universal as it depends on the actual model and the specific elastic potential that is used.
Remark: applications to nematic particles. Even though the focus of this article is mainly on the thin semi-flexible polymers, the model presented in this section is not restricted to this particular application and can be used to describe a behaviour of various physical and biological systems. The most prominent example might be elongated active swimmers like bacilli, which are a typical example of active nematic particles. Models for such systems where nematic particles have a constant speed, like the famous Vicsek model, [31][32][33] can be mimicked by simply adding a non-potential force in the form F = F 0 t.

Diffusion of stiff chain
In the previous section we have described the dynamics of single segments, which we intend to use as basic building blocks for the description of semi-flexible polymers. Therefore the typical length of these segments has to be much shorter than the persistence length of the polymer. Consequently, we present one of the possible bending models to reflect the semiflexibility of the polymer. Moreover in order to obtain a model for the full semi-flexible polymer we need to add suitable constraints on the independent movement of the segments. The presented approach of how to deal with such constraints is well known, 4,15 however the actual form is heavily model dependent.

Constraints
The basic constraint is that the positions of end points of neighbouring segments coincide at each instance of time, i.e. they have to move along the same trajectory, where by (i) we denote the index of the corresponding segment, see Fig. 3. Further on segments will be indexed from 0 to N À 1 as it is the convention used in C-like programming languages. Clearly these constraints are holonomic and as such they are in the Langevin formalism associated with Lagrange multipliers, which can be interpreted as virtual forces U (i) applied at the joints. In order to solve the equation of motion we need to determine these virtual forces from the constraints and include their contributions into the equation of motion.
Virtual forces contribute to the total force and the total effective torque in the Langevin equations of motion (14) like any other force where in order to avoid unnecessary verbosity we assemble the applied forces and effective ''torques'', Gaussian white noise and thermal effective ''torque'' (15) into the total force F (Á) tot respectively the total effective ''torque'' T (Á) tot . The last and first The constraint (21) is also valid for displacements, i.e. dx (i) À = dx (i+1) + , which provides us with the closure of eqn (22) and (23). Thus the virtual forces U (i) associated with constraints (21) are given as the solution of the set of linear equations Notice that the virtual forces, obtained as a solution of this set of equations, are in fact random variables that depend on the actual realization of all coupled white noises. Moreover if we can clearly see that the structure on the left-hand side corresponds to the matrix multiplication of the vector of virtual forces U (i) with the tridiagonal symmetric block matrix Moreover, the corresponding matrix is not necessarily positive definite. However, it is always invertible as any non-zero eigenvector corresponding to the zero eigenvalue would mean non-uniqueness in the resulting velocity profile and consequently in the non-uniqueness of positions. Such behaviour is deemed unphysical as one can expect that the problem is always well conditioned. From a computational point of view we can take advantage of algorithms like LDL Tdecomposition, which scale linearly in time with the polymer length, to find the desired virtual forces.

Bending models
The usual model describing the bending of polymers is the worm-like chain (WLC) model 34,35 with free energy where a is the bending stiffness which can be later directly associated with the persistence length l p , r(s) is the curvature radius parametrized by the length s and L is the total length of the polymer. One of the few measurable quantities that is analytically accessible from the WLC is the mean squared end-to-end distance 34 where the persistence length is defined as Depending on the ratio between the total length L of the polymer and the persistence length we can distinguish two regimes by observing the squared end-to-end distance. One is the rod-like regime L{l p where the thermal fluctuations are not strong enough to fully bend the polymer and hence the polymer will oscillate around the straight configuration, i.e. its mean squared end-to-end distance is close to the value for a straight polymer L 2 . The second regime is the collapsed regime L c l p . . .
where the ends of the polymer are locally freely diffusing. The simplest discretisation of the continuous WLC is the Kratky-Porod model 36,37 with interaction potential where t i is the normalized tangential vector of i-th segment.
However such a simple model has the disadvantage of enabling the segments to be aligned in an anti-parallel fashion. In the case of stiff semi-flexible polymers it is very unlikely that thermal fluctuations alone will ever cause such a perturbation, but in complex environments like the cytoskeleton we cannot exclude such possibility: external forces produced by molecular motors like Myosin-V can exceed thermal fluctuations by several orders of magnitude, 11 and hence in the extreme case it can cause an anti-parallel alignment of some neighbouring segments. This situation corresponds to two segments occupying the same physical space, which is unphysical. Moreover, it also corresponds to the situation where the constraints (24) introduced in the previous subsections are not uniquely solvable, thus making the simulation unstable. One can argue that forces exerted by molecular motors are mainly tangential to the polymer. However, in dense complex networks, micro-filaments of molecular motors like Myosin-II, that have heads on both ends of their polymer chain, 38 can be attached to two polymers at once. As two polymers generally are not parallel, the force exerted tangential to one polymer is translated to the other one as a force with a non-negligible perpendicular component.
Moreover, Myosin-II micro-filaments have a structure resembling a semi-flexible polymer. 38 Hence they represent another example of a system subjected to large forces and as such is suitable to be described by the model presented here and which is subjected to large forces, although in the scope of this article we restricted ourselves to more common examples of polymers like F-actin and DNA.
As our aim is to provide a model suitable for such violent environments we need to address this stability issue. An improved model is obtained by choosing a different discretisation scheme. Notably we assume a constant curvature between the centres of two neighbouring segments which is equal to the radius of the inscribed circle, while at the same time we assume the polymer length between the centre of two neighbouring segments to be constant and equal to l, see Fig. 4. This leads to the bending potential This particular choice of discretisation resembles the behaviour of the Kratky-Porod model for small angles, as we have shown in Fig. 5, and for angles close to p it diverges, thus preventing the anti-alignment.
It can be shown that the effective ''torque'' applied on a single segment and caused by a single joint (i,i + 1) is in general given by Fig. 4 Scheme for obtaining the improved bending model. The curvature is assumed to be constant between the centres of two neighbouring segments and is given as an inverted radius r of the inscribed circle, which is directly related to the angle between those segments' axes, t i Át i+1 = cos W.
The length of the corresponding polymer arc is always assumed to be l.  (26), case B is the quadratic approximation of the Kratky-Porod W 2 and case C is the improved bending model given by (27).

Soft Matter Paper
Open similarly it can be shown that the resulting forces are zero,

Simulation scheme
To summarize our model we present a step by step overview of the integration scheme of the model presented in this paper. Let us note that we employ the Euler-Maruyama integration scheme as a basic integration method.
(1) For each segment the random effective ''torque'' and force are generated according to eqn (14): where Dt is the used time step.
(2) For each segment the thermal effective ''torque'' T (i) thermal is evaluated according to the eqn (15).
(3) For each segment the elastic effective ''torque'' T (i) el is evaluated according to the eqn (8).
(4) For each joint the bending effective ''torques'' are evaluated according to eqn (28). (5) All other external forces applied to the segment, like the repulsion of polymers or forces exerted by molecular motors, not explicitly discussed in this paper, or the presence of optical tweezers, see Section 6, are evaluated and included where r (i) a is the projection of the point where the force is applied x a onto the segments axis Virtual forces per joint U (i) , see eqn (24), are evaluated based on the total forces and total effective ''torques'' per segment (7) Positions are updated according to eqn (22) and (23) by using the forward Euler(-Maruyama) scheme.
In order to be in the overdamped regime we need to make sure that the time step is bigger then the lower bound given by the relaxation time of canonical momenta, see eqn (51).

Time complexity
The most complex task in the basic simulation scheme is the determination of the virtual forces that are binding the polymer together (step 6). It was discussed in Section 3.1 that this task has a linear time complexity in the number of segments per polymer, hence the full simulation scheme overall scales linearly with the total number of segments presented in the simulation. Therefore, it belongs to the same complexity class of algorithms as the traditional bead-spring approach of the model presented by Montesi et al. 15 The main difference usually is the pre-factor which, when compared to the most simple bead-spring model, will be much larger in our case, because for each segment we need to generate two random numbers, as opposed to one in case of the bead-spring model (step 1), we need to perform the LDL T decomposition of the A ðiÞ AE (3 Â 3) matrices along with multiple matrix multiplications (step 6), etc. However the possibility to scale the length of the segments and thus decreasing the number of segments in the model can compensate the higher cost per single segment.
The key parameter of the system which determines the maximal possible advantage of our model with respect to the traditional bead-spring model is the ratio of the persistence length to the diameter of the polymer l p /d. This is caused by the fact that in the traditional bead-spring approach the beads need to be spaced much less than 2d apart in order to avoid the possibility of two polymers passing through each other, while the size of the segment is limited by the persistence length l p .
The second factor contributing to the efficiency of a given algorithm is how large the time step can be for a given set-up. The maximal allowed time step is determined by multiple factors, which can be split into two categories: the desired resolution of the simulation (both spatial and time) and the stability of the simulation. While the first category is determined solely by the considered set-up, the latter is determined only by physical properties and interactions in the system. Surprisingly, this means that the maximal time step also depends on the chosen segment length.
Some of these constraints are summarized for both the model presented here and the bead-spring model in Table 1, where the row labelled translations is the estimate of the maximal time step based on the requirement to prevent Table 1 Estimates for the upper bound of the time step Dt for different constraints each associated with some stability issue. Asymptotic behaviour of these estimates is show in the last column for large segment lengths l 0 c d

Constraint
Bead-spring New model

Model Exact Asymptotic
Translations translational thermal forces to cause two segments to pass through each other in single time step without colliding. From the stochastic nature of these forces, the additional parameter c s determining the confidence interval of the estimate, is needed. For further details we refer to Appendix D. The row labelled rotations is similar to that of translations but now with respect to thermal rotational forces. Diffusion is the estimate of the maximal time step for all thermal forces combined. This estimate is based on autocorrelations (16). g* is the effective friction matrix, see (56) for exact expressions. Stiffness in the third row corresponds to the upper bound on the time step. This stems from the requirement that the reaction of the elastic potential (7) to the perturbation cannot cause an even bigger perturbation. For bending in row four, a similar requirement applies but now for the action of the bending potential (27). From the asymptotic behaviour of the estimates on the upper bound of the time step we see that the maximal time step achievable, for an arbitrarily large segment length l 0 , is determined by the relaxation time of the longitudinal deformation t l . However in the typical case of stiff thin polymers these estimates are (much) higher than those limiting the bead-spring model mainly due to the increased friction, thus boosting the efficiency of the approach discussed here. In case pairwise interactions (e.g. Coulomb or Debye-Hückel electrostatic interaction) are present, the time to evaluate them depends quadratically on the number of elements, thus shifting the time complexity of the whole simulation up to a quadratic scaling in the worst case scenario. Hence purely by increasing the segment length, thus decreasing the number of elements in the simulation, can lead to a significant improvement, albeit sometimes at the cost of a more complex interaction potential.

Model performance
We demonstrate our simulation method by applying it on F-actin, which is a typical thin, d = 6 nm, stiff, l p = 17.6 mm, polymer. The common parameters of all simulations are listed in Table 2. The value relaxation time of longitudinal deformations for single segment t l was chosen arbitrarily as there, to the best knowledge of authors, are no experimental data available. The aim was to have the highest value possible in order to achieve large time steps while preserving the stability of simulation. After couple of trials, the value listed in Table 2  We can see that we are truly in the overdamped regime as the time step Dt is much bigger than any of the relaxation times. The mass of a single segment was assumed to be M = 1.26 MDa. All simulations presented here were created using the Mpacts framework. 39

Speed comparison
In the first set of tests we provide a practical comparison between the time complexity of the model proposed here and the bead-spring model with spacing of the beads l 0 = d. First we simulate the polymer with the contour length of 1.2 mm using the bead-spring model with 201 beads and using time step Dt = 10 ps for a total simulation time of T = 10 ms thus creating the baseline for comparison. Subsequently, we simulate the same polymer using the model introduced here with different segment lengths l 0 while maintaining the contour length and taking the same time step Dt = 10 ps and the same total simulation time T = 10 ms. The result of this comparison is shown in Fig. 6, where the speed-up is given as the ratio of the simulation duration † of our model and the bead-spring model. We can see that the maximal speed achieved was more than 10Â higher than for the bead-spring model, when the full polymer was represented by a single segment. From the same figure we can also obtain the information about the overhead of the algorithm presented in Section 4. We see that the simulation containing 120 segments is approximately as fast as the simulation of the bead-spring model with 201 beads thus being about 40% slower per element than the bead-spring model. Moreover the slower increase of the speed for longer segments, which deviates from the linear scaling, can be attributed to the overhead of the simulation toolbox for a low number of elements. Unfortunately this bias is unavoidable in this comparison, because increasing the number of elements for long segments to sufficient numbers, e.g. by increasing the number of non-interacting polymers, will lead to impractical times necessary for the simulation of the bead-spring model of the actin polymer.
As we have discussed in Section 4.1, the maximal possible time step also depends on the segment length. In order to illustrate this issue we have listed estimates for maximal time steps, see Table 1, for various segment lengths in Table 3, ‡ where we used the confidence interval 3s, i.e. c s = 3. There we can see that for short segment lengths the main limiting factor is the high bending rigidity. As the segment length increases the limitation caused by the bending rigidity quickly fades away due to its scaling as l 0 4 /ln(l 0 /d), see Table 1. Hence, in the intermediate range of segment lengths, diffusion becomes the main limiting factor for the maximal time step. Namely, thermal fluctuations of the segments' length combined with translational Brownian motion along the segments' axis become the main limiting factors in this range. For long segments, we reach the saturation point, where we cannot increase the time step over a certain threshold. In case of F-actin the threshold is given by the relaxation time for longitudinal deformations, Dt { 2t l . Hence the maximal speed-up theoretically achievable for the model presented in this paper by increasing the time step is more than 350Â, i.e. more than two orders of magnitude. This can significantly enhance the performance of the overall simulation as is demonstrated in Fig. 6, where after adjusting the time step according to Table 3, we have achieved the maximal simulation speed-up of more than 2000Â compared to the bead-spring model.

Diffusive properties of segments
Here we evaluate the rotational diffusion of free independent segments where translations of the segments are disabled, hence the equations of motion (14) simplify to It can be shown that the mean orientation of the segments relaxes towards equilibrium, i.e. to a uniform distribution of orientations, hti = 0, independent of the specific choice of the elastic potential (8) as where we can identify the rotational relaxation time This theoretical prediction is compared with the simulation of N = 728 independent segments, initially all oriented in the x-axis direction, hti t=0 = ê x , see Fig. 7. The count of segments N is sufficiently high in order to have a reasonably large statistical ensemble for evaluation. The rotational relaxation time is found to be t r = 2.685 AE 0.028 ms which is in good agreement with the theoretical prediction t = 2.676 ms. Moreover variances of the mean orientations' components converge to 1= ffiffi ffi 3 p , whose value is obtained from the uniform distribution of orientations.
The second test verifies the translational properties of the segments that are also undergoing rotational diffusion. Starting from the same initial configuration, i.e. hti t=0 = ê x , we let the particles evolve according to the complete equations of motion (14) with zero external force and torque. We shall observe the transition from inhomogeneous diffusion towards homogeneous diffusion as the orientation of the segments becomes random. Because the mean orientation dissipation is independent of translations we can obtain the mean square displacement over the full range of time Fig. 6 The log-log plot of the measured relative speed-up of the presented model for various segment lengths with respect to the bead-spring model based on the simulation of a single polymer with a contour length of L = 1.2 mm. The first data set uses the same time step Dt = 10 ps for all simulations, while the second data set increases the time step for each simulation according to the speed-up row of Table 3, thus preserving the minimal upper bound of the time step to the actual time step ratio. ‡ In the Table 3 we have used a more exact expression (59) for the maximal estimate on the time step based on the stability with respect to bending for the bead-spring model with l 0 = d.
where we have assumed an initial condition hti = ê x . In Fig. 8 we can see that the full time range solution agrees well with the simulation.

Equilibrium squared end-to-end distance of a single polymer
In a first simulation of the full polymer, we investigate the convergence to equilibrium of a free single polymer. The initial configuration is that of a straight polymer, 8i, j:t i = t j , with segment lengths l equal to the equilibrium segment length l 0 . In other words, the polymer's initial configuration corresponds to the configuration with the lowest possible energy or the configuration at T = 0 K. The polymer is then let to thermalise with the environment for 1 s. Let us note that for this test we artificially decreased the viscosity of the environment to Z = 1.71 Â 10 À4 pN ms nm À2 , i.e. ten times lower than the value listed in Table 2, in order to decrease the relaxation time. The mean squared end-to-end distance hR 2 i is estimated based on the ergodic average over the time interval of 1 s where the t k 's are evenly spaced snapshots of time with a period of 1 ms, hence t i = 1 s, t f = 2 s, N = 1000. Our simulation results are compared with theoretical prediction for the WLC (25), see Fig. 9. We can see that the results comply very well with the theoretical prediction for polymers shorter than the persistence length, L o l p . However, polymers longer than the persistence length depart from the theoretical prediction. The difference in standard deviation estimates suggests that the system is not yet relaxed to the thermal equilibrium. Indeed, from the plot of the relative difference between the theoretically predicted value (25) and the actual value or hR 2 i in time, Fig. 10, we estimate relaxation times which are larger than the overall simulation time for polymers longer than their persistence length.
This issue is not caused by the simulation method itself, but rather by the system we have chosen for this demonstration. For polymers much longer than their persistence length L c l p , the bending stiffness does not constrain possible configurations of the polymer, which means that the Rouse model represents a suitable approximation. On the other hand, polymers much shorter than the persistence length L { l p can be approximated by stiff rods. During the thermalisation, various sections of the polymers need to reorient themselves and thus the relaxation   (31) and full range (32) theoretical predictions. The mean square displacement is evaluated in all principal directions, i.e. parallel and perpendicular to the initial orientation t(t = 0) = ê x . Note that the long range and full range predictions will be asymptotically parallel for large times. Table 3 Upper bound on the time step (in ms) for various constraints and segment lengths with actual values for the F-actin model, see Table 2. The confidence interval used for diffusions is 3s, i.e. c s = 3. The minimal values for given segment lengths are denoted by the bold cursive. The speed-up refers to a relative speed-up with respect to the bead-spring model given as a ratio of the minimal upper bound for a given segment length and the minimal upper bound for the bead-spring model. For segments longer than 1.8 mm we have reached a maximal possible speed-up of 352Â time associated with the thermalisation will follow the behaviour of the rotational relaxation time, which oscillates between the Rouse model and a stiff rod depending on the contour length L in question. It is known that the rotational relaxation time of the Rouse model is proportional to the square of the contour length L 2 , cf. Doi and Edwards, 2 while for the stiff rod model the rotational relaxation time t r (29) is proportional to L 3 /ln L. Notice that the relaxation time scales at best quadratically L 2 with the contour length independently of the discretisation and the persistence length. Hence, if we choose a different system with a smaller persistence length, we may be able to overcome this issue by decreasing the contour length of the polymer while preserving the contour length to persistence length ratio L/l p . Indeed, if we choose DNA § instead of F-actin we are able to thermalise, within the simulation environment, polymers longer than the persistence length, see Fig. 11, and confirm a good agreement between the theoretical prediction for WLC and our simulations. We note here that the difficulties with simulations of the F-actin polymers longer than persistence length do not present a problem in studying real-life system, as the typical length of the F-actin polymer does not typically exceed 1 mm, cf. Biron and Moses, 46 and even in highly controlled environments, it rarely exceed its persistence length. 41,47,48

Single polymer under tension
Another simulation set-up is to have a single long polymer under tension, L = 23.04 mm. In particular the initial set-up of the simulation is the same as in the previous case, i.e. a straight polymer with segments lengths equal to equilibrium value l 0 and the viscosity also decreased to Z = 1.71 Â 10 À4 pN ms nm À2 . The polymer is then allowed to relax towards the thermal equilibrium for 10 s while constant anti-parallel forces are applied at both ends of the polymer, see Fig. 12. Estimates of the mean values of various equilibrium quantities are then again achieved by averaging over snapshots that are 1 ms apart over the time interval of 5 s from the initial time t i = 10 s, for reference see the ergodic average in the context of the mean squared end-to-end distance (33). The relevant quantity in this case is the projected end-to-end vector R in the direction of the force ê F , in this context called ''extension'', which for the WLC can be related to the Fig. 9 The mean squared end-to-end distance hR 2 i as a function of the contour length L of the polymer plotted on a logarithmic scale, obtained from the simulation of a single free polymer. Each point represents an ergodic average (33) over the timespan of 1 s. The simulated mean squared end-to-end distance fits the theoretical prediction, given by (25) for short lengths of polymers, well. The deviation from the theoretical prediction for polymers longer than the persistence length L 4 l p is caused by the increase in the relaxation time, as can be seen in Fig. 10. The vertical dotted line denotes the upper bound of the fitting range, which corresponds to the situation, where the contour length equals the persistence length. Fig. 10 Relaxation of the mean squared end-to-end distance towards its theoretical equilibrium value hR 2 i th on a semi-logarithmic scale. Each point represents an ergodic average (33) over the timespan of 50 ms. Estimated relaxation times from linear fits of the simulation data of a single polymer on the full interval are t r (23.04 mm) = 3.03 AE 0.35 s, t r (46.08 mm) = 5.51 AE 0.81 s which are exceeding the total simulation time. The shortest polymer's relaxation time estimate is overweighted by its error estimate, t r (11.52 mm) = 2.5 Â 10 4 AE 4.9 Â 10 7 s, which suggests that the polymer is already thermalised. Fig. 11 The mean squared end-to-end distance hR 2 i as a function of the contour length L of the polymer plotted on a logarithmic scale, obtained from the simulation of a single free DNA. Each point represents an ergodic average (33) over the timespan of 1 s. The simulated mean squared endto-end distance fits the theoretical prediction, given by (25)  magnitude of the applied force JFJ, by an approximate relation 1,5,35,45,49  : The above relation is only applicable to small applied forces where the bending dominates the behaviour, i.e. the polymer's behaviour is that of an entropic spring, especially for long polymers. For large applied forces, when the polymer is close to the straight configuration, the stiffness of individual segments prevails. In our case elastic properties of individual segments are given by the quadratic potential (7), so the asymptotic behaviour of the polymer corresponds to that of chained harmonic oscillators. Hence the force relation changes to where N is the number of segments in the polymer and k is the stiffness of a single segment. Theoretical predictions (34) and (35) are compared with the simulation data in Fig. 13 and 14. The error bars in Fig. 13 denote the fluctuations in the end-to-end distance of the polymer, which are caused by thermal forces bending the polymer. Bending becomes increasingly difficult when larger forces are applied and indeed we see diminishing fluctuations when the pulling force is increased. Using (34) we can determine the actual persistence length l est p = 16.38 AE 0.77 mm, which is in good agreement with the value from the set-up of the model l p = 17.2 mm. The estimate for the persistence length is obtained by the standard w 2 fitting method limited to the range hRÁê F i/L o 0.99, denoted by the dotted vertical line in Fig. 13. In order to demonstrate the transition between the two regimes of the polymer, we introduce an effective spring constant k eff of the polymer defined as Let us note that such an effective spring constant diverges at hRÁê F i = L, because we have chosen the straight configuration as the equilibrium value of such an effective spring, which also corresponds to polymer's configuration in the zero temperature limit, T -0 K. Although such choice is unphysical, it has its own merits in terms of well defined asymptotes for both small and large JFJ. The simulation results are presented in Fig. 14. For large forces the effective spring constant coincides with the effective spring constant for a chained harmonic oscillators (35), i.e. Fig. 12 Illustration of a polymer under the tension, where we apply forces F of equal magnitude F but opposite direction ê F at the ends of the polymer. R denotes the end-to-end vector and the upper grey contour shows the initial condition of the simulation, where the polymer is in a straight configuration and oriented in the direction of the applied force, i.e. R(t = 0) = ê F L.  Each point corresponds to ergodic average over the timespan of 5 s. In the small JFJ regime the actual exponent for the effective spring constant k eff p JFJ x was determined to be x = 1.534 AE 0.0074 which is in good agreement with the theoretical value x ¼ 3 2 , see (38). In the high JFJ regime the ratio between the theoretical value (37) and effective spring constant is k eff /k 0 = 1.0004 AE 0.0038. The crossover estimated from the intersection of the two linear fits is estimated to be JF cross J = 23.93 AE 0.55 pN. The dotted vertical lines denote the ranges for fitting of the parameters listed above. View Article Online which is also supported by the simulation data k eff /k 0 = 1.0004 AE 0.0038, where the estimate is based such data points for which JFJ 4 10 2 pN. For small forces, the effective spring constant will depend on the force. By making the use of the fact that in this regime the relation between the extension and the force (34) is dominated by its first term 35 we can recover the asymptotic spring constant for small forces From the simulation data based on points in the range of JFJ o 1 pN we obtain the exponent x = 1.534 AE 0.0074. Moreover, from the intersection of these asymptotes we can estimate the transition point to be at JF cross J = 23.93 AE 0.55 pN, which is comparable to the position of the peak in the plot, estimated to be around JFJ = 16.4 pN.

Relaxation dynamics towards thermal equilibrium
While in the previous set-ups we investigated the equilibrium properties of our model, in this section we aim to verify some aspects of the dynamics. The initial configuration of the system consist of 200 straight independent polymers of total length L = 23.04 mm with segment lengths equal to their equilibrium values l 0 . For each polymer one of its ends is trapped in a harmonic potential with spring constant k trap = k. The relaxation in equilibrium is then observed on different time-scales with different time steps in the range Dt A [10 À4 ,10 À2 ] ms. In Fig. 15 we plot the difference of the contour length from the initial configuration as a function of the total simulation time t on a log-log scale where ê t=0 is the normalized orientation of the polymer at the initial time t = 0. After the noisy period (region I in Fig. 15) we observe two different relaxation regimes. The exponent of those two regimes can be estimated to be x = 1.114 AE 0.010 in the first regime (region II) and x = 0.3436 AE 0.0030 in the second regime (region IV). We can see that polymer lengths first relax towards equilibrium almost linearly in time, while later the relaxation slows down and is close to the theoretical prediction x ¼ 1 3 for the WLC, 50,51 which is experimentally verified in case of DNA. 52 In order to understand the transition between those two relaxation regimes we investigate the mean energy stored in the bending (27) and elastic (7) potentials. To have some referential values we obtain analytical expressions for the mean value and variance of the bending and elastic potentials under the assumption that the associated degrees of freedom with those two potentials are decoupled and in thermal equilibrium. Furthermore we also consider neighbouring elements, i.e. nodes or segments, independent. The mean value of the bending energy (27) for a single joint is given by where E k (x) is the generalized exponential integral ¶ and l 0 is the segment's equilibrium length while its variance is We can see that the behaviour is determined by a single dimen- In our case x c 1, cf. Table 2, hence we expand the mean value and the variance in where we are close to the regime where the equipartition theorem is valid. Similarly, the mean value of the elastic energy (7) for a Fig. 15 Contour length difference D as a function of the simulation time, from simulations of 200 independent polymers, all initially in straight configuration of the total length L = 23.04 mm, with one end fixed in a harmonic potential with spring constant k trap = k, while the other end is free to diffuse. Several simulations with time steps chosen from the range Dt A [10 À4 ,10 À2 ] ms are evaluated in order to cover the relaxation in the broad range of time scales and also to demonstrate that the resulting behaviour is independent of a chosen time step. We distinguish two relaxation regimes, the first is the cooperative regime (region II) with the behaviour x = 1.114 AE 0.010, while the second is the uncoupled regime (region IV) with x = 0.3436 AE 0.0030, where respective fitting ranges for x are denoted by the vertical dotted lines. Error bars are omitted for clarity as for short times they exceed the full range of the plot. ¶ Generalized exponential integral is defined as E k ðxÞ ¼ Ð 1 1 t Àk expðÀxtÞdt. In case k = 0 the exponential integral reduces to E 0 (x) = e Àx /x. Asymptotic behaviour was provided for example by Olver. 53 single segment is where again we can define a dimensionless parameter y ¼ k' 0 2 2k B T determining the behaviour, which in our case is y c 1. Consequently the variance is given by Hence, we can again expand the exact results up to the leading order and obtain which again shows that we are close to the regime where the equipartition theorem is valid. The total mean energies for a single polymers are then obtained by multiplying with the number of nodes or segments in the polymer as we assume they are independent. Now we compare the theoretical equilibrium values with those obtained from the simulation during the relaxation. In Fig. 16 we observe that the elastic energy first relaxes towards thermal equilibrium (region I) with a relaxation time comparable to the relaxation time of a single segment t l . The unfreezing of transversal modes denoted by an increase of a mean bending energy hV b i initially causes an increase in the mean elastic energy over the expected equilibrium value (40) (region II). The reason for such behaviour in a straight configuration is that every bending event also causes stretching of the involved segments. At the same time, the straight configuration also prevents these perturbations to relax locally as a consequence of applied constraints (21). At later times, when the polymer departs from the straight configuration, the influence of constraints is weakened and transversal modes are decoupled from longitudinal modes, which can be seen by the return of the mean elastic energy to its equilibrium value (40) (region III). By comparing Fig. 15 with Fig. 16, we hypothesize that the increase in the power law exponent above the theoretical prediction x ¼ 1 3 is caused by a coupling of transversal and longitudinal modes leading to the ''cooperative'' regime. This goes against the assumption of their independence, hence invalidating the theoretical predictions in this regime provided by Hallatschek et al. 50,51 When those modes decouple further relaxation towards the equilibrium is governed solely by bending properties of the polymer (region IV) and hence reproduces the theoretical value of the power law exponent for the WLC. In order to verify theoretical predictions (40), (41) we compare the total mean elastic energy obtained from the simulation for Dt = 10 À2 ms, which is hE b i = 34.2 AE 6.2k B T. This is in reasonable agreement with the theoretical prediction hE b i eq = 32.0 AE 5.7k B T, where the difference may be caused by the fact that transversal and longitudinal modes never truly decouple and that individual segments are not independent.
In order to further verify this hypothesis we inspect the dependency of such behaviour on the relaxation time of the longitudinal deformations of a single segment t l , see Fig. 17. As can be seen in Fig. 17b the position of the peak in the mean elastic energy hV el i, i.e. transition between regions II and III in Fig. 15 and 16, is independent of the relaxation time t l and corresponds to the increase of the mean bending energy, cf. Fig. 17c. This suggests that the position of the peak is determined solely by the dynamics of the transversal modes. The power law exponent in the cooperative regime also depends on the relaxation time t l as can be seen in the Fig. 17a, namely x(t l = 10 À2 ms) = 1.178 AE 0.014, x(t l = 10 À1 ms) = 1.114 AE 0.010, x(t l = 10 ms) = 0.8941 AE 0.0078, while the behaviour in the decoupled regime is mostly unchanged x(t l = 10 À2 ms) = 0.3494 AE 0.0029, x(t l = 10 À1 ms) = 0.3436 AE 0.0030, x(t l = 10 ms) = 0.3313 AE 0.0032. Notice also the large change in the  behaviour when t l = 1 ms, which is close to the typical time-scale where transversal modes unfreeze, Fig. 17c.
On the other hand, if the relaxation time t l gets bigger than the characteristic time of the transversal modes, the behaviour of the system changes completely, see Fig. 18. As the longitudinal modes do not dissipate the energy fast enough they invalidate the assumption that these modes are independent and hence rise to a value much higher than the one estimated by (42), see Fig. 18c. The times where the mean elastic energy hV el i saturates correspond to the moment of transition in the dynamics of the relaxation, as can bee seen in Fig. 18a.
However, the power law exponent x no longer agrees with the theoretical prediction for the WLC x ¼ 1 3 . Moreover the saturation of longitudinal modes also denotes the rapid increase in the mean bending energy, Fig. 18c.

Possible experimental verification
The independence of the observed behaviour on the actual value of the relaxation time of longitudinal deformations t l , which is the only free parameter of the simulation that is not related to any available experimental value, suggests the Fig. 18 The behaviour of the system when the relaxation time of a single segment length t l is larger than the characteristic time of the transversal modes deviates, even in long time asymptotic, from theoretical prediction for the WLC, see (a), where the power law exponent x strongly depends on t l . The difference in the behaviour is even more pronounced in the time evolution of the mean bending energy, (b), which now depends on the relaxation time t l in contrast to Fig. 17c. Moreover, the mean value of the elastic energy no longer converges to its equilibrium value as can be seen in (c), which is in contrast with the behaviour for small t l , see (d). possibility of an experimental verification of this quantity. In order to verify the existence of the new ''cooperative'' regime presented in this work in the case of F-actin, one needs to achieve a time resolution less then 1 ms along with a space resolution of less than 0.1 nm at the same time. Even though experimental methods such as fluorescence microscopy, optical tweezers 54,55 and fast scanning atomic force microscopy 56,57 are capable of working with such a resolution, it still might be challenging as the relative extension D is obscured by the thermal noise. Hence, it will require a large amount of consecutive measurements to obtain a large enough statistical ensemble. Moreover has to be careful not to add large masses or volumes at the ends of the polymer, as this might drastically alter its behaviour.

Conclusion
The aim of this work was to develop a coarse grained model for thin long stiff semi-flexible polymers, applicable to an environment with large external forces, e.g. the cytoskeleton. The model is based on deformable cylinders as basic building blocks and exploits the assumption of the high stiffness and thinness, i.e. the Young's modulus is high compared to the bending stiffness, hence the polymers rather bend than deform, and the diameter of the polymer is much smaller than its typical length and also than its persistence length. We have derived from first principles (1) the equations of motion for the cylinders' end points in the overdamped limit (14), which inherently preserve autocorrelations of the noise applied on the end points (16). This is one of the key differences between our model and bead-rod models, 12,13 where the beads' displacements are uncorrelated similar to bead-spring models, [1][2][3][4] or to the second kind of bead-rod models, 4,15 where the displacements of beads on opposite sides of a single segment are correlated only in the direction of the segment's axis (17). Our description of the overdamped diffusion differs from standard approaches based on the balance of forces and moments 17 by the presence of thermal ''torque'' (15), a correction naturally arising from the intrinsically non-homogeneous diffusion of a cylindrical segment. Let us stress here that the presented model, already at the level of independent cylinders, is well suited for studying the behaviour of active nematic particles, e.g. active swimmers at low Reynolds number. 32,33 The full model of the polymer is achieved by chaining independent cylindrical segments by constraints (21), which are included in terms of virtual forces, obtained as a solution of the eqn (24), and which induced virtual effective ''torques'', and by prescribing the bending potential (27) representing the semi-flexibility, which is included in the equations of motion in terms of an effective ''torque'' (28). The advantage of the bending model presented in this paper, the prevention of antialignment of neighbouring segments in case large forces are applied, thus increasing the overall stability of the simulation, is another difference with traditional approaches. Namely in the cytoskeleton, myosin motors apply forces of magnitude JFJ p pN thus exceeding thermal fluctuations. 11,12 Furthermore the same polymer model can be as well used as a basis for models of molecular motors like Myosin-II, because of their polymer-like nature. 38 These aspects of our model in particular will be fully utilised in our future work.
The full model is summarized in Section 4 from where it can be seen that the model has a similar computational complexity to the traditional bead-spring model introduced by Montesi et al. 15 The main benefit of our approach is not only in decreasing the number of elements used in simulation, but also in the possibility to increase the time step, mainly owing to the increased friction for larger elements, thus improving the overall performance of the model even more, see the discussion in Section 4.1. In Section 5.1 we have provided a specific comparison of the model presented in this paper with a traditional bead-spring model, where we have achieved a simulation speed-up of more than 3 orders of magnitude, for a single polymer in the extreme case of very long segments. Although this approach truly show its strength for long segments, we have demonstrated that using the presented approach can be beneficial even for segments, which are short relative to the persistence length l p . This suggests the applicability of the model to wide range of (bio-)polymeric systems.
The proposed model very well recovers the equilibrium behaviour of the WLC as is presented in Section 4. Both the polymer's thermal equilibrium end-to-end distance without, see 5.3, and with an applied tension, see Section 5.4, are evaluated and compared with theoretical predictions. 34,35 Testing of the dynamics is achieved by analysing the relaxation from the fully elongated state towards thermal equilibrium, see Section 6, where we have described a new ''cooperative'' regime for short times, which suggests that on short time scales the longitudinal and transversal modes are coupled. This constitutes our second main result. For longer times, when transversal and longitudinal modes decouple, we recover the theoretically predicted behaviour. 50,51 The emergency of ''cooperative'' behaviour on short-time scales (t o 10 ms in case of actin filament) might manifest itself as a faster relaxation towards equilibrium after the applied stress on the polymer network is released. Consequently, the response to the action of a molecular motor might deviate from that of the ideal WLC, which can contribute to the dynamics of filopodia and lamellipodia. Hence, the methods studying the actin-myosin networks based on the rigid rod models 4,10,15 might have to be re-evaluated. This will also be the subject of further studies.

A Time evolution of macroscopic degrees of freedom
Here we derive eqn (4) from (1) by using the ansatz (2). First notice that Because the cylinder is symmetric with respect to point inversion p j = p(i) 3 x i À X = X À x j , all material properties are also invariant with respect to the same point inversion, i.e. m p(i) = m i and g p(i) = g i . By inserting (2) into (1) we obtain The equation of motion for the centre of mass is obtain by summing (44) over all i where we have used the symmetry of the cylinder, i.e. all odd terms in x i À X have been cancelled. The equation governing rotations is obtained by summation of the cross product between x i À X and dv i , i.e. P i x i À X ð ÞÂdv i , where as for dv i we use (44). By using the symmetries of a cylinder we obtain X The last equation, for stretching is then obtained by summation (x i À X)dv i , If we insert x i À X = r i t and (3) and denote we can finish all the remaining summations in the continuous limit, m i -0 + , Although the noises present in this set of equations Because of the symmetries of the cylinder, , we can conclude that they are in fact independent. The last step in order to obtain eqn (4) is noticing, that xÁt = 0.
The remaining set of eqn (5) follow. The time evolution of the length is obtained from where x N and x 1 are coordinates of end points of the cylinder, x N À x 1 = lt, and where we have also used (43). Similarly the orientation t is governed by

B Equivalence of Smoluchowski and Langevin equations
Here we show how to obtain a Smoluchowski equation from a set of coupled Langevin equations in the case of an inhomogeneous environment by using the Itô calculus. 18,25,58,59 We restrict ourselves to the set of Langevin equations of the following shape dX AE (t,W (a) ,W (b) ) = A AE (X + ,X À ) dt + B(X + ,X À )ÁdW (a) where coordinates X + and X À are random variables with explicit depend on time t and three dimensional Wiener processes W (a) , W (b) , and coefficients A À , A + , B, C are now arbitrary integrable functions of the coordinates. We express the joint probability distribution m(x + ,x À ) as the mean value over all the possible realisations of Wiener processes starting from arbitrary initial condition at time t = 0 m t (x + ,x À ) = hd(x + À X + (t,W (a) ,W (b) ))d(x À À X À (t,W (a) ,W (b) ))i.
The Smoluchowski equation then corresponds to the expression for time derivative of the joint probability distribution. In our case the only time dependent quantities in the mean value are random variables for the positions X AE , which depends on time explicitly as well as via the realization of Wiener processes W. By using an Itô chain rule for the derivative we obtain where r AE corresponds to the gradient with respect to x AE and analogously the r AE 2 -terms represent a Hessian matrix. The other derivatives are then given by a corresponding Langevin eqn (45), i.e. q t X AE = A AE (X + ,X À ), = W(a) X AE = B(X + ,X À ), = W(b) X AE = AEC(X + ,X À ).
In order to evaluate these expressions we exploit the fact that we have a Dirac d-function present in the expression, e.g.

C Overdamped limit
Here we review the main results of the procedure for time-scale separation described by Pešek 24 and apply it to obtain a Smoluchowski equation from a Fokker-Planck equation.

C.1 General framework
The procedure is based on a couple of assumptions. The first one is that we can decompose the time evolution of the joint probability distribution to one part describing the evolution of the fast degrees of freedom and another for the slow degrees of freedom. Furthermore it assumes that the time evolution is dominated by the dynamics of the fast degrees of freedom and that the evolution of slow degrees of freedom can be considered as a perturbation, where 1 e L F Ã denotes the part of the backward Kolmogorov generator corresponding to the time evolution of the fast degrees of freedom, where we explicitly introduce the degree of time scale separation e, and L S * is the part corresponding to the slow degrees of freedom. The last assumption is that given the particular configuration of slow degrees of freedom x S the backward Kolmogorov generator for fast degrees of freedom 1 e L F Ã has a unique steady state distribution r(x F |x S ), to which every initial configuration of fast degrees of freedom converges where G F (x) represents a phase space element of fast degrees of freedom. This implies that by taking the limit of infinite time scale separation, i.e. e -0, we reach a steady state for fast degrees of freedom as well, Under these conditions it can be shown 24 that the time evolution of the joint probability distribution up to the first power in e is fully determined by the marginal distribution of the slow degrees of freedom  60 and l i /e are non-zero eigenvalues of À 1 e L F Ã . Hence the dynamics for slow degrees of freedom is autonomous and Markovian up to the first order in degree of time scale separation e and can be described by @ t n t x S ð Þ % n tþDt x S ð Þ À n t x S ð Þ Dt for Dt much larger than the characteristic time of the fast degrees of freedom t F ¼ where the correction is given by eDm t x S ; x F ð Þ¼Àe 1 L F Ã L S Ã n t r ½ þO e 2 Dt; Dt 2 À Á ; (49) from where it also follows that it does not contribute to the marginal probability distribution ð dG F ðyÞ eDm t x S ; y ð Þ¼0: In the case of the overdamped limit we expect that velocities, or in phase space formalism canonical momenta, relax towards equilibrium much faster than positions, hence they will represent the fast degrees of freedom while the positions corresponds to the slow degrees of freedom. Following the general procedure, the first step is obtaining the canonical momenta p AE associated with the positions of the cylinders' end points x AE . Then we will provide a time evolution equation for the joint probability distribution m(x + ,x À ,p + ,p À ;t) where we will identify the time evolution of the fast degrees of freedom. And in the end we derive an autonomic time evolution eqn (48) for the marginal probability distribution n(x + ,x À ;t), for which we later find the corresponding Langevin equations. In the process we also obtain an estimate for relaxation time of the fast degrees of freedom t F as well as the correction to the joint probability distribution (49).

C.2 Canonical momentum
Under the assumption that the external forces applied on the segment don't depend on velocities, no presence of magnetic fields is assumed, the only part of the Langrangian that does depend on velocities is the kinetic energy given by where we again used the symmetries of the cylinder, i.e. all terms odd in x i À X are inherently zero. The kinetic energy can also equivalently be expressed in terms of the velocities of the end points As the kinetic energy is the only part of the Lagrangian that depends on velocities, we have all the information necessary in order to derive canonical momenta Notice that for pure translation both moments are the same and correspond to the moment of a point particle with the same mass and velocity, however for pure rotations they have an opposite sign and the same magnitude, corresponding to the angular momentum of the segment divided by its length. 8 It is a linear operator which maps an argument to its inverse whenever invertible and to 0 whenever it belongs to the kernel of 1 e L F Ã .

D.1 Spatial resolution of diffusion
The first constraint comes from the fact that we want to restrict the possibility of two segments to pass through each other in a single iteration without interacting. The main limiting dimension in our case it the width of the polymer, which leads to a constraint for the maximal displacement of the end-points in one single time step Dx AE where the factor 2 comes from the fact that both elements can be moving. The displacement for all other points along a single segment is given by the linear combination of displacements of its end-points. The displacement in general depends on all forces applied to the single segment, however let us restrict ourselves now only to diffusion. The diffusional motion can be divided into translations and rotations. Translations: the translational diffusion is described by eqn (14) dx AE ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2k B Tg trans À1 q Á dW ðaÞ ; which is discretised as q Á W ðaÞ ðDtÞ: As we can see, the displacement is determined by the realisation of the Wiener process over a finite time interval, which is unfortunately unbound, hence the constraint (53) in its naive version can in principle not be applied. Instead we determine the confidence interval in which the displacement Dx AE will not be bigger than the diameter of the segment d in terms of multiples of the variance of the Wiener process c s , i.e. 1s confidence interval corresponds to c s = 1, 2s corresponds to c s = 2, etc. Then the constraint (53) is replaced by Rotations: similarly to the translational diffusion, the rotational diffusion is described by