Elsevier

Applied Numerical Mathematics

Volume 156, October 2020, Pages 422-445
Applied Numerical Mathematics

A stable finite element method for low inertia undulatory locomotion in three dimensions

https://doi.org/10.1016/j.apnum.2020.05.009Get rights and content

Abstract

We present and analyse a numerical method for understanding the low-inertia dynamics of an open, inextensible viscoelastic rod - a long and thin three dimensional object - representing the body of a long, thin microswimmer. Our model allows for both elastic and viscous, bending and twisting deformations and describes the evolution of the midline curve of the rod as well as an orthonormal frame which fully determines the rod's three dimensional geometry. The numerical method is based on using a combination of piecewise linear and piecewise constant finite element functions based on a novel rewriting of the model equations. We derive a stability estimate for the semi-discrete scheme and show that at the fully discrete level that we have good control over the length element and preserve the frame orthonormality conditions up to machine precision. Numerical experiments demonstrate both the good properties of the method as well as the applicability of the method for simulating undulatory locomotion in the low-inertia regime.

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 x:[0,1]×[0,T]R3 and an orthonormal frame e0,e1,e2:[0,1]×[0,T]R3 up to some final time T>0. We call u[0,1] the material parameter. Since we do not allow shear deformations, we set e0τ:=xu/|xu| the unit tangent to the midline. The other two coordinates e1 and e2 form an orthonormal basis of the cross section of the rod to the midline. See Fig. 1. We can form the skew system:1|xu|eu0=αe1+βe2,1|xu|eu1=αe0+γe2,1|xu|eu2=βe0γe1, where α and β are smooth fields denoting the curvature of the midline in directions e1 and e2 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 α0,β0,γ0. The viscous contribution is proportional to the time derivatives αt,βt,γt.

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.

We see these three properties as allowing our method to provide a computational tool both for the understanding of viscoelastic rods and also for domain experts working on undulatory locomotion in the low inertia regime. Further work is required to link this model for undulatory locomotion with a more detailed model of the surrounding fluid which would give a more accurate model and allow a greater variety of behaviours to be captured. Where such links are to be used a balance must be struck between the computational efficiency of a one dimensional approach against the greater detail provided by a fully three dimensional model.

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 R3 over a time interval [0,T] for some 0<T<. The rod can be described by a (non-arc length) parameterization of the centre line x:[0,1]×[0,T]R3 and an oriented frame of reference Q:[0,1]×[0,T]SO(3). 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 ()u and ()t, respectively. Our

Finite element spaces

We take a partition of [0,1] by N points u1=0<u2<<uN=1. We call {u1,,uN} the mesh. We will use a combination of piecewise linear and piecewise constant functions. We introduce the spaces:Vh:={vhC([0,1]):vh|[ui,ui+1] is affine, for i=1,,N1}Qh:={qhL2(0,1):qh|[ui,ui+1] is constant, for i=1,,N1}. We will denote by Vh,0 the space of finite element functions in vhVh such that vh(0)=vh(1)=0. 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)

  • C.W. Wolgemuth et al.

    Dynamic supercoiling bifurcations of growing elastic filaments

    Phys. D, Nonlinear Phenom.

    (April 2004)
  • S.S. Antman

    Dynamical problems for geometrically exact theories of nonlinearly viscoelastic rods

    J. Nonlinear Sci.

    (January 1996)
  • J.W. Barrett et al.

    The approximation of planar curve evolutions by stable fully implicit finite element schemes that equidistribute

    Numer. Methods Partial Differ. Equ.

    (2011)
  • J.W. Barrett et al.

    Parametric approximation of isotropic and anisotropic elastic flow for closed and open curves

    Numer. Math.

    (October 2011)
  • S. Bartels

    A simple scheme for the approximation of the elastic flow of inextensible curves

    IMA J. Numer. Anal.

    (January 2013)
  • O.A. Bauchau et al.

    The vectorial parameterization of rotation

    Nonlinear Dyn.

    (April 2003)
  • H.C. Berg

    Motile behavior of bacteria

    Phys. Today

    (January 2000)
  • M. Bergou et al.

    Discrete viscous threads

    ACM Trans. Graph.

    (July 2010)
  • M. Bergou et al.

    Discrete elastic rods

    ACM Trans. Graph. (TOG)

    (2008)
  • S. Berri et al.

    Forward locomotion of the nematode C. elegans is achieved through modulation of a single gait

    HFSP J.

    (June 2009)
  • F. Bertails et al.

    Super-helices for predicting the dynamics of natural hair

    ACM Trans. Graph.

    (July 2006)
  • A. Bilbao et al.

    Roll maneuvers are essential for active reorientation of Caenorhabditis elegans in 3d media

    Proc. Natl. Acad. Sci.

    (April 2018)
  • J.J. Blum et al.

    Biophysics of flagellar motility

    Q. Rev. Biophys.

    (May 1979)
  • D. Bray

    Cell movements: from molecules to motility

    Garland Sci.

    (2001)
  • C. Brennen et al.

    Fluid mechanics of propulsion by cilia and flagella

    Annu. Rev. Fluid Mech.

    (January 1977)
  • R. Casati et al.

    Super space clothoids

    ACM Trans. Graph.

    (July 2013)
  • S. Childress

    Mechanics of swimming and flying

    (1981)
  • N. Cohen et al.

    A new computational method for a model of Caenorhabditis elegans locomotion: Insights into elasticity and locomotion performance

  • N.A. Croll

    The behaviour of nematodes: their activity, senses and responses

    (1970)
  • G. Daviet et al.

    A hybrid iterative solver for robustly capturing coulomb friction in hair dynamics

    ACM Trans. Graph.

    (December 2011)
  • T.A. Davis

    Algorithm 832

    ACM Trans. Math. Softw.

    (June 2004)
  • K. Deckelnick et al.

    Error analysis for the elastic flow of parametrized curves

    Math. Comput.

    (October 2008)
  • K. Deckelnick et al.

    Computation of geometric partial differential equations and mean curvature flow

    Acta Numer.

    (April 2005)
  • J.E. Denham et al.

    Signatures of proprioceptive control in Caenorhabditis elegans locomotion

    Philos. Trans. R. Soc. Lond. B, Biol. Sci.

    (September 2018)
  • E.H. Dill

    Kirchhoff's theory of rods

    Arch. Hist. Exact Sci.

    (1992)
  • D. Durville

    Numerical simulation of entangled materials mechanical properties

    J. Mater. Sci.

    (November 2005)
  • G. Dziuk

    An algorithm for evolutionary surfaces

    Numer. Math.

    (December 1990)
  • G. Dziuk et al.

    Evolution of elastic curves in Rn: Existence and computation

    SIAM J. Math. Anal.

    (January 2002)
  • C. Elliott et al.

    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.

    (2016)
  • M. Elmi et al.

    Determining the biomechanics of touch sensation in C. elegans

    Sci. Rep.

    (September 2017)
  • Cited by (0)

    This work was supported by a Leverhulme Trust Early Career Fellowship.

    View full text