Elastic rods are a ubiquitous model of semi-flexible biopolymers. In biophysics, the worm-like chain ((WLC) model [1] underpins many theoretical and numerical studies of semiflexible polymers. The WLC model is a linearization of the classical Kirchoff rod model [2] which is itself a limiting case where the product of the local curvature and filament thickness is everywhere small. Discretizing the continuum Kirchoff model leads to equations of motion that are stiff and difficult to integrate numerically; artificial dissipation is often needed for long-term stability. On the other hand, in discrete dynamical systems it is known that symplectic integration methods give superior long-term stability in comparison with either high-order explicit or implicit integration methods [3]; the most common symplectic integrator is the Verlet algorithm. Our algorithm for elastic filaments is based on a discretization of the Hamiltonian line integral, including shear and extensional degrees of freedom. The model systematically approximates the dynamics of an elastic filament, without the limitation to small deflections inherent in the WLC model. Since the nodal forces and torques follow from an exact differentiation of a potential function, the equations of motion are guaranteed to be Hamiltonian, although the potential function itself is only an approximation to the continuum limit. This is in contrast to finite-element methods, where the continuum equations of motion are discretized in space; in this case the Hamiltonian structure is not preserved, even if the total energy is conserved. Our numerical scheme is entirely stable and energy conserving even for large deformations, yet the algorithm is entirely explicit, with the same level of complexity as molecular dynamics. A report of this work will be published in the Journal of Chemical Physics [4].
Here we illustrate the method by a simulation of a straight filament, which is intially bent into a circle and then released. As the rod evolves from its initial configuration, flexural waves propagate along the filament, leading to a surprising variety of configurations; a sampling of the filament shapes is illustrated in figure below. Initially the ends move slowly, and the filament assumes a teardrop shape (t = 300 t0), followed by a hairpin (t = 600 t0) as the ends of the filament accelerate. The time unit t0 = d/c_l, where cl is the longitudinal wave speed. The inverted U shape (t = 900 t0) straightens out (t = 1200 t0), and then develops a ``double-minimum'' shape (t = 1500 t0). The center of the filament moves down to complete the inversion and the filament approximately retraces the sequence of shapes in reverse order, to arrive at the inverted configuration at roughly half the period of the main oscillation. However, the motion is not exactly periodic because of the strong coupling between the flexural modes. The interaction of flexural waves can lead to large local stresses, exceeding that of the initial configuration; for example at the top of the teardrop (t = 300 t0) and at the bends in the hairpin (t = 600 t0). It has been shown that flexural modes can cause unexpected fractures by this mechanism [5]. |
Filament shapes at different times: t = 300 t0 (solid), t = 600 t0 (long dashes), t = 900 t0 (dashes), t = 1200 t0 (dot dash), and t = 1500 t0 (dotted). The time scale t0 = d/c_l is the time for a longitudinal wave to cross the diameter of the filament, where cl is the longitudinal wave speed. A movie of the filament dynamics can be found here (1.2M MPEG-1). |