Weicheng
Huang
and
M. K.
Jawed
*

Department of Mechanical and Aerospace Engineering, University of California, Los Angeles, Los Angeles, California 90095, USA. E-mail: khalidjm@seas.ucla.edu

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.

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 works^{13} showed that Kirchhoff elastic rod^{28} – 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).

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.

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

(1) |

(2) |

(3a) |

(3b) |

τ_{i} = θ^{i} − θ^{i−1} + _{i}, | (4) |

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

(5a) |

(5b) |

(5c) |

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 = from time t_{k} to t_{k+1} (h = t_{k+1} − t_{k} is the time step size) by integrating the following 4N − 1 equations of motion:

(6) |

(7) |

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 u_{j} ≡ [_{4j},_{4j+1},_{4j+2}]^{T}; recall that is the 4N − 1 sized velocity vector of the DOFs. This nodal velocity u_{j} can be decomposed into two components: (i) a flow velocity (u_{f})_{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 (u_{h})_{j} that is induced by the motion of the head. For the first component, we use LSBT to relate the velocity (u_{f})_{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 is^{13}

(8) |

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

U_{f} = F, | (9) |

Knowing the velocity field, U_{f}, generated by the flagellum and the matrix , 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 F^{ext} 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 F^{ext} ≡ [f_{0},0,f_{1},0,f_{2},…,0,f_{N−1}]^{T}. The force on the head, f_{0}, is described next.

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, U_{h} ≡ [_{0},_{1},_{2}], and angular velocity, Ω_{h} (along _{b}). The viscous flow at the j-th node from the spherical head is given by Higdon's model,^{32,35}

(10) |

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πμbU_{h} (and torque −8πμb^{3}Ω_{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}

(11) |

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 t_{k}, we know the DOF vector q_{k}, its velocity _{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 U_{f} 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

f_{0} = −6πμbU_{h} + f_{h}, | (12) |

(13) |

To quantify the buckling angular velocity, _{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, , 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 = 1.02. Regardless of the size of the head, we observe that, at small values of , 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, _{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, _{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 (−X_{b} in the body-fixed frame) decreases. Since a fluid flow along −X_{b} 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 = 0 and = 1.02 – as a function of angular velocity, ω. The distance ‖s‖ increases approximately linearly with angular velocity when < _{b}, and the flagellated system covers a straight path in this phase. This relation between ‖s‖ and 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 ( < _{b}). On the other hand, the translational displacement ‖s‖ shows a nonlinear behavior when the angular velocity is beyond the threshold _{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, = 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, = , as a function of . As expected from the data on torque vs. angular velocity, the required power increases quadratically with angular velocity.

In Fig. 5(a), we explain the turning process: the flagellated system is moving along an initial direction, d_{0}, at time t = 0 and it has to change its swimming direction to t_{1}. Without any loss of generality, we define a world frame x − y − z based on the body-fixed frame X_{b} − Y_{b} − Z_{b} at t = 0, such that the x, y, and z-axes are along −X_{b}, −Y_{b}, and Z_{b}-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 d_{0} to t_{1} using two parameters: turning angle, θ, and out of plane angle, ϕ. Turning angle, θ, is the angle between the initial direction, d_{0}, and the target direction, t_{1}. Out of plane angle, ϕ, is the angle between the target, t_{1}, and its projection, t^{xy}_{1}, on the x–y 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 t_{c} influences the turning process, quantified by θ and ϕ. In Fig. 5(d), we plot the turning angle, θ, and out of plane angle, ϕ, between d_{0} and t_{1}, as functions of t_{c}. Each value of t_{c} corresponds to a pair of θ and ϕ angles. We observe that the turning angle, θ, is nonlinearly related to the running time, t_{c}. Within a certain range (0 ≲ _{c} ≲ 0.55), the longer the buckled phase at ω_{h}, the greater the turning angle, θ. Depending on the duration of the higher angular velocity, t_{c}, 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, d_{0}, and final direction, t_{1}, we can compute a turning angle, θ*, and an out of plane angle, ϕ**. The turning angle, θ*, corresponds to a running time, t_{c}* from Fig. 5(d). This combination of θ* and t_{c}* 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, t_{c}*. This (θ*,t_{c}*) 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 t_{c} 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 x–y plane (z = 0; x–y–z 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 x–y plane. As desired, the trajectory remains almost entirely in the x–y plane, with relatively small displacement along the z-axis, Δz/ΔL < 5%, during the buckled phase. The corresponding control signal – the angular velocity as a function of time – is presented in Fig. 6(c).

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. |

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.

- E. M. Purcell, Proc. Natl. Acad. Sci. U. S. A., 1997, 94, 11307–11311 CrossRef CAS PubMed .
- H. C. Berg and R. A. Anderson, Nature, 1973, 245, 380 CrossRef CAS PubMed .
- E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 096601 CrossRef .
- 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 .
- R. M. Macnab and M. K. Ornston, J. Mol. Biol., 1977, 112, 1–30 CrossRef CAS PubMed .
- K. Son, J. S. Guasto and R. Stocker, Nat. Phys., 2013, 9, 494 Search PubMed .
- 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 .
- Y. Man, W. Page, R. J. Poole and E. Lauga, Phys. Rev. Fluids, 2017, 2, 123101 CrossRef .
- F. T. Nguyen and M. D. Graham, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2018, 98, 042419 CrossRef CAS .
- N. C. Darnton and H. C. Berg, Biophys. J., 2007, 92, 2230–2236 CrossRef CAS PubMed .
- C. Calladine, Nature, 1975, 255, 121 CrossRef CAS PubMed .
- 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 .
- M. K. Jawed, N. Khouri, F. Da, E. Grinspun and P. M. Reis, Phys. Rev. Lett., 2015, 115, 168101 CrossRef CAS PubMed .
- R. Vogel and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2012, 35, 15 CrossRef CAS PubMed .
- E. E. Riley, D. Das and E. Lauga, Sci. Rep., 2018, 8, 10728 CrossRef PubMed .
- 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 .
- 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.
- R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone and J. Bibette, Nature, 2005, 437, 862 CrossRef CAS PubMed .
- P. Tierno, R. Golestanian, I. Pagonabarraga and F. Sagués, Phys. Rev. Lett., 2008, 101, 218304 CrossRef PubMed .
- 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 .
- A. Ghosh and P. Fischer, Nano Lett., 2009, 9, 2243–2245 CrossRef CAS PubMed .
- K. E. Peyer, E. Siringil, L. Zhang and B. J. Nelson, Bioinspiration Biomimetics, 2014, 9, 046014 CrossRef CAS PubMed .
- Y. Osada, H. Okuzaki and H. Hori, Nature, 1992, 355, 242 CrossRef CAS .
- D. Ahmed, T. Baasch, B. Jang, S. Pane, J. Dual and B. J. Nelson, Nano Lett., 2016, 16, 4968–4974 CrossRef CAS PubMed .
- M. Kaynak, A. Ozcelik, A. Nourhani, P. E. Lammert, V. H. Crespi and T. J. Huang, Lab Chip, 2017, 17, 395–400 RSC .
- S. Sánchez, L. Soler and J. Katuri, Angew. Chem., Int. Ed., 2015, 54, 1414–1444 CrossRef PubMed .
- B. Kim, D.-H. Kim, J. Jung and J.-O. Park, Smart Mater. Struct., 2005, 14, 1579 CrossRef CAS .
- G. Kirchhoff, J. Reine Angew. Math., 1859, 56, 285–313 Search PubMed .
- J. Lighthill, SIAM Rev., 1976, 18, 161–230 CrossRef .
- F. T. Nguyen and M. D. Graham, Biophys. J., 2017, 112, 1010–1022 CrossRef CAS PubMed .
- M. Jabbarzadeh and H. C. Fu, Phys. Rev. E, 2018, 97, 012402 CrossRef CAS PubMed .
- A. Thawani and M. S. Tirumkudulu, J. Fluid Mech., 2018, 835, 252–270 CrossRef .
- M. Bergou, M. Wardetzky, S. Robinson, B. Audoly and E. Grinspun, ACM Trans. Graph., 2008, 27, 63 CrossRef .
- M. Bergou, B. Audoly, E. Vouga, M. Wardetzky and E. Grinspun, ACM Transactions on Graphics (TOG), 2010, p. 116 Search PubMed .
- J. Higdon, J. Fluid Mech., 1979, 90, 685–711 CrossRef .
- M. Jawed and P. M. Reis, Phys. Rev. Fluids, 2017, 2, 034101 CrossRef .
- M. K. Jawed, A. Novelia and O. M. O'Reilly, A primer on the kinematics of discrete elastic rods, Springer, 2018 Search PubMed .
- A. T. Chwang and T. Y.-T. Wu, J. Fluid Mech., 1975, 67, 787–815 CrossRef .
- J. Blake and A. Chwang, J. Eng. Math., 1974, 8, 23–29 CrossRef .
- J. Martindale, M. Jabbarzadeh and H. Fu, Phys. Fluids, 2016, 28, 021901 CrossRef .
- B. Duričković, A. Goriely and J. H. Maddocks, Phys. Rev. Lett., 2013, 111, 108103 CrossRef PubMed .
- M. K. Jawed and P. M. Reis, Soft Matter, 2016, 12, 1898–1905 RSC .
- ESI† available: video showing a uniflagellar system following a 2D square trajectory.

## Footnote |

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

This journal is © The Royal Society of Chemistry 2020 |