2.1 Geometric model

The model bacterium consists of one rigid spherical cell body, one puller flagellum (dark slate grey) and one pusher flagellum (grey), as shown in figure 1. The cell body has centroid position denoted by $\boldsymbol {X}^{(B)}$ and orientation described by the basis $\{\boldsymbol {e}_1^{(B)},\boldsymbol {e}_2^{(B)},\boldsymbol {e}_3^{(B)}\}$. The pusher flagellum has position $\boldsymbol {X}^{(1)}$ and basis $\{\boldsymbol {e}_1^{(1)},\boldsymbol {e}_2^{(1)},\boldsymbol {e}_3^{(1)}\}$, and the puller flagellum has position $\boldsymbol {X}^{(2)}$ and basis $\{\boldsymbol {e}_1^{(2)},\boldsymbol {e}_2^{(2)},\boldsymbol {e}_3^{(2)}\}$. All of the aforementioned bases are right-handed and orthonormal. In this study, it is assumed that the two flagella have identical physical and elastic properties and their initial and rest configurations are right-handed helices with centrelines given by

(2.1)\begin{equation} \boldsymbol{\varLambda}^{(i)}(\xi) = \boldsymbol{X}^{(i)} + \xi\boldsymbol{e}_1^{(i)} + \varXi(\xi) \cos{(2{\rm \pi}\xi/p)}\boldsymbol{e}_2^{(i)} + \varXi(\xi)\sin{(2{\rm \pi}\xi/p)}\boldsymbol{e}_3^{(i)},\end{equation}

where $i=1,2$ for the pusher and the puller flagella, respectively; the variable $\xi$ parameterizes the distance along the axis of the helix with $0\leqslant \xi \leqslant L_F$, $\varXi (\xi )=a(1-\exp ({-(k_E\xi )^2}))$ is the helix amplitude function and $a$, $p$ and $k_E$ represent the maximum helix amplitude, the helix pitch and the amplitude growth rate (the amplitude grows from zero to $a$ over a region of length roughly ${2}/{k_E}$), respectively. The position and orientation of the puller flagellum is specified by two angles $\alpha$ and $\beta$ defined on the $\boldsymbol {e}_1^{(B)}$–$\boldsymbol {e}_3^{(B)}$ plane through the centre of the cell body. The pusher flagellum is placed symmetrically on the other side of the cell body with the same acute angles as the puller flagellum (figure 1). In all simulations, except one set of simulations where we specifically study the influence of the angle $\beta$, we assume that the rotors are normal to the cell membrane (i.e. $\beta =\alpha$). In studying the influence of $\beta$, the angle $\alpha$ is constant and $|\alpha -\beta |$ represents how much the rotors deviate from being normal to the cell membrane.

Figure 1. A schematic view of the model bacterium in which different bases and vectors are used to describe the position and orientation of the components. Here, the $\alpha$ and $\beta$ angles represent the position and orientation of the rotors on the cell body, defined with respect to $\boldsymbol {e}_1^{(B)}$. The internal moment between the $n$th and $(n+1)$th segments is denoted by $\boldsymbol {N}^{(i)n+{1}/{2}}$. Note that the thickness of the flagella in the figures does not reflect the actual thickness of the flagella in the model bacterium.

The experimental observations have shown that the cell body of MC-1 is approximately spherical (Reference Bente, Mohammadinejad, Charsooghi, Bachmann, Codutti, Lefèvre, Klumpp and FaivreBente et al., 2020; Reference Tokárová, Perumal, Nayak, Shum, Kašpar, Rajendran and NicolauTokárová et al., 2021). Based on the measurements for MC-1, the cell body diameter is $1.3\pm 0.1\,\mathrm {\mu }$m and the flagellum length is $3.3\pm 0.4\,\mathrm {\mu }$m (Reference Bente, Mohammadinejad, Charsooghi, Bachmann, Codutti, Lefèvre, Klumpp and FaivreBente et al., 2020). We are unaware of any study that measures the flexibility of the flagella or flagellar bundles in MC-1. Recalling that a filament in our model represents multiple flagella in a bundle, we use values of the rigidity approximately 1.5–11 times that of a flagellum in E. coli ($3.5\,{\rm pN}\,{\mathrm {\mu }}{\rm m}^2$ Reference Darnton and BergDarnton & Berg, 2007). Other parameters in this study, such as the motors’ torques, helical pitch and amplitude are chosen from the values given by Reference Bente, Mohammadinejad, Charsooghi, Bachmann, Codutti, Lefèvre, Klumpp and FaivreBente et al. (2020), Reference Mohammadinejad, Faivre and KlumppMohammadinejad et al. (2021) and Reference ShumShum (2019). Here, these parameters are non-dimensionalized by the averaged cell body radius $0.65\,\mathrm {\mu }$m (Reference Bente, Mohammadinejad, Charsooghi, Bachmann, Codutti, Lefèvre, Klumpp and FaivreBente et al., 2020), MC-1's motor torque which is approximately $12\,{\rm pN}\,\mathrm {\mu }$m (Reference Bente, Mohammadinejad, Charsooghi, Bachmann, Codutti, Lefèvre, Klumpp and FaivreBente et al., 2020) and the swimming fluid viscosity $\mu =10^{-3}$ Pa s.

2.2 Hydrodynamic interactions

In modelling of bacterial locomotion, the Reynolds number is very small and so the hydrodynamic interactions are governed by the incompressible Stokes equations

(2.2a,b)\begin{equation} {-}\boldsymbol{\nabla} p + \mu \Delta \boldsymbol{u} + \boldsymbol{F_b} = \boldsymbol{0},\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0,\end{equation}

where $\mu$ is the fluid viscosity, $p$ is the fluid pressure, $\boldsymbol {u}$ is the fluid velocity and $\boldsymbol {F_b}$ is the force per unit volume applied to the fluid by the immersed body. The model bacterium exerts a distribution of viscous stress $\boldsymbol {f}_{head}$ over the surface $S$ of the cell body (three-dimensional spherical cell body) and distributions of the viscous stress $\boldsymbol {f}_{fla}$ and viscous torque $\boldsymbol {n}$ along the centrelines of the flagella $\varGamma ^{(i)}$. We use the superscript $(i)$ to distinguish the pusher flagellum $(i=1)$ from the puller one $(i=2)$. Since both flagella are right-handed, the rotors in the pusher and puller flagella rotate in $-\boldsymbol {e}_1^{(1)}$ and $\boldsymbol {e}_1^{(2)}$ directions, respectively, to provide propulsion. In a Lagrangian description, the elastic filaments $\varGamma ^{(i)}$, which rotate and deform in time, and the cell body $S$ can be represented by a three-dimensional space curve $\boldsymbol {\gamma }(s,t)$ and surface $\,\boldsymbol {\varPsi }(\theta,\phi,t)$, respectively. The variables $s$, $\theta$ and $\phi$ are material coordinates along the filament (initialized as arclength) and the cell body surface, respectively, and $t$ is time. To ease the presentation, we let the variables $\boldsymbol {f}_{fla}(s,t)$, $\boldsymbol {n}(s,t)$ and $\boldsymbol {\gamma }(s,t)$ denote the respective quantities for both flagella with the understanding that the integral over $\varGamma ^{(i)}\ (i=1,2)$ involves the variables associated with the respective flagellum. The body force can be written as

(2.3)\begin{align} \boldsymbol{F_b}(\boldsymbol{x},t)& = \oint_{S} {\boldsymbol{f}_{head}(\theta,\phi,t) \delta(\boldsymbol{x}-\boldsymbol{\varPsi}(\theta,\phi,t))\,{\rm d} A} + \sum_{i = 1}^{2}{\int_{\varGamma^{(i)}} {\boldsymbol{f}_{fla}(s,t)\delta(\boldsymbol{x}-\boldsymbol{\gamma}(s,t))\,{\rm d} s}} \nonumber\\ &\quad + \sum_{i = 1}^{2}{\frac{1}{2}\boldsymbol{\nabla}\times\int_{\varGamma^{(i)}} {\boldsymbol{n}(s,t)\delta(\boldsymbol{x}-\boldsymbol{\gamma}(s,t))\,{\rm d} s}}. \end{align}

The evaluation point $\boldsymbol {x}$ can be anywhere in $\mathbb {R}^3$ including the model bacterium surface. To represent the flagellum thickness and avoid the singularity and also to enhance the stability of the solution, we use a regularized Stokes formulation for the flagella. Therefore, we replace the delta function ($\delta$) by a cutoff function $\psi _\epsilon$ in the second and third integrals. The radially symmetric cutoff function approximates the delta function in three dimensions. Following previous studies using the method of regularized stokeslets (Reference Cortez, Fauci and MedovikovCortez, Fauci, & Medovikov, 2005; Reference Olson, Lim and CortezOlson et al., 2013), we choose the cutoff function as

(2.4)\begin{equation} \psi_\epsilon(\boldsymbol{x}) = \frac{15\epsilon^4}{8{\rm \pi} ({\Vert \boldsymbol{x} \Vert}^2 + \epsilon^2)^{7/2}},\end{equation}

where we assume that $\epsilon ={d}/{2}$ to represent the filament radius. The solution for the flow field satisfying (2.2a,b) and (2.3) can be written as the boundary integral equations (Reference Cortez, Fauci and MedovikovCortez et al., 2005; Reference Olson, Lim and CortezOlson et al., 2013)

(2.5)\begin{align} \boldsymbol{u}(\boldsymbol{x},t)& = \oint_{S}{\boldsymbol{f}_{head}J_3 + [\,\boldsymbol{f}_{head}\boldsymbol{\cdot} (\boldsymbol{x}-\boldsymbol{\varPsi})](\boldsymbol{x}-\boldsymbol{\varPsi})J_4\,{\rm d} A} + \sum_{i = 1}^{2}{\int_{\varGamma^{(i)}}{\boldsymbol{f}_{fla}J_1 + [\,\boldsymbol{f}_{fla}\boldsymbol{\cdot} (\boldsymbol{x}-\boldsymbol{\gamma})](\boldsymbol{x}-\boldsymbol{\gamma})J_2\,{\rm d} s}} \nonumber\\ &\quad + \sum_{i = 1}^{2}{\frac{1}{2}\boldsymbol{\nabla}\times\int_{\varGamma^{(i)}}{[\boldsymbol{n} \times(\boldsymbol{x}-\boldsymbol{\gamma})]P_1\,{\rm d} s}}; \end{align}

in this regard, the angular velocity can be expressed as

(2.6)\begin{align} \boldsymbol{\omega}(\boldsymbol{x},t) = \frac{1}{2}\boldsymbol{\nabla}\times\boldsymbol{u}(\boldsymbol{x},t)& = \frac{1}{2} \oint_{S}{[\,\boldsymbol{f}_{head}\times(\boldsymbol{x}-\boldsymbol{\varPsi})]P_2\,{\rm d} A} + \sum_{i = 1}^{2}{\frac{1}{2}\int_{\varGamma^{(i)}}{[\,\boldsymbol{f}_{fla}\times (\boldsymbol{x}-\boldsymbol{\gamma})]P_1\,{\rm d} s}}\nonumber\\ &\quad + \sum_{i = 1}^{2}{\frac{1}{4}\int_{\varGamma^{(i)}}{[\boldsymbol{n}K_1 + \boldsymbol{n}\boldsymbol{\cdot} (\boldsymbol{x}-\boldsymbol{\gamma})(\boldsymbol{x}-\boldsymbol{\gamma})K_2]\,{\rm d} s}}, \end{align}

where $J_1,\ldots,J_4$, $P_1$, $P_2$, $K_1$ and $K_2$ are defined in terms of the distances between the points, the regularized parameter ($\epsilon$) and the viscosity, and are given in the supplementary material, § S1, available at https://doi.org/10.1017/flo.2022.34.

To evaluate these integrals, a finite number ($N_{head}$) of triangular elements is generated on the cell body surface and each flagellum is discretized into $N_{fl}$ connected equal-length straight rods. To yield a better accuracy, a tessellation of curved triangles is used to cover the cell body surface, where six nodes are required on the surface to construct an element. Three nodes are vertices and other three nodes are at the middle of the three edges (see figure 2a). Then, by following the scheme presented by Reference PozrikidisPozrikidis (2002), the surface of a curved triangle is mapped into a right-angle isosceles flat triangle. In the next step, Gauss–Legendre quadrature method with 12 Gauss points is implemented to evaluate the integrals over the standard triangles. The stokeslets (red points) are calculated by using cardinal interpolation functions and interpolating the nodal force densities in the evaluation points (black points) (Reference ShumShum, 2011). The integrals along segments of the flagella are evaluated by employing 8 Gauss points. In this regard, a second-order polynomial function is used as an interpolant to calculate the stokeslets along the segments. The evaluation points on the flagella are located at the middle and ends of each segment. Below, we let $N_{EPH}$ denote the number of evaluation points on the cell body and $N_{EPF}=2N_{fl}+1$ denote the number of evaluation points on each flagellum. By applying the presented scheme to (2.5) and (2.6) and satisfying the no-slip boundary condition on the model bacterium, a linear relationship between the nodal force and torque densities and the translational and angular velocities of the evaluation points is constructed. The relationship can be summarized in the form

(2.7)\begin{equation} \begin{bmatrix}\boldsymbol{u}_1\\ \vdots\\ \boldsymbol{u}_{N_{EPH} + 2N_{EPF}}\\ \boldsymbol{\omega}_1\\ \vdots\\ \boldsymbol{\omega}_{2N_{EPF}} \end{bmatrix} = \begin{bmatrix} {\boldsymbol{\mathsf{A}}}_1\\ {\boldsymbol{\mathsf{A}}}_2 \end{bmatrix} \begin{bmatrix} \boldsymbol{f}_1\\ \vdots\\ \boldsymbol{f}_{N_{EPH} + 2N_{EPF}}\\ \boldsymbol{n}_1\\ \vdots\\ \boldsymbol{n}_{2N_{EPF}} \end{bmatrix}, \end{equation}

where $\boldsymbol {u}_i$ is the translational velocity of the $i$th evaluation point on the cell body and the flagella, and $\boldsymbol {\omega }_j$ is the angular velocity of $j$th evaluation point on the flagella. Here, ${\boldsymbol{\mathsf{A}}}_1$ and ${\boldsymbol{\mathsf{A}}}_2$ are dense matrices with dimensions of $3(N_{EPH}+2N_{EPF}) \times 3(N_{EPH}+4N_{EPF})$ and $3(2N_{EPF}) \times 3(N_{EPH}+4N_{EPF})$ constructed by using the coefficients in (2.5) and (2.6), the mapping and the interpolation matrices.

Figure 2. (a) Distribution of stokeslets and rotlets over the curved triangular elements and along the connected straight rods. (b) Discretization of the flagellum into $N_{fl}$ segments. The triads locations and orientations are shown on two successive segments.

2.3 Kinematics

In the proposed model, the flagella complexes are driven by constant-torque motors. Let $\boldsymbol {U}^{(B)}$ and $\boldsymbol {\varOmega }^{(B)}$ denote the translational and rotational velocities of the cell body, respectively. Let $\boldsymbol {\omega }_s^{(i)1}$ denote the angular velocity vector of the first segment of the $i$th flagellum relative to the cell body and let $\boldsymbol {\omega }_s^{(i)n}$ denote the angular velocity of the $n$th segment of the $i$th flagellum with respect to the $(n-1)$th segment, for $n=2, 3,\ldots, N_{fl}$. Then, the overall instantaneous translational velocity of any given evaluation point $\boldsymbol {X}^E$ on the swimmer can be written as

(2.8) \begin{equation} \boldsymbol{U}(\boldsymbol{X}^E) = \begin{cases} \boldsymbol{U}^{(B)} + \boldsymbol{\varOmega}^{(B)}\times(\boldsymbol{X}^E - \boldsymbol{X}^{(B)}), & \boldsymbol{X}^E \textrm{ on cell body,}\\ \boldsymbol{U}^{(B)} + \boldsymbol{\varOmega}^{(B)}\times(\boldsymbol{X}^E - \boldsymbol{X}^{(B)}) + \displaystyle\sum_{n = 1}^{m}{\boldsymbol{\omega}_s^{(i)n}\times\boldsymbol{X}_{rel}^{(i)n}}, & \boldsymbol{X}^E \textrm{ on } m\textrm{th segment of $i$th flagellum,} \end{cases} \end{equation}

where

(2.9)\begin{equation} \boldsymbol{X}_{rel}^{(i)n} = \boldsymbol{X}^E-\boldsymbol{X}^{(i)}-\boldsymbol{\gamma}^{(i)n-{1}/{2}}\quad n = 1, 2,\ldots, N_{fl},\end{equation}

and $\boldsymbol {\gamma }^{(i)n-{1}/{2}}$ is the position vector of the $n$th joint with respect to the $i$th flagellum fixed frame (as shown in figure 1). The angular velocity of any given point $\boldsymbol {X}^E$ on the flagella can be written as

(2.10)\begin{equation} \boldsymbol{\omega}(\boldsymbol{X}^E) = \boldsymbol{\varOmega}^{(B)} + \sum_{n = 1}^{m}{\boldsymbol{\omega}_s^{(i)n}},\quad \boldsymbol{X}^E \textrm{ on $m$th segment of $i$th flagellum.} \end{equation}

Following (2.8) and (2.10), the translational velocities at the $N_{EPH}$ evaluation points on the cell body and the translational and rotational velocities at the $2N_{EPF}$ evaluation points on the flagella can be expressed in terms of $\boldsymbol {\omega }_s^{(i)n}$, $\boldsymbol {U}^{(B)}$ and $\boldsymbol {\varOmega }^{(B)}$ in the form

(2.11a,b)\begin{equation} \begin{bmatrix} \boldsymbol{u}_1\\ \vdots\\ \boldsymbol{u}_{N_{EPH} + 2N_{EPF}}\\ \end{bmatrix} = {\boldsymbol{\mathsf{A}}}_3 \begin{bmatrix} \boldsymbol{U}^{(B)}\\ \varpi \end{bmatrix},\quad \begin{bmatrix} \boldsymbol{\omega}_1\\ \vdots\\ \boldsymbol{\omega}_{2N_{EPF}}\\ \end{bmatrix} = {\boldsymbol{\mathsf{A}}}_4\varpi,\quad \varpi = \begin{bmatrix} \boldsymbol{\varOmega}^{(B)}\\ \boldsymbol{\omega}_s^{(1)1}\\ \vdots\\ \boldsymbol{\omega}_s^{(1)N_{fl}}\\ \boldsymbol{\omega}_s^{(2)1}\\ \vdots\\ \boldsymbol{\omega}_s^{(2)N_{fl}} \end{bmatrix}. \end{equation}

The matrices ${\boldsymbol{\mathsf{A}}}_3$ and ${\boldsymbol{\mathsf{A}}}_4$ are determined by using the position vectors employed in (2.8) and (2.10). Since the position vectors vary as the model bacterium makes the progress, these matrices are updated at each time step.

2.4 Elasticity

The hydrodynamic forces exerted on the flagella deform the flagella out of their static equilibrium configurations. To model the deformations, we use the standard Kirchhoff rod model and assume that the flagella are inextensible, unshearable and only allowed to bend and twist. The centrelines of the flagella at the initial/rest configurations are represented by the space curve $\boldsymbol {\gamma }(s,t)$. Right-handed orthonormal frames $\{\boldsymbol {D}^{(i)}_1(s,t),\boldsymbol {D}^{(i)}_2(s,t),\boldsymbol {D}^{(i)}_3(s,t)\}$ are introduced to describe the orientation of the material points in the cross-section of the flagella at $s$. To simplify the model, we assume that $\boldsymbol {D}^{(i)}_3(s,t)$ is always tangent to the curve $\boldsymbol {\gamma }(s,t)$ i.e. $\boldsymbol {D}^{(i)}_3(s,t)=\boldsymbol {\gamma }'(s,t)$. Based on the linear theory of elasticity, the internal moments $\boldsymbol {N}(s,t)$ transmitted along the flagella can be computed by (Reference Goriely and TaborGoriely & Tabor, 1997)

(2.12)\begin{equation} \boldsymbol{N}(s,t) = EI[(\kappa_1(s,t)-\hat{\kappa_1}(s))\boldsymbol{D}_1(s,t) + (\kappa_2(s,t)-\hat{\kappa_2}(s))\boldsymbol{D}_2(s,t) + \varUpsilon(\kappa_3(s,t)-\hat{\kappa_3}(s))\boldsymbol{D}_3(s,t)],\end{equation}

where $\boldsymbol {\kappa }(s,t)=(\kappa _1,\kappa _2,\kappa _3)$ is the twist vector at point $s$ and time $t$, $\boldsymbol {\hat {\kappa }}(s)$ is the rest twist vector and $\varUpsilon ={GJ}/{EI}$ is the ratio of twisting stiffness $GJ$ to bending stiffness $EI$. In this study, we assume that the flagella are homogeneous, isotropic and $\varUpsilon =1$ (Reference Park, Kim and LimPark et al., 2019a, Reference Park, Kim and Lim2019b).

Spatial discretization of each flagellum into $N_{fl}$ segments is done by introducing uniform intervals $\Delta {s}$ of the Lagrangian variable $s$. In our numerical scheme, as shown in figure 2(b), the triads $\boldsymbol {D}_{\hat {i}}^{(i)n}$ are placed at the middle of each segment with positions $\boldsymbol {\gamma }^{(i)n}$ and are updated over time as the segments rotate. To compute the twist vectors at the joints, an interpolation of two successive triads $\boldsymbol {D}_{\hat {i}}^{(i)n+{1}/{2}}$ is required. In this regard, the principal square root of the rotation matrix ${\boldsymbol{\mathsf{M}}}$ that maps the triad $\boldsymbol {D}_{\hat {i}}^{(i)n}$ to the triad $\boldsymbol {D}_{\hat {i}}^{(i)n+1}$ is used to interpolate the triads

(2.13a,b) \begin{equation} {\boldsymbol{\mathsf{M}}} = \sum_{{\hat{i}} = 1}^3{\boldsymbol{D}_{\hat{i}}^{(i)n + 1} (\boldsymbol{D}_{\hat{i}}^{(i)n})^{\rm T}},\quad \boldsymbol{D}_{\hat{i}}^{(i)n + {1}/{2}} = \sqrt{M}\boldsymbol{D}_{\hat{i}}^{(i)n}. \end{equation}

Following the scheme from Reference Lim, Ferent, Wang and PeskinLim et al. (2008), the discretized form of (2.12) is written as

(2.14a,b)\begin{equation} N_{\hat{i}}^{(i)n + {1}/{2}} = EI\left(\frac{\boldsymbol{D}_{\hat{j}}^{(i)n + 1}- \boldsymbol{D}_{\hat{j}}^{(i)n}}{\Delta{s}}\boldsymbol{\cdot}\boldsymbol{D}_{\hat{k}}^{(i)n + {1}/{2}}- \hat{\kappa}_{\hat{i}}^{(i)n + {1}/{2}}\right),\quad \boldsymbol{N}^{(i)n + {1}/{2}} = \sum_{{\hat{i}} = 1}^3{N_{\hat{i}}^{(i)n + {1}/{2}} \boldsymbol{D}_{\hat{i}}^{(i)n + {1}/{2}}},\end{equation}

where $({\hat {i}},{\hat {j}},{\hat {k}})$ is any cyclic permutation of $(1,2,3)$, $\boldsymbol {N}^{(i)n+{1}/{2}}$ is the internal moment transmitted from the $n$th to the $(n+1)$th segment of the $i$th flagellum, $n=0,1,\ldots,N_{fl}-1$ for both pusher and puller flagella and $\hat {\kappa }_{\hat {i}}^{(i)n+{1}/{2}}$ represents the twist vector components in the rest configuration. Also, $\boldsymbol {N}^{(i)({1}/{2})}$ is the internal moment transmitted from the rotor to the first segment of the $i$th flagellum. In the present scheme, the magnitude and direction of $\boldsymbol {N}^{(i)({1}/{2})}$ are determined by an iterative method to impose the motor torque and satisfy the Kirchhoff rod model, simultaneously. This method is further explained in supplementary material § S2.

2.5 Torque and force balance equations

Since bacteria swim at low Reynolds number, the inertial term is neglected and it is assumed that the total torques and forces acting on a bacterium complex are zero (Reference ShumShum, 2019; Reference Shum, Gaffney and SmithShum, Gaffney, & Smith, 2010) i.e.

(2.15)\begin{equation} \sum_{n = 1}^{N_{head}}{\int_{S_n}{\boldsymbol{f}_{head}\,{\rm d} A_n}} + \sum_{i = 1}^{2}{\sum_{n = 1}^{N_{fl}}{\int_{\varGamma^{(i)n}}{\boldsymbol{f}_{fla}\,{\rm d} s_n}}} = 0. \end{equation}

In our model, the total force balance equation (2.15) includes the integrals of viscous force densities over the triangular curved elements and along the straight segments of the flagella. By applying the Gauss–Legendre quadrature method, these integrals are expressed in terms of the nodal force densities at the evaluation points (i.e. $\boldsymbol {f}_1,\ldots, \boldsymbol {f}_{N_{EPH}+2N_{EPF}}$).

In the torque balance equation (2.16), the integrals represent the total viscous torques about the centre of the cell body. Like the force balance equation, the torque balance equation is also written in terms of the nodal force and torque densities at the evaluation points (i.e. $\boldsymbol {f}_1,\ldots, \boldsymbol {f}_{N_{EPH}+2N_{EPF}}, \boldsymbol {n}_1,\ldots, \boldsymbol {n}_{2N_{EPF}}$)

(2.16)\begin{align} \sum_{n = 1}^{N_{head}}{\int_{S_n}{(\boldsymbol{\varPsi}-\boldsymbol{X}^{(B)})\times \boldsymbol{f}_{head}\,{\rm d} A_n}} + \sum_{i = 1}^{2}{\sum_{n = 1}^{N_{fl}}{\int_{\varGamma^{(i)n}}{\boldsymbol{n}\,{\rm d} s_n}}} + \sum_{i = 1}^{2}{\sum_{n = 1}^{N_{fl}}{\int_{\varGamma^{(i)n}}{(\boldsymbol{X}^{(i)} + \boldsymbol{\gamma}-\boldsymbol{X}^{(B)})\times\boldsymbol{f}_{fla}\,{\rm d} s_n}}} = 0. \end{align}

To complete the system of equations, we balance the viscous torques about each joint of the flagellum chain with the transmitted internal moment, expressed as the equations

(2.17)\begin{equation} \sum_{n = m}^{N_{fl}}{\left(\int_{\varGamma^{(i)n}}{(\boldsymbol{\gamma}-\boldsymbol{\gamma}^{m-{1}/{2}})\times \boldsymbol{f}_{fla}\,{\rm d} s_n} + \int_{\varGamma^{(i)n}}{\boldsymbol{n}\,{\rm d} s_n}\right)} + \boldsymbol{N}^{(i)m-{1}/{2}} = \boldsymbol{0}, \end{equation}

where $m=1,\ldots,N_{fl}$ and $i=1,2$ for the pusher and the puller flagella, respectively. In fact, the torque balance equation (2.17) is written for all the joints and so $2N_{fl}$ equations are obtained in total. These equations are also written in terms of the nodal force and torque densities at the evaluation points on the flagella.

2.6 Overview

Before solving the equations, we apply the motors’ torques to their respective flagella. In our model, the orientation of the rotor is determined by the axial directions $\boldsymbol {e}_1^{(i)}$, $i=1,2$, which are fixed relative to the cell body frame. The transverse direction vectors $\boldsymbol {e}_2^{(i)}$ and $\boldsymbol {e}_3^{(i)}$ also rotate with the cell body and have an additional rotation about the $\boldsymbol {e}_1^{(i)}$ axis. At the joint connecting the flagellum to the rotor, the projection of the internal moment onto $\boldsymbol {e}_1^{(i)}$ is equal to the motor torque, i.e.

(2.18)\begin{equation} \boldsymbol{N}^{(i)({1}/{2})}\boldsymbol{\cdot}\boldsymbol{e}_1^{(i)} = T_i,\quad i = 1,2 . \end{equation}

In this equation, $\boldsymbol {e}_1^{(i)}$ and $T_i$ are known, and $\boldsymbol {N}^{(i)({1}/{2})}$ is determined by employing a sub-iterative method. In this method, the orientation of the rotor (specifically $\boldsymbol {e}_2$ and $\boldsymbol {e}_3$) is adjusted iteratively at each time step to satisfy the motor torque condition (2.18) and the Kirchhoff rod model (2.14a,b) (see supplementary material § S2). This method is used because the prescribed motor torque condition would generally not be satisfied if the rotor were updated with an explicit time-stepping scheme.

In this study, we characterize the flexibility of the flagella by a relative stiffness defined as

(2.19)\begin{equation} k_f = \frac{EI}{\bar{T}R} ,\end{equation}

where $E$ is the Young's modulus of the material, $I$ is the moment of inertia of the flagellum cross-section, $\bar {T} = ({T_1+T_2})/{2}$ is the averaged motor torque and $R$ is the radius of the cell body. By (2.19), the dimensionless relative stiffness value $k_f=1$ is achieved for the motor torque $\bar {T}=12\,pN\,\mathrm {\mu }$m, cell body radius $R=0.65\,\mathrm {\mu }$m and flexural rigidity $EI=7.8\,pN\,\mathrm {\mu }$m$^{2}$, which is 2.2 times the rigidity of an E. coli flagellum.

To sum up, substituting (2.11a,b) into (2.7) gives $3(N_{EPH}+4N_{EPF})$ linear equations in which the unknowns are the components of the nodal force and torque densities at the evaluation points, the components of the angular velocities of the segments ($3(2N_{fl})$ unknowns) and the components of the cell body's angular and translational velocities (${\rm six}$ unknowns). By adding (2.15), (2.16) and (2.17), a system of linear equations is constructed. In our study, linsolve solver in Matlab is used to evaluate the system of the equations and determine the unknowns.

We use quaternions to represent orientations and rotations between frames of reference. As a brief introduction, a quaternion $q$ has four components and is defined as the sum of a scalar part and a vector part

(2.20)\begin{equation} q = q_0 + \boldsymbol{q} = q_0 + q_1\hat{i} + q_2\hat{j} + q_3\hat{k}. \end{equation}

Only four of the nine components of a rotation matrix are independent and, in fact, the four components of a quaternion are sufficient to represent a rotation matrix. More details about the quaternion algebra and their relationships with rotation matrices are available in Reference Sarabandi and ThomasSarabandi and Thomas (2018); Reference ShepperdShepperd (1978) and Reference HornHorn (1987). If the connections between the frames are established by quaternions, the state of the microswimmer at each time step can be represented by a state vector $\boldsymbol {Q}$ defined as

(2.21)\begin{equation} \boldsymbol{Q} = \begin{bmatrix} \boldsymbol{X}^{(B)},q^{(B)},q^{(1)1},\ldots,q^{(1)N_{fl}},q^{(2)1},\ldots,q^{(2)N_{fl}},\end{bmatrix} \end{equation}

where $q^{(B)}$ represents the orientation of the cell body with respect to the global frame and $q^{(i)j}$, $i=1, 2,\ j=1, 2,\ldots, N_{fl}$ represent the orientations of the segments of the flagella with respect to the body frame. The state vector $\boldsymbol {Q}$ evolves over time according to the system of ordinary differential equations (ODEs)

(2.22a–c)\begin{equation} \dot{\boldsymbol{X}}^{(B)} = \boldsymbol{U}^{(B)},\quad \dot{q}^{(B)} = \tfrac{1}{2}{\boldsymbol{\mathsf{W}}}(q^{(B)})\boldsymbol{\varOmega}^{(B)},\quad \dot{q}^{(i)n} = \tfrac{1}{2}{\boldsymbol{\mathsf{W}}}(q^{(i)n})\sum_{j = 1}^{n}\boldsymbol{\omega}_s^{(i)j},\quad n = 1, 2,\ldots, N_{fl},\enspace i = 1,2,\end{equation}

where

(2.23)\begin{equation} {\boldsymbol{\mathsf{W}}}(q) = \begin{bmatrix} -q_1 & -q_2 & -q_3 \\ q_0 & q_3 & -q_2 \\ -q_3 & q_0 & q_1 \\ q_2 & -q_1 & q_0 \end{bmatrix}\end{equation}

and $\boldsymbol {U}^{(B)}$, $\boldsymbol {\varOmega }^{(B)}$ and $\boldsymbol {\omega }_s^{(i)n}$ have already been determined by solving the system of linear equations.

Combining a BEM with a Kirchhoff rod model to simulate the dynamics of a flexible filament in viscous fluid leads to a stiff set of ODEs that is computationally expensive to solve using general implicit schemes. Instead, we use an explicit multirate time integration scheme, as suggested by Reference Bouzarth and MinionBouzarth and Minion (2010) (without implementing their proposed spectral deferred corrections). In this approach, we update the nodal force densities and the angular and translational velocities of the cell body on coarse time steps while nodal force and torque densities and angular velocities of the flagellar segments are updated on finer time steps (100 fine time step per coarse time step).

This splitting procedure significantly decreases the computational cost because the non-stiff cell body portion is solved less frequently. Depending on the flagellar stiffness, the fine time step varies from $\delta t=5\times 10^{-5}$ to $\delta t=3\times 10^{-4}$ in our simulations. Comparing the multirate method with simply using a single (fine) time step in test simulations, we found a $55\,\%$ reduction in computational time and a $0.85\,\%$ difference in the computed net displacement.