Conformations, hydrodynamic interactions, and instabilities of sedimenting semiflexible filaments

The conformations and dynamics of semiflexible filaments subject to a homogeneous external (gravitational) field, e.g., in a centrifuge, are studied numerically and analytically. The competition between hydrodynamic drag and bending elasticity generates new shapes and dynamical features. We show that the shape of a semiflexible filament undergoes instabilities as the external field increases. We identify two transitions that correspond to the excitation of higher bending modes. In particular, for strong fields the filament stabilizes in a non-planar shape, resulting in a sideways drift or in helical trajectories. For two interacting filaments, we find the same transitions, with the important consequence that the new non-planar shapes have an effective hydrodynamic repulsion, in contrast to the planar shapes which attract themselves even when their osculating planes are rotated with respect to each other. For the case of planar filaments, we show analytically and numerically that the relative velocity is not necessarily due to a different drag of the individual filaments, but to the hydrodynamic interactions induced by their shape asymmetry.


Introduction
Semiflexible filaments are fundamental constituents of microbiological systems, where microtubules and actin filaments serve as scaffolds for cellular structures and as routes to sustain and guide cellular transport systems. 1 Microtubules are also the main structural elements of cilia and sperm flagella, where their relative displacement and deformation due to motor proteins gives rise to the flagellar beat and hydrodynamic propulsion. 2,3 Microtubules and flagella can be seen as elastic filaments interacting with their own flow field. The ability to visualize, assemble, and manipulate biological and artificial semiflexible polymers [4][5][6][7] poses new fundamental questions on the dynamics of filaments when elastic and hydrodynamic forces compete.
The dragging of stiff rods through a viscous fluid has been studied in detail. 8 A single rod does not reorient, but falls with its initial orientation. A more complex dynamical behavior can be expected and is indeed observed for semiflexible filaments when the curvature or stretching elasticity competes with the hydrodynamic interactions. [9][10][11] Single dragged semiflexible filaments bend into a shallow V-shape to balance the higher drag at both ends 12 and their end-to-end vector aligns perpendicularly to the external field. 11 For strong drag, higher modes have been found to be excited; this generates W-shapes initially, which then relax back into horseshoe-like U-shape. 12 Here, the dynamics seems to be constrained to the plane initially defined by the direction of the external field and the filament itself. However, these investigations address the problem from a deterministic point of view, and little attention has been paid to the dynamic stability of the resulting shapes. In all cases, the dragged and deformed semiflexible filament initially defines the settling plane, but the stability of the filament's planar shape has not been investigated as function of the external field or the relative position of possible neighboring filaments.
Here, we focus on the full three-dimensional shape of one, two, and three semiflexible filaments sedimenting in a homogeneous external field. We incorporate the hydrodynamics into the equations of motion for the filament shape via the Oseen tensor, valid in the limit of zero Reynolds number. As a result of our numerical and analytical analysis, we find that the deformations confined to a plane become unstable with respect to normal perturbations at a threshold value B 1 * of the strength B of the external field, which is smaller than the threshold B 2 * where initial, transient W-shapes become excited, see Fig. 1. Thus, with increasing strength of the external field, two instabilities and transitions to new sedimentation modes are predicted. The first transition is from a stable planar U-shape with little bending to a stationary horseshoe-like U-shape with out-of-plane bending. The second transition at stronger fields excites a metastable W-shape, also with out-of-plane bending, which then ''relaxes'' into a non-stationary asymmetric U-shape. As result, there exist two families of shapes, where the elastic forces are balanced by a conformation-dependent drag.
We consider next the interaction between two filaments in an external field. Indeed, while the dynamics of an isolated filament is an indispensable knowledge needed to understand the case of n 4 1 interacting filaments, many situations are characterized by elastic slender objects interacting via the generated flow field: cilia, 7,13 sperm, 14,15 and E. coli bundles 16,17 are probably the most relevant from a biological point of view.
It is known that the sedimentation behavior of colloids can be quite complex. The interaction of sedimenting particles has been studied in considerable detail for spherical colloids. 18,19 Two particles sediment together, but don't follow the direction of the external field, and move instead under an angle with respect to it. For more particles, many different dynamical behaviors can be found, in particular periodic motions where particles ''dance'' around each other. 19 For dragged semiflexible filaments, the dynamical behavior is even more complex. 10 In particular, we show that two filaments ( Fig. 1) attract each other, repel each other, or spin around the field depending on the intensity of the external field.
We focus here on the stability of the sedimentation plane for different field intensities and on the origin of the relative velocity. In particular, we want to see whether the velocity difference is due to different shapes or to the broken up-down symmetry. For even more filaments, the dynamics becomes unsteady at much weaker external field strength than expected from the two-filaments case.

Discrete model
In the simulation, filaments of length L = b(N À 1) are represented by N mass points of mass m connected by harmonic bonds of length b, with the potential where R i = r i+1 À r i is the vector connecting the consecutive points, i A {0,. . .,N}, and k b is the force constant. To account for filament stiffness, we introduce the bending potential with the bending rigidity k. 20 In addition, the mass points are exposed to an external gravitational field with the force on a particle with the unit vector e z along the z axis of the Cartesian coordinate system (cf. Fig. 1).
Inter-and intrafilament hydrodynamic interactions substantially influence the filament dynamics. Hence, we apply the equation of motion g 0 ( : for the overdamped dynamics, 21 where, u i is the background flow velocity at the position r i of the particle, g 0 = 3pZb is its friction coefficient for a fluid of viscosity Z, and F C i is the sum of all conservative forces. 10 We compute the background flow velocity explicitly via the Oseen tensor. 21 Thus, the final equations of motion are with the hydrodynamic tensor and r ij = r i À r j , r ij = |r ij |, and r ij = r ij /r ij . If needed, excluded-volume interactions are implemented via a short-range Yukawa potential between points of different filaments, which implies a minimal effective distance during the simulations. A low-amplitude white noise is added to avoid metastable states. The noise is not considered to be of thermal origin as it is chosen to be negligible compared to the other hydrodynamic and mechanical forces and barely influences the stationary settling dynamics for the considered external fields.

Parameters/methods
We set the hydrodynamic diameter of a point equal to the bond length, as in the Shish-Kabab model of ref. 10, 12 and 21, thereby fixing the aspect ratio to b/L. Lengths are measured in units of the hydrodynamic diameter b and time in units of g 0 b 3 /k. This choice eliminates the friction coefficient and the bending rigidity from the equations of motion. The force constant for the bonds is set to k b b 3 /k = 1, resulting in a maximum extension of AE0.6% of the length over the investigated range of parameters. In these units, the external field strength mg becomes G = mgb 2 /k. For convenience and an easier comparison with results of ref. 10, we characterize the external force by B = N 2 G, or B = mgL 2 /k. For B { 1, the bending rigidity dominates and the filament is essentially straight. We consider only filaments of length L = 30b in the following. For excluded volume interactions, the minimal effective distance is approximately 5b. The equations of motion are integrated with an adaptive time-stepping Velocity-Verlet algorithm. 22,23

Continuum model
For an analytical description of the filament dynamics, we adopt a continuum model. The equation of motion of the point r n (s,t) (ÀL/2 r s r L/2) along the contour of filament n is given by 24 where f m is the external force density and the index m indicates the various filaments. As before, the hydrodynamic tensor H(r n (s) À r m (s 0 )) comprises the Oseen tensor and the local friction. Explicitly, it reads as Here, Y(x) is the Heaviside function, g = 3pZ is the friction per unit length, and R = r n (s) À r m (s 0 ). 24,25 The force density f comprises bond, bending, and gravitational forces. In the limit of a rather stiff filament, it can be written as with the persistence length l p . 26,27 In the following, we will neglect the bond term, i.e., the term with the second derivative and focus on bending stiffness only. The expansion in terms of the eigenfunctions f n of the biharmonic operator, i.e., with suitable boundary conditions, 26,27 leads to the equations of motion of the mode amplitudes The matrix representation of the hydrodynamic tensor is The eigenfunctions f n (s) and relaxation times t n are well known. 5,24,27 For convenience, we summarize them in Appendix A. However, eqn (12) is nonlinear and thus cannot straightforwardly be solved. For the current analysis, the second and forth mode are most important; they are responsible for the V-and W-shape displayed in Fig. 1.
To characterize the numerically obtained filament conformations, we calculate the mode amplitudes v n n ðtÞ ¼ ð L=2 ÀL=2 ds r n ðs; tÞf n ðsÞ: The components of the vector v n n indicate the importance of the mode in the Cartesian directions. For example, the mode amplitude w n 2,z measures how much the filament bends along the z direction into a V-like shape.
As for the discrete model, we scale lengths by the filament diameter b and time by g 0 b 3 /k B Tl p . The latter is identical to the time scaling of the discrete model, since k = k B Tl p . This implies for the strength of the external force G = rgb 3 /k B Tl p = mgb 2 /k.

Deformation and dynamics of single filament
The filament is initially oriented along the x axis of the reference frame. After a certain time, the dragged filament reaches a stationary shape and velocity. Examples of conformational sequences for various field strengths are displayed in Fig. 1. We characterize the shapes via eqn (14) in terms of the mode amplitudes. In Fig. 2, the most important stationary amplitudes are presented. Below a critical field B 1 * C 1200, the filament shape is governed by planar modes (green and black lines), where w 2,z dominates and, thus, the characteristic V-shape appears.
In simulations restricted to a two-dimensional plane, or in three-dimensional simulations without noise, 12 the filament dynamics is localized in the xz plane and filaments bend into a planar W-shape for fields B 4 B 2 * E 1800. In contrast, in our three-dimensional simulations with weak noise, we find that the planer filament conformations are metastable for B 1 * o B o B 2 *, and also modes along the y axis are excited. We characterize the out-of-plane filament shape and dynamics by the mode amplitude w 2,> (t), where In the stationary state, an U-shaped and deck-chair-like conformation is assumed with out-of-plane bending (see Fig. 1). The filament orientation is fixed and w 2,> = w 2,y (blue line in Fig. 2). Since its shape is asymmetric, the filament drifts sideways while settling in the external field. When B \ B 2 *, the mode w 4,z becomes important at early times, leading to a temporary W-shape (Fig. 1). The trajectory for B C 3000, displayed in Fig. 1, shows the initial W, which later turns into an asymmetric U-shape, in which one arm is longer than the other. The appearing shape is stable; however, because of its asymmetry, the mode amplitude w 1,z is non-zero and the filament rotates around the z axis with frequency o, see Fig. 2 (orange line), which we determined via eqn (15). In Fig. 3, we characterize the helical trajectories by the pitch, radius, and rotation frequency (B 4 B 2 *). As the field increases the rotation frequency increases, and the helix becomes more tight because the radius decreases and the pitch shortens. The ratio between the pitch and the radius defines the helix angle a E 4p/9, constant for all B 4 B 2 *. Approaching the transition point from above, the rotation frequency vanishes, and both radius and pitch diverge because the trajectory straightens.
In contrast, in the deterministic dynamics of previous studies, 12 the W-shape was found to decay only into the stable and symmetric planar horseshoe shape.

Conformations and dynamics of two interacting filaments
3.2.1 Weak field -relative velocity. As shown in Section 3.1, the stationary shape of a single filament in weak fields B o B 1 * is of V-shape, which breaks the bottom-top symmetry. This is sufficient to generate an effective attraction between sedimenting filaments with the same shape. To characterize this interaction, we compute the relative velocity Dv between the centers of mass of two filaments of equal shape along the sedimentation direction. The filaments remain localized in the xz plane and are separated by a distance d. As shown in Fig. 4, the relative velocities exhibit a significant dependence on the filament separation. We especially find that Dv B d À2 for Fig. 2 Stationary mode amplitudes of a single semiflexible filament as function of the external field B. The shaded areas indicate the 66% confidence interval. When B o B 1 *, only planar modes are excited, and the filament stays in the plane defined by its initial orientation and the orientation of the applied field, here the xz plane. For B 4 B 1 *, an out-ofplane mode w 2,> is excited. For B 4 B 2 *, the out-of-plane component w 2,> , the bending component w 2,z saturates, and the amplitude w 4,z becomes important (visualized in Fig. 1). In black crosses indicate the maximum value of w 4,z before it decays. The resulting shape is asymmetric and spirals around the z axis with frequency |o|/o 2 (orange line), with o 2 the frequency of the second mode (eqn (15)).   (17), save for a common factor d E 11. The theory describes correctly the trend on d, and the trend on w 2,z holds up to w 2,z = 8 Â 10 À3 . distances larger than the filament length. The distance dependence can be understood by the theoretical model introduced in Section 2.3. For the considered filament shapes, r 1 (s,t) = (w 1,x (t)f 1 (s), 0, w 2,z f 2 (s)) T (16) and r 2 (s,t) = r 1 (s,t) + de z , the general expression for the velocity difference derived in Appendix B yields in the limit d -N. Evidently, the filaments attract each other due to the top-bottom asymmetry of their shapes. In the simulations, the filament shapes are determined initially by imposing the amplitude w 2,z , which is then kept fixed. The simulation results of Fig. 4 are in agreement with our theoretical prediction down to roughly the filament length. The d À2 power law is indeed a universal scaling, unaffected by the filaments shape and external field as evident from the theoretical considerations in Appendix B. The dependence of Dv on w 2,z (eqn (17)) is also verified for very small bending.
3.2.2 Weak field -stability. We now relax the imposed shape constraint and consider collective effects for two filaments, which are initially straight, oriented along the x axis, and displaced along the z axis by a distance d (cf. Fig. 5(a)). For easier comparison with ref. 10 and 12, we employ the dimensionless number D = (A upper À A lower )/(L/2) to quantify the bending asymmetry, where A upper,lower is the total z extension of the upper/lower filament. As indicated in Fig. 6, the filament curvature changes with time and the upper filament is bent stronger than the lower one. Fig. 6(a) and (b), show the curvature asymmetries D and the relative velocities for various external field strengths. D decreases with increasing distance d, indicating more similar shapes at larger distances. Hydrodynamic interactions lead to an attraction of the two filaments (v upper 4 v lower ), in agreement with the imposed-shape approximation studies of the last section. Indeed, the constant-shape approximation still gives the correct (L/d) 2 power-law dependence for d/L c 1, while the magnitude of the deformation, w (eff) 2,z , has to be fitted. When the filaments approach each other, the generated flow field depends on the details of their shapes that, in turn, depends on the external field, hence we expect a non-universal behavior. Note that in contrast to ref. 10, we find that the upper filament bends more than the lower filament (see also Fig. 5).
The planar configuration of a filament is also stable with respect to filament rotations around the field axis, see Fig. 5(b). Filaments that are initially displaced along the z axis (as in the previous case) and rotated with relative orientation angle y around the external field axis spin until the relative angle vanishes, as illustrated by Fig. 6(c). Also in this case, the upper filament drifts and rotates faster than the lower one, see Fig. 6(d). The relative velocity is essentially the same as in the planar case.
Thus, two filaments sedimenting in weak fields relax toward a stable planar configuration one behind the other. The shape of the filaments is dominated by the second mode, pointing downwards, as shown in Fig. 5. This mode dominates and it breaks the mirror symmetry of the hydrodynamic interactions even for filaments of the same shape. Note that, in contrast to the single filament case, the system does not reach a stationary state velocity or shape, since the upper filament is always faster than the lower filament until the filaments touch each other.

Strong field.
For strong fields, we consider two filaments, which are initially displaced by 6L along the field direction. We measure the shape eigenvalues when the distance is 5L, in the quasi-stationary regime, and find that the eigenvalues exhibit the same behavior as those of a single filament. This means that for B 4 B 1 * the dynamics of each filament is dominated by the local flow field and not by the interactions with the other filament. Indeed, we find no correlations between the orientations of the out-of-plane components of the two filaments for B 4 B 1 *: the two filaments can by chance bend out-of-plane and drift in arbitrary directions.
When B 4 B 2 *, the filaments undergo the same transitions as a single filament: each of them reaches the same stationary shape and rotation velocity as an isolated filament. We find no correlations between the rotation directions of the two filaments: some filaments spin in the same direction, others in The two filaments approach each other after initialization in a rotated configuration. Both filaments spin around the z axis, with the upper filament spinning faster (see Fig. 6(c and d)).
opposite directions, with no preference. This highlights the relevance of hydrodynamic interactions between two filaments for external fields weaker than B 1 *. At stronger fields, the non-planar configurations generate forces that compete with and dominate over hydrodynamic interactions among the filaments and the emergent behavior is, essentially, the same as that of an isolated filament.

Three filaments
Given the complex dynamics of two interacting filaments, it is interesting to consider also the collective behavior of several filaments. We find in simulations of systems with more than two filaments an intriguing collective dynamic behavior even for very weak fields (B C 60) and in the absence of noise.
We focus here on the case of three filaments, see Fig. 7. For most (randomly chosen) initial configurations, the nearest two filaments form a bundle that settles faster than the third filament that is then left behind. However, we find also some initial configurations where all three filaments attract each other and form a bundle. In this case, the relative positions are not stationary; instead, the filaments follow a periodic trajectory, see Fig. 7 (inset). In the inset of Fig. 7, we show also that the shapes of the three filaments are not stationary. The mode amplitude w 2,z of each filament changes periodically, with a constant phase shift between them.
Our results for one and two filaments indicate that triggering of a time-periodic bifurcation requires strong fields. However, the three-filaments results suggest that systems with more filaments display a very complex dynamics even for weak fields due to complex hydrodynamic interactions.

Discussion and conclusions
We have investigated the dynamics and stability of semiflexible filaments exposed to an external homogeneous field and interacting only via hydrodynamic fluid fields. Due to the competition between hydrodynamic interactions and bending stiffness, the appearing dynamical behavior is richer than for entropydominated polymers or interacting rods.
We have shown that, for weak fields B o B 1 *, co-planar configurations of two filaments are stable upon perturbations that rotate the shapes relative to each other around the field  axis. With simulations of fixed shape filaments, we have highlighted that a V-or U-shape is sufficient to break the hydrodynamic symmetry at low Reynolds numbers, leading to a relative velocity that scales with distance as (L/d) 2 . Hence, the difference in drag coefficients between filaments is not necessary to explain the faster settling velocity of the upper filament.
For external field strengths exceeding the critical value B 1 *, the hydrodynamic interactions bend the filament out of its principal plane. Simulations of a single bent filament show that the hydrodynamic forces balance the elastic force, stabilizing the out-ofplane shape. The resulting trajectory shows a drift in the direction of out-of-plane bending, superimposed to the settling motion. This is a novel result, not to be confused with the previously reported metastable W-state 12 that is excited when B 4 B 2 *. A careful analysis of the eigenmodes indicates that the decay of the metastable state does not, in general, lead to the reported planar horseshoe shape, but also excites an average rotation mode with respect to the field axis (w 1,z ) and our out-of-plane bending mode w 2,> . The filaments spin then around the field axis.
Finally, we have demonstrated that three filaments display an unexpected periodic dynamics even at field strengths far weaker than B 1 *. This is in contrast to the dynamics of a pair of filaments that either displays a monotonic dynamics that relaxes the attractive force (when B is weak) or a dynamics dominated by the single-filament (when B 4 B 1 *).
The interesting external fields B are in the range 10 1 t B t 10 4 . We can estimate these parameters for biopolymers like actin or microtubules. Actin has a persistence length of l p C 17 mm, an average length L C 20 mm, and the bending rigidity k C 60 Â 10 À3 pN mm 2 . 1 The external gravitational field, corrected for buoyancy, is about G E 10 À7 pN mm À1 , which implies B gravity C 10 À2 . Microtubules, on the other hand, are stiffer, longer, and heavier with l p B 1 mm, L B 100 mm { l p and G E 10 À6 pN mm À1 . 28 This yields the effective field strength B gravity C 10 À1 . An experimental test of our predictions is therefore within reach of modern centrifuges with accelerations of about 10 3 g.

A Eigenfunctions of a semiflexible filament
The eigenvalue eqn (11) with the boundary conditions yields the eigenfunctions þ sin z n s sin z n L=2 ; n 4 1; odd; (19) þ cos z n s cos z n L=2 ; n 4 1; even: (20) The wave numbers are approximately given by z n = (2nÀ 1)p/2L (n 4 1), and the corresponding relaxation times by More precise eigenfunctions are provided in ref. 24 and 27. The set of functions is complemented by the eigenfunction of the center-of-mass translation 27 and that of rotation of the rodlike object with the relaxation time

B Relative velocity of two filaments
We derive here an equation for the relative velocity between the centers of mass of two filaments. We restrict our analysis to the case of small bending amplitudes, that is equivalent to consider small external fields, and filaments of identical shape.
Since Ð L=2 ÀL=2 f n ðsÞds ¼ ffiffiffi ffi L p d n;0 for the exact eigenfunctions, the difference in the center-of-mass velocity Dv cm = v cm 1 À v cm 2 of two isolated filaments is given by ds@ t r 1 ðs; tÞ À r 2 ðs; tÞ Â Ã Substitution of eqn (12)  The first two terms on the right-hand side account for selfinteractions of the individual filaments, the other two terms for the hydrodynamic interactions between the filaments. We simplify our considerations by assuming identical shapes of the filaments, i.e., we set v n 1 = v n 2 : = v n . Moreover, for the constant external force the relation applies f n nG = f n 0G d 0n independent of the particular filament. Hence, its contribution vanishes, which yields We are primarily interested in the distance dependence of the relative center-of-mass velocity. Hence, we additionally neglect the dyadic term in the hydrodynamic tensor (8). Moreover, the local friction term vanishes in eqn (26), and the hydrodynamic tensor can be written as ð L=2 ÀL=2 f n ðsÞf 0 ðs 0 Þ jr n ðsÞ À r m ðs 0 Þj dsds 0 : This journal is © The Royal Society of Chemistry 2015 Using the eigenfunction expansion eqn (10), we obtain r n ðsÞ À r m ðs 0 Þ ¼ Dr nm cm þ X 1 n¼1 w n f n ðsÞ À f n ðs 0 Þ ð Þ ¼ Dr nm cm þ Nðs; s 0 Þ: (28) With this definition, we obtain for DH 12 0n = H 21 0n À H 12 0n DH 12 0n ¼ 1 8pZ Thus, the relative velocity decreases quadratically with the distance between the filaments. There is evidently no velocity difference when Dr 21 cm is perpendicular to N(s,s 0 ). In particular, there is no force between two specifically aligned rods as long as their director w 1 is perpendicular to Dr 21 cm .