Amoeboid propulsion of active solid bodies, vesicles and droplets: a comparison

Abstract

We present a unified discussion of three types of near-spherical amoeboid microswimmers, driven by periodic, axially symmetric, achiral deformations (swim strokes): a solid deformable body, a vesicle with incompressible fluid membrane, and a droplet. Minimal models are used, which characterize the swimmer type only by boundary conditions. We calculate the swimming velocities, the dissipated power and the Lighthill efficiencies within a second order perturbation expansion in the small deformation amplitudes. For solid bodies, we reproduce older results by Lighthill and Blake, for vesicles and for droplets we add new results. The unified approach allows for a detailed comparison between the three types of microswimmers. We present such comparisons for swim strokes made up of spherical harmonics of adjacent orders l and l + 1, as well as for a manifold of swim strokes, made up of spherical harmonics up to order l = 4, which respect volume- and surface-incompressibility. This manifold is two-dimensional, which allows to present swimming velocities and efficiencies in compact graphical form. In a race in which each swimmer can choose the stroke that maximizes its speed, the droplet always comes in first, the vesicle comes in second, while the particle finishes third. However, if the three swimmers perform the same stroke, other order of rankings become possible. The maximum of the total efficiency of a droplet is greater than that of a vesicle if the internal dissipation is small. The efficiency of the solid body turns out to be typically two orders of magnitude smaller than that of vesicles and droplets. Optimizing the Lighthill efficiency and optimizing the swimming velocity result in different optimal swim strokes.

---

1 Introduction

Self-propelling microorganisms have developed a wide variety of mechanisms to exert driving forces on different environments.^1,2 The crawling of cells on a solid substrate, for example, is accomplished by specific adhesion sites. This mechanism has been extensively studied since the work of Abercrombie.³ It has also been shown that self-propulsion in confinements such as solid micro-channels is possible without adhesion sites.^1,4 If microorganisms move in fluid environments, the transmitted forces are dominated by viscous drag, as the Reynolds number of the generated flow is very small. Biological cells create propelling tractions by time-asymmetric sequences of morphological changes, known as swim strokes, either through the beating of short (cilia) or longer (flagella) active filaments, or by shape deformations of the entire cell (amoeboid motion).

A number of eukaryotic cells were previously thought to need an underlying substrate to move (crawl) using shape deformations, but growing evidence now suggests that many types of cells are capable of migrating without focal adhesion, enabling them to swim just as effectively as they crawl. Among the organisms which have been found to move by shape changes, are not just amoebae,⁵ but also neutrophils,⁶ mutants of Dicties^6,7 and Euglenids.⁸ In fact, some of these organisms can change between different modes of swimming strategies depending on the environment. Neutrophils and Dicties were shown to perform chemotaxis in solution and, under appropriate environmental conditions, crawl on solid surfaces as well.^6,7

Inspired by the swimming of microorganisms, many efforts have been made to construct self-propelled or field-propelled synthetic microrobots, which can serve many different purposes in applications. These attempts start from a solid body, a vesicle, or a droplet, driven by a variety of active mechanisms.^9–17

In the present work, we focus on amoeboid motion which is understood in the following as motion caused by actively controlled shape deformations, irrespective of the origins of the deformations. Amoeboid motion differs from other phoretic or chemical propulsion mechanisms, which can operate without shape changes. Theoretical approaches¹⁸ help to clarify the fundamental hydromechanical principles, underlying this type of propulsion. The problem can be decomposed into two parts: the generation of deformations by a variety of active drives^{15,16,19–23} and the propulsion and fluid flows due to these deformations. Here, we concentrate on the second part of the problem, which builds on the work of Taylor,²⁴ who computed the swimming speed of an inextensible sheet performing wavelike oscillations. Lighthill²⁵ extended Taylor's work to a sphere with small surface deformations (squirmer) and computed its swimming velocity and energetic efficiency. Later, Lighthill's work was corrected and extended by Blake,²⁶ and it has been analyzed in great detail in the following years.^27–31 Energy dissipation and efficiency have been calculated and optimized with respect to possible shape deformations.^32,33 A two-dimensional variant of the model has been applied to eukaryotic cells,³⁴ which were experimentally observed to swim with only viscous traction.

As a bridge between the first and second part of the problem, models that start from active tractions localized at the surface^35,36 have also been considered. Deformations of droplets driven by Marangoni stresses,^37,38 by electric fields,³⁹ or by polar molecules⁴⁰ have been analyzed. Farutin et al.³⁶ model an amoeboid swimmer as a vesicle driven by active membrane forces that cause shape deformations and self-propulsion. More strongly coarse-grained models⁴¹ discuss the motion of a deformable particle in terms of coupled equations for the flow and a quadrupolar order parameter, whereas more detailed computational models^18,42–44 aim to mimic specific biological processes such as pseudopod formation, including biochemical reactions in the membrane.

The intention of the present work is not to discuss detailed models of specific biological or artificial amoeboid swimmers, but to study universal properties of this mode of swimming. Our objectives are twofold: first, we aim to understand amoeboid motion of active liquid droplets, which are driven by normal active tractions. We compute the propulsion velocity for given dynamic deformations. Second, we compare these results to those of vesicles and deformable solid bodies, which are considered a standard model for swimming micro-organisms. In order to achieve these goals, we derive a unified approach to minimal models of single, nearly spherical amoeboid microswimmers of three different types: (i) a deformable solid, (ii) a vesicle with a fluid, incompressible membrane, and (iii) a droplet. The swimmers move in a Newtonian fluid, and the vesicle and droplet are filled with another Newtonian fluid. For all three types, we calculate the swimming velocity and the efficiency due to given swim strokes using the same analytical perturbation method applicable for small swim stroke amplitudes.^45–47 Although we do not treat it in detail, it will be obvious that the method also gives the fluid flow fields. It reproduces the results of Lighthill²⁵ and Blake²⁶ for a squirmer, and it complements the results of Farutin et al.,³⁶ who consider given membrane tractions, which produce the swim strokes. The results for active droplets are new to the best of our knowledge. Subsequently, we compare the swimming velocities and efficiencies for given swim strokes between all three types of amoeboid swimmers. The deformations of the sphere during a swim stroke are constrained by volume incompressibility. To include the vesicle in the comparisons, they also have to obey the surface incompressibility constraint.

The models of amoeboid swimmers are introduced in the next section. In our discussion, the shape deformations are considered as given and there is no need to specify the details of the shape control mechanisms in terms of active tractions, external fields or other means. The perturbative solution for the swimming velocities up to second order is presented in Section 3, and the dissipation due to viscous fluid flow is calculated in Section 4. In Section 5, we compare the velocities and the hydrodynamic efficiencies of the three types of particles for periodic swimming strokes. In the simplest case, the average swimming velocity does not depend on the time course of the deformations. Beyond this case, we construct a two-parameter manifold of simple harmonic deformations. Within this manifold, we compare the velocities and efficiencies for all possible strokes and identify the corresponding optimal swim strokes. We also discuss the dependence of these quantities on the viscosities of the ambient and internal fluids. Our main results are summarized in Section 6. Some more technical points of the calculations are delegated to the Appendices A–F.

2 Models

We study the self-propulsion resulting from time-varying shape deformations of a microorganism or an artificial microswimmer, which can be modeled in one of three ways: as a deformable solid body, as a vesicle, or as a fluid droplet. These systems, which are all referred to as “particles” in the following, are submerged in a Newtonian fluid and are neutrally buoyant. The unit of mass is chosen such that the mass density ρ of the materials is 1. Each particle occupies a volume V with a smooth boundary ∂V. In the absence of particles, the ambient fluid is at rest in the laboratory frame. For low Reynolds number, the outer flow field obeys the Stokes equation

∇·σ = η∇²v − ∇p = 0, (1)

supplemented by the incompressibility condition ∇·v = 0. The viscosity of the external fluid is denoted by η, and its stress tensor σ is given by the Cartesian components σ_ij = −pδ_ij + η(∂_iv_j + ∂_jv_i) = −pδ_ij + σ^visc_ij, with the pressure p determined from incompressibility.

The particles are driven by small, time-dependent deformations of nearly spherical shapes. Such deformations, caused by active forces, have been calculated for vesicles³⁶ and droplets.³⁵ Here we take a different point of view: we consider the time-dependent deformations as given and compute the corresponding propulsion of the particle. The deviations of the shape of a particle from a reference sphere, which has the same volume as the particle, are described by a function f (see Fig. 1), which determines the positions r_s of surface points via the relationship

r_s = R(1 + f)e_r. (2)

Fig. 1 Illustration of setup and notation. The (uniaxial and achiral) shape changes (dr_s/dt = ḟ(θ,t)e_r) are assumed to be radial with respect to the center of a reference sphere of radius R (dashed line), which has the same volume as the particle. The interior can be a solid or a Newtonian fluid. For a droplet, the interface is a geometric surface, whereas for a vesicle it represents an incompressible fluid membrane.

We use R as the unit of length. The driving, which we will study in this work, is given by shape changes ḟe_r (the dot denotes the time derivative), which are restricted to radial directions with respect to the center of the reference sphere. We therefore call this point the center of deformation in the following. The radial drive can be studied analytically for all three types of particles by using a perturbation expansion, as will be shown below. For simplicity, we also restrict our considerations to uniaxial (and achiral) deformations f(θ,t). This assumption is easily lifted, but the uniaxial case is sufficient to quantitatively compare the properties of amoeboid translational motion of the different particles. Our aim is to calculate the propulsion velocities U and the dissipated powers Ė, which result from given deformations f(t). We can expand the function f(θ,t) in Legendre polynomials

and characterize the deformations by the amplitudes f_l. The sum does not include the term l = 1 that would correspond to rigid body motion. For all three types of particles, the f_l are restricted by the condition that the particle volume remains constant, as the interior fluid or solid is taken to be incompressible. The time dependence of the amplitudes f_l(t) has to be restricted to slow dynamics. More precisely we require the product of the Reynolds number and the Strouhal number to be small so that the partial time derivative in the Navier–Stokes equation can be neglected.⁴⁸

Solid bodies, vesicles, and droplets are modeled in the simplest possible way: we ignore the rheology of the interface and instead formulate appropriate boundary conditions which are presented below. The interior of vesicles and droplets consists of a Newtonian fluid with the same density as the ambient fluid, but in general with a different viscosity λη. The inner flow will contribute to the energy that is dissipated during the motion. For the solid, we do not include inner dissipation mechanisms because this would go beyond the scope of a hydrodynamical model and leaves many possibilities. As a consequence, the dissipated energy calculated from the outer flow is only a lower bound to the total dissipation of the solid particle.

All calculations will be performed in a frame, in which the center of deformation is at rest (CODF), because it is in this frame that the shapes f(θ,t) are defined. In the CODF, the flow field v(r,t) obeys the boundary condition

at infinity. The swimming velocity U is defined here as the velocity of the center of deformation and it points in the z-direction, due to the uniaxial symmetry, i.e.U = Ue_z. We could have chosen any other point of a deformable particle to define its velocity. From the perspective of dynamics, the center of mass is particularly emphasized. In the CODF, it moves with a non-vanishing velocity v_cm that is easily calculated from the deformations f(t),and does not depend on the type of particle. The velocity of the center of mass U_cm in the laboratory frame is then given by U_cm = U + v_cm.

The models for the particles are now completed by specifying the boundary conditions, which uniquely determine the flow field.

2.1 Deformable solid body

For a solid particle, we use no-slip boundary conditions on its surface ∂V, implying

Together with eqn (4), this determines the outer flow uniquely.

2.2 Vesicle

In the simplest model of a fluid membrane, we consider only the local inextensibility, which means that the surface divergence of the flow vanishes, that is ∇_s·v(r_s) = 0 for every point r_s on ∂V. We can rewrite this constraint in terms of viscous tractions t^vis_j = n_iσ^visc_ij,

∇_s·v(r) = t^vis(r)·n = 0, r ∈ ∂V. (7)

where n is the outward normal to the surface.

The shape evolution of the vesicle is described by a level set function H(r,θ,t) = r − r_s(θ,t) = 0 for all t, which implies the kinematic equation

The boundary condition (7) together with the kinematic eqn (8) and the far-field asymptotics eqn (4) uniquely determine the outer flow field v and the swimming velocity U for given f_l and ḟ_l. To calculate the dissipated energy, one also needs the flow and the stress inside the vesicle, which will be denoted by V and Σ. The internal flow obeys the same boundary conditions on the membrane as the flow of the surrounding fluid.

In contrast to particles and droplets, the deformations f_l of vesicles are further restricted because their surface area

with z = cos θ, remains constant. Instead of A, we will use the excess area Δ, defined by the relation A = 4πR²(1 + Δ), to express this restriction.⁴⁶

2.3 Droplet

The droplet consists of an incompressible Newtonian fluid with viscosity λη. The internal and external fluids are completely immiscible and are therefore separated by a sharp interface characterized by a constant surface tension. We use lower case symbols for the exterior and upper case symbols for the interior of the droplet. The inner flow field (which is again indicated by V) must remain finite in the center of deformation.

The boundary conditions require continuity of the flow across the interface

v(r) = V(r), r ∈ ∂V. (10)

In the following, we also assume continuity of the tangential tractions

(δ_ij − n_in_j)(t_j(r) − T_j(r)) = 0, r ∈ ∂V. (11)

This excludes inhomogeneous surface tensions or any other tangential forces at the interface and, therefore, restricts the class of possible physical mechanisms generating the deformations to active normal forces. Prominent examples of such mechanisms are droplet deformation by electric fields^23,39,49,50 and deformation by polar molecules with strong anchoring.⁴⁰

The balance of normal forces at the interface of a droplet with homogeneous surface tension [small gamma, Greek, circumflex] and normal active tractions t_act = t_actn is given by

n·(t(r) − T(r)) = [small gamma, Greek, circumflex]∇·n + t_act. (12)

This relation is not needed to determine the flow at given f(θ,t). However, it can be used to calculate the active tractions responsible for the given deformations if desired.³⁵ Using the expansion in small deformations, ∇·n = 2(1 − εf) − ε∇²f + [scr O, script letter O](ε²), we can write eqn (12) in linear order in ε as

t_act = t_r − T_r + [small gamma, Greek, circumflex](2f + ∇²f). (13)

This equation either determines the deformations for given active tractions or yields the active tractions for given deformations as worked out in detail in ref. 35, but is not needed here.

Just as for vesicles, the shape evolution of the droplet is described by a level set function and eqn (4), (8), (10) and (11) determine a unique solution for the outer and inner flow fields.

3 Perturbation expansion

To set up a perturbation expansion in small deformations of a spherical particle, we introduce an order counting parameter ε, which multiplies the amplitudes f_l and is set to 1 at the end of the calculations. It turns out that the expansion proceeds in powers of both deformation f and its time derivative ḟ; therefore, the expansion requires f and ḟ to be of the same order of smallness. Up to the second order, the constraints of constant volume and constant surface take the form⁴⁶

The volume constraint eqn (14) will be used to eliminate f₀ at [scr O, script letter O](ε²).

For uniaxial and achiral systems, we expand the general external solution of Stokes equations³⁵ in Legendre polynomials P_l(θ) and their derivatives :

In a perturbative approach, the coefficients a_l and b_l of the flow field are expanded in powers of ε,

a_l = εa⁽¹⁾_l + ε²a⁽²⁾_l + [scr O, script letter O](ε³), (18)

b_l = εb⁽¹⁾_l + ε²b⁽²⁾_l + [scr O, script letter O](ε³) (19)

giving rise to a corresponding expansion of the flow field: v = εv⁽¹⁾ + ε²v⁽²⁾ + [scr O, script letter O](ε³).

The perturbative analysis is carried out to the second order, since this is the minimum order that results in a non-zero swimming velocity U. In the following, we present the results for deformable solid particles, vesicles, and droplets. All three use the above expansion eqn (18) in the appropriate boundary conditions. The computations are easiest for the solid particle; we therefore give the explicit steps of the calculation for this case in the next subsection and delegate some details of the calculations for vesicle and droplet to the appendices.

3.1 Deformable solid body

As a first step we expand the boundary condition eqn (6) in powers of f,

v(1 + f,θ) = v(1,θ) + f∂_rv(1,θ) + … = ḟe_r. (20)

To lowest order this equation reads

with the solutionfor l ≥ 2. The boundary condition in second order becomes

Substitution of the first order solution into the above equation yields

These equations determine the flow field as a function of the prescribed deformations to [scr O, script letter O](ε²). In order to obtain the coefficients {a⁽²⁾_l,b⁽²⁾_l} explicitly, one has to project the above equations onto P_l, .

Projection of eqn (25) onto P₀ yields

Here we have used the orthogonality of the Legendre polynomials and the constraint of constant volume, , which guarantees that the projection onto P₀ vanishes. We are mainly interested in the propulsion velocity of the particle, and hence project eqn (25) onto P₁ and eqn (26) onto , the latter denoting the associated Legendre polymomial P¹₁. The calculation requires integrals of 3 associated Legendre polynomials, which are given in Appendix C. As a result of this procedure, we obtain 2 equations for a⁽²⁾₁, b⁽²⁾₁ and U:

The particle has to be autonomous, that is, the total force has to vanish, which implies that the flow has to decay faster than 1/r for large r and hence b⁽²⁾₁ = 0.

The solution of this system is a linear combination of the inhomogeneities (right hand sides) and therefore consists of terms proportional to f_lḟ_l+1 and to f_l+1ḟ_l. This property also applies to the other two types of particles, so that we can always give the swimming velocity in the form

For the solid body, the coefficients R_l and S_l are obtained from eqn (28) and take on the form

Note that the swimming velocity of deformable solid bodies is independent of the viscosity η, and therefore, it is a purely geometric quantity.²⁹ These results are not new, but have been derived previously, first by Lighthill²⁵ and later corrected by Blake.²⁶ We have presented them here in detail to demonstrate our strategy, which we now apply to both the droplet and the vesicle.

3.2 Vesicle

For the vesicle, the calculation steps are nearly identical to those used for the solid particle; it is only the form of the boundary conditions that differs. To first order, the boundary conditions (eqn (7) and (8)) readand the solution is

To second order we get

v⁽²⁾_r = ḟ₀ + f′v⁽¹⁾_θ, (34)

σ⁽²⁾_rr = 2f′σ_rθ⁽¹⁾ − f∂_rσ⁽¹⁾_rr. (35)

Substituting the first-order solution from eqn (17) and (33) and the results of Appendix B, yields

Using both the constant volume and the constant area constraint, the projection of the above two equations onto P₀ can be shown to vanish. The projections onto P₁ lead to two equations for a⁽²⁾₁, b⁽²⁾₁ and U. We require the system to be force free, implying b⁽²⁾₁ = 0, so that the resulting equations determine a⁽²⁾₁ and U in second order in the deformation. The swimming velocity takes the form of eqn (29) with

Just as for the solid particle, the swimming velocity U does not depend on the viscosity η. The velocity U differs from the result given in ref. 36 by a total time derivative, which may arise from a different choice of reference point within the deformable vesicle.

3.3 Droplet

The analysis of the droplet requires both inner and outer solutions of the Stokes equation. The outer solutions v are taken from eqn (16) and (17), the inner solutions V are given by

The boundary conditions in first order read

v⁽¹⁾_θ = V⁽¹⁾_θ, (43)

σ⁽¹⁾_θr = Σ⁽¹⁾_θr. (44)

The required stress components of the inner and outer flow are calculated in the Appendix B.

In second order the eqn (8), (10) and (11) take on the form:

v⁽²⁾_r + f∂_rv⁽¹⁾_r − f′v⁽¹⁾_θ = ḟ₀,

V⁽²⁾_r + f∂_rV⁽¹⁾_r − f′V⁽¹⁾_θ = ḟ₀,

v⁽²⁾_θ + f∂_rv⁽¹⁾_θ = V⁽²⁾_θ + f∂_rV⁽¹⁾_θ,

σ⁽²⁾_θr + f∂_rσ⁽¹⁾_θr + f′(σ⁽¹⁾_rr − σ⁽¹⁾_θθ) = Σ⁽²⁾_θr + f∂_rΣ⁽¹⁾_θr + f′(Σ⁽¹⁾_rr − Σ⁽¹⁾_θθ). (45)

We substitute the results of the first-order solution for the velocities and stress components in the above equations, which must be solved for the coefficients a⁽²⁾₁, c⁽²⁾₁, d̂⁽²⁾₁ and U. The resulting R_l and S_l in eqn (29) (and therefore also U^drop) are rational functions of the form

The second-order polynomials r⁽ⁿ⁾_l(λ) and s⁽ⁿ⁾_l(λ) can be found in Appendix D. Unlike a solid particle or a vesicle, the swimming velocity of a droplet depends on the viscosity contrast λ. Note that the l-components of the velocity U have well-defined limits for λ → 0 and for λ → ∞.

3.4 Center of mass velocity

The calculated swimming velocities U of the particles refer to the motion of the center of deformation in the laboratory frame. We can transform them to center-of-mass velocities in the laboratory frame by adding to U the second order result of eqn (5), given by

4 Energy dissipation

The dissipated power in an incompressible Stokes flow in a volume V is given by

Pressure does not contribute to this expression, so σ_ij can be replaced by viscous stress σ^visc_ij. In spherical coordinates the integral takes on the form

To evaluate the leading order, we need the flow field only in [scr O, script letter O](ε) and the integration is over the reference sphere. For both outer and inner flow fields, the two terms are evaluated in Appendix E.

Our model does not take into account interior dissipation of the solid particle. The power due to the ambient flow,

is therefore only a lower bound on the total dissipation.

For the vesicle, the dissipated power due to the outer flow takes on the form

and the inner flow contributes

The l-components of outer and inner dissipation of the droplet,

are given explicitly in Appendix E.

5 Comparisons

We will now use the general results of eqn (29) and (51)–(54) to discuss and compare the swimming velocities, the dissipation and the energetic efficiencies of the three types of particles, when they are driven by periodic deformations with period T.

As a result of Section 3, the average swimming velocities

in second order can be represented as a sum of terms Ū_l, which contain the time averages of Φ_l(t) = ḟ_l+1f_l. Note that total time derivatives, such as v_cm, do not change Ū. Therefore, the average swimming velocity is a physical quantity that is independent of the chosen body fixed reference frame. Differences in swimming velocities between the particle types arise solely from the numerical values of R_l − S_l. These values are constants for the solid and the vesicle (see Table 1), whereas for droplets they decrease with increasing viscosity contrast λ, such that the bubble (λ = 0) is the fastest.

Table 1 The coefficients R_l and S_l, which determine the swimming velocities (see eqn (29)) for l = 2 and 3. These coefficients are used in Section 5.4. In the third column, μ = (λ + 1)(3λ + 2)

	Solid	Vesicle	Droplet
R 2	−9/35	−47/35	−3(219λ^2 + 545λ + 286)/(245μ)
S 2	−6/35	−32/35	6(−47λ^2 − 61λ + 8)/(175μ)
R 3	−16/63	−96/63	−8(23λ^2 + 57λ + 31)/(63μ)
S 3	−5/63	−65/63	(−103λ^2 + 13λ + 160)/(147μ)

---

In the following, we discuss two special examples in detail to illustrate our general results: (1) periodic driving with only two adjacent amplitudes f_lf_l+1 and (2) simple harmonic motion with three adjacent amplitudes, which allow for optimization of strokes either with respect to speed or efficiency. Whenever necessary for de-dimensionalization, we use T/2π as a unit of time in the following subsections.

5.1 Swimming velocity ((l, l + 1)-strokes)

In general, the time average [capital Phi, Greek, macron]_l depends on the detailed functional form of the deformation amplitudes f_l(t) and f_l+1(t). For droplets and deformable solids these amplitudes are unconstrained, but for vesicles they have to obey the area constraint eqn (15). If only one (l, l + 1)-pair of amplitudes is nonvanishing, this restricts the manifold of possible deformations to an ellipse in the (f_l,f_l+1)-plane, and T[capital Phi, Greek, macron]_l becomes the area of the ellipse,^25,36 that iswith β_l = 2(2l + 1)/[(l + 2)(l − 1)]. This shows that Ū does not depend on the time course of the deformation amplitudes. For l = 2 we recover the result given in ref. 36.

In Fig. 2, we show the speed |Ū_l| of the droplet for two (l, l + 1)-pairs as a function of λ in comparison to those of a vesicle. All components U_l decrease with λ and have finite and non-vanishing limits for small and large λ as stated above. The l-dependence of the swimming speeds of a vesicle and a solid body are shown in the inset of Fig. 2.

Fig. 2 Average swimming speed Ū_l (divided by excess area Δ), generated from (f_l, f_l+1)-pairs, of a droplet (curved line) and a vesicle (horizontal lines) versus viscosity contrast λ for two pairs. Dotted lines correspond to l = 2, 3 and dash-dotted lines to l = 3, 4. The inset shows the λ-independent Ū_l/Δ for a vesicle and a solid body versus l.

In general, all |Ū_l| of droplets decrease with increasing λ. For small λ, the swimming speeds also decrease with increasing l, while for λ ≳ 1.3, they begin to increase with increasing l. Note also that unlike a vesicle the values of |Ū_l| for a solid body increase for l = 2, 3 before they start to decrease significantly from l = 5.

5.2 Dissipation

In leading order, each deformation amplitude f_l(t) contributes to the dissipated energy – even if only a single driving amplitude is present and the particle is not propelled. We make use of eqn (51)–(54) to compute the energy that is dissipated in a period T. For our model of a deformable solid, Dⁱⁿ = 0, so that Ē_l is independent of λ. For a vesicle it is linear in the viscosity contrast λ, and for a drop it is the quotient of two second order polynomials in λ, multiplied by a factor λ. Due to the absence of internal dissipation, the total dissipated energy of the solid particle will always be the smallest for large λ. In Fig. 3 we compare D_l = D^out_l + Dⁱⁿ_l of vesicle and drop (with the same internal and external viscosities) to that of the solid particle for l = 2 and 3.

Fig. 3 Dissipation coefficients D_lvs. viscosity contrast λ. Solid lines: drop, dashed lines: vesicle, dotted lines: solid particle. Upper lines refer to l = 2, lower lines to l = 3.

For all particle types, the dissipation coefficients D_l become smaller with increasing l. The dissipated energy is always lower for vesicles than for droplets. For droplets, D₂ (D₃) reaches the value of the solid particle at λ ≈ 2.9 (λ ≈ 3.4). For vesicles the crossover appears at larger values (λ ≈ 20 for l = 2 and 22 for l = 3).

5.3 Simple harmonic deformations with area constraint

For deformations consisting of more than two adjacent amplitudes f_l, f_l+1, the average speed Ū depends on the time course of the constrained swim strokes. In the following, we restrict the time dependence to simple harmonic motion, for which a more detailed discussion is possible. The deformations consist of a sum of terms f_l(t)P_l(cos θ) with

f_l(t) = F_l(e^iα_l+iωt + e^{−iα_l−iωt}). (57)

Such strokes lead to particle velocities

when inserted in eqn (29). Note that one of the phases α_l can be arbitrarily chosen by appropriately adjusting the origin of the time axis. The time-independent term is the mean swimming velocity Ū. Fourier components with different frequencies ω_l for different l-components do not contribute to the time average Ū, so we can consider a single ω without loss of generality. The choice of time unit 2πT corresponds to ω = 1 for the harmonic deformations.

In the following, we discuss the simple non-trivial case of three l components, l = 2, 3 and 4. A swim stroke is then characterized by the parameters F₂, F₃, F₄, α₂, α₃ (α₄ is set to zero by choosing a suitable initial time). These five parameters must satisfy three equations, derived by inserting eqn (57) into eqn (15) and separating the time independent and time dependent parts. The resulting two-parameter manifold can be conveniently described by the l = 2 parameters F₂, α₂, as detailed in Appendix F. For each pair (F₂, α₂) there are two swim strokes, corresponding to motion in either the (+z)- or the (−z)-direction with identical speed. Consistent strokes can be found for every α₂, while the range of F₂ is limited to 0 ≤ F₂ ≤ F_2,max (in Appendix F we show that ). Now we can compare the trajectories, the swimming velocities and the dissipated power of a solid body, a vesicle, and a drop (all three of equal volume), which execute all possible area-conserving swim strokes composed of l = 2, 3 and 4 components.

5.4 Swimming velocities with area constraint

For a comparison of average swimming velocities with constrained (l = 2, 3, 4)-swim strokes, we use the entries of Table 1. An example of trajectories for a fixed stroke is shown in Fig. 4. It illustrates that the speeds and oscillation amplitudes vary between particles, and it is even possible for the same swim stroke to propel different particles in opposite directions.

Fig. 4 Example of trajectories: α₂/2π = 0.952, F₂/F_2,max = 0.78, λ = 1. Solid line: droplet, dashed line:vesicle, dotted line:solid particle. The dotted straight line Ūt highlights the motion of the solid particle in the (−z)-direction.

The results for the average swimming velocities are summarized in Fig. 5. The parameter space of all possible constrained swim strokes lies in the (F₂,α₂)-plane. In the contour plots, the maxima are marked with black dots. Slices along parts of the black lines are shown in Fig. 6. The white line indicates zero propulsion velocity.

Fig. 5 Contour plots of swimming velocity Ū over the entire parameter space (F₂, α₂) at λ = 1: (a) vesicle: the maximum vesicle speed (black dot) is 0.372 at F₂/F_2,max = 0.6224, (b) droplet: the maximum drop speed is 1.001 at F₂/F_2,max = 0.5081, (c) solid particle: the maximum particle speed is 0.113 at F₂/F_2,max = 0.4089. The slices marked by (black) solid lines are shown in Fig. 6.

Fig. 6 Swimming speed along the solid lines shown in Fig. 5 near the zeros of |U|. (a) At fixed α₂/2π = 0.91, (b) at fixed F₂/F_2,max = 0.95.

The general form of the dependence of Ū on the parameters F₂ and α₂ is the same for all three types of particles and is determined by eqn (58). The maxima of Ū are always at α₂ = π and the patterns are symmetric with respect to the line α₂ = π. The strokes leading to maximum swimming speed differ for the three types of particles, the F₂ value being largest for vesicles and smallest for solid particles. Note that the swimming velocity reverses its direction across the white line. This line is thus a one-dimensional manifold of swim strokes, which do not propel the particle at all.

If each active swimmer is free to select the swim stroke that maximizes its velocity, the droplet wins for all λ, followed by the vesicle and the solid particle in the third position. Note that the maximum speed of the solid particle is a factor ≈10 smaller than that of the drop. However, if the three swimmers perform the same stroke, their order of ranking can change with the parameters. This can be seen in Fig. 6, which shows the swimming speed |U| versus F₂/F_2,max (a) and versus α₂ (b) for λ = 1 along the cuts shown in Fig. 5 in the vicinity of the zeros of Ū. At each crossing of two lines in the graph, a pair of particles changes their order in the speed ranking.

If λ is increased, the average speed of a droplet will decrease. Nevertheless, the maximum speed always remains higher than the maximum speed of the vesicle and the solid particle, as shown in Fig. 7.

Fig. 7 Maximum speed vs. λ. The ranking (drop fastest, followed by vesicle, followed by solid particle) remains for all values of λ.

5.5 Lighthill efficiency with area constraint

A given swim stroke, parametrized by (F₂,α₂), propels different particles with different velocities, but also needs different amounts of dissipated power. A frequently used measure of the energetic efficiency of the particle is due to Lighthill.²⁵ It is defined as the ratio of the power dissipated by a solid sphere (with the same volume as the particle) moving at the particle's mean swimming velocity Ū to the power dissipated by the particle itself, that is,

In Fig. 8, we show contour plots of the Lighthill efficiencies for a vesicle, a droplet and a solid particle at λ = 1. The lines of vanishing Ū in Fig. 5 show up here as lines of zero efficiency. In the vicinity of these zeros, the efficiencies change their order as shown in Fig. 9 for a drop and a vesicle. The efficiency of the solid particle does not vanish at the zeros of drop and vesicle but is too small to be visible on the scale of the plots. However, it is typically two orders of magnitude smaller than the efficiencies of droplets and vesicles. The small ε can already be inferred from the Fig. 3 and 5, which show that the maximum speeds of the drop and the solid particle differ by a factor of ≈10, while their dissipation coefficients D_l differ only by factors between 1 and ≈2.

Fig. 8 Contour plots of Lighthill efficiency over the entire parameter space (F₂, α₂) at λ = 1 for (a) a vesicle, (b) a droplet, and (c) a solid. Dissipation from the ambient fluid and the interior (none in case of the solid) is taken into account. The maximum efficiency of a vesicle is 2.565 at F₂/F_2,max = 0.217, the maximum efficiency of a droplet is 3.200 at F₂/F_2,max = 0.364 and the maximum efficiency of a solid is 0.0280 at F₂/F_2,max = 0.195.

Fig. 9 Lighthill efficiencies along the dashed lines of Fig. 8 near the zeros of ε. (a) At fixed α₂/2π = 0.91, (b) at fixed F₂/F_2,max = 0.95. The insets show the full range of ε.

For small internal viscosities, ε_max of the droplet is higher than that of the vesicle, but it decreases with increasing λ as shown in Fig. 10. When λ ≈ 1.35 the maximum efficiencies become equal, and with more viscous internal fluids the vesicle is the most efficient, although it moves slower than the droplet (as can be seen from Fig. 7). Consequently, optimizing for speed and optimizing for efficiency yield different outcomes.

Fig. 10 Total maximum efficiency vs. λ.

6 Summary

In this work, we have presented a unified discussion of swimming velocities, hydrodynamic dissipation and efficiencies of three types of near-spherical amoeboid microswimmers, driven by axially symmetric achiral swim strokes consisting of radial deformations: a solid deformable body, a vesicle with incompressible fluid membrane, and a droplet. The intention of the present work is to study universal properties of this mode of swimming, which can be found in microorganisms as well as in artificial microswimmers, irrespective of the physical mechanisms generating the deformations. Our approach uses (scalar) spherical harmonics to represent surface deformations and a system of general solutions of the Stokes equation based on vector spherical harmonics. The swimming velocity and the hydrodynamic dissipation have been calculated to second order of small deformation amplitudes and their time derivatives. In our discussion, there is no need to consider the details of the physical mechanism generating the deformation if one is not interested in the additional dissipation that is caused by the mechanism. The calculated hydrodynamic dissipation is always present and represents a lower bound to the total dissipation. The restriction to axial symmetry is sufficient to compare translational motion between the three types of swimmer. The minimal models describe the types of swimmers by appropriate boundary conditions. For solids and vesicles, these conditions couple the interior material to the surrounding Newtonian fluid only via surface deformations, whereas for a droplet, the interior and the ambient flows are additionally coupled by the condition of continuity of tangential stress at the interface. As a consequence, the velocities of solids and vesicles do not depend on the viscosity ratio λ, in contrast to droplets.

To compare all three types of microswimmer, the swim strokes have to obey both the constraints of volume- and of surface-incompressibility. First, we studied swim strokes, which consist of pairs f_l, f_l+1 of deformations, varying periodically in time. In this case, the average swimming velocity does not depend on the detailed time course of deformations. In general, all |Ū_l| of droplets resulting from f_l, f_l+1 decrease with increasing λ. For small λ, the swimming speeds also decrease with increasing l, while for λ ≳ 1.3, they begin to increase with increasing l. Note also that unlike a vesicle the values of |Ū_l| for a solid body increase for l = 2, 3 before they start to decrease significantly from l = 5. In a next step, we have constructed a two-dimensional manifold of simple harmonic strokes, which contains l modes up to l = 4 and can be parametrized by the l = 2 parameters F₂ and α₂. This allows us to present the results for the swimming velocities U and the efficiencies ε in a clear and complete form as contour plots.

From these plots, the results of races between the three microswimmers are easily obtained. When each swimmer can choose the stroke that maximizes its speed, the droplet always comes first regardless of λ, although its velocity decreases as λ increases. The vesicle comes second, while the solid particle comes third. However, if the three swimmers perform the same stroke, their order of ranking can change with the parameters.

Optimizing the Lighthill efficiency and optimizing the swimming velocity result in different optimal swim strokes and rankings. The maximum total efficiency (based on the dissipation in the interior and surrounding flow) of a droplet is greater than that of a vesicle only if the dissipation ratio λ is small. Beyond λ ≈ 1.35 the maximum efficiency of the droplet falls below that of the vesicle. In all cases, the efficiency of the solid is typically two orders of magnitude smaller than that of vesicles and droplets, except in the vicinity of zeros of the efficiency. The efficiency of the solid particle will drop further if internal viscosity mechanisms are included in the model of the solid.

Experimental studies of the dependence of propulsion velocities on the time sequence of deformations require the control over particle shapes. For clean and surfactant-covered droplets, this becomes possible through the application of time-dependent electric fields.⁵⁰ For a typical experiment, a glass-walled box with electrodes embedded on opposite sides is filled with the ambient liquid. A comparison with our model, which restricts the forces at the interface to radial directions, requires negligible electric currents in the ambient fluid because these would generate tangential electro-hydrodynamical tractions. Droplets are dispersed using a micro-pipette, high voltage is applied via the electrodes, and the droplet dynamics is captured by cameras. To compare our results of hydrodynamic dissipation and efficiencies with experimental results, one must be able to split the total dissipation into the internal and hydrodynamic parts or provide a driving mechanism with negligible internal dissipation. For magnetic field-driven microswimmers, the Lighthill efficiency has been directly measured.⁵¹

Our work can be extended in several directions. First, the restriction to axially symmetric swim strokes can be lifted by using deformation amplitudes f(t,θ,ϕ), which depend on both polar angle θ and azimuthal angle ϕ, and expanding f in spherical harmonics. The calculational strategies presented above can be used if the set of vector spherical harmonics is enlarged to include rotational motion.²⁰ Furthermore, it is possible to study deformations other than radial (such as normal or tangential) by starting from appropriate variants of eqn (2). This brings our model closer to physical mechanisms like Marangoni flows or deformations by an active cytosceleton, which include tangential deformations. Second, the complete internal and ambient flow can be calculated, by projecting the boundary conditions onto higher l-components. Third, surface or body forces (for example Marangoni forces), which generate the deformations f, can be included. This makes the model more difficult because the relation between forces and deformations becomes a differential equation and introduces new relaxation time scales. For vesicles with well-separated relaxation time scales, an approach to this problem can be found in ref. 36. Finally, the internal materials of the particles may be changed to complex (for example visco-elastic) fluids or elastic solids.

Author contributions

Both authors contributed equally.

Data availability

All data is contained in the article.

Conflicts of interest

There are no conflicts to declare.

Appendices

A First-order solution for the droplet

We insert the expansions introduced in eqn (16), (17), (40) and (41) in the first-order boundary conditions eqn (42)–(44), use eqn (68) for the stress, and project on fixed l components. The resulting 4 × 4 linear system takes on the form

−(l + 1)(a⁽¹⁾_l − b⁽¹⁾_l) = ḟ_l, (60)

l(c⁽¹⁾_l + d⁽¹⁾_l) = ḟ_l, (61)

The four equations are easily reduced to two, because interior and exterior of the droplet are decoupled in the first two of them, implying

Substituting these expressions into eqn (62) yields the following equation:

Together with eqn (63) this fixes the coefficients c⁽¹⁾_l and a⁽¹⁾_l:

The remaining coefficients b⁽¹⁾_l and d⁽¹⁾_l are obtained from eqn (64).

B Stress components

The first-order stress components for the exterior flow, which are used in the boundary conditions of the vesicle and the droplet are calculated from

Substituting the flow velocity from eqn (16) and (17), we find for the off-diagonal component:

Concerning diagonal components, we only need the deviatoric stress Δσ = σ_rr − σ_θθ, which is independent of the pressure. To first order it needs to be evaluated at r = 1 and we get

eqn (68) apply equally well to the interior flow field and yield

for the corresponding stress components.

C Integrals

In this appendix, we summarise the integrals over products of associated Legendre polynomials which are needed in the main text, when projecting onto P₁:

D Droplet 2nd order

We start from the four eqn (45), multiply the first two by P₁ and the remaining 2 by P¹₁, and integrate over z = cos θ to project onto the l = 1 mode. The resulting 4 × 4 system of linear equations takes on the form

−4a − 2U = 3I_r⁺, (77)

2c + 2d − 2U = 3I_r⁻, (78)

4a − 4c − 8d = 3I_θ, (79)

−8η⁺a − 4η⁻d = I_σ, (80)

with

Here, v⁺ = v and v⁻ = V. The term Δσ represents the deviatoric stress, defined as σ_rr − σ_θθ, while the angular brackets signify the discontinuity across the surface of the unit sphere, that is, [A(r,θ)] = A(1 − ε,θ) − A(1 + ε,θ) for ε → 0₊. The velocity and the stress fields in the inhomogeneities are obtained from the first-order solutions given by eqn (64), (66) and (67) and (69)–(72). The swimming velocity U is obtained from the system (77)–(80) in the form

To evaluate the four inhomogeneities I^±_r, I_θ, I_σ we insert the expansions eqn (16), (17), (40) and (41) and (69)–(72) in eqn (81)–(83) and perform the integrations with the help of eqn (73)–(76). As an example, consider I^±_r. It consists of a double sum of terms arising from , from and from . It takes on the form

We insert eqn (73) and (74) and perform the summation over l′. The remaining l-series can be written in the form

with

B_l = l(l + 2)A_l. (88)

The other terms are treated in the same way. I_θ contains terms and gives the following result:

The inhomogeneity I_σ contains terms from f∂_rσ_θr. The contributions from f′Δσ produce terms (denoted by f′Δσ⁽⁰⁾) and terms (denoted by f′Δσ⁽²⁾), which can be read from eqn (70) and (72). The resulting l-series takes on the form

All first-order quantities are proportional to ḟ_l, so that the result for U can be written in the form of eqn (29) in the main text.

The explicit calculation of the terms R^drop_l, S^drop_l and Ū_l is straightforward but tedious, if not done with the help of computer algebra. The polynomials r⁽ⁿ⁾_l(λ), s⁽ⁿ⁾_l(λ) in eqn (46) and (47) take on the form

and

These expressions complete the analytical solution of the drop velocity.

E Energy dissipation

For the outer flow fields eqn (16) and (17) and the inner flow fields eqn (40) and (41), the integral eqn (50) can be evaluated by inserting the first-order solutions. For outer flow one gets

The energy dissipation for the inner space can equally well be computed and results in

The flow inside the vesicle can be calculated in analogy to the outside flow and yields

After inserting the first-order solutions into these expressions, one gets the eqn (51)–(54). For the droplet, the dissipation coefficients take on the explicit form

and

F Swim strokes with fixed surface area

With the designationwe find constraint equations from eqn (15) and (57), which take on the simple form

We restrict the deformations to the three l components l = 2, 3 and 4. For F₄ ≠ 0, we eliminate g₄ = 1 − g₂ − g₃ in eqn (102) and write the remaining equations in the form

2ĝ₃ sin α₃ = 1 − 2g₂ sin² α₂, (103)

2ĝ₃ cos α₃ = −2g₂ sin α₂ cos α₂, (104)

with ĝ₃ = g₃ sin α₃. First, we square and then add the two eqn (103) and (104), and so we can express ĝ₃ in terms of g₂ and α₂. Then we can obtain sin α₃ and cos α₃ from eqn (103) and (104) and finally, g₃ = ĝ₃/sin α₃. Note that the solutions of eqn (103) and (104) remain unchanged for α₂ → α₂ + π. They can be written in the formwith C = 1 + 4(g₂² − g₂)sin² α₂. The phase α₃ is determined only up to an angle π, due to the two roots ±C^1/2. In the following, we discuss only the +-branch. Changing α₃ → α₃ + π changes the sign of the swimming velocity U. Thus, it should always be remembered that for every motion discussed in the following, there is a corresponding motion in the opposite direction. The values of g₂ must be restricted to 0 < g₂ < 1/2 (corresponding to ) to ensure that g₃ and g₄ are positive and |cos α₃| < 1.

The special case F₄ = 0 has only 3 parameters F₂, F₃, α₂. In this case g₂, g₃ and α₂ are completely determined to g₂ = g₃ = 1/2, α₂ = π/2, 3π/2, and no adjustable parameters remain.

References

1. E. K. Paluch, I. M. Aspalter and M. Sixt, Annu. Rev. Cell Dev. Biol., 2016, 32, 469.
2. H. G. Othmer, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2019, 9, 1.
3. M. Abercrombie, J. E. Heaysman and S. M. Pegrum, Exp. Cell Res., 1970, 62, 389.
4. M. Bergert, A. Erzberger, R. A. Desai, I. M. Aspalter, A. C. Oates, G. Charras, G. Salbreux and E. K. Paluch, Nat. Cell Biol., 2015, 17, 524.
5. H. Ebata, Y. Nishigami, H. Fujiwara, S. Kidoaki and M. Ichikawa, Amoeboid movement utilizes the shape coupled bifurcation of an active droplet to boost ballistic motion, arXiv, 2024, preprint, arXiv:2402.03759 DOI:10.48550/arXiv.2402.03759.
6. N. P. Barry and M. S. Bretscher, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 11376.
7. A. J. Bae and E. Bodenschatz, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, E165.
8. M. Arroyo, L. Heltai, D. Millan and A. DeSimone, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 17874.
9. S. Michelin, Annu. Rev. Fluid Mech., 2023, 55, 77–101.
10. H. Boukellal, O. Campas, J.-F. Joanny, J. Prost and C. Sykes, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2004, 69, 061906.
11. E. A. Shah and K. Keren, eLife, 2014, 3, e01433.
12. M. Weiss, J. Frohnmayer, L. Benk, B. Haller, J. J. T. Heitkamp, M. Börsch, R. Lira, R. Dimova, R. Lipowsky, E. Bodenschatz, J. Baret, T. Vidakovic-Koch, K. Sundmacher, I. Platzman and J. Spatz, Nat. Mater., 2017, 17, 89–96.
13. S. Palagi and P. Fischer, Nat. Rev. Mater., 2018, 3, 113–124.
14. A. Mateos-Maroto, F. Ortega, R. G. Rubio, J.-F. Berret and F. Martínez-Pedrero, Particle, 2019, 36, 1900239.
15. C. Abaurrea-Velasco, T. Auth and G. Gompper, New J. Phys., 2019, 21, 123024.
16. H. R. Vutukuri, M. Hoore, C. Abaurrea-Velasco, L. van Buren, A. Dutto, T. Auth, D. A. Fedosov, G. Gompper and J. Vermant, Nature, 2020, 586, 52–56.
17. R. Golestanian, Active Matter and Nonequilibrium Statistical Physics: Lecture Notes of the Les Houches Summer School, Oxford University Press, 2022, ch. 8, vol. 112.
18. A. Callen-Jones, Front. Cell Dev. Biol., 2022, 10, 1000071.
19. R. Kree and A. Zippelius, Eur. Phys. J. E:Soft Matter Biol. Phys., 2021, 44, 6.
20. R. Kree and A. Zippelius, Eur. Phys. J. E:Soft Matter Biol. Phys., 2022, 45, 15.
21. R. Kree, L. Rückert and A. Zippelius, Phys. Rev. Fluids, 2021, 6, 034201.
22. A. R. Sprenger, V. A. Shaik, A. M. Ardekani, M. Lisicki, A. J. T. M. Mathijssen, F. Guzmán-Lastra, A. M. M. Hartmut Löwen and A. Daddi-Moussa-Ider, Eur. Phys. J. E: Soft Matter Biol. Phys., 2020, 43, 58.
23. G. Taylor, Proc. R. Soc. London, Ser. A, 1966, 291, 159–166.
24. G. Taylor, Proc. R. Soc. London, Ser. A, 1951, 209, 447.
25. M. J. Lighthill, Commun. Pure Appl. Math., 1952, 5, 109.
26. J. R. Blake, J. Fluid Mech., 1971, 46, 199.
27. E. Lauga and T. Powers, Rep. Prog. Phys., 2009, 72, 096601.
28. O. S. Pak and E. Lauga, J. Eng. Math., 2014, 88, 1–28.
29. A. Shapere and F. Wilczek, J. Fluid Mech., 1989, 198, 557.
30. B. U. Felderhof and R. B. Jones, Phys. A, 1994, 202, 94.
31. B. U. Felderhof and R. B. Jones, Phys. A, 1994, 202, 119.
32. A. Shapere and F. Wilczek, J. Fluid Mech., 1989, 198, 587.
33. B. U. Felderhof and R. B. Jones, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2014, 90, 023008.
34. Q. Wang and H. G. Othmer, J. Math. Biol., 2018, 76, 1699.
35. R. Kree and A. Zippelius, Phys. Fluids, 2023, 35, 042002.
36. A. Farutin, S. Rafai, D. K. Dysthe, A. P. Duperrey, P. Peyla and C. Misbah, Phys. Rev. Lett., 2013, 111, 228102.
37. M. Morozov and S. Michelin, J. Fluid Mech., 2018, 860, 711–738.
38. N. Yoshinaga, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2014, 89, 012913.
39. S. Torza, R. Cox and S. Mason, Proc. R. Soc. London, Ser. A, 1971, 269, 295–319.
40. L. Giomi and A. DeSimone, Phys. Rev. Lett., 2014, 112, 147802.
41. T. Ohta, J. Phys. Soc. Jpn., 2017, 86, 072001.
42. E. J. Campbell and P. Bagchi, Phys. Fluids, 2017, 29, 101902.
43. Q. Wang and H. G. Othmer, J. Math. Biol., 2016, 72, 1893.
44. A. Moure and H. Gomez, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2016, 94, 042423.
45. G. I. Taylor, Proc. R. Soc. London, Ser. A, 1934, 146, 501–523.
46. U. Seifert, Eur. Phys. J. B, 1999, 8, 405–415.
47. J. M. Rallison, J. Fluid Mech., 1980, 98, 625–633.
48. S. Kim and S. J. Karrila, Microhydrodynamics: Principles and Selected Applications, Dover Publ., 2005.
49. D. Das and D. Saintillan, J. Fluid Mech., 2017, 810, 225–253.
50. M. S. Abbasi, H. Farooq, H. Ali, A. H. Kazim, R. Nazir, A. Shabbir, S. Cho, R. Song and J. Lee, Materials, 2020, 13, 2984.
51. C. Calero, J. García-Torres, A. Ortiz-Ambriz, F. Sagués, I. Pagonabarraga and P. Tierno, Nanoscale, 2019, 11, 18723–18729.

---