Phase Dependent Forcing and Synchronization in the three-sphere model of Chlamydomonas
The green alga Chlamydomonas swims with synchronized beating of its two flagella, and is experimentally observed to exhibit run-and-tumble behaviour similar to bacteria. Recently we studied a simple hydrodynamic three-sphere model of Chlamydomonas with a phase dependent driving force which can produce run-and-tumble behaviour when intrinsic noise is added, due to the non-linear mechanics of the system. Here, we consider the noiseless case and explore numerically the parameter space in the driving force profiles, which determine whether or not the synchronized state evolves from a given initial condition, as well as the stability of the synchronized state. We find that phase dependent forcing, or a beat pattern, is necessary for stable synchronization in the geometry we work with.
pacs:05.45.45.Xt, 87.16.Qp, 47.63.-b
Microorganisms swim in the low Reynolds number regime where viscous forces dominate, inertia is negligible and the familiar propulsion methods of larger organisms become ineffective [1, 2, 3, 4]. Fluid flow is governed by the Stokes equation, which is time reversible. A necessary condition on a periodic swimming stroke in order to achieve net propulsion is that it is non-time reversible . Inspired by sperm cells, which achieve propulsion by propagation of bending waves through their flagellum, Taylor demonstrated that propulsion is possible in a viscous environment by studying the propagation of waves on an infinite sheet . Purcell showed that a swimmer needs at least two compact degrees of freedom to break the time reversal symmetry and achieve net propulsion .
Many microorganisms swim using flagella ; there are two fundamentally different types of flagella: bacterial flagella and eukaryotic flagella (or cilia). Eukaryotic flagella form bends when microtubules on one side of the flagella ‘walk’ or ‘slide’ along the microtubules on the other side . The propagation of bends allows the flagella to form beat patterns that can break the time reversal symmetry. For example, the first half of an individual cilium’s beat cycle, called the power stroke, has the cilium sticking out and pushing the fluid, while the second half, called the recovery stroke, has the cilium bent as it returns to its original position .
Our understanding of propulsion at low Reynolds number has been developed by theoretical model microswimmers. Lighthill demonstrated a model that can achieve net propulsion by studying periodic shape deformations of a nearly spherical swimmer, showing that the swimming velocity is at most of the order of the square of the amplitude of the deformations . Purcell’s three-link swimmer was studied by Becker et al., who determined the swimming direction and velocity for different angle amplitudes and relative link lengths [1, 11]. A useful one-dimensional model is the linear three-sphere swimmer, where three beads are connected by two rods that change length with a non-reciprocal pattern . Dreyfus et al. studied a rotational analogue of the three-sphere swimmer . Avron et al. presented a more efficient swimmer consisting of a pair of bladders which exchange their volume and vary the distance between them . There have been several experimental realisations of artificial low Reynolds number swimmers [15, 16, 17].
When two sperms swim close to each other, their tails beat in synchrony  and Taylor studied this using his waving sheet model with hydrodynamic interactions . Coordinated beating of flagella or cilia is important for a range of processes including motility, efficient pumping of fluid and symmetry breaking in developing embryos [18, 20, 19, 21, 22]. Theoretical and experimental models have been studied to show that synchronization can occur through hydrodynamic interaction and that it is relevant to bacterial swimming and pumping by arrays [24, 25, 26, 23, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44]. Flagellar synchronization is observed in Chlamydomonas, a unicellular green alga that swims using two flagella that beat with a breaststroke pattern [45, 46, 47]. The cell has diameter and swims with velocity so the Reynolds number is and inertia is negligible. During normal swimming, the flagella beat in synchrony. These periods of synchrony are interrupted by periods of asynchronous beating and during these asynchronies, there is a large change in the cell’s orientation [48, 49]. This is analogous to run-and-tumble behaviour observed in bacteria.
Simple models have helped us understand better the intricacies of low Reynolds number swimming and hydrodynamic synchronization, and a recent development has been to combine these two effects in the context of a simple three-sphere model for the swimming of Chlamydomonas [50, 51, 52]. This simple model captures some of the important features of Chlamydomonas, namely, the ability to swim, the exact role of hydrodynamic interactions [50, 52], the existence of stable synchronized states and an emergent run-and-tumble behaviour which is observed when we add white noise to the driving force . Here, we consider the model without added noise and explore its phase diagram and full parameter space to see when the model evolves into the synchronized state. We also investigate the stability of the synchronized state under various conditions, and the different types of behaviour that can be obtained form the model. These studies have revealed a number of intriguing features.
2 The Model
We arrange three beads in the plane, each of radius , on a frictionless scaffold, as shown in figure 1.
We refer to the left, right and back beads with the subscripts ’’, ’’ and ’’, respectively. Let be the origin of the cell reference frame with respect to a lab frame. The cell axes make an angle with the lab axes . The left and right beads model the flagella and move on circular trajectories in the cell frame of radius in opposite directions and with phases and ; the back bead models the cell body and is fixed with respect to the cell frame. The positions and velocities of the beads are
where the unit vectors and are in the normal and tangent directions of the circular trajectory of . The left and right beads are driven by tangential forces and respectively. Normal forces and are exerted by the beads in order to be constrained to the circular trajectories. The force on the back bead is such that the swimmer is force free and torque free:
where for and for .
The forces and velocities are related through hydrodynamic interactions between the beads:
where is the friction coefficient of each bead ( is viscosity of the ambient fluid). In the limit when is small compared with all other length scales, the hydrodynamic interaction is described by the Oseen tensor with .
The phase difference evolves according to
We solve equation 8 numerically for a choice of stroke pattern (driving forces) , . The long term behaviour of depends on the stroke pattern and in many cases the initial condition. We compute and to leading order in , since we do not require hydrodynamic interactions for synchronization , but we include the next order hydrodynamic term when computing the velocities, since we need this second order affect to achieve a net swimming velocity.
3 Swimming velocity in the synchronized state
First we consider the synchronized state where and , so that . We do not worry about the stability of the synchronized state, which we consider in the next section, and assume that the swimmer stays in this state. Since the Reynolds number is low, we need to ask ‘does the model achieve net propulsion?’ If hydrodynamic interactions are not included, then the cell just moves forwards and backwards and there is no net motion. However, if we include hydrodynamic interactions, which vary in strength around the cycle, then the symmetry in the swimming stroke is broken and net propulsion is achieved.
The magnitude and direction of the net swimming velocity depends on the ratios and . Figure 2 shows the net swimming velocity in the direction for a range of and constant driving force . For , swimming is in the positive direction, otherwise the cell swims in the negative direction. Polotzek and Friedrich in reference  give the following explanation of why the cell may swim in either direction. The instantaneous velocity is the ratio of the force , which has to be applied to the back bead to prevent it from moving, and the friction coefficient associated with towing the swimmer in the direction. The force oscillates during a stroke cycle and the hydrodynamic interactions, which reduce the magnitude of , are strongest when the beads are closest together. On the other hand, the friction coefficient is largest when the beads are furthest apart and smallest when the beads are close together. The geometry of the swimmer determines which effect dominates and therefore whether the net swimming is in the positive or negative direction.
Herein we fix the values and , which results in forwards swimming for .
The the net velocity is also affected by the force profile and we consider driving forces such that , where is a fixed average force. The net velocity can be written as , where , is the period, and we can write . In the synchronized state the force dependence cancels in the ratio , so the force only enters the net velocity expression through the term. In order to maximise the net velocity, we must minimise the period , where we can write . Minimising with the constraint tells us that a constant force profile maximises the net velocity. Clearly, increasing increases the net velocity. However as we shall see in the next section, the synchronized state is not stable when we choose a constant force profile. Friedrich et al. showed that a constant driving force can give a stable synchronized state if we change the direction of rotation of the beads, so [50, 52], which is equivalent to changing the sign of , but here we choose to work with and .
4 Synchronization and stability
We consider force profiles of the forms and , where and . The definition of the synchronized state that we use here is zero (or integer multiple of ) phase difference between the two flagella, i.e. when , . Initially we tried to analyst the the synchronization stability by linear stability analysis, but we are unable to do this because we cannot perform a valid Taylor expansion when , where . We work numerically to avoid this Taylor expansion; for further details see the appendix.
We identify five main types of stability of the synchronized state by looking at the evolution of from different initial conditions for a number of equally spaced in the range : (i) All the initial conditions evolve into the synchronized state for all (the synchronized state is stable). (ii) Some choices of evolve into the synchronized state and others choices evolve into an oscillating state, but there is a larger number of that lead to synchronization than the number of that leads to oscillations. (iii) Some choices of evolve into the synchronized state and others choices evolve into an oscillating state; the numbers of that lead to each type of behaviour are similar. (iv) Some choices of evolve into the synchronized state and others choices evolve into an oscillating state, but there is a larger number of that lead to oscillations than the number of that leads to synchronization. (v) All choices of evolve into the oscillating state, (the synchronized state is unstable).
Although the choice of initial condition is arbitrary, we want to know how likely it is that a small perturbation from the synchronized state will decay back to synchronization, or whether it is likely to evolve into an oscillating state. This choice of is suitable for this purpose. For type (v) stability, if we start in the synchronized state, then it is likely that some numerical noise will kick into an oscillating state. It is possible for a small amount of noise to kick the synchronized state into an oscillating state for types (ii), (iii), (iv), with low probability for type (ii), then increasing probability for type (iii) and then type (iv).
4.1 Equal beat patterns
First we consider the case where . For each choice of coefficient and initial condition, evolves either to an integer multiple of and remains at this value (synchronization); or it reaches a state where it oscillates about sinusoidally; or it reaches a periodic state near zero, but never reaching zero. Figure 3 shows examples of these three cases and the corresponding orientation .
When oscillates about , then the orientation oscillates about some fixed value. The cycle averaged motion is in a straight line, but the cell jiggles from side to side as well as backwards and forwards as it moves along. When and the oscillations are near zero, there is a net drift in the orientation so the net motion of the cell is along a curved trajectory.
For many choices of , the initial condition determines whether evolves into the synchronized state or the oscillating state. Figure 4 shows the the dependence of evolution on initial condition for . A 63 69 grid is shown where each square represents an initial condition . A black square represents an initial condition for which after a sufficiently long time; a white square represents an initial condition for which continues to oscillate periodically as .
For , all initial conditions lead to an oscillating state and the synchronized state is unstable (type (v)). The black squares in figure 4 are initially in the synchronized state. Many squares along the line are white because a small amount of numerical noise drives the system away from the synchronized state. For , initial conditions close to lead to the synchronized state, but initial conditions far from lead to an oscillating state. The synchronized state is stable (type (i)). For , most initial conditions close to lead to the synchronized state, but a few initial conditions close to lead to an oscillating state and we have type (ii) stability; if starts close to the synchronized state, it is likely to evolve into the synchronized state and it will stay in the synchronized state if there is no noise, but it is also possible for the cell to start close to the synchronized state and move away into an oscillating state. In an oscillating state the cell can still swim, but there will be more side to side movement. There appear to be a few white squares on the line , however this is because the grid does not lie exactly on the line, the grid contains points and and this small deviation from the synchronized state is enough for evolution into an oscillating state for a few choices of .
Similar behaviour is observed for other harmonics, but with different ranges of coefficients giving the different stability types for the synchronized state. For example, figure 5 shows the 4th harmonic for three different coefficients of cosine and initial condition . Figure 5(a) shows that oscillate about for , oscillations are close to zero for (it is interesting to note the 4 peaks in every cycle), and evolves into the synchronized state for .
Figure 6 shows the stability of the synchronized state for the first 10 harmonics for a discrete range of coefficients when we choose the cosine term, i.e. .
We see that there are more type (i) and (ii) force profiles for negative coefficients than for positive, showing that we are more likely to end up in the synchronized state when the coefficient is negative than when the coefficient is positive, for an arbitrary choice of harmonic. We see that for there is no type (i) behaviour, so we cannot guarantee that we will reach the synchronized states for these higher harmonics. The unstable type (v) band around gets narrower as increases, so we only need weak phase dependence to reach a synchronized state for higher harmonics, but the region of initial conditions which does lead to the synchronized state is very small. Even after reaching the synchronized state, it is likely that noise will drive away from the synchronized state. Usually the stability moves towards the lower types of stability (more stable types) as increases for each harmonic, although there are some exceptions. For example, when we can start in a type (i) region, then as increases we move into a type (ii) region. We also see a type (iii) region surrounded by type (iv) regions on both sides for and .
For each type of stability shown in figure 6, we select an arbitrary force profile and show the full initial condition phase diagram in figure 7. In the type (i) stable case, we see that oscillating states can still evolve (see also figure 4), but the initial conditions for an oscillating state are not close to the synchronized case. In a few cases when we replace the cosine with sine, all initial conditions lead to the synchronized state, for example, the force profile for , except for a very thin dotted white curve through the middle of the phase diagram indicating an unstable oscillating state. However, if we look at this oscillating state for long enough, after some time the numerical noise causes it to move into the synchronized state.
We see from the bright green band down the center of figure 6 that we need some form of phase dependence in order to achieve synchronization. If we choose a constant driving force then evolves into an oscillating state for all initial conditions. The synchronized state is unstable, so even if we start in the synchronized state, a small amount of numerical noise can drive the system into the oscillating state. Friedrich et al. showed that synchronization can occur with constant forcing when the direction of rotation of the beads is reversed (equivalent to and reversing swimming direction) but here we focus on the case and .
The inspiration for our run-and-tumble model in reference  came from the initial condition phase diagrams. To see run-and-tumble we want to start in a stable synchronized state, then allow noise to move us temporarily into a white region of the phase diagram, before moving back into a black region. In reference , we allowed noise to vary the coefficients , so that a black square in the noiseless phase diagram can change to a white square when the instantaneous effect of the noise changes the value of , then changes back to a black square when the instantaneous effect of the noise is smaller. This allows us to start in the synchronized state (with fluctuations due to the noise), then move away from the synchronized state when the instantaneous noise is large and we are at a suitable point in the phase diagram, then move back into the fluctuating synchronized state when the instantaneous noise is small. In reference , we chose to work with the first harmonic and use , where is the noise term. We chose the value 0.7 because it is in the stability type (i) region, but lies close to the type (ii) region. In the fluctuating synchronized state, when we move through the values of which are surrounded by white squares in the type (ii) phase diagram, there is the possibility to move into an oscillating state, but there are still plenty of black squares surrounding the line , so we can have long periods in the run phase.
The noise causes fluctuations of the position in the phase diagram, and this could also be a cause of run-and-tumble behaviour. For example, consider the phase diagram in figure 7(a). If we are in the synchronized state and noise is small enough such that fluctuations in are within the black region, then tumbles will not occur. If the noise is larger, so fluctuations move into the white region, then the cell could begin to move into the oscillating state, and after a few oscillations noise could kick the oscillations into the black region and the cell would move back towards a synchronized state. Elsewhere, we will consider the effects of adding noise to and , without any noise in the coefficient, to see if we obtain run-and-tumble behaviour this way and compare the statistics to the run-and-tumble obtained when noise is added to the coefficient.
4.2 Mismatched coefficients
If we start in the synchronized state, then
which is non-zero for , so the system does not stay in the synchronized state. We focus on the first harmonic and consider the case and where (and we have dropped the upper index on the coefficient). When and , synchronization is frustrated and oscillates about , , shown in figure 8(a) for . The orientation of the cell drifts, shown in figure 8(b) so the cell swims along a curved trajectory. For , the centre of oscillations drifts away from the synchronized state, but remains close to . Swapping the signs of the coefficients swaps the direction of the orientation drift.
When the coefficients have the same sign but different magnitudes, oscillates about and the orientation oscillates about some fixed value. Figure 8(c) shows the evolution of for . This type of behaviour is also seen for some choices of equal coefficients, for example, in figure 3(e), (f) where . Choices of coefficients which give this type of behaviour can be used to model a mutant of Chlamydomonas which swims with antiphase synchrony . When the model swims with antiphase beating, the beat frequency is higher than when it swims with in phase beating, which has been observed in real Chlamydomonas cells .
When the coefficients have opposite signs and different magnitudes, there are two main types of behaviour that occur. For and or and , then oscillates periodically and the corresponding orientation drifts in the negative direction in the former case and in the positive direction in the latter case. This is shown in figure 8(e), (f) for . For and or and , then oscillations in phase difference, , drift in the negative (positive) direction and the orientation drifts in the positive (negative) direction in the former (latter) case. Figure 9 shows examples of this this case where one bead completes more cycles than the other bead. There is a difference in the mean angular velocity of the two beads without choosing a different for each bead.
4.3 Combinations of harmonics
So far we have considered force profiles with only one harmonic term. Now we consider driving forces with contributions from two harmonics. For simplicity we choose equal profiles for the left and right beads, , and of the form , and with the ’s chosen such that for all real . Figure 10 shows the stability of the synchronized state for and . Each grid square represents a choice of driving force with coefficients , and the colour represents the stability of the synchronized state for that particular driving force.
We see that there are large regions for which the synchronized state is type (i) stable.
This simple mechanical model is able to evolve into a stable synchronized state for certain choices of parameters in the driving force when the initial condition is within some region of the synchronized state. We do not need hydrodynamic interactions to achieve stable synchronization; we include hydrodynamic friction on each bead, with force free and torque free conditions and a phase dependent driving force, which can be constructed with a suitable combination of harmonic terms. For many choices of force profile, some initial conditions allow the model to evolve into the synchronized state, while other initial conditions that are very close to the synchronized state lead to an oscillating state. There are some force profiles, including constant forcing, where there are no initial conditions that evolve into the synchronized state, and if the system starts in the synchronized state when the driving forces are equal, even a small amount of numerical noise can drive the system into an oscillating state. There are different types of periodic behaviour for different choices of parameters in the driving force; often the phase difference oscillates about , but sometimes the oscillations can occur close to zero with multiple peaks per cycle.
When the parameter in the driving force is different for the left and right beads, we can get periodic oscillating states about a range of values, or we can get a drifting oscillating state, where one bead has a higher average angular velocity than the other. When the coefficients have equal magnitude and opposite sign, this can lead to oscillations about the synchronized state. This frustrated synchronization is interesting when we add intrinsic noise to the driving force, because then the behaviour of is very similar for both opposite coefficients and equal coefficients, although the behaviour in the orientation is different in the two cases.
The nonlinear mechanics of the system make it difficult to study analytically and it is not easy to predict the parameter ranges which give stable synchronization. Our numerical results have highlighted some of the main types of behaviour of the model.
An important feature is that it is necessary to have some sort of phase dependent driving force in order to have a stable type (i) or type (ii) synchronized state. When the phase dependence is only weak, then the synchronized state is unstable, which we see from figure 6 when the coefficient is small. The value of the coefficient at which the phase dependence becomes strong enough to give synchronization depends on the harmonic, whether we choose a positive or negative coefficient, and whether we choose sine or cosine. These latter choices are equivalent to adding a constant phase , or in the harmonic term.
This simple mechanical model shows a wide range of behaviour when we vary the parameter in the driving forces. The variety of stabilities suggests possibilities for developing run-and-tumble models, where noise can be used to to jump between regions of a phase diagram that lead to synchronization or oscillations, or jump between phase diagrams as we did in reference .
We note that linear stability analysis can not be used to probe the stability of the synchronization, as we cannot perform a valid Taylor expansion when , where . For equal driving force profiles that we linearize, , , if we were to Taylor expand, then the linearized expression for would be
which has a singularity at . The apparent singularity actually occurs at and at in the full expression, but the choice of constraining force ensures this zero in the denominator is canceled by the numerator. However, when we expand in and shift the singularity so that it occurs at , then the numerator is no longer zero at this point. The reason we have this zero in the denominator is the following: The torque free condition (4) is
Along with equations (1-3,5-7), we use (11) to solve for the constraining forces , . However, at , is multiplied by a term which vanishes, so the torque free condition can be satisfied without specifying . We over-constrain the system when we divide by zero and specify at . Geometrically, corresponds to the phase where is parallel to . Our numerical analysis of the full expression avoids this singularity.
- Purcell E M 1977 Am. J. Phys. 45 3
- Happel J and Brenner H 1973 Low Reynolds Number Hydrodynamics (Dordrecht; London: Kluwer Academic)
- Lauga E and Powers T R 2009 Rep. Prog. Phys. 72 096601
- Childress S 1981 Mechanics of Swimming and Flying (Cambridge: Cambridge University Press)
- Shapere A and Wilczek F 1987 Phys. Rev. Lett. 58 2051
- Taylor G I 1951 Proc. R. Soc. London, Ser. A 209 447
- Lighthill M J 1976 SIAM Rev. 18 161
- Alberts B, Johnson A, Lewis J, Raff M, Roberts K and Walter P 2008 Molecular Biology of the Cell (New York; Abingdon: Garland Science) 5th ed
- Blake J R and Sleigh M A, 1974 Biol. Rev. 49 85
- Lighthill M J 1952 Commun. Pure Appl. Math. 5 109
- Becker L E, Koehler S A and Stone H A 2002 J. Fluid Mech. 490 15 (2002)
- Najafi A and Golestanian R 2004 Phys. Rev. E 69 062901
- Dreyfus R, Baudry J and Stone H A 2005 Eur. Phys. J. B 47 161
- Avron J E, Kenneth O and Oaknin D H 2005 New. J. Phys. 7 234
- Dreyfus R, Baudry J, Roper M L, Fermigier M, Stone H A and Bibette J 2005 Nature 437 862
- Tierno P, Golestanian R, Pagonabarraga I and Sagués F 2008 Phys. Rev. Lett. 101 218304
- Leoni M, Kotar J, Bassetti B, Cicuta P and Cosentino Lagomarsino M 2009 Soft Matter 5 427
- Gray J 1928 Ciliary Movement (Cambridge: Cambridge University Press)
- Berg H C 2004 E. coli in Motion (New York: Springer-Verlag)
- Gueron S and Levit-Gurevich K 1999 Proc. Natl. Acad. Sci. U.S.A. 96 12240
- Osterman N and Vilfan A 2011 Proc. Natl. Acad. Sci. U.S.A. 108 15727
- Nonaka S, Yoshiba S, Watanabe D, Ikeuchi S, Goto T, Marshall W F and Hamada H 2005 PLoS Biol. 3 1467
- Golestanian R, Yeomans J M and Uchida N 2011 Soft Matter 7 3074
- Cosentino Lagomarsino M, Bassetti B and Jona P 2002 Eur. Phys. J. B 26 81
- Cosentino Lagomarsino M, Jona P and Bassetti B 2003 Phys. Rev. E 68 021908
- Kim M, Bird J, Van Parys A J, Breuer K S and Powers T R 2003 Proc. Natl. Acad. Sci. U.S.A. 100 15481
- Kim M and Powers T R 2004 Phys. Rev. E 69 061910
- Reichert M and Stark H 2005 Eur. Phys. J. E 17 493
- Kim Y W and Netz R R 2006 Phys. Rev. Lett. 96 158101
- Vilfan A and Jülicher F 2006 Phys. Rev. Lett. 96 058102
- Lenz P and Ryskin A 2006 Phys. Biol. 3 285
- Guirao B and Joanny J F 2007 Biophys. J. 92 1900
- Niedermayer T, Eckhardt B and Lenz P 2008 Chaos 18 037128
- Qian B, Jiang H, Gagnon D A, Breuer K S and Powers T R 2009 Phys. Rev. E 80 061919
- Elfring G J and Lauga E 2009 Phys. Rev. Lett. 103 088101
- Uchida N and Golestanian R 2010 Phys. Rev. Lett. 104 178103
- Uchida N and Golestanian R 2010 Europhys. Lett. 89 50011
- Kotar J, Leoni M, Bassetti B, Cosentino Lagomarsino M and Cicuta P 2010 Proc. Natl. Acad. Sci. U.S.A. 107 7669
- Golestanian R and Uchida N 2011 Phys. Rev. Lett. 106 058104
- Wollina C and Stark H 2011 Eur. Phys. J. E 34 42
- Bruot N, Kotar J, de Lillo F, Cosentino Lagomarsino M and Cicuta P 2012 Phys. Rev. Lett. 109 164103
- Di Leonardo R, Búzás A, Kelemen L, Vizsnyiczai G, Oroszi L and Ormos P 2012 Phys. Rev. Lett. 109 034104
- Lhermerout R, Bruot N, Cicuta G M, Kotar J and Cicuta P 2012 New J. Phys. 14 105023
- Golestanian R and Uchida N 2012 Eur. Phys. J. E 35 135
- Ringo D L 1967 J. Cell. Biol. 33 543
- Rüffer U and Nultsch W 1985 Cell Mot. 5 251
- Rüffer U and Nultsch W 1998 Cell Mot. 41 297
- Polin M, Tuval I, Drescher K, Gollub J P and Goldstein R E 2009 Science 325 487
- Goldstein R E, Polin M and Tuval I 2009 Phys. Rev. Lett. 103 168103
- Friedrich B M and Jülicher F 2012 Phys. Rev. Lett. 109 138102
- Bennett R R and Golestanian R 2013 Phys. Rev. Lett. 110 148102
- Polotzek K and Friedrich B M 2012 arXiv:1211.598 to appear in New J. Phys.
- Oseen C V 1927 Neue Methoden und Ergebnisse in der Hydrodynamik (Leipzig 1927).
- Goldstein R E, private communication.