Numerical exploration on buckling instability for directional control in flagellar propulsion

Weicheng Huang and M. K. Jawed *
Department of Mechanical and Aerospace Engineering, University of California, Los Angeles, Los Angeles, California 90095, USA. E-mail:

Received 12th September 2019 , Accepted 9th December 2019

First published on 10th December 2019

We report a numerical method to control the swimming direction by exploiting buckling instability in uniflagellar bacteria and bio-inspired soft robots. Our model system is comprised of a spherical rigid head and a helical elastic flagellum. The rotation of the flagellum in low Reynolds environment generates a propulsive force that allows the system to swim in fluid. The locomotion is an intricate interplay between the elasticity of the flagellum, the hydrodynamic loading, and the flow generated by the moving head. We use the Discrete Elastic Rods algorithm to capture the geometrically nonlinear deformation in the flagellum, Lighthills Slender Body Theory to simulate the hydrodynamics, and Higdons model for the spherical head in motion within viscous fluid. This flagellated system follows a straight path if the angular velocity of the flagellum is below a critical threshold. Buckling ensues in the flagellum beyond this threshold angular velocity and the system takes a nonlinear trajectory. We consider the angular velocity as the control parameter and solve the inverse problem of computing the angular velocity, that varies with time, given a desired nonlinear trajectory. Our results indicate that bacteria can exploit buckling in flagellum to precisely control their swimming direction.

1 Introduction

A primary mode of locomotion in bacteria is the rotation of a single helical flagellum in a fluid medium.1–3 A flagellum – typically several microns in length and a few nanometers in cross-sectional radius1 – is effectively an elastic rod, and the fluid flow – at the length and time scales of bacterial locomotion – is in the Stokes regime (Reynolds number ≪ 1).4 A non-trivial coupling between the geometrically nonlinear deformation in the flagellum and the hydrodynamics of the low Reynolds flow leads to a number of hallmarks of flagellar propulsion, e.g. tumbling,5 turning,6 bundle formation,7–9 and polymorphic transformations.10,11 Previous investigations involving both experiments6,12 and simulations13,14 show that the hydrodynamic force can lead to buckling in the flagellum or in the soft connector between the flagellum and the cell body.15,16 In this article, we assume that the connector and the flagellum are made of the same material; this assumption simplifies the system and makes the design of bio-inspired robot easier to fabricate. We seek to understand how buckling of soft flagellum can be exploited to precisely steer the swimming direction. We are also motivated by revolutionary advances in flagellated micro-robots17 and the potential application of the steering by buckling mechanism in the control of such robots. Bacteria-inspired flagellated robots are typically controlled by external magnetic field,18–22 electric field,23 acoustic excitation waves,24,25 and chemically powered propulsion.26 As a result, their manufacturing process often involves advanced specialized materials, e.g. ionic polymer–metal composites.27 Buckling, on the other hand, can be induced in any elastic material simply by changing a single scalar variable – the rotational speed of the flagellum.13

The propulsion of a uniflagellar system (bacteria or robot) can be divided into three components: (1) elasticity of the flagellum, (2) hydrodynamic flow by the flagellum, and (3) hydrodynamics of the head. Prior works13 showed that Kirchhoff elastic rod28 – as a model for the flagellum – and Lighthills slender body theory (LSBT)29 – as the hydrodynamic force model – can accurately describe the dynamics of a flexible flagellum in a viscous medium, including the buckling instability.6,14,30,31 The role of a head, the flow generated by its motion, and its coupling with the flow generated by the flagellum, were not explored in these studies. Thawani et al.32 incorporated the effect of the head on the locomotion of a model uniflagellar bacterium; however, the flagellum was assumed to be rigid. We seek to bridge this gap and incorporate all the three aforementioned components to demonstrate that the buckling instability can be used by the uniflagellar system to follow a prescribed 3D trajectory (Fig. 1).

image file: c9sm01843c-f1.tif
Fig. 1 Buckling for direction change. (a) At angular velocity ω = 10 rpm, bacteria soft robot can move along a straight line. (b) By increasing angular velocity to ω = 60 rpm and last for certain time, flagella robot can change its orientation.

In this numerical study, we consider a model system comprised of a rigid spherical head and a rotating flexible flagellum in a viscous medium at low Reynolds number. The flagellum is treated as an elastic rod and modeled using the Discrete Elastic Rods (DER)33,34 algorithm – a numerical method developed in the computer graphics community for visually dramatic simulation of hair, fur, and other filamentary structures. The hydrodynamic force applied by the surrounding medium on the flagellum is modeled using LSBT that can capture the long range hydrodynamic interaction between flows generated by distant parts of the structure. The hydrodynamics of the spherical head is included in the framework following Higdons model.35 These three ingredients – DER, LSBT, and Higdons model – in concert form a numerical framework that seamlessly captures the geometrically nonlinear deformation in the flagellum. Exploiting the computational efficiency of the framework, DER in particular, we systematically explore parameter space to quantify the change in swimming direction resulting from buckling instability. We demonstrate that the swimming direction can be controlled simply by changing a single scalar input – the angular velocity of the flagellum. While the numerical tool is general enough to include a hook – a soft joint between the head and the flagellum,2 we do not include it in our numerical exploration and, instead, focus on the most simple system of a head and a single flagellum that can follow any desired trajectory. This study indicates that uniflagellar bacteria can apply buckling in flagella to turn. It also provides a blueprint for self-contained soft robots, particularly suited for the micro-meter scale due to the assumption of low Reynolds number.

This paper is organized as follows. In Section 2, the numerical framework is described. This tool is used in Section 3 to illustrate how the uniflagellar system can follow a prescribed trajectory. Potential challenges in developing bio-inspired robots using this numerical investigation are discussed in Section 4. Concluding remarks and avenues for future research are presented in Section 5.

2 Numerical method of flagella locomotion

The simulation of fluid–structure interaction of flagellar robot combines three components: (i) DER that simulates the flexible flagellum as a linear elastic rod,33,34 (ii) LSBT that models the force exerted by the viscous flow on the flagellum,29 and (iii) Higdons model35 for flow generated by the spherical head. This section is organized as follows. The geometry of the system is described in Section 2.1. DER model and LSBT are summarized in Sections 2.2 and 2.3, respectively. Next, in Section 2.4, the spherical rigid head is incorporated into the framework. Finally, the complete framework is employed to quantify the onset of buckling in the flagellum in Section 2.5.

2.1 Geometry of flagella soft robot

The system has two parts shown schematically in Fig. 2(a) – a rigid spherical head and a soft helical flagellum. The physical parameters are: helix pitch λ = 32.6 mm, helix radius R = 6.04 mm, axial rod length L = 13 cm, (resulting in a contour length l = 20 cm), rod radius r0 = 1 mm (and, therefore, second moment of inertia I = πr04/4, polar moment of inertia J = πr04/2, and cross sectional area A = πr02), Young's modulus E = 10 MPa, shear modulus G = E/3 (assuming incompressible material), and material density ρr = 1.273 g cm−3. The radius of head is b = 1.2 cm. The viscosity of the fluid is μ = 2.7 Pa s, and its density ρm is the same as the rod density ρr. Overall, the system is neutrally buoyant. The Reynolds number in the current study is ρrωRr0/μ ≤ 4 × 10−2, i.e., always in the Stokes limit.36 These physical parameters were chosen partly based on prior experimental works.13,36
image file: c9sm01843c-f2.tif
Fig. 2 (a) Geometry of flagella-inspired soft robot. (b) Schematic of the relevant quantities. Bending curvature is related to the misalignment between two edges and twist is associated with rotating angle between two frames.

2.2 Discrete elastic rods

We use the DER method33,34 to model the geometrically nonlinear dynamics of the soft filament on uniflagellar robot. In the discrete setting of DER shown schematically in Fig. 2(a), the structure of this bio-inspired soft robot is discretized into N nodes: x0, x1,…,xN−1, which correspond to N − 1 edge vectors: e0,…,eN−2 such that ei = xi+1xi and i = 0,…,N − 2. We use subscripts for node-based quantities and superscripts for edge-based quantities. The head is located at x0 and the edge connecting the head and helical filament is e0. The length of this edge is equal to the radius of the head. In the body-fixed XbYbZb Cartesian frame in Fig. 2(a), the Xb axis is along −e0. The second edge from x1 to x2 is e1 = −R[x with combining circumflex]b + b, where R is the radius of the helix and [x with combining circumflex]b (and ŷb) is the unit vector along the Xb (and Yb) axis. All the remaining nodes (x2,…,xN−1) in their undeformed state at time t = 0 fall on a helical shape and the length of each edge is 2δ, where δ is the natural cutoff length described in Section 2.3. Each edge, ei, has an orthonormal adapted reference frame {di1,di2,ti} and a material frame {mi1,mi2,ti}; both of the frames share the tangent ti = ei/‖ei‖ as one of the directors. The reference frame is updated at each time step through parallel transport in time, and, referring to Fig. 2(b), the material frame can be obtained from the time-parallel reference frame using a scalar twist angle θi. See ref. 33, 34 and 37 for a detailed exposition of the DER algorithm. The nodal coordinates, xi, and the twist angles, θi, constitute the 4N − 1 sized degrees of freedom (DOF) vector, q = [x0,θ0,x1,θ1,…,xN−2,θN−2,xN−1]T, that completely describes the configuration of the robot. Hereafter, superscript T denotes transpose operator. Based on this kinematic representation, in the remainder of this section, we discuss the formulation of elastic strains, energies, and the time stepping procedure of the DER algorithm.

The strains of a deformed Kirchhoffs rod are comprised of three parts: stretching, bending, and twisting. Stretching strain associated with the i-th edge, ei, is

image file: c9sm01843c-t1.tif(1)
where double bar (‖·‖) denotes the magnitude of a vector. Hereafter, quantities with an overbar indicate evaluation in the undeformed state, e.g.ēi‖ is the undeformed length of the i-th edge. Bending strain is captured by the curvature binormal which measures the misalignment between two consecutive edges at a node xi,
image file: c9sm01843c-t2.tif(2)
and its norm is ‖(κb)i‖= 2tan(ϕi/2). The material curvatures are given by the inner products between the curvature binormal and material frame vectors,
image file: c9sm01843c-t3.tif(3a)
image file: c9sm01843c-t4.tif(3b)
The twisting strain at the i-th node, in the discrete setting of DER, is measured using the discrete twist
τi = θiθi−1 + [m with combining low line]i,(4)
where [m with combining low line]i is the reference twist associated with the reference frame.33

The total stretching, bending and twisting energies of the helical filament are related to the quadratic form of strains,

image file: c9sm01843c-t5.tif(5a)
image file: c9sm01843c-t6.tif(5b)
image file: c9sm01843c-t7.tif(5c)
where Δli = (‖ēi‖ + ‖ēi−1‖)/2 is the Voronoi length associated with the i-th node.

In the time stepping scheme of the DER method, implicit Euler integration is used to update the DOF vector q and its velocity (time derivative of DOF) v = [q with combining dot above] from time tk to tk+1 (h = tk+1tk is the time step size) by integrating the following 4N − 1 equations of motion:

image file: c9sm01843c-t8.tif(6)
where subscript k + 1 (and k) denotes the evaluation of the quantity at time tk+1 (and tk), image file: c9sm01843c-t9.tif is the internal forces (associated with nodal positions) and internal moments (associated with the twist angles) vector, Fextk+1 is the external force vector (e.g. hydrodynamic force discussed in the subsequent sections), and image file: c9sm01843c-t10.tif is the diagonal mass matrix comprised of the lumped masses. The Jacobian associated with eqn (6) used for Newtons iteration is
image file: c9sm01843c-t11.tif(7)
where i and j are integers between 0 to 4N − 2. Since the gradient of the external hydrodynamic force cannot be analytically evaluated, this term is neglected in eqn (7). As such, the hydrodynamic drag is treated explicitly in this time marching scheme. Importantly, the Jacobian [Doublestruck J] is a banded matrix and the time complexity of this algorithm is O(N).34 This computationally efficient algorithm can be treated as a data generator for the use of trajectory design discussed in Section 3.

2.3 Lighthill slender body theory

We use LSBT to model the viscous drag experienced by a slender filament in motion within a low Reynolds environment and implement this as the external force Fext in eqn (6). Previous studies have already successfully coupled DER algorithm and LSBT to show the buckling instability of helical rods rotating in viscous fluid and compared the numerical results against experiments.13,36 In this section, we briefly review the LSBT-DER implementation, without considering the effects of moving head, which will be detailed in Section 2.4.

The velocity at each node on the rod is equal to the fluid velocity at that point (no-slip boundary condition). The velocity at the j-th node on the flagellum is uj ≡ [[q with combining dot above]4j,[q with combining dot above]4j+1,[q with combining dot above]4j+2]T; recall that [q with combining dot above] is the 4N − 1 sized velocity vector of the DOFs. This nodal velocity uj can be decomposed into two components: (i) a flow velocity (uf)j that is generated by the force exerted by the flagellum onto the fluid (equal and opposite to the hydrodynamic force on the flagellum), and (ii) another flow velocity (uh)j that is induced by the motion of the head. For the first component, we use LSBT to relate the velocity (uf)j and the hydrodynamic force on the flagellum.

Due to the linearity of the Stokes equations for low Reynolds number flow,38,39 LSBT assumes a series of Stokeslets and dipoles along the centerline of rod, and builds a relationship between the velocity field contributed by the flagellum with the hydrodynamic force applied on it. In the discrete setting, that relation is13

image file: c9sm01843c-t12.tif(8)
where fi ≡ [Fext4i,Fext4i+1,Fext4i+2]T is the external force at the i-th node, (fj) = ([Doublestruck I]tjtj)fj is the projection of fj in the plane perpendicular to the tangent tj at the j-th node, [Doublestruck I] is the identity tensor, ⊗ denotes the tensor product, rij is the position vector from the i-th to the j-th node, [r with combining circumflex]ij is the unit vector along rij, and image file: c9sm01843c-t13.tif is the natural cutoff length (r0 is the radius of the circular cross-section of the rod and e is the Napier's constant). The node-based tangent, tj, is computed from the average of the preceding and following edge vectors such that image file: c9sm01843c-t14.tif, where tj (and tj−1) is the unit vector along the j-th (and j − 1-th) edge. While image file: c9sm01843c-t15.tif has been commonly used in prior works,13,32 recent studies suggest that this value can be tuned to achieve the most accurate solution.40 However, in the model setup used in this paper, the buckling threshold does not strongly depend on the value of δ. This discrete formulation of LSBT requires that the length of each edge be 2δ, resulting in a total node number of N = 122 in our model system.

Using eqn (8), a linear system can be formed to describe the relationship between the velocity and the hydrodynamic forces:

Uf = [Doublestruck A]F,(9)
where the velocity vector, Uf ≡ [(uf)1,(uf)2,…(uf)N−1], and the hydrodynamic force vector, F ≡ [f1,f2,…,fN−1], both have size 3(N − 1). The matrix [Doublestruck A], with size 3(N − 1) × 3(N − 1), depends on the geometry of the flagellum and can be evaluated from eqn (8). Of the total N nodes, the first one is the head and will be treated in the next section.

Knowing the velocity field, Uf, generated by the flagellum and the matrix [Doublestruck A], the viscous drag force, F, along the centerline can be computed using eqn (9). This force is then used as the external force in the equations of motion (eqn (6)). To avoid numerical issues while solving this inverse problem in eqn (9), we assume that the force varies smoothly along the arc-length of the rod. Details of this algorithm can be found in a previous work.36 Also note that the (4N − 1)-sized external force vector Fext in eqn (6) is simply a rearranged form of the 3(N − 1)-sized force vector F in eqn (9); the former can be expressed as Fext ≡ [f0,0,f1,0,f2,…,0,fN−1]T. The force on the head, f0, is described next.

2.4 Interaction between head and flagellum

LSBT provides a formulation for the hydrodynamic flow generated by the flagellum and the viscous drag force. However, the flow generated by the head and the hydrodynamic forces on it have not been considered. Here, we discuss the interaction between the rigid spherical head and the soft flagellum in viscous fluid environment.

We first formulate the flow generated by the moving spherical head with radius b. In Fig. 3, we consider a rigid head with a translation velocity, Uh ≡ [[q with combining dot above]0,[q with combining dot above]1,[q with combining dot above]2], and angular velocity, Ωh (along [x with combining circumflex]b). The viscous flow at the j-th node from the spherical head is given by Higdon's model,32,35

image file: c9sm01843c-t16.tif(10)
where (rh)j is the position vector of the j-th node relative to the center of the head, (rh)j = ||(rh)j||, and × notation is the cross product. As a reminder, when combining eqn (8) and (10), the actual velocity at the j-th node is uj = (uf)j + (uh)j (no-slip boundary condition).

image file: c9sm01843c-f3.tif
Fig. 3 Interaction between head and filament and filament itself in viscous fluid.

Next we model the viscous drag experienced by the moving head. The hydrodynamic force applied on the head is also comprised of two parts: (i) drag force −6πμbUh (and torque −8πμb3Ωh) caused by its own motion, and (ii) the force from the Stokeslets on the flagellum. Because of the long range hydrodynamic interaction in viscous environment, the Stokeslets on the flagellum also contribute force and torque to the spherical head,32,35

image file: c9sm01843c-t17.tif(11)
where (rh)i is the position vector of the i-th node relative to the center of the head and its norm is (rh)i = ||(rh)i||.

Note that, without considering the method of image,35 the presence of flow from the Stokeslets on the flagellum does not satisfy the no-slip boundary on the surface of the spherical head. Method of image assumes that there is an imaginary Stokeslet located inside the head, such that the sum of the flow exerted by the real Stokeslet and the fake Stokeslet can ensure the no-slip boundary on the spherical surface.35 However, the error can be small when compared with experimental data if the fake Stokeslets are ignored.32

At this point, we have all the pieces to build a framework to capture the coupling of DER, LSBT, and moving head. At each time step tk, we know the DOF vector qk, its velocity [q with combining dot above]k, and the angular velocity of the head (Ωh)k. We first evaluate the flow generated by the head based on eqn (10). Then the flow from the flagellum Uf can be obtained by subtracting the flow due to the head from the velocity of each node. Next, LSBT (eqn (9)) is used to compute the external force F applied on the filament. The external force on the head can be obtained from

f0 = −6πμbUh + fh,(12)
where fh is given by eqn (11). With the external force on both the spherical head and the elastic flagellum, we can compute Fext ≡ [f0,0,f1,0,f2,…,0,fN]T and use Newtons method to solve the equations of motion in eqn (6). This gives us the updated DOF vector qk+1 and the velocity [q with combining dot above]k+1 = (qk+1qk)/h. Finally, the angular velocity of the head (Ωh)k+1 is updated from the torque balance of the entire structure,
image file: c9sm01843c-t18.tif(13)
These quantities (qk+1, [q with combining dot above]k+1, (Ωh)k+1) evaluated at time tk+1 are used in the next iteration to move from t = tk+1 to t = tk+2. The time step size used in the current work is h = 1 ms following a convergence study.

2.5 Buckling instability of soft filament

We use this numerical framework to systematically investigate the deformation and locomotion of the uniflagellar bacteria and robots in viscous fluid. We apply a Dirichlet boundary condition to specify the twisting angle on the first edge, θ0, to perform the rotation of the soft filament at a prescribed angular velocity ω – the control parameter in this study. All other nodes and edges are free and evolve based on the balance between the elastic and external forces. The motor speed is the angular velocity of the flagellum relative to the angular velocity of the head, Ωh. The rotation of the flagellum generates a propulsive force that moves the entire structure – the head and the flagellum – forward. At sufficiently small angular velocity in Fig. 4(a), the motion of the head, i.e. the position of the first node, traces a linear path, and the helical structural remains stable. However, the soft filament undergoes a buckling instability when the angular velocity ω exceeds a critical value and the swimming trajectory becomes highly nonlinear. A representative configuration of the flagellum upon buckling is shown in Fig. 4(b). A time-scale – μL4/EI – is obtained from the balance between bending force and viscous force,36 which can be used to normalize the angular velocity and time. Hereafter, overbar represents non-dimensional quantities, e.g. [small omega, Greek, macron] = ωμL4/(EI) and [t with combining macron] = tEI/(μL4) are normalized angular velocity and normalized time, respectively.
image file: c9sm01843c-f4.tif
Fig. 4 Configurations at normalized time t = 100 s ([t with combining macron] = 1.02): (a) stable configuration (with [small omega, Greek, macron] = 102.82); and (b) buckling configuration (with [small omega, Greek, macron] = 616.92). Relation between normalized angular velocity and (c) normalized length, L′/L at [t with combining macron] = 1.02; (d) translational displacement, ‖s‖, over [t with combining macron] = 1.02; (e) normalized external applied torque, [T with combining macron] = TL/EI; and (f) normalized external power [P with combining macron] = [T with combining macron][small omega, Greek, macron].

To quantify the buckling angular velocity, [small omega, Greek, macron]b, for application in trajectory design in the next section (Section 3), we perform a parameter sweep along angular velocity to systematically study the mechanical response of the uniflagellar system in viscous fluid. In Fig. 4(c), we plot the normalized length, L′/L, as a function of the normalized angular velocity, [small omega, Greek, macron], at different values of the head radius, where L′ is the Euclidean distance between the first node and the last node, and L is the end-to-end distance in the undeformed configuration. In this parameter sweep, a flagellated system in rest position starts moving at a prescribed angular velocity and keeps that angular velocity constant for 100 s (1.02 in normalized time) – a time sufficiently longer than the transient dynamics. The normalized length of the flagellum is measured at [t with combining macron] = 1.02. Regardless of the size of the head, we observe that, at small values of [small omega, Greek, macron], the soft filament slightly stretches so that L′/L is above 1; this behavior is reminiscent of the observations in optical tweezers experiments with DNA molecule.41 There exists a critical value of angular velocity, [small omega, Greek, macron]b, beyond which the helical structure undergoes buckling instability and deforms into a shape substantially different from its natural helical shape. The buckling direction is deterministic and the rotational symmetry is broken by the connection between the flagellum and the head. This buckling velocity, [small omega, Greek, macron]b, is strongly influenced by the radius of the head, b: compared with the case where the effect of head is ignored (b = 0.0), different values of the head size, b = {0.4, 0.8, 1.2} cm, reduces ωb by 19%, 32%, and 40%, respectively. A physical explanation behind the trend of decreasing ωb with increasing head size is provided next. If the angular velocity of the flagellum, ω, remains fixed, a larger head requires more propulsive force to overcome the drag and, therefore, the overall swimming velocity of the flagellated system decreases. As a result, the relative fluid velocity between the medium and the system along the axis of the flagellum (−Xb in the body-fixed frame) decreases. Since a fluid flow along −Xb direction resists buckling,42 a larger head brings about the buckling instability sooner, i.e. at a lower value of angular velocity.

We turn to another measure of the motion to explore the effect of instability in Fig. 4(d) and use the same data to plot the traversed displacement ‖s‖ – the Euclidean distance between the position of the head at [t with combining macron] = 0 and [t with combining macron] = 1.02 – as a function of angular velocity, ω. The distance ‖s‖ increases approximately linearly with angular velocity when [small omega, Greek, macron] < [small omega, Greek, macron]b, and the flagellated system covers a straight path in this phase. This relation between ‖s‖ and [small omega, Greek, macron] would be exactly linear in case of a rigid flagellum. Due to the flexibility of the flagellum, it deforms and slightly deviates from the linear displacement – angular velocity relationship. However, the distance covered, and as such the speed of the whole structure, depends on the radius of the head, b. Compared with the case without a head (b = 0), the traversed distance drops by 16%, 32%, and 43% at b = {0.4, 0.8, 1.2} cm, respectively, when swimming at a constant angular velocity in this regime ([small omega, Greek, macron] < [small omega, Greek, macron]b). On the other hand, the translational displacement ‖s‖ shows a nonlinear behavior when the angular velocity is beyond the threshold [small omega, Greek, macron]b. This is due to the structural instability in the helical filament and the resulting nonlinear 3D trajectory.

Considering the potential application of this study in design of soft robots, we turn to the torque and power required by the flagellated system. Fig. 4(e) shows the normalized external torque, [T with combining macron] = TL/EI, applied by the head to maintain the rotation of the flagellum as a function of normalized angular velocity. We observe that the torque increases almost linearly with the angular velocity. Moreover, there is negligible variation in torque with the size of the head, even though the swimming speed varies with the head radius. In Fig. 4(f), we plot the normalized power, [P with combining macron] = [T with combining macron][small omega, Greek, macron], as a function of [small omega, Greek, macron]. As expected from the data on torque vs. angular velocity, the required power increases quadratically with angular velocity.

3 Trajectory design

Due to the existence of a critical angular velocity for buckling, [small omega, Greek, macron]b, the flagellated system follows a straight line at [small omega, Greek, macron] < [small omega, Greek, macron]b and goes into a complex 3D trajectory at [small omega, Greek, macron] > [small omega, Greek, macron]b. Based on this observation, we present a control policy for directional control of uniflagellar robots and bacteria. This far, the numerical framework computes the trajectory of the system given the physical parameters and the angular velocity, ω(t), as a function of time. We now turn to the inverse problem of computing the angular velocity, ω(t), – a single scalar input – given a prescribed 3D trajectory. The prescribed trajectory can be divided into a finite number of discrete line segments and, therefore, the path following process can be reduced to a simpler problem comprised of three parts: (i) initially the flagellated system goes along a straight line with an initial direction, d0, by rotating its flagellum at a low angular velocity ωl (where ωl < ωb), and the traversed distance depends linearly on the running time of this lower angular velocity; (ii) when there is a misalignment between its initial direction d0 and the target direction t1 (see Fig. 5(a)), its angular velocity increases to ωh (where ωh > ωb) and stays at that value for tc seconds such that the helical filament undergoes buckling instability and the orientation changes; (iii) finally the angular velocity switches from ωh to ωl, the helical structure reshapes to the unbuckled state, and the flagellated system swims along a straight line with a new orientation, t1. By repeating these three phases, the uniflagellar system can cover a given 3D trajectory. The lower angular velocity for straight path used in this study – chosen rather arbitrarily – is ωl = 10 rpm ([small omega, Greek, macron]l = 102.82) and the higher angular velocity for directional change is ωh = 60 rpm ([small omega, Greek, macron]h = 616.92). The head size considered in this section is b = 1.2 cm and, from Fig. 4, the critical normalized angular velocity for buckling is ωb ≈ 52 rpm ([small omega, Greek, macron]b ≈ 534.66). The outstanding problem is to compute the running time, tc, at high angular velocity, ωh, given the desired change in swimming direction.
image file: c9sm01843c-f5.tif
Fig. 5 (a) Definition of turning angle, θ, and out of plane angle, ϕ, in 3D trajectory. (b) Illustration of different turning angle, θ, on xy plane (ϕ = 0). (c) Illustration of different out of plane angle, ϕ, with a fixed turning angle, θ = 90°. (d) Turning angle, θ, and out of plane angle, ϕ, as functions of running time, tc, of higher angular velocity, ωh. (e) Schematic of body-fixed frame of rotating helix.

In Fig. 5(a), we explain the turning process: the flagellated system is moving along an initial direction, d0, at time t = 0 and it has to change its swimming direction to t1. Without any loss of generality, we define a world frame xyz based on the body-fixed frame XbYbZb at t = 0, such that the x, y, and z-axes are along −Xb, −Yb, and Zb-axes. This world frame remains fixed with time and the body-fixed frame moves and rotates with the flagellum. We can define the turn from d0 to t1 using two parameters: turning angle, θ, and out of plane angle, ϕ. Turning angle, θ, is the angle between the initial direction, d0, and the target direction, t1. Out of plane angle, ϕ, is the angle between the target, t1, and its projection, txy1, on the xy plane. Fig. 5(b) shows the initial and final directions, corresponding to different values of θ, at a fixed out of plane angle ϕ = 0°. Next, to illustrate the out of plane angle, ϕ, Fig. 5(c) shows different out of plane angles at a fixed turning angle, θ = 90°.

We first investigate how tc influences the turning process, quantified by θ and ϕ. In Fig. 5(d), we plot the turning angle, θ, and out of plane angle, ϕ, between d0 and t1, as functions of tc. Each value of tc corresponds to a pair of θ and ϕ angles. We observe that the turning angle, θ, is nonlinearly related to the running time, tc. Within a certain range (0 ≲ [t with combining macron]c ≲ 0.55), the longer the buckled phase at ωh, the greater the turning angle, θ. Depending on the duration of the higher angular velocity, tc, the system can make a turn from 0° to a maximum of 154°. For a turn requiring an even greater turning angle, we can divide it into two turns so that each turn necessitates θ < 154°.

For a given initial direction, d0, and final direction, t1, we can compute a turning angle, θ*, and an out of plane angle, ϕ**. The turning angle, θ*, corresponds to a running time, tc* from Fig. 5(d). This combination of θ* and tc* then corresponds to a specific value of out of plane angle, ϕ*. However, the values of ϕ* and ϕ** are not necessarily equal and this may imply that the problem is overdetermined. Referring to Fig. 5(e), we can resolve this issue from the simple observation that the desired out of plane angle, ϕ**, can be varied by up to 2π with a single rotation of the flagellum and the distance covered during a single rotation at the lower angular velocity, ωl, is only 0.37 cm. This is negligible compared with the length of the flagellum L = 13 cm. A methodological explanation of the turning process, as implemented in our simulation, is provided next.

When the uniflagellar system, with the flagellum rotating at a low angular velocity of ωl, needs to turn, it first computes the turning angle, θ*, and the running time, tc*. This (θ*,tc*) pair corresponds to an admissible value of the out of plane angle, ϕ*, based on Fig. 5(d). The flagellated system will keep rotating at ωl, i.e. hovering, until the required out of plane angle is equal to the admissible value of ϕ*. The flagellum needs to make less than one turn during hovering and moves less than 0.37 cm. Following this step, the angular velocity of the flagellum will be increased to ωh for tc seconds and the swimming direction will change due to buckling. Subsequently, the angular velocity is switched back to ωl to return to straight line motion. The sequence of steps discussed above has to be repeated for each turn in the trajectory. With this path planning strategy, the uniflagellar system can follow an arbitrary 3D trajectory with only a single control parameter – angular velocity ω – by utilizing the buckling instability of soft filament.

An example implementation is presented in Fig. 6 and Movie S1 (ESI)43 where the system covers a square trajectory on xy plane (z = 0; xyz is the space-fixed frame) with each side of length ΔL = 0.25 m, approximately equal to only two body lengths. Initially (t = 0), the flagellated system is located at z = 0 along the x-axis. Fig. 6(a) presents snapshots of the deformed structure at five points along the path. In Fig. 6(b), the target trajectory (dashed line) and the achieved trajectory (solid line) are shown on the xy plane. As desired, the trajectory remains almost entirely in the xy plane, with relatively small displacement along the z-axis, ΔzL < 5%, during the buckled phase. The corresponding control signal – the angular velocity as a function of time – is presented in Fig. 6(c).

image file: c9sm01843c-f6.tif
Fig. 6 (a) A planar square trajectory obtained by buckling instability. Also see Movie S1 (ESI).43 (b) The comparison between the prescribed trajectory (dashed line) and real path (solid line). (c) The angular velocity ω as a function of time for the achieving of the desired trajectory.

4 Discussion

A natural next step is to apply this numerical investigation to design and fabricate bio-inspired swimming robots. An important consideration in the design is the buckling threshold of the flagellum. The framework to follow a prescribed trajectory relies on instigating the instability. This threshold value, ωb, scales approximately linearly with EI/(μL4) (under the assumption that the role of the cross-sectional radius on the generated flow can be ignored) and decreases with the size of the head. The geometry of the helix, represented by the non-dimensional parameters λ/L and R/λ, also significantly influences this threshold.13 The data presented in this paper as well as prior works13 can be used to get a good estimate of ωb. The main challenge is likely going to be the imperfection sensitivity of this instability. Manufacturing errors associated with fabrication and measurement errors associated with physical parameter estimation render unlikely a perfect match between experiments and simulations. To overcome these uncertainties, a recommended approach is to construct a plot similar to the one presented in Fig. 5(d) using the experimental robot and use that data to formulate the control strategy.

5 Conclusion

We developed a computational framework to simulate the geometrically nonlinear deformation of soft filament and bio-locomotion of uniflagellar soft robot rotating in low Reynolds fluid environment and solved the inverse problem of following a desired trajectory simply using one control parameter – angular velocity ω. The computational framework fully accounts for three components: elasticity of the flagellum, long-range hydrodynamic forces on the flagellum, and the flow due to the head. We quantified the effect of the size of the head on the instability of the flagellum and the motion of the entire system. Our finding indicates a strong role of the head size on flagellar propulsion. We then studied the application of this instability to control the swimming direction of uniflagellar bacteria and robots. Supported by the robustness and efficiency of the numerical tools, we performed systematic parameter sweeps to quantify how the turning angles change with the time series of the angular velocity. The results are used to solve the problem of making a prescribed turn by varying the angular velocity with time. A sequence of such turns can then be employed to follow any prescribed trajectory in three dimensions. The results were presented in non-dimensional form to emphasize the scale independence of the phenomenon; turning by buckling may be used in the microscopic scale of bacteria and can be equally applicable to centimeter-sized swimming robots.

Our findings can lead to design and control of bioinspired soft robots that applies the same mechanism for swimming. Such robots will require a single control input, leading to simplified control scheme. The structure of these robots will be comprised of motor embedded in the head and one flexible flagellum, leading to simplified structural design. The material of the flagellum is linear elastic and, as such, a wide variety of polymeric materials can be used. The trajectory design in the current study has some limitations, i.e. after each turn, we allow the flagellum to relax to its unbuckled helical state and only then we can consider another turn. The numerical tool, in future, can be used in conjunction with machine learning to develop model-based control policies that are more generic. Importantly, we used a long range hydrodynamic force model and our framework may be extended to study multi-flagellated systems. In these more complex systems, the numerical simulation can be used as the source of data for machine learning-based data-driven models. We hope that our numerical results can inspire future work on all these fronts to motivate fundamental understanding of the biophysics of microorganisms and support the design of advanced functional soft robots.

Conflicts of interest

There are no conflicts to declare.

Notes and references

  1. E. M. Purcell, Proc. Natl. Acad. Sci. U. S. A., 1997, 94, 11307–11311 CrossRef CAS PubMed .
  2. H. C. Berg and R. A. Anderson, Nature, 1973, 245, 380 CrossRef CAS PubMed .
  3. E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 096601 CrossRef .
  4. M. Kim, J. C. Bird, A. J. Van Parys, K. S. Breuer and T. R. Powers, Proc. Natl. Acad. Sci. U. S. A., 2003, 100, 15481–15485 CrossRef CAS PubMed .
  5. R. M. Macnab and M. K. Ornston, J. Mol. Biol., 1977, 112, 1–30 CrossRef CAS PubMed .
  6. K. Son, J. S. Guasto and R. Stocker, Nat. Phys., 2013, 9, 494 Search PubMed .
  7. M. T. Brown, B. C. Steel, C. Silvestrin, D. A. Wilkinson, N. J. Delalez, C. N. Lumb, B. Obara, J. P. Armitage and R. M. Berry, J. Bacteriol., 2012, JB-00209 Search PubMed .
  8. Y. Man, W. Page, R. J. Poole and E. Lauga, Phys. Rev. Fluids, 2017, 2, 123101 CrossRef .
  9. F. T. Nguyen and M. D. Graham, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2018, 98, 042419 CrossRef CAS .
  10. N. C. Darnton and H. C. Berg, Biophys. J., 2007, 92, 2230–2236 CrossRef CAS PubMed .
  11. C. Calladine, Nature, 1975, 255, 121 CrossRef CAS PubMed .
  12. M. A. Constantino, M. Jabbarzadeh, H. C. Fu, Z. Shen, J. G. Fox, F. Haesebrouck, S. K. Linden and R. Bansil, Sci. Rep., 2018, 8, 14415 CrossRef PubMed .
  13. M. K. Jawed, N. Khouri, F. Da, E. Grinspun and P. M. Reis, Phys. Rev. Lett., 2015, 115, 168101 CrossRef CAS PubMed .
  14. R. Vogel and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2012, 35, 15 CrossRef CAS PubMed .
  15. E. E. Riley, D. Das and E. Lauga, Sci. Rep., 2018, 8, 10728 CrossRef PubMed .
  16. I. Spöring, V. A. Martinez, C. Hotz, J. Schwarz-Linek, K. L. Grady, J. M. Nava-Sedeño, T. Vissers, H. M. Singer, M. Rohde and C. Bourquin, et al. , PLoS Biol., 2018, 16, e2006989 CrossRef PubMed .
  17. J. Edd, S. Payen, B. Rubinsky, M. L. Stoller and M. Sitti, Intelligent Robots and Systems, 2003. (IROS 2003). Proceedings. 2003 IEEE/RSJ International Conference on, 2003, pp. 2583–2588.
  18. R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone and J. Bibette, Nature, 2005, 437, 862 CrossRef CAS PubMed .
  19. P. Tierno, R. Golestanian, I. Pagonabarraga and F. Sagués, Phys. Rev. Lett., 2008, 101, 218304 CrossRef PubMed .
  20. L. Zhang, J. J. Abbott, L. Dong, K. E. Peyer, B. E. Kratochvil, H. Zhang, C. Bergeles and B. J. Nelson, Nano Lett., 2009, 9, 3663–3667 CrossRef CAS PubMed .
  21. A. Ghosh and P. Fischer, Nano Lett., 2009, 9, 2243–2245 CrossRef CAS PubMed .
  22. K. E. Peyer, E. Siringil, L. Zhang and B. J. Nelson, Bioinspiration Biomimetics, 2014, 9, 046014 CrossRef CAS PubMed .
  23. Y. Osada, H. Okuzaki and H. Hori, Nature, 1992, 355, 242 CrossRef CAS .
  24. D. Ahmed, T. Baasch, B. Jang, S. Pane, J. Dual and B. J. Nelson, Nano Lett., 2016, 16, 4968–4974 CrossRef CAS PubMed .
  25. M. Kaynak, A. Ozcelik, A. Nourhani, P. E. Lammert, V. H. Crespi and T. J. Huang, Lab Chip, 2017, 17, 395–400 RSC .
  26. S. Sánchez, L. Soler and J. Katuri, Angew. Chem., Int. Ed., 2015, 54, 1414–1444 CrossRef PubMed .
  27. B. Kim, D.-H. Kim, J. Jung and J.-O. Park, Smart Mater. Struct., 2005, 14, 1579 CrossRef CAS .
  28. G. Kirchhoff, J. Reine Angew. Math., 1859, 56, 285–313 Search PubMed .
  29. J. Lighthill, SIAM Rev., 1976, 18, 161–230 CrossRef .
  30. F. T. Nguyen and M. D. Graham, Biophys. J., 2017, 112, 1010–1022 CrossRef CAS PubMed .
  31. M. Jabbarzadeh and H. C. Fu, Phys. Rev. E, 2018, 97, 012402 CrossRef CAS PubMed .
  32. A. Thawani and M. S. Tirumkudulu, J. Fluid Mech., 2018, 835, 252–270 CrossRef .
  33. M. Bergou, M. Wardetzky, S. Robinson, B. Audoly and E. Grinspun, ACM Trans. Graph., 2008, 27, 63 CrossRef .
  34. M. Bergou, B. Audoly, E. Vouga, M. Wardetzky and E. Grinspun, ACM Transactions on Graphics (TOG), 2010, p. 116 Search PubMed .
  35. J. Higdon, J. Fluid Mech., 1979, 90, 685–711 CrossRef .
  36. M. Jawed and P. M. Reis, Phys. Rev. Fluids, 2017, 2, 034101 CrossRef .
  37. M. K. Jawed, A. Novelia and O. M. O'Reilly, A primer on the kinematics of discrete elastic rods, Springer, 2018 Search PubMed .
  38. A. T. Chwang and T. Y.-T. Wu, J. Fluid Mech., 1975, 67, 787–815 CrossRef .
  39. J. Blake and A. Chwang, J. Eng. Math., 1974, 8, 23–29 CrossRef .
  40. J. Martindale, M. Jabbarzadeh and H. Fu, Phys. Fluids, 2016, 28, 021901 CrossRef .
  41. B. Duričković, A. Goriely and J. H. Maddocks, Phys. Rev. Lett., 2013, 111, 108103 CrossRef PubMed .
  42. M. K. Jawed and P. M. Reis, Soft Matter, 2016, 12, 1898–1905 RSC .
  43. ESI available: video showing a uniflagellar system following a 2D square trajectory.


Electronic supplementary information (ESI) available. See DOI: 10.1039/c9sm01843c

This journal is © The Royal Society of Chemistry 2020