Phase Dependent Forcing and Synchronization in the three-sphere model of Chlamydomonas

The green alga {\it 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 {\it 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.


Introduction
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 [5]. 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 [6]. Purcell showed that a swimmer needs at least two compact degrees of freedom to break the time reversal symmetry and achieve net propulsion [1].
Many microorganisms swim using flagella [7]; 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 [8]. 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 [9].
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 [10]. 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 nonreciprocal pattern [12]. Dreyfus et al. studied a rotational analogue of the three-sphere swimmer [13]. Avron et al. presented a more efficient swimmer consisting of a pair of bladders which exchange their volume and vary the distance between them [14]. 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 [18] and Taylor studied this using his waving sheet model with hydrodynamic interactions [6]. 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][19][20][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 . 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 ∼ 10µm and swims with velocity 100µm/s so the Reynolds number is Re ∼ 10 −3 and inertia is negligible. Figure 1. A three sphere model of Chlamydomonas. The left and right beads represent the flagella and move on circular trajectories in the cell frame. The back bead represents the cell body. The scaffold is shown to define the reference frame of the cell which has origin R 0 in a lab frame, but it does not interact with the fluid. The green underlay is a schematic of a Chlamydomonas cell.
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 [51]. 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.

The Model
We arrange three beads in the x − y plane, each of radius a, on a frictionless scaffold, as shown in figure 1. We refer to the left, right and back beads with the subscripts 'l', 'r' and 'b', respectively. Let R 0 be the origin of the cell reference frame with respect to a lab frame. The cell axesx,ŷ make an angle θ(t) with the lab axesX,Ŷ. The left and right beads model the flagella and move on circular trajectories in the cell frame of radius b in opposite directions and with phases φ l and φ r ; 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 vectorsn i andt i are in the normal and tangent directions of the circular trajectory of R i . The left and right beads are driven by tangential forces F t l and F t r respectively. Normal forces F n l and F n r 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 F i = F t it i + F n in i for i = l, r and T j = R j × F j for j = b, l, r. The forces and velocities are related through hydrodynamic interactions between the beads:Ṙ where ξ = 6πηa is the friction coefficient of each bead (η is viscosity of the ambient fluid).
In the limit when a is small compared with all other length scales, the hydrodynamic interaction is described by the Oseen tensor G ij = 1 8πη|r ij | (I+r ijrij ) with r ij = r i −r j [53]. The phase difference δ = φ r − φ l evolves according tȯ where φ = (φ l + φ r )/2 anḋ We solve equation 8 numerically for a choice of stroke pattern (driving forces) F t i (φ i ), i = l, r. The long term behaviour of δ depends on the stroke pattern and in many cases the initial condition. We computeδ,θ andφ to leading order in a/L, since we do not require hydrodynamic interactions for synchronization [50], 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.

Swimming velocity in the synchronized state
First we consider the synchronized state where F t l (φ) = F t r (φ) and δ = 0, so that φ l = φ r = φ. 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 H/b and L/b. Figure 2 shows the net swimming velocity in theŷ direction for a range of (H/b, L/b) and constant driving force F t (φ) = F 0 . For L 2.25|H|, swimming is in the positive direction, otherwise the cell swims in the negative direction. Polotzek and Friedrich in reference [52] give the following explanation of why the cell may swim in either direction. The instantaneous velocity v = f /g is the ratio of the force −fŷ, which has to be applied to the back bead to prevent it from moving, and the friction coefficient g associated with towing the swimmer in theŷ direction. The force f oscillates during a stroke cycle and the hydrodynamic interactions, which reduce the magnitude of f , 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 a/L = 1/33, H/b = 1/5 and L/b = 2, which results in forwards swimming for F 0 > 0.
The the net velocity is also affected by the force profile and we consider driving forces F t (φ) such that 2π 0 dφF t (φ) = 2πF 0 , where F 0 is a fixed average force. The net velocityv can be written asv = 1/T In the synchronized state the force dependence cancels in the ratioṘ b /φ, so the force only enters the net velocity expression through the 1/T term. In order to maximise the net velocity, we must minimise the period Minimising T with the constraint 2π 0 dφF t (φ) = 2πF 0 tells us that a constant force profile F t (φ) = F 0 maximises the net velocity. Clearly, increasing F 0 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 F 0 < 0 [50,52], which is equivalent to changing the sign of H, but here we choose to work with F 0 > 0 and H > 0. 3.

Synchronization and stability
We consider force profiles of the forms The definition of the synchronized state that we use here is zero (or integer multiple of 2π) phase difference between the two flagella, i.e. when δ = 2πn, n ∈ Z. 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 φ = φ s , where H cos φ s + L sin φ s = 0. 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 δ(φ 0 ) = 0.1 for a number of φ 0 equally spaced in the range φ 0 ∈ [0, 2π): (i) All the initial conditions δ(φ 0 ) = 0.1 evolve into the synchronized state for all φ 0 (the synchronized state is stable). (ii) Some choices of φ 0 evolve into the synchronized state and others choices evolve into an oscillating state, but there is a larger number of φ 0 that lead to synchronization than the number of φ 0 that leads to oscillations. (iii) Some choices of φ 0 evolve into the synchronized state and others choices evolve into an oscillating state; the numbers of φ 0 that lead to each type of behaviour are similar. (iv) Some choices of φ 0 evolve into the synchronized state and others choices evolve into an oscillating state, but there is a larger number of φ 0 that show a small part of the plot in more detail. The bottom axis φ/2π shows the number of cycles which increases monotonically with time. Bottom row: corresponding θ(φ).
lead to oscillations than the number of φ 0 that leads to synchronization. (v) All choices of φ 0 evolve into the oscillating state, (the synchronized state is unstable).
Although the choice of initial condition δ 0 = 0.1 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 δ 0 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).

Equal beat patterns
First we consider the case where F t l (φ) = F t r (φ) = F t (φ). For each choice of coefficient and initial condition, δ evolves either to an integer multiple of 2π 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 b (1) = 0.4 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 F t (φ), the initial condition determines whether δ evolves into the synchronized state or the oscillating state. Figure 4 shows the the dependence of δ φ 0 2π evolution on initial condition for F t (φ) = F 0 (1 + a (1) cos φ). A 63 × 69 grid is shown where each square represents an initial condition (φ 0 , δ 0 ). A black square represents an initial condition for which δ → 0 after a sufficiently long time; a white square represents an initial condition for which δ continues to oscillate periodically as t → ∞. For a (1) = 0.4, 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 δ = 0 are white because a small amount of numerical noise drives the system away from the synchronized state. For a (1) = 0.6, initial conditions close to δ = 0 lead to the synchronized state, but initial conditions far from δ = 0 lead to an oscillating state. The synchronized state is stable (type (i)). For a (1) = 0.8, most initial conditions close to δ = 0 lead to the synchronized state, but a few initial conditions close to δ = 0 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 δ 0 = 2π, however this is because the grid does not lie exactly on the 2π line, the grid contains points δ = 6.2 and δ = 6.3 and this small deviation from the synchronized state δ = 2π is enough for evolution into an oscillating state for a few choices of φ 0 .
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 δ(φ = 0) = 0.1. Figure 5(a) shows that δ oscillate about π for a (4) = 0.4, oscillations are close to zero for a (4) = 0.8 (it is interesting to note the 4 peaks in every cycle), and δ evolves into the synchronized state for a (4) = 0.9. Figure 6 shows n a (n) Figure 6. Stabilities for range of coefficients and harmonics. Each horizontal band represents a harmonic, which is broken up into 49 blocks and each block represents a coefficient a (n) . Orange blocks represent stability (i); dark blue blocks represent stability (ii); pale blue blocks represent stability (iii); pale green blocks represent stability (iv); and bright green blocks represent stability (v).
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. F t i (φ i ) = F 0 (1 + a (n) cos nφ i ). 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 n > 2 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 a (n) = 0 gets narrower as n 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 |a (n) | increases for each harmonic, although there are some exceptions. For example, when n = 1 we can start in a type (i) region, then as |a (1) | increases we move into a type (ii) region. We also see a type (iii) region surrounded by type (iv) regions on both sides for n = 3 and a (3) > 0. 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 F t (φ) = F 0 (1 + b (1) sin φ) for b (1) 0.6, 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 H → −H and reversing swimming direction) but here we focus on the case F t i > 0 and H > 0. The inspiration for our run-and-tumble model in reference [51] 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 [51], we allowed noise to vary the coefficients a (1) i , 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 a (1) i , 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 [51], we chose to work with the first harmonic and use a (1) = 0.7 + ζ, 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 δ = 0, 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.

Mismatched coefficients
If we start in the synchronized state, then  which is non-zero for F t l = F t r , so the system does not stay in the synchronized state. We focus on the first harmonic and consider the case F t l = F 0 (1 + a l cos φ l ) and F t r = F 0 (1 + a r cos φ r ) where a l = a r (and we have dropped the upper index on the coefficient). When a l = −a r and |a i | 0.6, synchronization is frustrated and δ oscillates about 2πn, n ∈ Z, shown in figure 8(a) for −a l = a r = 0.7. The orientation of the cell drifts, shown in figure 8(b) so the cell swims along a curved trajectory. For |a i | < 0.6, the centre of δ oscillations drifts away from the synchronized state, but remains close to 2πn. 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 (2n + 1)π and the orientation oscillates about some fixed value. Figure 8(c) shows the evolution of δ for a l = 0.7, a r = 0.8. This type of behaviour is also seen for some choices of equal coefficients, for example, in figure 3 Choices of coefficients which give this type of behaviour can be used to model a mutant of Chlamydomonas which swims with antiphase synchrony [54]. 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 [54].
When the coefficients have opposite signs and different magnitudes, there are two main types of behaviour that occur. For |a l | < |a r | and sgn(a l ) = −1 or |a r | < |a l | and sgn(a r ) = −1, 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 a l = −0.7, a r = 0.8. For |a l | < |a r | and sgn(a l ) = +1 or |a r | < |a l | and sgn(a r ) = +1, 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 F 0 for each bead.

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, F t l (φ) = F t r (φ) = F t (φ), and of the form F t (φ) = F 0 (1 + a (m) cos mφ + a (n) cos nφ), m = n and with the a (i) 's chosen such that F t (φ) > 0 for all real φ. Figure 10 shows the stability of the synchronized state for m = 1, n = 2 and m = 1, n = 3. Each grid square represents a choice of driving force with coefficients (a (1) , a (j) ), j = 2, 3 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.

Conclusion
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 ). The colour of each square represents the following stability: Orange represents type (i); dark blue represents type (ii); pale blue represents type (iii); pale green represents type (iv); bright green represents type (v).
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 0, π or ±π/2 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 [51]. the phase where R i − R b is parallel ton i . Our numerical analysis of the full expression avoids this singularity.