A stable finite element method for low inertia undulatory locomotion in three dimensions☆
Introduction
The dynamics of active microswimmers has long captured the interest of physicists, mathematicians and engineers, not to mention researchers in various biological fields. Undulation, passing bending waves down the body, is a common strategy used across many orders of magnitude [40] ranging from bacteria [57], [61] and spermatozoa [37] to larger fish and whales [38]. It is especially common in the low inertia regime, where viscous forces of the surrounding media dominate over inertial forces and the scallop theorem prohibits self-propelled locomotion for time-reversal symmetric sequences of body postures [68], [79]. For additional reading, we refer the reader to a large body of excellent reviews on animal locomotion [19], [44], [58], fluid dynamics at low Reynolds number [16], [18], [37], [60], [61], [68], [83], and the biophysics and biology of cell motility [7], [8], [14], [15], [47], [48]. There are many computational studies in this area which have focussed on solving a fully detailed three dimensional formulation capturing many aspects of the problem including the dynamics and its interaction with the body (e.g. [34], [65], [66], [73], [74], [76]). These large scale computational studies provide many interesting results giving insights on both complex fluids and animal locomotion (see, e.g., the review [58] and references therein). An alternative is to choose a reduced, lower dimensional (in our case one dimensional) description of the body and represent all forces on this reduced description. This allows orders of magnitude faster simulations by considering only a one dimensional problem which significantly reduces the numbers of computational degrees of freedom.
In this work we develop, analyse and show experiments using a novel computational approach to the evolution of an open viscoelastic rod, representing the body of a long, thin microswimmer - a three dimensional object with one axis much longer than the other two. We assume the body is embedded in three dimensions and can undergo bending and twisting deformations. The resulting problem can be described as a system of one dimensional partial differential equations for a midline curve and an orthonormal frame which describes the conformation of cross sections to the midline. The model is a natural generalisation of the locomotion model [46] to three spatial dimensions. The elastic terms we consider arise from the classical Kirchhoff-Love model for an elastic rod [52] which are combined in a simple linear Voigt viscoelastic model [1], [53]. Our model avoids considering shearing or extensional deformations from the full Cosserat or Timoshenko models. In general this nonlinear system can not be solved analytically and requires the use of computer methods. The reduction to a one dimensional object allows for significant reduction in the complexity of the resulting mathematical model in particular with respect to the computational effort required to solve the problem.
Our new method tackles three key challenges. First when considering active locomotion one must be able to map directly between anatomical detail of the organism under consideration and the geometry defined by the model. For example, we should be able to clearly identify where the muscles are located and in which direction do they apply force. Our approach to capture this results in equations for the midline of the rod coupled to equations for the orientation of the cross section of the rod (see the discussion in [56] and Section 4.2 for more details). These equations must be carefully coupled to ensure that we have an accurate and robust representation of the geometry which allows us to apply the biological forces appropriately. Second we have a moving geometry which is a priori unknown. Our scheme captures this with a parameterization which is equivalent to a moving mesh. As is typical in this type of problem we must ensure that the moving mesh does not become too distorted (see, for example, [4], [33] and references therein). Finally, in many biological systems various parameters and problem data may be unknown or poorly characterized. Any computational method should be both cheap enough to run so that parameterization studies may be run and also robust to input parameters so that a wide variety of behaviours can be observed. See [49], [75] for example. We will demonstrate through analysis of our method and numerical experiments that we can address all three challenges.
Similar models to those considered in this paper arise in many areas of natural science and engineering. For example, similar models have been considered for elastic ribbons and filaments [9], [10], [39], tangled hair and fibres [12], [24], [30], [81], plants [42], [43], and woven cloth [50]. A historical overview of the model used in this paper is given in [29].
We are interested in the locomotion of low inertia microswimmers. We explore the applicability of our method through the case study of the 1 mm long nematode Caenorhabditis elegans. This animal has become a model organism studied by physicists, computer scientists and engineers partly due to its simple, undulating, periodic gait but also its experimental tractability and simple neuroanatomy [21]. In its natural environment, C. elegans grows mostly in rotten vegetation: a highly complex, heterogeneous three dimensional environment. It is cultured on the surface of agar gels for extensive use in laboratories. On agar, C. elegans move by propagating bending waves from head to tail to generate forward thrust [11], [23], [36], [45], [59], [67] where the bending arises from body wall muscles which line sides of the body. A two dimensional version of our viscoelastic model [46], developed originally for capturing worms or snakes moving across land, has previously been applied to C. elegans locomotion both on surfaces and in a range of both Newtonian and non-Newtonian media [20], [28], [36], [77], [78]. However, previous formulations have either linearized the equations, which means although postures are recovered trajectories are not, or viscosity is neglected, which limits the applicability of the model in less resistive environments. Recent experiments [13], [72] have shown that C. elegans achieves a different gait in three dimensions. Here we demonstrate that our computational tool is capable of accurately and robustly capturing such behaviours. The simulation results are meant to be indicative of the capabilities of the method rather than demonstrating any properties of the underlying model which is left to future work. The numerical method used in this section builds on the unpublished work [20] (see also [28]).
In our model, the rod is described through a smooth, time-dependent parameterization of the midline and an orthonormal frame up to some final time . We call the material parameter. Since we do not allow shear deformations, we set the unit tangent to the midline. The other two coordinates and form an orthonormal basis of the cross section of the rod to the midline. See Fig. 1. We can form the skew system: where α and β are smooth fields denoting the curvature of the midline in directions and and γ is a smooth field denoting the twist of the orthonormal frame about the midline. At the core of the model is a moment which is the sum of elastic and viscous contributions. The elastic terms are proportional to differences between the fields and some desired values . The viscous contribution is proportional to the time derivatives .
Inertial terms are ignored and external forces are simply modelled as linear drag terms [51]. This model can be seen as a quasi-static or low inertia approximation of the full rod dynamics. Although this very simple approach to modelling external forces for undulatory locomotion has been used previously in undulatory locomotion models (see the references above), a more accurate model would include non-local terms [51], [60] or capture the external fluid using dynamic equations (e.g. [22], [62], [64], [70]). Any model for the surrounding fluid must be coupled to a model of the body. The focus in this work is to accurately capture the body mechanics combined with a simple model of surrounding fluid. Study of a more detailed model of the environment coupled to the model of body mechanics presented here is left to future work.
Our model can be seen as a simplification of the model presented in [54], [55] or a three dimensional version of [46]. Full details of the model are given in Section 2.
Our approach is to discretize an appropriate formulation of the continuous equations directly. This allows us to use the structure of the equations to recover a robust numerical scheme through careful discretization choices. The discretization extends the approach of [63] to the case of non-constant twist and open curves.
Key to our numerical method is a well suited formulation of the continuous partial differential equation system. We start from the balance laws for linear and angular momentum and the linear viscoelastic constitutive law. This is combined with geometric constraints so that the solution variables are the position of the midline, line tension (a Lagrange multiplier for enforcing the length constraint), the curvature of the midline, the twist and angular velocity of the frame and two auxiliary variables which describe the bending and twisting moments. The continuous system of equations is discretized in space by a mixed finite element method where we use a mix of piecewise linear and piecewise constant approximation spaces. We use a Lagrange multiplier to enforce inextensibility (i.e. that the length of the curve is locally fixed) in a similar approach to [5]. Finally, we discretize in time using a semi-implicit method which results in a linear system of equations to solve at each time step, inspired by the approach in [31], [32], [63] (see also [4]), together with an update formula for the frame. The use of lower order finite element approximations, along with mass lumping [80], allows us to derive identities for various geometric quantities at the mesh vertices. The frame is updated using a Rodrigues formula where the rotation is specified by the frame's angular moment which can be derived from the solution variables. We use the angular momentum and Rodrigues formula as an alternative to using Euler angles [71] or quaternions [54], [55]. Although direct use of the Rodrigues formula is not recommended for numerical applications in general [6], we will show that our formulation avoids problems for any angles. More details of the numerical scheme are given in Section 3.
We will demonstrate three key properties of our scheme:
- •
a semi-discrete stability result (coupled with computational evidence of fully discrete stability) which shows that we recover a discrete Lyapunov functional for our scheme;
- •
control over the length element which ensures that vertices in the moving mesh do not collide;
- •
preservation of the frame orthogonality conditions to ensure that the frame really does remain orthonormal over time.
In our previous unpublished work [20], we presented a similar scheme (also used in [28]). This paper builds on that work extending the scheme to three dimensions including twisting as well as bending contributions and generalising the material law to include viscous terms.
There are several very successful methods from discrete differential geometry which solve either the problem we consider or a generalisation. In this approach the rod is first discretized and then equations are derived to evolve those discrete quantities. We mention in particular the approach of [10] (extended to viscous threads in [2], [9]) which uses a discrete set of vertices (in our notations a piecewise linear curve) with a material frame on each edge between the vertices (piecewise constant in our notation) to represent a rod which was stored using a reduced curve-angle formulation [56]. The approach uses curvature and twist defined as integrated quantities based at vertices to define a Kirchhoff-Love energy and then applies discrete parallel transport and variation of holonomy to derive the update equations. This work was generalized to Cosserat rods by [41] who revert away from the reduced curvature-angle formulation storing the full frame. The super-helix and super-clothoid approaches [12], [17] use a piecewise constant or piecewise linear approximation of generalized curvature of a rod and then recover the geometry of the rod (position of midline and frame) using analytic expressions. The scheme results in a smooth curve with well defined curvatures and twist.
In the context of these schemes our model can be seen as using a set of vertices with a frame at each vertex to define the discrete geometry of the rod. The equations are straight discretization of the continuous scheme using a finite element approach and discrete equations to define curvature and relate twist with tangential angular velocity. Our choice to discretize continuum equations allows us to use standard finite element tools to both analyse and implement the method in a single unified framework. This approach is well suited to the parabolic nature of the equations we consider. It is our particular choice of both geometric discretization and direct discretization of the equations that allows us to demonstrate the properties of our scheme.
In Section 2, we present the continuous model we use and the discretization is shown in Section 3. Numerical experiments to demonstrate the efficacy of the method are shown in Section 4. The restriction of our scheme to a two dimensional problem is given in Appendix A.
Section snippets
Geometry
We consider a smooth, inextensible, unshearable rod embedded in over a time interval for some . The rod can be described by a (non-arc length) parameterization of the centre line and an oriented frame of reference . For a discussion of rod representations see [56]. We will write derivatives with respect to the first coordinate, the material coordinate u, and the second coordinate, time t, with subscripts and , respectively. Our
Finite element spaces
We take a partition of by N points . We call the mesh. We will use a combination of piecewise linear and piecewise constant functions. We introduce the spaces: We will denote by the space of finite element functions in such that . We will use the subscript h, defined to be the maximum mesh spacing, to denote discrete quantities.
Results
We provide three test cases for our numerical scheme. In the first we relax a straight rod to one with prescribed curvatures and twist and in the other two we demonstrate the applicability of the method to C. elegans locomotion.
Discussion
We have presented a new finite element scheme for viscoelastic rods suitable for simulating undulatory locomotion in three dimensions. We have shown analytically that the semi-discrete problem preserves the energy structure of the continuous problem and that the scheme preserves geometric constraints exactly (up to machine precision in practical computations). Further our numerical experiments show, first, that the analytic properties are realized in a practical fully discrete scheme and that
Declaration of Competing Interest
None
Acknowledgements
The author would like to thank Netta Cohen and Felix Salfelder for discussions which have improved the writing of the manuscript and the computer implementation of the method. I would also like to thank the anonymous reviewers for their careful consideration and helpful comments.
References (83)
- et al.
A discrete geometric approach for simulating the dynamics of thin viscous threads
J. Comput. Phys.
(November 2013) - et al.
Modelling biological and bio-inspired swimming at microscopic scales: Recent results and perspectives
Comput. Fluids
(January 2019) - et al.
Nematode locomotion: dissecting the neuronal–environmental loop
Curr. Opin. Neurobiol.
(April 2014) Regularized stokeslet segments
J. Comput. Phys.
(December 2018)- et al.
Numerical aspects in the dynamic simulation of geometrically exact rods
Appl. Numer. Math.
(October 2012) - et al.
Locomotion control of Caenorhabditis elegans through confinement
Biophys. J.
(June 2012) - et al.
Application of smoothed particle hydrodynamics to modeling mechanisms of biological tissue
Adv. Eng. Softw.
(August 2016) - et al.
Direct measurements of drag forces in C. elegans crawling locomotion
Biophys. J.
(October 2014) - et al.
A fully three-dimensional model of the interaction of driven elastic filaments in a stokes flow with applications to sperm motility
J. Biomech.
(June 2015) - et al.
Material Properties of Caenorhabditis elegans Swimming at Low Reynolds Number
Biophys. J.
(2010)
Dynamic supercoiling bifurcations of growing elastic filaments
Phys. D, Nonlinear Phenom.
Dynamical problems for geometrically exact theories of nonlinearly viscoelastic rods
J. Nonlinear Sci.
The approximation of planar curve evolutions by stable fully implicit finite element schemes that equidistribute
Numer. Methods Partial Differ. Equ.
Parametric approximation of isotropic and anisotropic elastic flow for closed and open curves
Numer. Math.
A simple scheme for the approximation of the elastic flow of inextensible curves
IMA J. Numer. Anal.
The vectorial parameterization of rotation
Nonlinear Dyn.
Motile behavior of bacteria
Phys. Today
Discrete viscous threads
ACM Trans. Graph.
Discrete elastic rods
ACM Trans. Graph. (TOG)
Forward locomotion of the nematode C. elegans is achieved through modulation of a single gait
HFSP J.
Super-helices for predicting the dynamics of natural hair
ACM Trans. Graph.
Roll maneuvers are essential for active reorientation of Caenorhabditis elegans in 3d media
Proc. Natl. Acad. Sci.
Biophysics of flagellar motility
Q. Rev. Biophys.
Cell movements: from molecules to motility
Garland Sci.
Fluid mechanics of propulsion by cilia and flagella
Annu. Rev. Fluid Mech.
Super space clothoids
ACM Trans. Graph.
Mechanics of swimming and flying
A new computational method for a model of Caenorhabditis elegans locomotion: Insights into elasticity and locomotion performance
The behaviour of nematodes: their activity, senses and responses
A hybrid iterative solver for robustly capturing coulomb friction in hair dynamics
ACM Trans. Graph.
Algorithm 832
ACM Trans. Math. Softw.
Error analysis for the elastic flow of parametrized curves
Math. Comput.
Computation of geometric partial differential equations and mean curvature flow
Acta Numer.
Signatures of proprioceptive control in Caenorhabditis elegans locomotion
Philos. Trans. R. Soc. Lond. B, Biol. Sci.
Kirchhoff's theory of rods
Arch. Hist. Exact Sci.
Numerical simulation of entangled materials mechanical properties
J. Mater. Sci.
An algorithm for evolutionary surfaces
Numer. Math.
Evolution of elastic curves in : Existence and computation
SIAM J. Math. Anal.
On algorithms with good mesh properties for problems with moving boundaries based on the harmonic map heat flow and the DeTurck trick
SMAI J. Comput. Math.
Determining the biomechanics of touch sensation in C. elegans
Sci. Rep.
Cited by (0)
- ☆
This work was supported by a Leverhulme Trust Early Career Fellowship.