Batchelor Prize Lecture Fluid dynamics at the scale of the cell

Abstract

The world of cellular biology provides us with many fascinating fluid dynamical phenomena that lie at the heart of physiology, development, evolution and ecology. Advances in imaging, micromanipulation and microfluidics over the past decade have made possible high-precision measurements of such flows, providing tests of microhydrodynamic theories and revealing a wealth of new phenomena calling out for explanation. Here I summarize progress in four areas within the field of ‘active matter’: cytoplasmic streaming in plant cells, synchronization of eukaryotic flagella, interactions between swimming cells and surfaces and collective behaviour in suspensions of microswimmers. Throughout, I emphasize open problems in which fluid dynamical methods are key ingredients in an interdisciplinary approach to the mysteries of life.

1 Introduction

The noted baseball player and amateur philosopher Yogi Berra is alleged to have said ‘you can observe a lot by watching’ (Berra 2010). Perhaps nowhere is this more true than in cellular biophysics, where with appropriate tools one can witness the noisy orchestrated to-and-fro that is life at the smallest scales. The existence of complex fluid flows around, and even inside, individual cells was recognized in some of the earliest microscopic studies of the living world. From the letters that Antoni van Leeuwenhoek wrote to the Royal Society in the late 17th and early 18th centuries (Hoole 1800) we learned about swimming spermatozoa, bacteria and protists, and in Bonaventura Corti’s celebrated treatise (Corti 1774) we find the first description of the persistent fluid flow inside large eukaryotic cells – the phenomenon we now call ‘cytoplasmic streaming’.

While these discoveries were made over two centuries ago, it is only in the last decade that there have been sufficiently precise measurements of the flows involved that it is possible to test fluid dynamical theories of their origins. As is often the case in science, it is not so much the sudden appearance of a single new technology that has made such developments possible, but rather it is the gradual development of a range of experimental techniques, and more importantly their synergistic use informed by theory, that is key. Such is the case in the work I will describe in this paper, in which a combination of micromanipulation, high-speed imaging, image processing and tracking methods, along with applications of ideas from statistical physics and nonlinear dynamics, has led to advances in our understanding of fluid dynamics at the scale of the cell.

A word or two is in order about scientific motivation. The work described below falls within the rapidly expanding field of ‘active matter’, the name given to systems in which energy in injected at the smallest scales (e.g. through biological processes involving molecular motors, biofilaments, cells, etc.) and produces larger-scale collective behaviour, most often in a fluid medium (Ramaswamy 2010). These systems raise many fundamental questions in non-equilibrium statistical mechanics and a significant body of work has addressed such physical issues as the implications of a lack of reversibility between microscopic transformations (Cates 2012), the multipolar decomposition of flow fields around microswimmers (Ghose & Adhikari 2014), the interpretation of the pressure field (Takatori & Brady 2014) and the rheological properties of the cytoplasm (Moeendarbary et al. 2013).

Here I want to emphasize a different perspective, one that is grounded in issues that arise in evolutionary biology. We might say that if the mechanical approach is about ‘how’ then the biology questions are (ultimately) about ‘why’. Of course, these approaches are not in conflict, for even with the biology-based problems it is the quantitative description that we seek in the end. And often one has to answer ‘how’ before ‘why’ in any event. But the ‘why’ questions tend to be much richer, even if inherently harder to answer. The danger of developing ‘just-so’ stories like those of Kipling (1902) will be mitigated if we take a hypothesis-based approach.

As an example of a ‘why’ question, we can look at uni- and multicellular green algae ranging from Chlamydomonas to Volvox shown in figure 1 (Goldstein 2015). More will be said about them in subsequent sections, but here I note that all of these algae rotate slowly around a body-fixed axis as they swim through the action of multiple eukaryotic flagella. We can ask how this swimming motion is achieved by studying the low Reynolds number problem of actuated elastic filaments attached to a spherical body, a problem whose history dates back to work of Lighthill (1952) and Blake (1971), but we can also ask why the organisms have evolved to rotate in the first place. The answer to the latter question turns out to involve phototaxis, directed movement toward light. Detailed study (Drescher, Goldstein & Tuval 2010b ) shows that there is a rather precise tuning of the organism’s rotational period and the internal relaxation time of the flagellar response to the periodic illumination by light of a directional photosensor on the cell body. As these organisms do not have a central nervous system, this tuning turns out to be a mechanism to achieve organism-scale synchrony of otherwise autonomous units. Thus this particular ‘why’ question is really one of control theory embedded in a fluid dynamical context.

This example touches on one particular aspect of the overarching question that has guided our research for the past decade: the origins of multicellularity. This issue of why and how the simplest unicellular organisms evolved to become multicellular and to exhibit cellular differentiation is one of the most fundamental issues in biology (Maynard-Smith & Szathmáry 1995), and has been of particular interest since the work of the great biologist August Weismann (1892). When one asks about the driving forces behind increased size it is important to recognize (Short et al. 2006; Solari et al. 2006; Solari, Kessler & Goldstein 2007, 2013) that one of the most important features of life is the exchange of materials (nutrients and waste products) with the surroundings and the physical processes of diffusion, reaction, mixing, motility and buoyancy are all involved. So the issue of multicellularity is deeply connected to problems of fluid mechanics!

Figure 1. The Volvocine green algae. (a) Chlamydomonas reinhardtii, (b) Gonium pectorale, (c) Eudorina elegans, (d) Pleodorina californica, (e) Volvox carteri and (f) Volvox aureus. C. reinhardtii (g) and V. carteri (h) held on glass micropipettes. Adapted from Goldstein (2015).

Finally, if we are to study active matter in vivo then there is the question of the appropriate organisms to use. Here we learn from the biology community that often it is best to use what are termed ‘model organisms’, those with specific features which allow for a comprehensive approach, and often provide the basis for extrapolation to other organisms. In plant science one often studies the weed Arabidopsis thaliana, not because it is an important crop but because it has many features that make it amenable to a holistic study. It is easy to grow and modify by genetic engineering, there are many mutants available, etc. In the early 1960s Sidney Brenner proposed the nematode Caenorhabditis elegans as a model for animal development because of its simplicity, ease of growth and genetic analysis (Brenner 1988). By analogy, we now understand (Goldstein 2015) that the green algae are extremely well suited to study problems in biological fluid dynamics because of their ease of growth, ready availability in nature and understanding from diverse perspectives (Kirk 1997). One should also not overlook the role that geometry can play in making organisms suitable for study. Recall that the ‘giant squid axon‘ used by Hodgkin & Huxley (1952) is not from a giant squid but is itself very large (up to 1 mm in diameter), making it straightforward to introduce electrodes and micropipettes into the axonal interior. Likewise, two of the great advantages of the Volvocine algae is that they are large enough to be captured on micropipettes (figure 1 g,h) and they are to a good approximation spherical in shape and axisymmetric in flagellar actuation – the spherical (green) ‘cows’ of the microscopic world – so theory is that much easier.

In the following sections I focus on experimental and theoretical progress on four distinct aspects of active matter; cytoplasmic streaming in eukaryotic cells, synchronization of eukaryotic flagella, interactions of microswimmers with surfaces and collective behaviour in concentrated suspensions of swimmers. The discussion in each section is chronological, so that the logical progression of our investigations is made clear.

Figure 2. Discovery of cytoplasmic streaming. Illustrations from the monograph (a) by Bonaventura Corti (1774); (b) the entire aquatic plant Chara, (c) result of osmotically shocking the cell, in which the vacuole shrinks to reveal the vacuolar membrane (tonoplast) and (d) a depiction of the continuous circulation of the cytoplasm from one end of each internodal cell to the other and back. Images courtesy of the British Library. (e) Digital image of Chara.

2 Cytoplasmic streaming

In 1774 the Italian Bonaventura Corti, who was ‘Maestro’ to the great fluid dynamicist Giovanni Battista Venturi, published an account (Corti 1774) of the wonders he found examining certain aquatic plants under the microscope (figure 2). He studied the freshwater species Chara which can grow several tens of centimetres in height, anchored in the mud with rhizoids. The plant itself has a chain of cylindrical cells each of which can be up to 10 cm long, with branches sprouting at the nodes. These are among the largest single cells to be found in nature. Corti discovered that there is a persistent circulation of the fluid contents of the cells, especially within the internodal ones. Rediscovered several decades later by Treviranus (1811), but properly attributed to Corti by Venturi in his historical work on optics (Venturi 1814), the phenomenon was rediscovered several other times over the centuries (Hughes 1959). It was only in the 1960s that the underlying molecular mechanism that drives streaming was understood (Shimmen 2007). That process is the motion of multitudes of cargo-carrying motor proteins moving along filamentary tracks arrayed along the interior cell wall, entraining fluid as they move. In plants the motors generally belong to a class of myosins moving along bundled actin filaments, while in animals streaming is frequently driven by kinesins moving along microtubules.

Plant myosins are among the fastest known in nature, and the myosin XI found in Chara can move at speeds $U\sim 100~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$ (Shimmen & Yokota 2004), more than an order of magnitude faster than kinesins. The fluid flows created by these moving proteins can be as fast as the motors, and with the characteristic size $L$ of these cells, say the diameter, approaching 1 mm, and $D\sim 10^{-5}~\text{cm}^{2}~\text{s}^{-1}$ a molecular-scale diffusion constant, we see that the Péclet number $Pe=UL/D$ can be ${\sim}10^{2}$ and considerably larger for higher molecular weight species. This is ‘life at high Péclet number’ (but still ‘life at low Reynolds number’, as so elegantly elucidated by Purcell (1977)). If we look at longitudinal motion along the internodal cell, now with $L\sim 10$  cm, the diffusion time scale $L^{2}/D$ would be prohibitively long ( $10^{7}$  s or 4 months), whereas the advective time $L/U\sim 10^{3}$  s or 16 min. This need for active transport when cellular sizes exceed even a few mm was explained eloquently by Haldane (1926) in his famous essay ‘On Being the Right Size’.

Streaming is now known to take place in a wide variety of organisms, from aquatic and terrestrial plants to animals (Verchot-Lubicz & Goldstein 2010; Goldstein & van de Meent 2015), particularly in certain stages of embryonic development. In the latter it has been of recent interest in the process of oogenesis in Drosophila, where the cytoplasmic flows can take on a variety of different forms as development proceeds, and modern fluid mechanical measurement methods are now being used (Ganguly et al. 2012; Khuc Trong et al. 2015). The one overarching question in all of these systems is: what is the biological function of cytoplasmic streaming, i.e. why does the organism do it? One extreme view is that the flows per se have no useful function, but are simply an inevitable by-product of nanoscale transport (Pickard 1974). Counter to this is the view that streaming can promote mixing in the cell and thereby help maintain homeostasis (Hochachka 1999). This debate prompted us to investigate the flows quantitatively and to consider their implications for transport and metabolism within cells.

Figure 3. Structure of internodal cells of Chara corallina and flow measurements. (a) The entire plant. (b) Bidirectional helical boundary conditions on the flow distinguish two indifferent zones ( $\text{IZ}_{\pm }$ ). (c) Cross-section showing theoretical longitudinal velocity profile (colours) and in-plane vortical flows. (d) Schematic cross-section of internodal cell. (e) Experimental set-up for magnetic resonance velocimetry (MRV). (f) MRV measurements of longitudinal flow in three successive regions along cell. (g) Average of the measurements in (f), rotated so indifferent zones coincide. (h) Theoretical velocity profile corresponding to (g). Reproduced with permission from van de Meent, Tuval & Goldstein (2008) and van de Meent et al. (2010).

The internodal cells of Chara have often been used in plant science for studies of electrophysiology. Its microstructure is shown in figure 3(d); bundled actin filaments attached to chloroplasts at the inside of the cell spiral along the cell in helical paths. A thin ( $10{-}20~\unicode[STIX]{x03BC}\text{m}$ ) layer of cytoplasm with a viscosity several times that of water fills the space between the inner cell wall and the vacuolar membrane (also called the tonoplast), inside of which is the vacuolar fluid, water containing many dissolved chemical species. The myosin molecules move along actin filaments in a particular direction (‘minus’ to ‘plus’ ends, as distinguished by the protein microstructure), and the filamentary arrangement is such that motors walk in opposite directions on two twisted hemicylinders, giving rise to a ‘barber-pole’ pattern of motor motion, and hence of fluid flow in the cytoplasm. Each junction between the two hemicylinders is termed an ‘indifferent zone’ (IZ) and is characterized by high shear; the streaming speed changes by up to $200~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$ over a distance of two chloroplasts ( ${\sim}20~\unicode[STIX]{x03BC}\text{m}$ ), so the shear rate $\dot{\unicode[STIX]{x1D6FE}}\sim 10~\text{s}^{-1}$ . This is pipe flow driven not by a pressure gradient but by a distribution of forces or velocities at the wall, not unlike the ‘squirmer’ model (Lighthill 1952) of microorganisms in which the action of beating flagella is recast as a prescribed velocity field on the organism’s surface.

The simplest model of the flow is one in which the helical winding of the indifferent zones is ignored and the velocity is prescribed at the interior cell wall to be piecewise constant ( $+U$ on one hemicylinder, $-U$ on the other), and the thin cytoplasm and tonoplast are ignored. Pickard (1972) found the solution to the Stokes equation in this case, in the absence of any external forcing, as a longitudinal velocity field

(2.1) $$\begin{eqnarray}u(r,\unicode[STIX]{x1D703})=\frac{2U}{\unicode[STIX]{x03C0}}\arctan \left(\frac{2(r/R)\sin \unicode[STIX]{x1D703}}{1-(r/R)^{2}}\right),\end{eqnarray}$$

where $R$ is the cylinder radius and the IZs are at $\unicode[STIX]{x1D703}=0$ and $\unicode[STIX]{x03C0}$ . This field has integrable singularities at the two indifferent zones. Under the assumptions (i) that the dominant contribution to the viscous dissipation rate per unit length inside the cell is localized at the indifferent zones, and therefore scales as $U^{2}$ independent of $R$ and (ii) the power input by the molecular motors scales with the cellular circumference, and hence as $R$ , one obtains the scaling $U\sim R^{1/2}$ that has been observed (Pickard 1974).

One might imagine that generalizing this calculation to the case of helical indifferent zones would lead only to a small quantitative change in the basic flow properties, but something more interesting happens because the cell is chiral (Goldstein, Tuval & van de Meent 2008; van de Meent et al. 2008); the lack of mirror symmetry leads to a distinction between the two indifferent zones, and in a plane orthogonal to the long axis of the cell there is a pair of vortices (figure 3 c). Although these bear a superficial resemblance to Dean vortices, which appear in curved pipes at finite $Re$ , these vortices are present when $Re=0$ . The effect of the transverse flow associated with these vortices is to create a kind of micromixer (van de Meent et al. 2008) in which fluid particles are shuttled from one indifferent zone to another as they progress along the cellular axis. This phenomenon is similar to that which occurs in the microfluidic ‘herringbone micromixer’ (Stroock et al. 2002) in which diagonal grooves in the channel deflect the flow from its primary downstream direction, but in the case of Chara the flow is not chaotic. Computations (van de Meent et al. 2008) show that there can be enhanced uptake of nutrients from the surrounding medium and more rapid mixing within the vacuole, providing support for the notion mentioned above (Hochachka 1999) that homeostasis is enhanced by streaming.

Early measurements (Kamiya & Kuroda 1956; Mustacich & Ware 1976) of the cross-sectional flow using endogenous particles as tracers and a laser Doppler technique revealed an S-shaped longitudinal velocity profile along one particular diameter, interpolating between the upward and downward wall-driven flows. But these were not capable of revealing the entire flow profile and testing even the most basic aspects of the global flow, such as whether the idea of piecewise constant forcing was a sensible result from multitudes of stochastic molecular motors moving along an array of filaments. To our surprise, we learned that the technique of magnetic resonance velocimetry could achieve sufficient spatial resolution ( ${\sim}10~\unicode[STIX]{x03BC}\text{m}$ ) to attempt such a measurement, provided the flow was steady on the time scale needed for averaging (several hours), which it is. The result of these measurements (van de Meent et al. 2010) is shown in figure 3(f–h), where we see striking agreement with a theory that incorporates piecewise constant helical forcing at the wall and complete shear transmission by the tonoplast. This is perhaps the first direct measurement of the velocity field inside a single living cell (albeit an enormous one).

These measurements and the earlier ones offer no indication that the tonoplast plays any role in hindering the transmission of shear from within the cytoplasm into the vacuolar fluid. Pickard (1972) wondered how this could be, saying ‘What manner of membrane is this?’ If we take a larger view of the internodal cell we recognize, as Corti saw, that the vacuolar membrane has the topology of a sphere, and in that sense comprises one enormous lipid vesicle. The so-called giant (synthetic) lipid vesicles studied in the laboratory are perhaps at most $100~\unicode[STIX]{x03BC}\text{m}$ in diameter, but here we have one that is 10 cm long and 1 mm in diameter! A very clever experiment (Vézy, Massiera & Viallat 2007) on lipid vesicles adhered to the bottom of a microfluidic chamber showed the appearance of a vortical flow pattern induced in the membrane by the action of a shear flow, and suggested to us that one could study the problem of shear transmission through the membrane by monitoring the entire flow fields in such a set-up. We used the same basic geometry as Vézy et al. (2007), shown in figure 4(a), but introduced microspheres into the fluid within and external to the vesicle and used a multicomponent lipid membrane in a regime of two-phase coexistence: either involving small domains of a ‘liquid disordered’ ( $L_{d}$ ) phase in a background of ‘liquid ordered’ ( $L_{o}$ ) phase or solid-like ‘gel’ domains in a background of $L_{d}$ phase (figure 4 b–d), made visible with fluorescent dyes. Tracking the microspheres and the domains under shear, it was possible, with confocal imaging, to map out stacks of two-dimensional slices of the flow field inside, outside and within these vesicles (Honerkamp-Smith et al. 2013), as seen in figure 4(e,f). Using the constraint of incompressibility of the fluid flow within the vesicle, the two-dimensional slices of the velocity field can be used to reconstruct the full three-dimensional flow within the vesicle, as shown in figure 5(a).

Figure 4. Microfluidic set-up to study a sheared vesicle. (a) A hemispherical vesicle (grey) is adhered to the base of a microfluidic chamber and subjected to shear flow. Confocal imaging of $L_{o}$ phase vesicle with $L_{d}$ domains (b) and tracking of (c) gel domains in $L_{d}$ background and (d) $L_{d}$ domains in $L_{o}$ background at gel apex. (e) Two-dimensional flow fields from particle image velocimetry (PIV) at heights $z/R=0.26,0.47,0.71$ above chamber base. (f) Confocal images at same heights, showing fluorescent microspheres (green) inside vesicle (red). (g) Theoretical velocity fields corresponding to positions shown in (e). Reproduced from Honerkamp-Smith et al. (2013).

Figure 5. Velocity fields inside and outside a sheared vesicle. (a) Experimental and (b) theoretical interior flows are in excellent agreement. Exterior flows (c) exhibit closed streamlines, consistent with the circulation of phase-separated domains seen from confocal imaging (d). Arrows indicate flow directions. Reproduced from Honerkamp-Smith et al. (2013).

A calculation (Woodhouse & Goldstein 2012a ) of the flows inside, within and outside a lipid vesicle adhered to a no-slip surface incorporates the key feature that a lipid membrane is laterally incompressible, so the flow patterns are different from those of a drop of one fluid surrounded by a second, immiscible one, where fluid can exchange freely between the surface and the bulk. The essential result in shown in figures 4(g) and 5(b), where we see excellent agreement with the experimental results. The constraint of lateral incompressibility of the membrane leads to a significant reduction in the membrane velocity relative to the value in the shear flow that would be there in the absence of the vesicle, even when the viscosity is negligible. But in addition, careful analysis of the radius dependence to that reduction allows for a direct determination of the membrane viscosity $\unicode[STIX]{x1D707}_{m}$ , the measurement of which has historically been based on tracking the fluctuating trajectories of embedded objects and subject to considerable uncertainty. The concept of a membrane viscosity can be understood from the essential force balance on a planar membrane $\unicode[STIX]{x1D6E4}$ through the two-dimensional Stokes equation for the membrane velocity $\boldsymbol{u}^{m}$ forced by the velocity gradients from either side,

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D707}_{m}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}^{m}+2[\unicode[STIX]{x1D707}_{+}\boldsymbol{e}_{\Vert }^{+}-\unicode[STIX]{x1D707}_{-}\boldsymbol{e}_{\Vert }^{-}]_{\unicode[STIX]{x1D6E4}}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F1},\end{eqnarray}$$

where $\boldsymbol{e}_{\Vert }^{\pm }$ are the in-plane normal rates of strain on either side of the membrane, $\unicode[STIX]{x1D707}_{\pm }$ the associated viscosities there and $\unicode[STIX]{x1D6F1}$ is the membrane pressure field. Experiments show membrane viscosities at most of order $\unicode[STIX]{x1D707}_{m}\sim 15~\unicode[STIX]{x03BC}\text{P}~\text{cm}$ , so if we take $\unicode[STIX]{x1D707}_{\pm }$ to be the viscosity $\unicode[STIX]{x1D707}_{w}$ of water, the Saffman–Delbrück length $\ell =\unicode[STIX]{x1D707}_{m}/\unicode[STIX]{x1D707}_{w}\sim 15~\unicode[STIX]{x03BC}\text{m}$ . On length scales large compared with $\ell$ the viscous dissipation is dominated by the surrounding fluid (Saffman & Delbrück 1975). This result suggests why the measured velocity profiles in Chara are so insensitive to the presence of the tonoplast.

One of the most interesting aspects of the cytoplasmic flow profile in the Characean algae is its bidirectional helical geometry. How does it get established in the developing organism? Tantalizing experimental evidence (Foissner & Wasteneys 2000) points to a self-organization process; when the actin and microtubule cytoskeletons are disrupted with chemical agents so that streaming stops, and then those agents are removed, streaming reconstitutes itself with a shifted IZ. As discussed further below in the context of a model for self-organization, an important ingredient in this process would appear to be flow-induced alignment and even distortion of the filaments along which motors move. Key insight into this kind of process was provided by the analysis of Young & Shelley (2007) of the buckling dynamics of an elastic filament at a hyperbolic stagnation point. As discussed by Batchelor (1970), if a filament of diameter $a$ and length $L$ (aspect ratio $\unicode[STIX]{x1D700}=a/L$ ) lies along the $x$ -axis, it experiences tension with the profile

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D70E}(x)=\frac{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}\dot{\unicode[STIX]{x1D6FE}}}{\ln (1/\unicode[STIX]{x1D700}^{2}e)}(L^{2}/4-x^{2}),\end{eqnarray}$$

vanishing at the free ends. If the filament is aligned to be under compression (with $\dot{\unicode[STIX]{x1D6FE}}<0$ ) it can deform through a process analogous to Euler buckling, but here with a distributed compressive force. In a long-wavelength approximation one can associate to the filament configuration $h(x)$ the energy functional

(2.4) $$\begin{eqnarray}\mathscr{E}[h]=\frac{1}{2}\int _{-L/2}^{L/2}\text{d}x\{Ah_{xx}^{2}+\unicode[STIX]{x1D70E}(x)h_{x}^{2}\},\end{eqnarray}$$

where $A$ is the filament bending modulus. Beyond a critical tension (or flow velocity) there is a bifurcation to a buckled state, which then proceeds to rotate and stretch out along the extensional direction. Experiments on macroscopic fibres (Wandersman et al. 2010) confirmed the details of this process, but left unexplored the important role of noise in a microscale realization of the instability. This experiment used the cross-channel microfluidic set-up in figure 6 within which were held individual fluorescently labelled actin filaments (Kantsler & Goldstein 2012). Figure 6(e–g) show the buckled filament shapes as the dimensionless control parameter

Figure 6. Buckling instability of single actin filaments in a cross-channel flow. (a) Microfluidic set-up, in which an actin filament is held at the intersection of rectangular channels by adjusting the pressure difference $\unicode[STIX]{x0394}P$ . (b) Velocity profile in a horizontal plane midway between bottom and top of chamber. (c) Contrast-enhanced image of an actin filament. (d) The function $h(s,t)$ describes the deviation from its mean orientation. (e)–(g) Montage of buckled filaments at $\unicode[STIX]{x1D6F4}=-0.55,-1.9,-47$ and signs of the shear rate, where $-(+)$ indicates compressional direction is along $x(y)$ . Times indicated are rescaled by $\dot{\unicode[STIX]{x1D6FE}}$ . (h) Bifurcation plot showing end-to-end distance $\mathscr{L}$ relative to the contour length $L$ as a function of $\unicode[STIX]{x1D6F4}$ . Red dashed lines shows theoretical bifurcation points for modes indicated, and red arrows show values of $|\unicode[STIX]{x1D6F4}|$ associated with images in (e)–(g). Grey band represents noise floor. Reproduced from Kantsler & Goldstein (2012).

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D6F4}=\frac{2\unicode[STIX]{x1D707}\dot{\unicode[STIX]{x1D6FE}}L^{4}}{\unicode[STIX]{x03C0}^{3}A\ln (1/\unicode[STIX]{x1D700}^{2}d)}\end{eqnarray}$$

is varied. We see clearly the first buckling mode and higher-order ones. A bifurcation plot (figure 6 h) using an order parameter based on the end-to-end distance $\mathscr{L}$ of the filament shows excellent agreement with the predicted threshold of the instability and clear indications that the transition is rounded by thermal fluctuations. This is to be expected because these experiments involve filaments for which $L$ is comparable to the persistence length $L_{p}=A/k_{B}T$ , so thermal energy can induce curvature ${\sim}L^{-1}$ . A theory of this rounded transition has recently been presented (Manikantan & Saintillan 2015).

Figure 7. Self-organisation of cytoplasmic streaming. Ingredients in mathematical model are: (a) preferential orientation away from direction $\boldsymbol{d}$ to minimize curvature energy, (b) alignment with flow is stable with $+$ end downstream, but unstable if it is upstream (c) and (d) filament friction against cell wall competes with drag from cargo vesicles. (e) Six panels show time evolution of streaming, where colour indicates longitudinal component of the vector $\boldsymbol{P}$ (purple for $P_{z}>0$ , green for $P_{z}<0$ , brightness increases with magnitude), and indifferent zones are white lines separating up- and downwardly directed streaming. Arrows indicate streamlines, with flow from the thin end to thick end. Reproduced from Woodhouse & Goldstein (2013).

Returning to the description of self-organisation of streaming, the experimental evidence points to several key processes at work. First, filaments tend to bundle together through the action of proteins, with a rate constant $\unicode[STIX]{x1D6FC}_{p}$ . Second, as the filaments are elastic they will tend to reorient away from regions of high cell wall curvature, as represented by a disfavoured direction $\boldsymbol{d}$ as in figure 7(a), here taken to be in the circumferential direction and with coupling $\unicode[STIX]{x1D705}$ . Third, while bare filaments will tend to orient nematically in the flow, the fact that motors walk along filaments in a polar manner implies that a filament oriented with the plus end downstream is stable in the presence of shear flow (figure 7 b), but it is unstable, for example to buckling, if oriented the opposite way (figure 7 c). Therefore there will be a tendency toward polar order, with coupling $\unicode[STIX]{x1D6FC}_{u}$ . Finally, a filament orthogonal to the flow will be reoriented by the flow through a hydrodynamic torque created in part by friction against the cell wall itself (figure 7 d), and this frictionally restricted advection is associated with a parameter $\unicode[STIX]{x1D716}$ .

A continuum description of the process of filament self-organisation in a thin shell near the cell wall (Woodhouse & Goldstein 2013) is based on the theory of ‘active suspensions’ (Saintillan & Shelley 2008b ; Woodhouse & Goldstein 2012b ), in which the fundamental quantity is the angular distribution function $\unicode[STIX]{x1D6F9}(\boldsymbol{x},\boldsymbol{p})$ of microfilaments with unit direction $\boldsymbol{p}$ . Moments of $\unicode[STIX]{x1D6F9}$ yield averages in the usual manner, so that the local filament concentration is $c(\boldsymbol{x})=\int _{\boldsymbol{p}}\unicode[STIX]{x1D6F9}(\boldsymbol{x})\,\text{d}\boldsymbol{p}$ and the average orientation director is $\boldsymbol{P}(\boldsymbol{x})=c(\boldsymbol{x})^{-1}\int _{\boldsymbol{p}}\boldsymbol{p}\unicode[STIX]{x1D6F9}(\boldsymbol{x})\,\text{d}\boldsymbol{p}$ . From the Smoluchowski equation obeyed by $\unicode[STIX]{x1D6F9}(\boldsymbol{x},\boldsymbol{p})$ one obtains dynamics of these moments under suitable closure approximations. Here, the concentration field can be taken to be uniform, leaving only the dynamics of the director field. In a form made dimensionless with the cell radius $R$ it reads

(2.6) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\boldsymbol{P}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D716}\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{P}=d^{(s)}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{P}-d^{(r)}\boldsymbol{P}+(\mathbb{I}-\boldsymbol{P}\boldsymbol{P})\boldsymbol{\cdot }[\unicode[STIX]{x1D716}(\unicode[STIX]{x1D735}\boldsymbol{u})\boldsymbol{\cdot }\boldsymbol{P}+\unicode[STIX]{x1D6FC}_{p}\boldsymbol{P}+\unicode[STIX]{x1D6FC}_{u}\boldsymbol{u}-\unicode[STIX]{x1D705}(\boldsymbol{P}\boldsymbol{\cdot }\boldsymbol{d})\boldsymbol{d}],\end{eqnarray}$$

where $d^{(s)}$ and $d^{(r)}$ are spatial and rotational diffusion constants. The fluid flow $\boldsymbol{u}$ is driven by motors walking along the filaments, corresponding to some force $\boldsymbol{F}$ ,

(2.7) $$\begin{eqnarray}-\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+\unicode[STIX]{x1D708}\boldsymbol{u}+\unicode[STIX]{x1D735}\unicode[STIX]{x1D6F1}=\boldsymbol{F},\end{eqnarray}$$

where $\unicode[STIX]{x1D708}$ represents the effects of wall friction in this thin-layer approximation. As bundling is known to be required for streaming to occur, the force $\boldsymbol{F}$ is nonlinear in the director field $\boldsymbol{P}$ and at lowest order is proportional to $|\boldsymbol{P}|^{2}\boldsymbol{P}$ . The model is closed by requiring incompressibility and vanishing net fluid flow in the longitudinal direction of a cylindrical shell computational domain to mimic the closed ends of the cell. The essential output of this model is indeed a process of self-organisation, as shown in figure 7. An initially disordered director field and induced flow evolve to become close to piecewise constant, pointing in opposite directions on the two hemicylinders. This process is driven by the synergy between the individual physical effects embodied in the model, each of which in isolation would not be sufficient to drive the self-organisation observed. The phenomenological model outlined assumes that filaments are non-self-advective, but a freely suspended filament with a cargo-carrying molecular motor acts on the fluid as a force dipole (a pusher stresslet) (Woodhouse & Goldstein 2012b ), and so it is not surprising that there can be interesting connections between cytoplasmic streaming and suspensions of microswimmers. This is discussed further in § 5.

A final point to mention is that there is a wealth of fascinating fluid–structure interactions in cytoplasmic streaming in higher organisms. A prime example is found in the oocytes of Drosophila, where molecular motors transport gene products to specific places in order to establish the body plan. Unlike in plants, the cytoskeleton is, during much of the process, highly disordered, as are the flows. As the cargo-carrying motors stochastically attach and detach from the filaments, diffusing and advected when detached, their motion along the network is analogous to passive scalar transport in turbulent flows (Ganguly et al. 2012) and is only beginning to be studied quantitatively (Ganguly et al. 2012; Khuc Trong, Guck & Goldstein 2012; Khuc Trong et al. 2015). In addition, at certain stages of oogenesis there are drastic rearrangements in the architecture of the microtubule cytoskeleton, presumably driven by back reaction of the motor-induced flows on the filaments themselves (Serbus et al. 2005). A quantitative understanding of this feedback is one of the important open problems in cytoskeletal dynamics. A second unsolved problem (in Chara) is understanding how a collection of cargo-carrying molecular motors moving stochastically along filaments along the cell periphery produce the very regular fluid flows observed, and the important role of the boundary conditions near those filaments (Wolff, Marenduzzo & Cates 2012).

3 Synchronisation of eukaryotic flagella

Since at least the early 20th century it has been known that groups of beating cilia and flagella can exhibit synchronous behaviour. This is particularly clear in the case of multiciliated tissues and organisms, where the phenomenon of ‘metachronism’ is seen: cilia beating at the same frequency but with a definite phase lag between neighbours, giving rise to patterns like the so-called Mexican wave in a stadium or, as Gray put it, ‘…the waves which pass over a field of corn when exposed to the breeze’ (Gray 1928). Synchronous beating is also observed in freely swimming cells such as spermatozoa, as in the studies of Rothschild (1949), and collectively these observations led Taylor (1951) to develop his ‘waving sheet’ model in which the problem of synchrony was abstracted to two laterally infinite sheets with a prescribed sinusoidal wave of deformation (amplitude $b$ , wavelength $\unicode[STIX]{x1D706}$ ) placed at some separation $h$ . Solving the zero Reynolds number fluid mechanical problem to second order in the assumed small parameter $b/\unicode[STIX]{x1D706}$ , Taylor found the intuitive result that the rate of energy dissipation was minimized when the two sheets execute in-phase motion, an effect increasing as $h$ decreases. Although the analysis does not yield an equation of motion for the phase shift that would allow an arbitrary initial condition to evolve to this state, it does suggest that if a principle of least dissipation holds then synchrony follows. Recent work (Elfring & Lauga 2011) has shown that relaxing the assumption of a fixed waveform leads to such a dynamical equation for which the synchronized state is an attractor.

The idea that hydrodynamic interactions could be at the root of flagellar synchrony has been a dominant theoretical idea ever since Taylor’s work, but for much of that time it remained untested experimentally. But starting in the mid-1980s studies of flagellar beating advanced significantly through the advent of high-speed imaging technology. One of the most important series of experimental studies was that of Rüffer & Nultsch (1985, 1987, 1998), who studied micropipette-held cells of Chlamydomonas using high-speed imaging and light-table tracings of each video frame to analyse time series containing perhaps 20–30 beats. Chlamydomonas typically displays a synchronized breast stroke, but it is known that the two flagella, termed cis and trans for their proximity to the eyespot (figure 8), have differential response to internal calcium levels and in isolation display distinct beating frequencies. Rüffer & Nultsch (1985) found that while the synchronized breast stroke was seen in 85 % of the cells examined, 10 % of the time series displayed a brief episode in which one flagellum transiently beat faster than the other, losing synchrony, only to recover it a few beats later. And 5 % of the time the two flagella displayed distinct beat frequencies that could differ by ${\sim}20\,\%$ . Their interpretation of these three modes was that they represented different subpopulations of cells.

Figure 8. The eukaryotic flagellum. (a) A Chlamydomonas cell body (green), eyespot (orange), basal bodies (B1 and B2) of the trans and cis flagella and filaments connecting them. (b) Cross-sectional structure of the flagellum, with the ‘ $9+2$ ’ arrangement of microtubule doublets. (c) Key molecular motors connecting to the microtubules. Reproduced from Wan, Leptos & Goldstein (2014).

Figure 9. Stochastic dynamics of Chlamydomonas flagella. (a,b) Top images are processed video frames during noisy synchronous beating (a) and drift (b). Middle image in (a) shows small interrogation windows for left (L) and right (R) flagella. Lower panels in (a,b) show time series of the signals $X_{L},X_{R}$ in those windows. (c) Phase difference $\unicode[STIX]{x1D6E5}$ versus time, during which appear periods of noisy synchrony, phase slips and drifts. Inset shows instantaneous frequencies of flagella before, during and after a drift period. Adapted from Polin et al. (2009).

Our initial investigations of this phenomenon (Goldstein, Polin & Tuval 2009; Polin et al. 2009) sought to utilize modern high-speed imaging technology, micropipette manipulation and dynamical systems methods not only to measure the synchronization dynamics with greater spatio-temporal precision, but also to interpret the resulting time series of beating in the context of emerging models of synchrony. While it is now routine to extract detailed beating waveforms from videos of single cells, it was conceptually simpler first to implement a Poincaré section method to study synchronization. As shown in figure 9(a), for each video frame the pixel intensities $X_{L},X_{R}$ are monitored in small interrogation windows on either side of the cell body and used to determine when a complete cycle has elapsed, thereby yielding a single scalar phase angle $\unicode[STIX]{x1D703}_{i}$ as a function of time for each flagellum $i$ . From this, one determines the normalized phase difference shown in figure 9(c):

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D6E5}=\frac{\unicode[STIX]{x1D703}_{1}-\unicode[STIX]{x1D703}_{2}}{2\unicode[STIX]{x03C0}}.\end{eqnarray}$$

A first discovery is that the beating dynamics is intrinsically noisy, even in the ostensibly synchronous state. The noise amplitude is orders of magnitude larger than that which arises from thermal fluctuations of a filament with the elasticity of the axoneme (Goldstein et al. 2009). We deduce therefore that this is internal biochemical noise, likely associated with calcium levels in the cell. Second, by studying time series over much longer than possible with the techniques available to Rüffer and Nultsch, we see that the three distinct behaviours they found (synchrony, transient asynchrony and decoupled motion), and which they attributed to distinct subpopulations of cells, are actually exhibited by a single cell, which stochastically switches from one state to another. Remarkably, the fraction of the population they found in each of those three behaviours is, within experimental error, the same as seen in our time series averages. This is the situation of ergodicity that is assumed to hold in statistical physics; that an instantaneous ensemble average yields the same result as temporal averages over a single member. It is not clear a priori that this should be the case in biology, but it is true here.

One way to extract information from the times series $\unicode[STIX]{x1D6E5}(t)$ is to consider the minimal model for coupled oscillators (Adler 1946), developed in the study of electronic circuits. The original Adler model is deterministic; the simplest stochastic generalization is

(3.2) $$\begin{eqnarray}\dot{\unicode[STIX]{x1D6E5}}=\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D708}-2\unicode[STIX]{x03C0}\unicode[STIX]{x1D716}\sin (2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6E5})+\unicode[STIX]{x1D709}(t).\end{eqnarray}$$

Here $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D708}=(\unicode[STIX]{x1D714}_{1}-\unicode[STIX]{x1D714}_{2})/2\unicode[STIX]{x03C0}$ is the scaled difference of intrinsic frequencies $\unicode[STIX]{x1D714}_{i}$ of the two oscillators, and $\unicode[STIX]{x1D716}$ is the coupling strength. The nonlinearity must respect the periodicity in $\unicode[STIX]{x1D6E5}$ and the sinusoidal coupling is the lowest-order term. The noise $\unicode[STIX]{x1D709}(t)$ is assumed to be the standard form, with $\langle \unicode[STIX]{x1D709}(t)\unicode[STIX]{x1D709}(t^{\prime })\rangle =2T_{eff}\unicode[STIX]{x1D6FF}(t-t^{\prime })$ , where $T_{eff}$ is an effective temperature. When $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D708}=0$ one sees immediately that if $\unicode[STIX]{x1D716}>0$ then the synchronized state $\unicode[STIX]{x1D6E5}=0$ is linearly stable, as are all equivalent states $\unicode[STIX]{x1D6E5}=2\unicode[STIX]{x03C0}n$ for $n=1,2,\ldots ,$ whereas the antiphase states $\unicode[STIX]{x1D6E5}=(2m+1)\unicode[STIX]{x03C0}$ for $m=0,1,2,\ldots$ are unstable. When $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D708}\neq 0$ the fixed points shift until a critical value $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D708}^{\ast }=\unicode[STIX]{x1D716}$ , at which point all fixed points are lost.

Figure 10. Time-averaged flow field around freely swimming Chlamydomonas and bead-spring model. (a) Experimental flow field, with streamlines (red) and schematic of cell body (green). (b) Theoretical flow field from three-Stokeslet model. (c) Decay of velocity field in various directions from cell. (d,e) Mechanism of elastohydrodynamic synchronisation in bead-spring model, as discussed in text. Reproduced from Drescher et al. (2010a ).

The noisy Adler model can be made more intuitive by observing that (3.2) can be interpreted as the overdamped dynamics of a Brownian microparticle in a potential $V(\unicode[STIX]{x1D6E5})$ , which generates a force given by the first two terms: $V(\unicode[STIX]{x1D6E5})=-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D708}\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D716}\cos (2\unicode[STIX]{x03C0}\unicode[STIX]{x1D6E5})$ . This is the ‘tilted washboard’ potential known widely in condensed matter physics. It is now possible to provide a heuristic interpretation of the observations first made by Rüffer and Nultsch: periods of noisy synchrony correspond to Brownian motion of the microsphere in one particular minimum of $V(\unicode[STIX]{x1D6E5})$ , periods of transient asynchrony (‘slips’) are thermally assisted hops left or right to an adjacent minimum, and ‘drifts’ are associated with a detuning $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D708}$ exceeding the critical value, so the potential has no local minima and $\unicode[STIX]{x1D6E5}$ evolves monotonically without bound. The three unknowns in the problem, $\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D708},\unicode[STIX]{x1D716}$ and $T_{eff}$ , can be found from experimental time series by using the mapping above. From the principles of statistical physics the ratio of left and right hopping probabilities satisfies

(3.3) $$\begin{eqnarray}\frac{p_{+}}{p_{-}}=\text{e}^{\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D708}/T_{eff}}.\end{eqnarray}$$

When the detuning is small the fluctuating synchrony is just the noisy dynamics of a particle in a quadratic potential, for which the autocorrelation function decays as

(3.4) $$\begin{eqnarray}R(t)=T_{\text{eff}}\unicode[STIX]{x1D70F}_{ac}\exp (-t/\unicode[STIX]{x1D70F}_{ac}),\quad \text{with}~\unicode[STIX]{x1D70F}_{ac}=\frac{1}{2\unicode[STIX]{x03C0}\sqrt{(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D716})^{2}-\unicode[STIX]{x1D6FF}\unicode[STIX]{x1D708}^{2}}}.\end{eqnarray}$$

Before discussing the results of extracting model parameters from the data, I note that there are many microscopic models that yield the Adler model; an intuitive class is based on representing a beating flagellum as a microsphere executing a periodic orbit. While this may seem a bit fanciful, given the complex beating waveform, there are several lines of reasoning suggesting that this is actually reasonable. First, the fluid flow at a given point due to some distant flagellum can be constructed from a superposition of Stokeslet contributions from the force distribution along the filament and the mean of that distribution is equivalent to a single microsphere on some orbit. Second, direct measurements of the flow around freely swimming Chlamydomonas cells have shown (Drescher et al. 2010a ) that the time-averaged flow (figure 10 a,c) is remarkably well captured by a superposition of three Stokeslets, one from the cell body pushing forward and one each at the approximate midpoints of the flagella, pushing backward, with overall force cancellation as required for a freely swimming organism. This extremely simple model captures the topology of the flow, with a stagnation point in front of the body and the swirls on either side and flow in from the anterior and posterior and outwards at the sides, as expected for a puller stresslet.

In the particularly insightful model of Niedermayer, Eckhardt & Lenz (2008), a beating flagellum is represented by a microsphere of radius $a$ and drag coefficient $\unicode[STIX]{x1D701}$ driven around an orbit by some tangential internal force $F$ that represents the action of molecular motors. A radial spring exerts a force $-\unicode[STIX]{x1D706}(R-R_{0})$ that tends to return the orbit to a radius $R_{0}$ . Two such spheres whose orbits are separated by a distance $d$ interact through the Stokeslet fields each creates, so the dynamics of the radius $R_{1}$ and angle $\unicode[STIX]{x1D719}_{1}$ of particle 1 are

(3.5a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D701}(R_{1}\dot{\unicode[STIX]{x1D719}}_{1}-\boldsymbol{e}_{\unicode[STIX]{x1D719}_{1}}\boldsymbol{\cdot }\boldsymbol{v}_{12})=F_{1}\quad \text{and}\quad \unicode[STIX]{x1D701}({\dot{R}}_{1}-\boldsymbol{e}_{R_{1}}\boldsymbol{\cdot }\boldsymbol{v}_{12})=-\unicode[STIX]{x1D706}(R_{1}-R_{0}),\end{eqnarray}$$

where $\boldsymbol{e}_{\unicode[STIX]{x1D719}_{1}}$ and $\boldsymbol{e}_{R_{1}}$ are unit vectors in the azimuthal and radial directions and $\boldsymbol{v}_{12}$ is the fluid velocity at the location of particle 1 due to the motion of 2, modelled as a Stokeslet. This model naturally leads to in-phase synchrony when the spheres orbit in the same sense. The physical mechanism is illustrated in figure 10(d,e) and relies upon the fact that the internal driving force $F$ is assumed constant. Therefore the drag force felt by the sphere, $\unicode[STIX]{x1D701}R\dot{\unicode[STIX]{x1D719}}$ , must also be constant. If particle 1 lags behind 2, its flow field will deflect 2 outward to a larger radius, thus slowing its angular velocity down, so 1 catches up with 2. If 1 leads 2, its flow field pushes 2 to a smaller radius and hence larger angular velocity, and 2 catches up with 1. Similar considerations hold from the opposite perspective. Clearly, the synchronous state is stable in this elastohydrodynamic mechanism.

The deterministic Adler equation for the difference in phase angles follows from this model under the assumptions that the radial dynamics is fast compared to azimuthal motion and that the orbits are well separated ( $d\gg R_{0}$ ). Interpreting the spring constant as that of a flagellum with length $\ell$ and bending modulus $A$ , one finds a coupling constant

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D716}\sim \frac{a}{d}\unicode[STIX]{x1D70F}_{r}\unicode[STIX]{x1D708}_{1}\unicode[STIX]{x1D708}_{2},\end{eqnarray}$$

with the relaxation time $\unicode[STIX]{x1D70F}_{r}\sim 3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}a\ell ^{3}/A$ , where $\unicode[STIX]{x1D707}$ is the fluid viscosity and $\unicode[STIX]{x1D708}_{1,2}$ are the intrinsic frequencies of the two oscillators. The results of analysing a very large amount of data is that the experimentally determined magnitude of the coupling constant is remarkably close to that predicted by the model for realistic parameter values. This was the first evidence that the coupling strength was consistent with a hydrodynamic mechanism. Subsequent analysis (Goldstein et al. 2009) showed that the time evolution of a slip is also well captured by the dynamics, showing the Adler model to be a viable candidate for the phase dynamics of coupled flagella.

When the two flagella fall out of synchrony during a drift state it is natural to expect the cell to turn, as one flagellum is beating faster than the other. The fact that an individual cell stochastically switches back and forth between periods of synchrony and asynchrony would then imply that the long-time trajectory is some kind of random walk. A population of these would then behave diffusively. A simple experiment was able to verify this prediction (Polin et al. 2009): a dilute suspension of cells was gently centrifuged to concentrate them at the chamber bottom, after which their random swimming motion diffusively spread out the concentration field. Careful analysis of the concentration profiles found consistency with Fick’s law, yielding a diffusion constant $D\sim 10^{-3}~\text{cm}^{2}~\text{s}^{-1}$ . As the translational diffusion constant arising from random walks can be interpreted as $D\sim u^{2}\unicode[STIX]{x1D70F}$ , where $u$ is a typical swimming speed and $\unicode[STIX]{x1D70F}$ the mean time between turns, the known value $u\sim 100~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$ for Chlamydomonas shows that $\unicode[STIX]{x1D70F}\sim 10$  s. It was just at this time that we had developed a three-dimensional tracking system (Drescher, Leptos & Goldstein 2009a ) in which swimming organisms as small as Chlamydomonas and as large as Volvox could be studied free from hydrodynamic effects from nearby walls. The device is a dual-view darkfield imaging system that uses red light-emitting diode ring lights for illumination and two charge-coupled device (CCD) cameras mounted on long-working-distance microscopes for long-distance observation of a centimetre-sized inner chamber suspended in a heat bath to prevent convective motions associated with differential heating. With this tracking system we confirmed that Chlamydomonas does indeed execute a kind of random walk, with the probability distribution function of duration of turnings agreeing precisely with that for the duration of drift events for pipette-held cells. Moreover, the distribution of free-flight times between turns was exponential with a time constant of ${\sim}11$  s, in agreement with the time scale inferred from the translational diffusion constant.

These results establish that the green alga Chlamydomonas executes a eukaryotic equivalent of the ‘run-and-tumble’ locomotion of peritrichously flagellated bacteria such as E. coli. Here it might be called ‘run-and-turn’. But the question of why it does this remain largely unanswered. In the bacterial case we understand the random motion as a way of searching space and performing chemotaxis (by biasing the probability of tumbles in response to chemical species), but there is no evidence yet to suggest whether this holds for Chlamydomonas. It has been suggested (Stocker & Durham 2009) that this might be a means of avoiding predators, but again there is no systematic study of this idea. As to how, we know that the two flagella respond differently to calcium, and it would appear likely that the transition between synchrony and asynchrony is driven by changes in the internal calcium levels. It remains an important experimental goal to achieve simultaneous measurements of flagellar synchrony and calcium levels within the cell.

In their vegetative state, wild-type Chlamydomonas flagella have a fixed natural length, precluding a test of the very strong length dependence to the coupling (3.6) predicted by the elastohydrodynamic mechanism. To test this, we took advantage of the fact that the volvocine algae have the ability to regrow their flagella after losing them through the action of mechanical forces or a change in pH. This process of regrowth takes over an hour, which is very long compared to the time it takes to make a measurement of synchrony. Hence, short experiments can be made during regrowth to capture the length dependence of the coupling constant. The results (Goldstein, Polin & Tuval 2011) show that whereas when the flagella are but a few microns long there is no visible synchrony, plateaus in $\unicode[STIX]{x1D6E5}(t)$ begin to appear more and more as time proceeds and the flagella lengthen. The scaling of the coupling constant with length is fully consistent with the elastohydrodynamic model.

As in most studies of biophysical phenomena, it is of interest to examine mutants for what their behaviour can reveal of the underlying mechanisms. In this regard, a class of mutants known as ptx has been of continuing interest to biologists because of their defective phototaxis and possible connections to human ciliopathies. Phototaxis in Chlamydomonas requires the two flagella to respond differently to internal calcium levels – displaying ‘flagellar dominance’ – with the cis flagellum dominant at low calcium levels and the trans dominant at high calcium levels. This differential response, coupled with rotation of the cell around a body-fixed axis as it swims, allows phototactic steering. In the mutants, this difference is lost and the two flagella behave nearly identically, frustrating phototaxis. Rüffer & Nultsch (1997) had observed in this mutant a different type of synchrony which, upon modern analysis, turns out to be precise antiphase (AP) synchrony (Leptos et al. 2013). As shown in figure 11(a,b) at very high temporal resolution (2000 fps), these cells exhibit stochastic transitions between the in-phase (IP) breaststroke and the AP mode, which is like the human freestyle swimming stroke, corresponding to half-integer jumps in $\unicode[STIX]{x1D6E5}$ . Using a fixed reference point on each flagellum to trace out the angle $\unicode[STIX]{x1D703}$ with respect to the cell midplane, the phase portraits in figure 11(c) show how accurately the single-cell limit cycles hug the upward and downward diagonals in IP and AP states, respectively. Averaging over many cells shows that the AP state is indeed precisely AP.

Figure 11. Antiphase synchronisation in the mutant ptx1. (a) Time series of phase difference $\unicode[STIX]{x1D6E5}$ and (b) instantaneous beat frequency (b). (c) Single-cell and (d) multicell average phase portraits in the IP and AP states. Reproduced from Leptos et al. (2013).

The existence of this antiphase state called for a re-examination of the elastohydrodynamic mechanism, and bead-spring models in particular. Note that while the mechanism shown in figure 10 supposes the two spheres to orbit in the same direction (appropriate to the ciliary arrangement in tissues), the two flagella of Chlamydomonas have their power strokes oriented in opposition, so that the appropriate sphere model has counter-rotating spheres. Reconsidering this situation reveals that antiphase synchrony is expected. Thus, the mutant behaves normally and the wild type is the outlier! This led us to conjecture (Leptos et al. 2013) that wild-type synchrony is controlled by a process internal to the cell, possibly associated with filamentary connections between the bases of the flagella.

We return to this point shortly, but note that the difficulty of disentangling the competing effects of hydrodynamics and internal couplings can be eliminated by examining the flagella on two physically separated cells whose only interaction is through the intervening fluid. This was accomplished (Brumley et al. 2014) by holding two somatic cells of Volvox on two separately movable micropipettes whose separation $d$ and orientation could be controlled. The two flagella of a given somatic cell beat in parallel, unlike in Chlamydomonas, and frequently beat together as a single unit, but one can also half-deflagellate each so indeed there are only two flagella interacting. In this way, the power strokes of the flagella could be oriented to be in the same or opposite directions, corresponding to the ciliary or Chlamydomonas-like orientations. The results of this experiment are clear (figure 12); flagella closer together than their length strongly phase lock, and as the separation increases the phase difference drifts at an ever increasing rate (figure 12 c). When synchronised, the fluctuations depend strongly on the cell–cell separation (figure 12 d), increasing with $d$ . This is easily understood within the stochastic Adler model (3.2) in the simplified case of identical oscillators, where by the interpretation of motion in a washboard potential the fluctuating dynamics around some synchronous state $\unicode[STIX]{x1D6E5}_{0}$ is that of a particle in a potential $V\sim \unicode[STIX]{x1D716}(\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D6E5}_{0})^{2}+\cdots \,$ , and by equipartition $\langle (\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D6E5}_{0})^{2}\rangle \sim T_{eff}/\unicode[STIX]{x1D716}$ , so with $\unicode[STIX]{x1D716}\sim 1/d$ we have a scaling collapse of the distribution of $(\unicode[STIX]{x1D6E5}-\unicode[STIX]{x1D6E5}_{0})/\sqrt{d}$ , as in figure 12(e). Furthermore, direct determination of the coupling strength shows consistency with a decay of $1/e$ , as expected for a Stokeslet interaction. When the power strokes are parallel the synchronization is in phase, and when they are opposed the synchronous state is antiphase, as predicted by the elastohydrodynamic model. These results show that hydrodynamic interactions are sufficient to synchronize flagella beating (but they do not show that such interactions are necessary).

Figure 12. Experimental test of hydrodynamically driven synchronization. (a,b) Set-up, with two cells held on micropipettes whose separation and orientation can be adjusted. (c) Phase difference versus time at various intercell separations. (d) Fluctuations about synchrony at various separations. (e) Probability distribution functions of fluctuations and scaling collapse, where $L=d/\ell$ , consistent with hydrodynamic coupling. Reproduced from Brumley et al. (2014).

Analysis of the beating dynamics in these experiments revealed many interesting features of flagellar waveforms, including evidence for the orbital compliance that lies at the heart of the elastohydrodynamic mechanism of synchrony. In addition, time-resolved PIV measurements of the flow field around a beating flagellum allowed a test of the orbiting microsphere model discussed earlier. The result, shown in figure 13 is that indeed the flow field can be accurately captured by a Stokeslet moving on an orbit, but not surprisingly there are significant variations in the magnitude through the beat period. And indeed, as pointed out by Uchida & Golestanian (2011, 2012), beating synchrony can arise not only from orbital compliance with fixed internal forcing, but also by variations in internal forcing along fixed orbits. Elegant experimental studies with colloidal oscillators (Brout & Cicuta 2016) probe in detail the competition between these effects.

Figure 13. Flagellar beating and orbiting Stokeslets. (a) Results of fitting the measured instantaneous velocity field around a single flagellum to that of a Stokeslet. (b) The fitted Stokeslet position (red dots) throughout the average flagellar beat cycle. (c) Stokeslet magnitude as a function of time during the beat period $T$ . Reproduced from Brumley et al. (2014).

Now it is necessary to consider some biology and return to ptx1. It has been known for decades (Ringo 1967) that there is an elaborate network of filaments clustered around the basal bodies (BBs) of the flagella in green algae, some of them directly linking the BBs together, and there is evidence for contractile behaviour of those filaments. We conjectured that the difference between the mutant and wild type had its origin in these filaments, and that modifications to those filaments might influence synchrony in a way that would be revealing of the mechanism. Interestingly, there is another class of mutants known as vfl for ‘variable flagellar number’ in which those filaments are missing and as a consequence there may be more than two flagella. The study by Quaranta, Aubin-Tam & Tam (2015) showed no evidence of synchrony in these mutants. A more extensive examination (Wan & Goldstein 2016) showed that when the flagella of these mutants are close enough together for strong hydrodynamic interactions then synchrony is observed, with the symmetries found in the measurements of separated Volvox somatic cells described above. Moreover, nature has provided us with not only biflagellates, but unicellular quadri-, octo- and even hexadecaflagellates, all with known and often elaborate networks of filamentary connections between the flagellar BBs. Systematic study (Wan & Goldstein 2016) of the beating patterns of these organisms show that often the symmetries of the BB connections are mirrored in the patterns of flagellar actuation. An example is shown in figure 14 in the case of quadriflagellates. The observed beat patterns mimic those of quadrupeds which can gallop, trot and pronk, and are amenable to a group-theoretic classification. These results strongly suggest that beating synchrony arises from a combination of direct hydrodynamic interactions and the effects of intracellular elastic couplings. Amusingly, the latter contribution is reminiscent of that which Huygens deduced led to the synchrony of nearby pendulum clocks. A complete theory of the hydrodynamic/elastic competition in flagellar synchrony remains to be developed.

Figure 14. Quadriflagellates. In different species the flagella are arranged in one of two possible configurations (types I, II). Observed patterns of actuation include (a) the trot, (b) pronk, (c) rotary and (d) transverse gallops, as indicated by coloured arrows in (a–d). In each panel, the coloured bars depict time series of beat phases. Reproduced from Wan & Goldstein (2016).

I close this section with a brief discussion of additional developments in the study of flagellar synchrony. The initial quantifications of noise in the beating of Chlamydomonas flagella (Goldstein et al. 2009; Polin et al. 2009) were not well resolved in time, but more recent precision studies (Wan & Goldstein 2014; Wan et al. 2014) have begun to reveal that there are important and largely unexplained systematics both within a beat period and on much larger time scales. For example, the noise amplitude peaks within a cycle during the transition between the power and recovery strokes. Is this a direct reflection of the statistics of a finite number of motors acting in concert? Just as it is known that the human heartbeat exhibits a broad distribution of inter-beat intervals, so too does Chlamydomonas flagellar beating. Here, there are long-range correlations in the beat intervals extending to hundreds of beats (tens of seconds). The origin of this is a mystery. Finally, we discovered that Volvox carteri displays robust metachronal waves (Brumley et al. 2012, 2015). This finding is significant in the search for a mechanistic understanding of metachronal waves because the spacing between the somatic cells of Volvox is comparable to or larger than the flagellar length, and there are no direct intercellular connections between the somatic cells. In this case we expect these flagella to be in a more weakly coupled regime and hence more amenable to theoretical interpretation than are the tightly packed cilia of Paramecium, which has been the standard organism of study in this area for some time. The emerging view from bead-spring models of synchrony is that metachronism arises from the proximity to the no-slip cell surface, but there has been no definitive test of this idea.

4 Surface interactions of microswimmers

Swimming microorganisms are frequently found near surfaces, whether it is spermatozoa in Fallopian tubes or bacteria in the soil. There has been a longstanding interest in the hydrodynamic interactions between swimmers and surfaces. For example, it is now well established that bacteria swimming through the action of helical flagella swim in circles near a solid surface due to the chirality of the flagella (Lauga et al. 2006), and that peritrichously flagellated bacteria can switch the pole at which they bundle their flagella when they encounter an obstacle (Cisneros et al. 2006). Observation of large-scale collective behaviour in suspensions of motile bacteria (Dombrowski et al. 2004) discussed in § 5 led naturally to the question of whether organisms with a more symmetric distribution of force generation on their surface would do the same. As bacteria push the fluid back with their flagella and exert a forward-directed force on the fluid by their body, they define a pusher stresslet, whereas an organism like Chlamydomonas is a puller. The high symmetry of Volvox suggests that at leading order it is neither a puller nor a pusher. Casual observations on suspensions of Volvox indicated that there was a tendency toward pairwise interactions, and further investigation in the dual-view system tracking system led to the discovery of what we term ‘hydrodynamic bound states’ (Drescher et al. 2009b ). These are formed when two colonies that have swum up to the top of a glass chamber are attracted together until they nearly touch, and then orbit about their centre as shown in figure 15. Volvox colonies all spin in the same sense around their body-fixed axis (clockwise when viewed from above), and the pair also orbit in a clockwise sense, so the adjacent colony surfaces are moving in opposition (unlike meshing gears).

Figure 15. Waltzing Volvox. (a) Images using dual-view apparatus: (i) overlaid temporal sequence showing orbiting motion in bound state, (ii) side view and (iii) top view of bound state with streamlines from PIV, (iv) cluster of colonies. (b) Calculational geometry of two downward-pointing Stokeslets and (c) two nearby corotating colonies. (d) Infalling trajectories compared to model. Inset shows orbital frequency as a function of colony rotation frequency and theoretical fit. Adapted from Goldstein (2015).

In searching to explain the apparent attractive force bringing them together we were reminded of the work of Squires (2001), who discussed how two particles hovering near an upper no-slip surface and acted upon by gravity would experience just such an attraction. The levitation was imagined to arise from e.g. electrophoretic effects, and the problem was conceptualized as the coupled dynamics of two Stokeslets a distance $h$ below the surface, separated by a lateral distance $r(t)$ as in figure 15(b). Placed in the context of the dynamics of two Volvox colonies, the coupled dynamics of their locations $\boldsymbol{x}_{i}$ would be

(4.1) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\boldsymbol{x}}_{i}=\boldsymbol{u}(\boldsymbol{x}_{i})+\boldsymbol{v}_{i}, & \displaystyle\end{eqnarray}$$

(4.2) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\boldsymbol{p}}_{i}=\frac{1}{\unicode[STIX]{x1D70F}}\boldsymbol{p}_{i}\times (\hat{\boldsymbol{z}}\times \boldsymbol{p}_{i})+\frac{1}{2}(\unicode[STIX]{x1D735}\times \boldsymbol{u})\times \boldsymbol{p}_{i}, & \displaystyle\end{eqnarray}$$

where $\boldsymbol{u}(\boldsymbol{x}_{i})$ is the fluid velocity at particle $i$ , say due to the other, $\boldsymbol{p}_{i}$ is the direction of its axis, $\boldsymbol{v}_{i}$ is its intrinsic swimming velocity, $\hat{\boldsymbol{z}}$ is the vertical direction and $\unicode[STIX]{x1D70F}=6\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}_{c}g\ell$ is the time scale for viscous relaxation of the axis due to gravity, with $\unicode[STIX]{x1D70C}_{c}$ the colony density and $\ell$ the displacement between the centres of gravity and geometry of the colony. When the colonies’ axes remain vertical and they touch the upper surface, the Stokeslet displacement $h=R$ . Using Blake’s result for a Stokeslet near a wall to deduce the in-plane velocity $u$ , one finds (Squires 2001) that the scaled separation $x(t)=r(t)/h$ obeys

(4.3) $$\begin{eqnarray}\frac{\text{d}x}{\text{d}t}=-\frac{3F}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}h^{2}}\frac{x}{(x^{2}+4)^{5/4}},\end{eqnarray}$$

where $F$ is the magnitude of the gravitational force. By studying the upswimming speed of flagellated colonies and the settling velocity of deflagellated ones, along with the orientational relaxation time, it is possible to determine experimentally all the parameters of this model. The measured infalling trajectories are in excellent agreement with theory (figure 15 d). With a bit more work it is possible to develop a lubrication theory for the rotational motion of the colonies and explain the observed orbital frequency $\unicode[STIX]{x1D6FA}$ as a function of rotational frequency of the colony $\unicode[STIX]{x1D714}$ , as in figure 15(c,d).

Figure 16. Flow field around a colony of Volvox carteri. (a) Residual magnitude, velocity vector field and streamlines after subtraction of fitted Stokeslet. (b) Velocity magnitude as a function of distance from centre of colony. Data (blue circles) averaged over 19 colonies (black dots), fitted average Stokeslet (dashed green line). Deviations from a pure Stokeslet appear for $r/R<5$ and are captured by the addition of a source doublet and a stresslet (red solid line). (c) Vertical section of the flow field through the colony centre. Forward–backward asymmetry arises from stresslet component. Reproduced from Drescher et al. (2010a ).

The degree to which the data on infalling trajectories were quantitatively explained by the advected-Stokeslets model was at first somewhat confusing: why should the dynamics of a Volvox colony be so accurately captured by a Stokeslet flow field? For sure the Stokeslet field, falling off as $1/r$ , should dominate sufficiently far from the organism, but the lack of precise uniformity of the surface somatic cell density and the presence of the colony surface itself should induce contributions that would be appreciable at small separations. This led to an investigation (Drescher et al. 2010b ) of the flow field around freely swimming colonies. For Volvox, the experimental set-up consisted of a chamber with a dilute suspension of colonies and fluorescent microspheres, illuminated with a laser light sheet and observed from the side by a lightweight CCD camera on a long-working-distance microscope. The camera itself was mounted on a dual-axis computer-controlled motorised microscope stage actuated so as to keep a given colony in the centre of the field of view as it swims upward against gravity. Because of the steadiness of the colonial motion it was possible to average over as many as 50 000 video frames to obtain extremely accurate PIV maps of the velocity field around the colony. The result is shown in figure 16, where we see that the flow field is indeed dominated by the Stokeslet contribution down to remarkably small distances from the colony surface, with small contributions arising from a source doublet and a stresslet. This provides ex post facto justification for the use of Squires’ result for the lateral advection of Stokeslets in explaining the infalling trajectories of colonies forming a hydrodynamic bound state.

Figure 17. Flows around freely swimming E. coli. Average flow field far from surfaces (a–d) and near chamber surface (d–h). Streamlines indicate the local direction of flow, and the logarithmic colour scheme indicates flow-speed magnitudes. (a) Measured flow field in swimming plane; inset shows anterior–posterior asymmetry close to the cell. (b) Best-fit force dipole flow and (c) residual of fit. (d) Radial decay of the flow speed from centre of cell body in different directions. (e) As in (a), but for bacteria swimming $2~\unicode[STIX]{x03BC}\text{m}$ from chamber bottom. (f) Best-fit force dipole flow and (g) residuals of fit. (h) As in (d), but now showing the $r^{-3}$ decay of the flow speed. Reproduced from Drescher et al. (2011).

Measurements of the flow fields around individual freely swimming cells were then extended to the case of bacteria using a different method (Drescher et al. 2010b , 2011). Here, rather than using a tracking microscope, we adopted a fixed (horizontal) field of view on an inverted microscope, with a typically narrow focal plane width, and monitored swimming cells and advected tracer particles. Using only those trajectories that stayed within the focal plane for some prescribed time period, a very large number of short trajectories of E. coli cells was collected, and the associated PIV flow fields were reoriented based on the local cell swimming direction to yield a PIV map in the frame of motion of the cell. This was done for cells swimming far from the chamber surface as well as close to it, providing a test of the predicted effects of a nearby no-slip surface. From general considerations we expect the velocity field in bulk to fall off as $u\sim F\ell /\unicode[STIX]{x1D707}r^{2}$ , where for a force-free swimmer there is a force dipole produced by a forward force $F$ from the head and a backward force $-F$ from the tail, separated by some microscopic distance $\ell$ . Figure 17 presents the data as PIV maps of the velocity fields and as plots of the radial falloff of the velocity, confirming the expected $r^{-2}$ behaviour. From the fits we obtain $\ell =1.9~\unicode[STIX]{x03BC}\text{m}$ and $F=0.42$  pN, consistent both with the geometry of the bacterium and calculations of the force from a rotating helical flagellum. For cells swimming $2~\unicode[STIX]{x03BC}\text{m}$ from the chamber surface, the PIV map shows the expected closed streamlines and a radial decay of the velocity as $r^{-3}$ , as expected for a stresslet parallel to a no-slip surface.

Armed with these results, one can investigate the competition between reorienting hydrodynamic interactions between pairs of cells and the decorrelating effects of rotational diffusion. The stresslet flow field falls off as $u\sim A/r^{2}$ and its vorticity decays as $\unicode[STIX]{x1D6FA}\sim A/r^{3}$ . Two bacteria swimming at speed $V_{0}$ and of size $a$ may be thought to interact on a time scale $\unicode[STIX]{x1D70F}\sim a/V_{0}$ . The mean squared angular displacement of one cell relative to the other during that time, due to hydrodynamic interactions, will be $\langle \unicode[STIX]{x0394}\unicode[STIX]{x1D719}^{2}\rangle _{H}\sim A^{2}\unicode[STIX]{x1D70F}^{2}/r^{6}$ . Rotational Brownian motion provides a competing mechanism for angular decorrelation, with $\langle \unicode[STIX]{x0394}\unicode[STIX]{x1D719}^{2}\rangle _{B}=4D_{r}\unicode[STIX]{x1D70F}$ , where $D_{r}$ is the rotational diffusion constant. Equating these two yields a length (the ‘hydrodynamic horizon’)

(4.4) $$\begin{eqnarray}r_{H}\sim \left(\frac{A^{2}\unicode[STIX]{x1D70F}}{D_{r}}\right)^{1/6}\end{eqnarray}$$

beyond which hydrodynamic interactions are irrelevant. The small exponent in this relation shows that this length is very insensitive to the details of the microscopic parameters of the calculation. Estimates using experimental values of $D_{r}$ and $A$ for E. coli, yield $r_{H}$ of just a few microns, or a few body diameters. This suggests that in dilute suspensions long-range hydrodynamic interactions between swimmers are unimportant and that when the fluid flows do matter it is on the scale at which short-range effects such as steric repulsion and direct flagellar interactions are also important. It is only when the concentration in a suspension is sufficiently high that the mean distance between cells becomes comparable to $r_{H}$ that hydrodynamic interactions start to dominate. Estimates suggest this requires volume fractions in the range of a few per cent, and indeed this is roughly the point at which collective behaviour begins to be seen, as discussed in § 5. The problem of cell–surface interactions with competing hydrodynamic effects and rotational Brownian motion leads to a non-trivial Kramers-like problem for escape from the surface that has only begun to be explored (Drescher et al. 2011).

It has been known for many years that swimming cells accumulate near surfaces. One of the prominent early observations was that of Rothschild (1963) who found this in the case of sperm cells, and more recent work on bacteria discovered a similar phenomenon (Berke et al. 2008) and attributed the effect to long-range hydrodynamic interactions between the swimmers (viewed as pusher stresslets) and the no-slip wall. An opposing point of view advocated by Li & Tang (2009), but actually first mentioned by Rothschild (1963), is that the accumulation is associated with what has come to be called ‘inelastic scattering’ of cells by surfaces. That is, rather than the kind of elastic reflection one might see with a ball bouncing off the surface, with equal incident and final angles with respect to the surface normal, here nearly all incident angles result in the same outgoing orientation, nearly parallel to the surface. The mechanism for this involves specific and non-hydrodynamic interactions of the flagella with the surface, and are thus beyond the reach of the simplest models based on singularities.

Following Berra’s advice, we sought to understand this issue in more detail by high-speed imaging of Chlamydomonas and sperm interacting with the chamber walls of microfluidic devices (Kantsler et al. 2013). Here I will focus on the key results from the algal experiments, shown in figure 18. This can be seen most clearly from the frames of high-speed videos of cells approaching and then interacting with surfaces. In the case of the wild type (figure 18 a), we see clearly that the first contact made by the alga with the surface is through one of its flagella, and as beating continues the cell body reorients to be just slightly away from parallel to the surface, and then it swims away. Heuristically, the small final angle is set by the larger spread of the beating waveform relative to the cell body size. This can be verified by studying mutants with longer flagella (the strain CC-289 lf3-2) and shorter flagella (CC-2347 shf1). In addition, the amusingly named mutant mbo1 (for ‘move backward only’) beats its flagella in a sinusoidal waveform at the anterior of the body and pushes itself in the opposite direction to the wild type. When it encounters a surface it is driven into it at an angle and remains trapped for some time. Taken together, these results show that ciliary contact interactions dominate the surface interactions of green algae. However, a more detailed analysis of the trajectories in a high-throughput experiment does show clear evidence for a competition between steric and hydrodynamic interactions as a function of scattering parameters (Contino et al. 2015). Results of this type can be used to design microfluidic devices that can rectify swimmers’ motion in interesting ways (Kantsler et al. 2013).

Figure 18. Chlamydomonas surface scattering. (a) Upper panels: images from high-speed video (with time indicated) showing a wild-type cell encountering a surface and reorienting through contact interactions. Flagella are highlighted in lower panels. (b) As in (a), but for the ‘move backward only’ mutant. Scale bar is $20~\unicode[STIX]{x03BC}\text{m}$ . (c) Conditional probability distributions for wild type, mbo1, short- and long-flagella mutants. (d) Scattering probability distributions shows how cilia length and swimming gait determine details of inelastic scattering. (e) Illustration of the scattering and trapping mechanisms. Reproduced from Goldstein (2015).

I close this section by noting that the strong interaction of swimming microorganisms with surfaces is thought to underlie the process of rheotaxis, in which they swim upstream against a bulk flow. Of particular interest recently has been bacterial rheotaxis (Marcos et al. 2012), where the chirality of the helical flagella plays an important role in reorientation of cells along the walls bounding the flow. Analogous phenomena have been seen in sperm rheotaxis (Kantsler et al. 2014). In all of these contexts, it remains an open problem to observe the development of the spatially non-uniform density distributions from an initially uniform suspension as a means of confirming the underlying mechanism.

Figure 19. Collective behaviour in a bacterial suspension. PIV map of velocity field in a suspension of B. subtilis, showing vortices on scales large compared to individual bacteria. Adapted from Dombrowski et al. (2004).

5 Collective behaviour in microswimmer suspensions

In this section I discuss the nature of collective behaviour in suspensions of motile bacteria. This has developed into an intensely studied subject which has been reviewed extensively elsewhere (Ramaswamy 2010; Marchetti et al. 2013), so my purpose here is to trace briefly the key historical developments and recent discoveries that connect with issues discussed in previous sections, particularly the work on cytoplasmic streaming. Motivated by collective behaviour exhibited by organisms as diverse as schooling fish and herding wildebeest, models began to emerge in the mid-1990s to describe the possible orientational self-organization of motile entities. The original discrete model, due to Vicsek et al. (1995), was a computational algorithm by which a moving particle updated its velocity at each time step to be the mean of its nearest neighbours, while continuum models (Toner & Tu 1995) were generalisations of spin models known from equilibrium statistical physics that incorporate alignment interactions to include advective nonlinearities. Both models predicted that there could be long-range order in such systems, a state in which all the organisms move in some randomly selected direction. In the two-dimensional case in particular it was noted that such long-range order does not violate the Mermin & Wagner (1966) theorem (which forbids long-range order in two-dimensional systems with short-range interactions and a continuous symmetry) precisely because of the non-equilibrium aspects of the problem. Later, Simha & Ramaswamy (2002) noted that for the case of swimming (pusher) microorganisms in suspension, hydrodynamic interactions between them would lead to a long-wavelength instability, frustrating such order.

Experiments by Wu & Libchaber (2000) on E. coli in soap films saw hints of the lack of such order, and later we discovered (Dombrowski et al. 2004; Tuval et al. 2005) suspensions of swimming Bacillus subtilis exhibit a dynamical state resembling turbulence, with transient recurring vortices and jets of collective swimming on length scales large compared to the individual cells, as shown in figure 19. This was experimental verification of the prediction of an intrinsic instability in pusher suspensions. Note that the collective swimming speed and vortex size are large enough that the Péclet number in these systems can also approach or even exceed unity, implying that flow-induced transport of chemical species can outcompete diffusion, even though at the single-cell level this is not the case. Subsequent experimental work on thin films (Sokolov et al. 2007) has shown analogous behaviour, and further investigations of quasi-two-dimensional systems (Cisneros et al. 2007, 2011; Wensink et al. 2012; Dunkel et al. 2013b ) have explored in more detail the nature of this mesoscale ‘bacterial turbulence’, and in particular the statistics of velocity distributions and correlation functions.

The key features necessary to describe the collective state were identified through numerical studies of suspensions of self-propelled rods (Saintillan & Shelley 2007), and a detailed kinetic theory (Saintillan & Shelley 2008a , b ). As in the description of self-organised streaming in § 2, the starting point for this analysis is the time evolution of the distribution function $\unicode[STIX]{x1D6F9}(\boldsymbol{x},\boldsymbol{p},t)$ for particle position $\boldsymbol{x}$ and orientation vector $\boldsymbol{p}$ , $\unicode[STIX]{x1D6F9}_{t}=-\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }(\dot{\boldsymbol{x}}\unicode[STIX]{x1D6F9})-\unicode[STIX]{x1D735}_{\boldsymbol{p}}\boldsymbol{\cdot }(\dot{\boldsymbol{p}}\unicode[STIX]{x1D6F9})$ , where the single-particle time evolution equations are

(5.1) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\boldsymbol{x}}=U_{0}\boldsymbol{p}+\boldsymbol{u}-D\unicode[STIX]{x1D735}_{\boldsymbol{x}}\left(\ln \unicode[STIX]{x1D6F9}\right), & \displaystyle\end{eqnarray}$$

(5.2) $$\begin{eqnarray}\displaystyle & \displaystyle \dot{\boldsymbol{p}}=(\boldsymbol{I}-\boldsymbol{p}\boldsymbol{p})\boldsymbol{\cdot }[(\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D640}+\unicode[STIX]{x1D652})\boldsymbol{\cdot }\boldsymbol{p}-d_{r}\unicode[STIX]{x1D735}_{\boldsymbol{p}}(\ln \unicode[STIX]{x1D6F9})]. & \displaystyle\end{eqnarray}$$

Here, $U_{0}$ is the particle’s swimming speed, $\boldsymbol{u}$ is the local fluid velocity, $D$ and $d_{r}$ are the translational and rotational diffusion constants, $\unicode[STIX]{x1D640}$ and $\unicode[STIX]{x1D652}$ are the rate of strain and vorticity tensors and $\unicode[STIX]{x1D6FE}$ is a shape parameter ${\sim}1$ for rods. These equations are supplemented by that for the fluid velocity,

(5.3) $$\begin{eqnarray}-\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}_{\boldsymbol{x}}^{2}\boldsymbol{u}+\unicode[STIX]{x1D735}_{\boldsymbol{ x}}q=\unicode[STIX]{x1D735}_{\boldsymbol{x}}\boldsymbol{\cdot }\unicode[STIX]{x1D72E},\end{eqnarray}$$

with $q$ the pressure and

(5.4) $$\begin{eqnarray}\unicode[STIX]{x1D72E}=\unicode[STIX]{x1D70E}_{0}\int \text{d}\boldsymbol{p}\unicode[STIX]{x1D6F9}(\boldsymbol{x},\boldsymbol{p},t)\left(\boldsymbol{p}\boldsymbol{p}-\frac{\boldsymbol{I}}{3}\right)\end{eqnarray}$$

the active contribution to the stress tensor (Batchelor 1971). Linear stability analysis around a fully ordered state shows that in the case of pushers (but not pullers) there is a finite-wavelength instability, and numerical studies reveal that the full model generates transient recurrent vortices and jets on large length scales, as seen in experiments. Many other models for this general phenomenon have been developed, both computational (Lushi, Goldstein & Shelley 2012) and continuum (Dunkel et al. 2013a , b ), the latter incorporating steric interactions between cells and exploring the role of particle shape (Wensink et al. 2014). Computational models of discrete particles in particular have been able to study carefully the role of long-range hydrodynamic interactions in generating the collective state.

In § 2 we discussed the apparent self-organisation that may underlie cytoplasmic streaming in plant cells, and noted that the hydrodynamic interactions between filaments laden with motors would tend to be like those of ‘pusher’ microswimmers. This idea led us to suggest, on the basis of a model of active matter similar to that above, that a suspension of stresslets would, if confined to sufficiently small domains, exhibit self-organisation even without self-locomotion (Woodhouse & Goldstein 2012b ). The model predicts that in a suspension of concentration $c_{0}$ there is a critical stresslet strength

(5.5) $$\begin{eqnarray}\unicode[STIX]{x1D70E}^{\ast }\simeq \frac{16\unicode[STIX]{x1D707}}{c_{0}}\left(9.33\frac{D}{L^{2}}+d_{r}\right),\end{eqnarray}$$

above which a domain of diameter $L$ exhibits a spontaneously circulating spiral vortex.

A test of this predicted self-organisation was made with motile bacteria by a method in which a suspension is turned into an emulsion by mixing it with surfactant and mineral oil (Wioland et al. 2013). The result is a suspension of droplets of concentrated bacteria in a matrix of oil. A single such droplet confined between two coverslips takes on a pancake shape as in figure 20, with a lateral size in the range of $15{-}150~\unicode[STIX]{x03BC}\text{m}$ . The surfactant-covered oil–water interface is approximately no slip. The main experimental result is a direct confirmation of the prediction of spontaneous circulation in the form of a spiral vortex below a critical radius. Roughly speaking, when the drop diameter is smaller than the typical vortex size observed in the bulk collective state then the confinement can stabilise it. In addition to the predicted global behaviour, this experiment had a surprising and unanticipated feature: a counter-rotating boundary layer at the drop periphery (‘edge current’), perhaps one or two cells wide. It was initially unclear the direction in which the cells were swimming in this layer, relative to the bulk, but numerical studies (Lushi, Wioland & Goldstein 2014) suggested the counterintuitive result that the cells within the bulk vortex actually swim upstream against the ‘backwash’ from the boundary layer. This prediction was verified by use of a two-colour fluorescent labelling technique that allows one to visualise the cell body separately from the flagella, and thereby arrive at an unambiguous determination of the cell orientations in the domain.

Figure 20. Confinement of a bacterial suspension. (a) A drop of bacterial suspension in a matrix of mineral oil. (b) Indication of bulk circulation in a drop, viewed from below. (c) PIV flow field in drop, with enlargement of boundary region showing edge current. Reproduced from Wioland et al. (2013). (d,e) Method for determining cellular orientation inside drop. Brightfield image with cell membrane (red at time $t=0$ and blue at $t=0.2$  s) and flagella (green at $t=0.1$  s) (d) Cell at the oil–water interface both points and moves to the top left corner. (e) Cell in bulk is pointing to the top left corner while moving overall in the opposite direction. (Scale bars: drop images, $10~\unicode[STIX]{x03BC}\text{m}$ ; individual bacteria, $5~\unicode[STIX]{x03BC}\text{m}$ .) Reproduced from Lushi et al. (2014).

The existence of these vortices suggested an investigation of the coupling between them, first as a pair of confining circular regions connected by a thin neck, and then extended to arrays of coupled domains (Wioland et al. 2016b ): a ‘bacterial vortex lattice’. Experiments showed that when the neck between two adjacent regions was small the edge currents did not penetrate into the constriction, but instead bypass it (figure 21 c) so that the neighbouring domains circulate in opposite directions, with antiparallel vorticity. In the language of statistical physics, this is a case of antiferromagnetic coupling. When the gap is wider, the edge currents trace along the continuous boundary and drive adjacent vortices to circulate in the same direction, and thus be ferromagnetically coupled, with parallel vorticity vectors. Intriguingly, the distribution of vorticity in these lattices can be understood using lattice field theory methods well known in statistical physics. In very recent work, we studied these same bacterial suspensions confined to annuli whose width was in the range $10{-}100~\unicode[STIX]{x03BC}\text{m}$ and up to several millimetres long (Wioland, Lushi & Goldstein 2016a ). Here, if the transverse dimension is sufficiently small there would be spontaneous unidirectional motion around the ‘racetrack’, but when the channel width exceeded a size comparable to the intrinsic vortex scale ( ${\sim}70~\unicode[STIX]{x03BC}\text{m}$ ), the mean flow would rapidly vanish.

Figure 21. Coupled bacterial vortices. (a) Microfluidic lattice in the regime of narrow connecting gaps promoting antiferromagnetic order. (b) Close-up of nearby domains and PIV velocity field. (c) Schematic of edge and bulk flows near a gap. (d) Lattice model of vortices and pillars in colour coded in antiferromagnetic state, with couplings indicated by solid and dashed lines. (e–h) The ferromagnetic case at large gaps. Reproduced from Wioland et al. (2016b ).

I end this section by pointing out that since the work of Wu & Libchaber (2000) and others (Angelani, Leonardo & Ruocco 2009) it has been apparent that a suspension of swimmers acts as a kind of heat bath in which the organisms play the role of molecules, driving fluctuations of suspended tracer particles. Taking advantage of the high-speed swimming of the green algae and their large size relative to typical tracer particles, it is possible to resolve the dynamics of those collisions and test this idea directly. The results show that there can be markedly non-Gaussian statistics of those tracers, despite particle displacements obeying a diffusive scaling (Leptos et al. 2009). Subsequent work (Rushkin, Kantsler & Goldstein 2010) showed that the applicability of the central limit theorem in the random superposition of flow fields from point singularities depends on the order of the singularity, with Stokeslets and stresslets on either side of the marginal power law. Detailed theories for the probability distributions are now emerging (Thiffeault 2015; Zaid & Mizuno 2016).

6 Outlook

As many people have said: predictions are difficult, especially about the future. Nevertheless I will offer some thoughts about research questions whose solution may well depend on fluid mechanical ideas and methods. The first is the broad class of problems in the general area of symbiosis between different kingdoms of life, such as between bacteria and algae. This has come to the fore through the discovery (Croft et al. 2005) that certain species of algae acquire the needed vitamin $\text{B}_{12}$ through a symbiotic relationship with bacteria (see also Helliwell et al. 2011). How do these species find each other and stay close enough for the symbiosis from one generation to the next? What is the interplay between diffusion, chemotaxis and advection in these natural contexts? More generally, how do we measure experimentally and quantify theoretically the advantages conferred by symbiotic relationships in space and time? Some ideas about this have been proposed using methods from theoretical economics (Peaudecerf & Goldstein 2015), but much more work is needed. A second area of interest is within the broad field of developmental biology, where the mechanics and dynamics of tissues have come to be recognized as having hydrodynamic aspects (Ranft et al. 2010). Examples of systems in which this has begun to be explored include gastrulation in Drosophila (He et al. 2014), where large-scale morphological transformations are driven by cell shape changes. Similar issues arise in the mechanics of embryonic inversion of Volvox (Haas & Goldstein 2015; Höhn et al. 2015). Finally, moving beyond evolutionary transitions from single- to multicellular organisms we encounter the development of animals. Here, there is growing interest in the choanoflagellates, the closest living animal ancestors. Not only can they exist in unicellular and colonial forms, and display fascinating motility (Kirkegaard, Marron & Goldstein 2016), but they are closely related to choanocytes, the building blocks of sponges, which occupy an important place in the evolutionary history of the animals (Butterfield 2015).