## Abstract

Despite evidence for a hydrodynamic origin of flagellar synchronization between different eukaryotic cells, recent experiments have shown that in single multi-flagellated organisms, coordination hinges instead on direct basal body connections. The mechanism by which these connections lead to coordination, however, is currently not understood. Here, we focus on the model biflagellate *Chlamydomonas reinhardtii*, and propose a minimal model for the synchronization of its two flagella as a result of both hydrodynamic and direct mechanical coupling. A spectrum of different types of coordination can be selected, depending on small changes in the stiffness of intracellular couplings. These include prolonged in-phase and anti-phase synchronization, as well as a range of multi-stable states induced by spontaneous symmetry breaking of the system. Linking synchrony to intracellular stiffness could lead to the use of flagellar dynamics as a probe for the mechanical state of the cell.

## 1. Introduction

Cilia and flagella are structurally identical, whip-like cellular organelles employed by most eukaryotes for tasks ranging from sensing and locomotion of single cells [1], to directing embryonic development [2] and driving cerebrospinal fluid flow [3] in animals. Originally observed in 1677 by the Dutch pioneer Antonie van Leeuwenhoek, groups of motile cilia and flagella have a seemingly spontaneous tendency to coordinate their beating motion and generate large-scale patterns known as metachronal waves [4]. Coordination has often been proposed to provide an evolutionary advantage by improving transport and feeding efficiency [5–10], although estimates of the magnitude of this effect are notoriously difficult. Despite the uncertainty on its biological role, however, the universality of flagellar coordination is an empirical fact, and it suggests the existence of a correspondingly general mechanism for synchronization. Mechanical forces, transmitted either by the surrounding fluid or internally through the cells, have often been proposed as responsible for this coordination [11–15]. Understanding how synchronization emerges could therefore highlight novel and potentially subtle roles played by physical forces in cell biology. Here, we develop a minimal model that links small changes in the mechanical properties of cells with the dynamics of their protruding flagella. In turn, this approach could lead to coordination being used as a probe to measure the internal mechanical state of a cell.

Reports of coordinated motion in nearby swimming sperm [16,17] hint at the importance of hydrodynamic coupling. Hydrodynamic-led coordination of self-sustained oscillators, mimicking the active motion of cilia and flagella [18,19], has been extensively investigated theoretically [12,20–22], numerically [23–25], and experimentally with colloidal rotors [26,27] and rowers [28,29]. Despite the peculiar constraints of low-Reynolds number hydrodynamics, these studies suggest that hydrodynamic interactions can lead to ciliary coordination when coupled to either a phase-dependent driving force [12,21], or axonemal elasticity [30]. Indeed, hydrodynamic-mediated synchronization has been confirmed experimentally between pairs of eukaryotic flagella from different somatic cells from the green alga *Volvox carteri* [14].

Inter-cellular coordination of flagella, however, does not necessarily imply intracellular coordination, and therefore it is not *a priori* clear whether hydrodynamic coupling is also responsible for the synchronization observed in individual multi-flagellated cells. Experimental studies of flagellar coordination within a single cell have focused mainly on the biflagellate green alga *Chlamydomonas reinhardtii* (CR) [31,32], whose flagella are usually locked in a characteristic in-phase breaststroke motion. Early studies of flagellar coordination in CR [33–35] were recently refined [11,36] and extended [37] using microfluidic devices and high-speed imaging of flagella [38]. These pointed at a fundamentally hydrodynamic origin for the observed synchronization, either through direct coupling or via a mechanism based on cell-body rocking [13]. However, a series of elegant novel experiments in CR and other flagellates challenged this view convincingly, showing instead that coordination requires the intracellular striated fibres that join flagellar basal bodies [39,40]. For the CR mutant *vfl3*, with impaired mechanical coupling and a variable number of flagella, synchronization is completely disrupted [39] except for sporadic periods of synchrony in cells with two close flagella oriented in the same direction (see fig. S1 in [40]). Even though the precise mechanism by which direct connections affect flagellar coordination remains to be clarified [15], the spontaneous transitions between extended in-phase (IP) and anti-phase (AP) beating in the CR mutant *ptx1* [41] already suggest that multiple synchronization states should be achievable through changes in the fibres’ mechanical properties within the physiological range.

Here we propose a minimal model for flagellar dynamics for CR which can sustain both stable IP and stable AP states even in the absence of hydrodynamic coupling. Within this framework, the phase dynamics are determined principally by the mechanical state of the basal body fibres [42], with both types of coordination possible within a physiological range of fibre stiffnesses [43]. The inclusion of hydrodynamic coupling leads to the emergence of a region in parameter space where non-trivial states can emerge as a result of spontaneous symmetry breaking through pitchfork bifurcations of limit cycles.

## 2. Minimal model and leading order dynamics

Figure 1*a* summarizes our minimal model of flagella coupled through hydrodynamic and basal body interactions. Following previous theoretical work [20,44], and experimental measurements of flagellar flow fields and waveform elasticity [14], the two flagella of *Chlamydomonas* are represented here by two spheres of radius *a*, immersed in a three-dimensional fluid of viscosity *μ*, and driven around circular orbits of radii () by constant tangential driving forces . Springs of stiffness λ resist radial excursions from the equilibrium length *R*_{0}; and the magnitude of the internal driving forces guarantees that, when isolated, the *i*th oscillator will rotate at the intrinsic angular speed . The orbits are centred along the *x*-axis, a distance *l* apart, and lie along the *xy* plane. Polar and radial coordinates define the oscillators’ instantaneous positions around the centres of their respective orbits. We will refer to as the external rotor radii.

Both hydrodynamic and direct elastic interactions couple these minimal cilia. Hydrodynamic coupling is mediated by the fluid disturbance generated by each sphere’s motion, modelled here as the flow from a point force. These interactions affect the instantaneous angular speeds both directly, through a hydrodynamic torque, and indirectly by modifying the orbits’ radii. For counter-rotating oscillators like those describing *Chlamydomonas* flagella, the resulting effective coupling will promote AP synchronization [15,41,44,45] (figure 1*b*). In addition to external flagellar interactions, considerable evidence [39,40] suggests that flagellar dynamics are strongly influenced by direct intracellular mechanical coupling, through striated fibres joining the basal bodies [33,46] that can lead the system to IP synchrony (figure 1*c*). Intracellular connections are modelled here by introducing, for each oscillator, an auxiliary arm of fixed length at an angle *θ* ahead of the rotating sphere. The endpoints of these arms (red spheres in figure 1*a*) are coupled via an *anisotropic* elastic medium acting as elastic springs of stiffnesses and equilibrium lengths in the *x* and *y* directions, respectively. This is intended to represent the intrinsically anisotropic structure of the fibre bundles connecting the basal bodies [42]. The equations of motion, derived in electronic supplementary material, S1, follow from the requirement of zero net force and torque on each oscillator, and in the limit can be approximated as (see electronic supplementary material, S1)
2.1where ( as ); is the viscous drag coefficient of the rotating sphere; and
2.2In order to model the configuration typical of *Chlamydomonas*, we will focus here on identical but counter-rotating oscillators, . Parameter values are given in table 1 unless otherwise specified.

## 3. Fibres-only coupling

Despite the apparent simplicity, this minimal system displays rich dynamics, and it is therefore convenient to analyse its behaviour following steps of increasing complexity. Let us begin by considering the case in which hydrodynamic coupling is completely neglected. The hydrodynamic drag is still necessary for each model cilium (through ) in order to balance the driving force, but there is no direct hydrodynamic coupling between the spheres. In this case, , and recasting the angular dynamics in terms of phase sum and difference , we obtain
3.1Requiring that the maximal torque exerted on each oscillator by the internal springs is always smaller than that from the driving force () guarantees that the cilia will always be beating (), and defines a physiologically plausible range for *k*’s, which in our case is . At the same time, unless , the system will monotonically converge to either IP () or AP () synchrony, depending on whether *k*_{y} is larger or smaller than . Figure 2*a*–*c* shows the convergence to either state, for a set of initial conditions and increasing internal stiffnesses. When , the phase difference evolves much faster than the sum, and follows the approximate dynamics
3.2as shown in figure 2*d*. The instantaneous and average phase speed profile can be solved analytically (see electronic supplementary material, S3). In the case of , the model predicts a lower average phase speed in IP than in AP which is in qualitative agreement with experimental observations [41]. Recent studies argue that hydrodynamics plays a negligible role in flagellar synchronization for single cells [39,40]. Without hydrodynamics, our model predicts that the effective interflagellar coupling should be given by . Using the known value for *Chlamydomonas*, [11,41], we obtain . The Young’s modulus for the bundle of striated fibres can then be estimated as , where *L* and *A* are the length and cross-sectional area of the bundle, respectively. Taking and results in the value [42,47]. This is a biologically plausible estimate, midway between the elastic modulus of relaxed skeletal muscle () and elastin () [43].

## 4. Stiff flagella hydrodynamics

We begin now to include the effect of hydrodynamic interactions under the assumption of artificially stiff flagella, implemented in this section by increasing the radial spring stiffness to (10× the value in table 1). Increasing λ reduces the typical response time of the radial coordinate () and, when this is much smaller than the angular timescale (), it allows us to simplify the dynamics by assuming an instantaneous radial response [44]. Then, to first order in the small parameter *ρ*, equations (2.1) imply that will obey
4.1once the dynamics have been averaged over the fast variable . The coefficient of above provides an intuitive understanding of the roles of basal body coupling and elasto-hydrodynamic interactions in determining the synchronization state of the system, which is useful as a general rule of thumb to assess the relative importance of these phenomena even beyond the ‘stiff flagella’ case. Within this approximation, hydrodynamics appears to simply shift the location of the transition between AP and IP states, determined by the steady state time average , from to . This is indeed confirmed by simulations of the full system for large *l* (see figure 4*b*).

A closer inspection, however, reveals more intriguing dynamics which become particularly evident for small interflagellar separations. Figure 3*a*–*c* shows the steady state dynamics of *σ* for , as *k _{y}* is swept across the transition. Between AP and IP coordination (figure 3

*a*,

*c*) there is a distinct region of

*k*values for which the system synchronizes in a non-trivial intermediate state (figure 3

_{y}*b*). This is accompanied by a permanent difference in the average values of the oscillators’ radii (figure 3

*e*), with the asymmetry depending on which of the equally probable signs of is chosen by the system. Figure 3

*g*shows the full positive branch of as

*k*is swept between AP and IP values (simulations: green solid line). This can be compared to the leading order behaviour with and without hydrodynamics (black dashed and solid lines); and the one predicted by refining equation (4.1) to next-to-leading order in (black dotted line; see electronic supplementary material, S2) 4.2The semi-quantitative agreement between the simulated and predicted dependence of the steady state extends also to the

_{y}*k*-dependence of the time-averaged radii (figure 3

_{y}*h*), which follows in the same approximation 4.3

Figure 4*a* shows that the agreement extends across the full range of separations , with particularly accurate estimates for the values of *k _{y}* marking the beginning and end of the transition (figure 4

*b*). These results suggest that the simple expression in equation (4.2) captures the essential features of the dynamics and can therefore be used to analyse the nature of the transition. For small deviations in

*k*around the transition from AP, equation (4.2) can be approximated as 4.4which therefore suggests that the emergence of non-trivial coordination follows a supercritical pitchfork bifurcation [48]. A similar argument leads to an equivalent conclusion for the bifurcation from IP as

_{y}*k*is decreased. AP and IP domains are therefore bounded by a pitchfork bifurcation of limit cycles.

_{y}Within the intermediate regime, the system becomes naturally bistable through a spontaneous symmetry breaking from a state where both oscillators follow the same limit cycle to one where they sustain a stable difference in their average oscillation amplitudes (see figure 3*e*). Transitions between homogeneous and inhomogeneous oscillation states, and bistability, have only recently been discovered in pairs of coupled limit cycle oscillators, also as a consequence of non-equilibrium symmetry breaking pitchfork bifurcations [49,50]. Here, we discover their emergence in a simple model of hydrodynamic- and basal body-coupled flagella. Intermediate equilibrium states appear when the internal elastic interaction promoting IP coordination is approximately compensated by the leading order hydrodynamic coupling favouring AP, amplifying the importance of higher-order hydrodynamic effects. In this parameter range, the oscillators are able to sustain a stable difference in their average amplitude (external rotor radii). Interestingly, the permanent difference in the average radii of the two oscillators after the bifurcation could be easily interpreted by an observer as a difference in intrinsic frequency (since ). The rotors would then appear intrinsically different despite in fact being identical. Eventually, a sufficient increase of the internal stiffness can overcome the antagonistic effect of hydrodynamic interactions at any given separation *l*, and drive the system to a stable IP state.

## 5. The full model

We conclude by looking at the full system with realistic parameters throughout (see table 1). In this case, radial and phase dynamics have comparable timescales () and the radii cannot be considered as approximately slaved to the phases anymore. Together with the sizeable radial deformations ( here and for *Chlamydomonas*-like ) this results in a complex interplay between radial and phase variables and implies the need to consider the full system of governing equations (see electronic supplementary material, S1). These will be explored here through numerical simulations only.

Figure 5*a*,*b* shows a representative set of curves for a *k _{y}* sweep with and . Electronic supplementary material movies 1–5 show the dynamics for (i)–(v) respectively. Similarly to the case of stiffer flagella, low and high values of

*k*correspond, respectively, to AP (i) and IP (v) synchronization, and in each of these states the oscillators follow the same dynamics as one another. The average phase speed, however, is observed to be lower in IP than in AP (a difference of between and ), in qualitative agreement with experimental observations of CR mutants [41] (see electronic supplementary material, S3). In between, there is a range of

_{y}*k*values for which the system synchronizes around intermediate values of , with the two oscillators following again different limit cycles (iii,iv). A new state, however, appears as

_{y}*k*approaches the symmetry breaking transition from the AP side (ii). Despite corresponding formally to AP (), the system displays symmetric excursions in the relative phase difference which are long lived and not much smaller than

_{y}*π*. In this state, the oscillators spend most of their time at values of

*σ*far from 0, and the null average is only guaranteed by the symmetry of the dynamics. Although the system oscillates here by about , amplitudes can be easily obtained just by increasing

*a*(electronic supplementary material, figure S2). In this condition, the system will not appear synchronized in AP at all, but will rather be continuously alternating between IP at

*π*and IP at . Figure 5

*c*,

*e*shows that this situation is typical for all of the separations displaying a discontinuous, rather than continuous, transition out of the AP state (here all ). In the IP case, instead, the transition maintains its continuous nature throughout, and in fact the bifurcation point is still well predicted by the first order expression from equation (4.1) (see inset).

From the AP side, discontinuous transitions are always preceded by a region of *k _{y}* values where the system displays large excursions in

*σ*, which act as a predictor of the impending discontinuity [51]. Figure 5

*e*shows the variance, , of the time-dependent signal, , averaged over 7 s ( beats). The large excursions evident in figure 5

*a*(ii) manifest in the peak of figure 5

*e*.

The presence of the discontinuity in and the preceding large fluctuations depend on both the separation *l* and *k _{x}*, as shown in figure 6 for a parameter sweep. For the realistic interflagellar separation (), figure 6

*a*shows that the region of discontinuous transition out of AP, marked by the large standard deviation of , exists only for . Above this value, the bifurcation changes its nature and becomes continuous but sharp. At the slightly larger separation of , the width of the extended transition zone observed previously for is reduced (see figures 6

*b*, 5

*c*). Further increasing the separation to reduces hydrodynamic forces by an order of magnitude compared to the case, and for all the

*k*values explored, the system follows a sharp continuous transition from AP to IP as

_{x}*k*is increased.

_{y}Although exploring in detail the nature of these bifurcations is beyond the scope of the present work, clear similarities with the case of stiff flagella suggest strongly that the qualitative nature of the continuous transition is the same in the two cases. We expect therefore the continuous transitions to be supercritical pitchfork bifurcations of limit cycles, inducing the observed symmetry breaking in the system (figure 5*a*(iii–iv)). For , the emergence of a discontinuity in implies that decreasing *k _{x}* can change the nature of the transition. Looking closely at the discontinuous case, we find that there is an extended region of overlap between the and the intermediate branches (see electronic supplementary material, S5). This is typical of a catastrophe-like transition which, given the symmetry of the system, is likely to be a subcritical pitchfork bifurcation.

Coexistence of three different states, all of which are locally stable for the dynamics, means that the system displays multi-stability: presence of noise might then induce the system to jump between these locally stable states and therefore alternate between periods of AP synchronization an other non-trivial types of coordination, with transitions dictated by escape rate arguments [52–54].

## 6. Conclusion

While mechanisms for hydrodynamic-led synchronization of cilia and flagella have been extensively studied [12,14,20,21,23–27,55], the impact of direct intracellular connections on flagellar dynamics is only starting to be recognized [15,39,40]. Here we have extended a simple and popular minimal model for the hydrodynamically interacting flagella pair of *Chlamydomonas* to account for intracellular mechanical coupling. The clearly anisotropic ultrastructure of striated fibres [42] is mirrored in the use of a non-isotropic elastic interaction between the oscillators () and results in a phase–phase coupling that can promote *by itself* either IP or AP synchronization, within a biologically plausible range of Young’s moduli. Transitions would then result simply through changes in the relative magnitude of *k _{x}* and

*k*. Given that intracellular calcium can control the contraction of striated fibres in

_{y}*Chlamydomonas*[56], we hypothesize that the transitions in coordination observed experimentally could be the result of localized apical variations in cytoplasmic [57–59]. Natural extensions of this model to amplitude–phase coupling do not influence the leading order coordination dynamics (see electronic supplementary material, S1, for brief discussion) and have been omitted here. IP coordination has recently been proposed to result from a different nonlinear interplay between hydrodynamic and intracellular mechanical coupling [15], with AP due to either one of them operating in isolation. However, several experimental observations, from the absence of phase-locking in mutants lacking striated fibres [39,40], to the complex synchronization observed in multi-flagellated algae [40], suggest that hydrodynamics plays in fact only a minimal role in this system. The model introduced here can sustain both AP and IP states without the need for external coupling, through a mechanism potentially under the direct control of the cell. Unequal tightening of the different fibres joining the basal bodies of cells with more than two flagella could then lead to the complex synchronization patterns observed experimentally [40]. Yet, subtle hydrodynamic effects can exist and need to be investigated through dedicated experiments blocking fluid-mediated coupling between flagella. These are currently underway.

With coordination of cilia and flagella within single cells being sensitive to the direct intracellular coupling between the filaments, we believe that better understanding its emergence will eventually enable synchronization to be used as a new and sensitive probe for the intracellular mechanical state of a cell.

## Data accessibility

Additional data are included in the electronic supplementary material.

## Authors' contributions

M.P. and D.R.B. designed the research; M.P., D.R.B., and R.C. developed the initial model; D.R.B. and Y.L. developed the full model and performed the simulations; M.P., D.R.B. and Y.L. analysed the data and wrote the manuscript.

## Competing interests

We declare we have no competing interests.

## Funding

This work was supported by: a Discovery Early Career Researcher Award (DECRA) from the Australian Research Council (D.R.B.); The University of Melbourne through an Albert Shimmins International Fellowship (M.P. and D.R.B.), and a Mathematics Vacation Scholarship (Y.L.).

## Acknowledgements

M.P. gratefully acknowledges the hospitality of the School of Mathematics and Statistics from The University of Melbourne, where part of this work was done.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.4238486.

- Received June 16, 2018.
- Accepted September 11, 2018.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.