Long-range interactions, wobbles, and phase defects in chains of model cilia

Eukaryotic cilia and flagella are chemo-mechanical oscillators capable of generating long-range coordinated motions known as metachronal waves. Pair synchronization is a fundamental requirement for these collective dynamics, but it is generally not sufficient for collective phase-locking, chiefly due to the effect of long-range interactions. Here we explore experimentally and numerically a minimal model for a ciliated surface: hydrodynamically coupled oscillators rotating above a no-slip plane. Increasing their distance from the wall profoundly affects the global dynamics, due to variations in hydrodynamic interaction range. The array undergoes a transition from a traveling wave to either a steady chevron pattern or one punctuated by periodic phase defects. Within the transition between these regimes the system displays behavior reminiscent of chimera states.

The ability of ensembles of oscillators to achieve collective motions is fundamental in biological processes ranging from the initiation of heartbeats to the motility of microorganisms.The emergent properties of coupled oscillators can vary dramatically depending on the intrinsic properties of the oscillators and the nature of the coupling between them [1].Flashing fireflies equally and instantaneously coupled to one another [2] can support very different behaviors to chemical micro-oscillators which are coupled only locally, and subject to time delays [3].
Eukaryotic cilia and flagella are chemo-mechanical oscillators which generate a variety of collective motions that can be quantified with high-speed microscopy in microfluidic environments [4][5][6].The molecular biology of these internally driven filaments is virtually identical in green algae [5], protists [7] and humans [8], and the flows they generate fulfill crucial roles in development, motility, sensing and transport.When close together, the mutual interaction between their oscillatory flow fields can cause them to beat in synchrony [9], and larger ensembles of flagella demonstrate striking collective motions in the form of metachronal waves (MWs) [10][11][12], akin to the 'Mexican wave' propagating around a packed stadium.Many surrogate models for flagellar dynamics have been proposed [12][13][14][15][16][17][18][19][20][21][22], typically with a set geometry which fixes the range and coupling between oscillators.
Here we relax this condition and study a linear array of colloidal oscillators [23] driven in circular trajectories at a controllable height above a no-slip wall.The oscillator couplings can be modified continuously from being primarily through nearest neighbors to a regime involving significant long-range interactions.As a function of rotor properties, a traveling wave found at small heights becomes either a chevron pattern or is punctuated by phase defects at large ones.The transition is not a gradual morphing between the two profiles, but rather a process involving generation and propagation of defects along the strip, where phase-locked and non-phase-locked subgroups of oscillators can coexist.A behavior arising from long-range interactions whose amplitude is modulated by the distance from the wall [16], these dynamics are reminiscent of chimera states, in which oscillators split into phase-locked and desynchronized clusters [24,25].
In our experiments, silica colloids of radius a = 1.74 µm (BangsLab, USA) suspended in a water-glycerol solution of viscosity µ = 6 mPas within a 150 µm-thick sample, are captured and driven by feedback-controlled time-shared (20 kHz) optical tweezers (OTs) based on acousto-optical deflection of a 1064 nm-wavelength diode-pumped solid-state laser (CrystaLaser IRCL-2W-1064) as previously described [26,27].The OTs describe a planar array of circular trajectories (Fig. 1a) of radius R = 1.59 µm and center-to-center separation = 9.19 µm, a distance h above the sample bottom, with 4.2(4) µm ≤ h ≤ 51.7(4) µm.This configuration, which reflects the capabilities and limitations of our OT setup, is similar to arrays of nodal cilia, but differs from the more common situation where the ciliary beating plane is perpendicular to an organism's surface.The OTs are arranged so that a colloidal particle on the ith trajectory (i ∈ {0, . . ., N − 1}) experiences a radial harmonic potential with spring constant λ = 2.06 ± 0.06 pN/µm resisting excursions from the prescribed radius, and a constant tangential force of magnitude F i = F dr D i−0.5 leading to rotation.Before each experiment we calibrate F dr 2.23 pN (see [28]; typical variation ±2%).D = 1 is used to break left-right symmetry along the chain and induce a stable traveling wave for small h [11].For the detuning adopted here, D = 1.01, the period of individual oscillators varies between τ ∼ 0.5 s and ∼ 1 s across the explored range of h.The oscillators are imaged under a Nikon inverted Eclipse Ti-E with a 60× Nikon Plan Apo VC water immersion objective (NA = 1.20), and recorded for up to 1200 s using an AVT Marlin F131B CMOS camera set at 229 fps.The particle's positions are measured using a profile image correlation with subpixel resolution, and used to track the rotor phases {φ i (t)} (Fig. 1a).Isolated oscillators rotate with a height-dependent angular velocity ω i = F i /Rζ 0 ζ w , where ζ 0 = 6πµa is the sphere's bulk drag coefficient, and ζ w (h) = 1 + 9 16 a h + O(a 3 /h 3 ) accounts for the presence of the wall [29].Experimental results are compared with deterministic hydrodynamic simulations based on the Oseen approximation to inter-particle coupling [11].Unavoidable delays in the OT's feedback response introduce a mismatch between experiments and simulations which for a pair of oscillators is corrected by increasing the simulation value of λ by a factor κ relative to the experimental one.We estimate κ = 2.21 in agreement with previous reports [27].
Consider first two rotors separated by a distance .For rotors with instantaneous positions {x i } and velocities {v i }, the hydrodynamic drag on the ith rotor is given by − where F ext j is the net external force acting on the jth sphere and G(x i , x j ) is the Green's function in the presence of the no-slip wall.For identical rotors (detuning D = 1), hydrodynamic coupling eventually leads to synchrony provided λ < ∞, by perturbing the angular velocities of the two rotors so that the leading and lagging rotors become slower and faster respectively [12,13].The timescale for synchronization depends on the spring constant λ and the strength of hydrodynamic interactions between rotors, which is itself a function of height h and spacing .The dynamics become richer if a discrepancy between the rotor's intrinsic frequencies is introduced (D = 1), for then the coupling must be sufficiently strong to overcome the natural tendency for the rotor's phase difference χ = φ 1 − φ 0 to drift.
Bifurcation plots in Fig. 1b show, for different h, the average phase drift between two oscillators as a function of D. The behavior is typical of a saddle-node bifurcation: the oscillators phase-lock until D reaches a critical value D * (h) and then drift with a monotonically increasing speed.D * (h) increases with h, reflecting the strengthening of inter-rotor hydrodynamic coupling with increasing distance from the wall.The phase-locking behavior is summarized in Fig. 1c, where the experimental synchronization boundary is based on a threshold of 5 slips in the whole experiment ( χav = 0.131 rad/s).As h is increased, the rotor pair moves deeper into the synchronized region: the coupling between the two strengthens due to lower hydrodynamic screening from the wall, leading to an enhanced stability of the synchronized state.This is reproduced by simulations (Fig. 1c) up to small shift in D, likely the result of a finite value of a/d and experimental noise.In the limit a, R , the evolution of the phase difference χ = φ 1 − φ 0 can be derived by a generalization of previous arguments [13,28].As phase-locking is slow compared to the rotation period, we average over this fast timescale and find where, and β = 2h/ .From Eq. ( 1), the average phase drift χav for non-phase-locked states reads (2) Given the functional form of the frequency detuning, F i = F dr D i−0.5 , Eq. ( 2) can be solved explicitly to yield the critical detuning D * (h) (solid line in Fig. 1c).The theoretical and numerical solutions for the boundary in Fig. 1c slightly under-and over-estimate the data, respectively, owing to neglect of temporal variations in the inter-particle spacing and the finite size of the beads, respectively.Both also neglect thermal fluctuations.
We now turn to the dynamics of a linear array of 6 rotors, with the ith rotor centered at x = (il, 0, h).This is the longest controllable chain with our active-feedbackbased OTs.The dynamics are studied experimentally as a function of h, but numerical simulations allow wider ex- ploration of parameters, including changes in the radial stiffness λ, which governs the coupling strength [9,[11][12][13]27] as in Eq. ( 1).In both experiments and simulations we introduce a mild frequency bias D = 1.01, typical also of Volvox colonies [12], which breaks the translational symmetry and induces a MW for h 10 µm.At all heights studied, this value of D is deep within the synchronized region of parameter space for two rotors.Figure 2a shows that at h = 6.7(4) µm the rotors phase-lock in a stable MW whose direction is set by the frequency bias.With increasing h, defects (phase slips) emerge, giving rise to a net drift in the cumulative phase difference between rotors at opposite ends of the chain.Phase defects always propagate in the direction of the fastest oscillator.At these intermediate heights, the phase profile also displays "wobbles" -perturbations to the MW that are not accompanied by a phase defect.Numerical results shown in Fig. 2b capture the traveling wave at h = 5 µm, the presence of defects and their propagation direction, and wobbles at larger heights.At the largest height, h = 50 µm, defects no longer propagate through the chain, and rotors 3-5 remain phase-locked.
The phase dynamics of wobbles and defects are shown in Fig. 3a for h = 11.7(4)µm.The first 25 seconds of the time-series show fluctuations in φ i − φ 0 (wobbles), even while the system is frequency-locked.Fluctuations start at the first oscillator pair and travel unidirectionally along the chain (Fig. 3b) with a preserved, soliton-like signature [30].Occasionally, they terminate within the chain with a slip (Fig. 3a); these are the phase defects observed in kymographs.Both wobbles and defects are characterized by initial excursions of amplitude W and recurrence time τ (Fig. 3c) which depend on h (Fig. 3d,e).The typical time τ ∼ 10 T (where T 1 s is the average period) depends less strongly on h than does W , which shows a pronounced growth (Fig. 3d) mirroring the increased probability that a wobble will terminate in a slip within the chain, causing a defect (Fig. 3f).Although their position can vary, defects tend to cluster, in this case at the middle of the chain (position i = 2), as seen also in simulations of longer chains [28].
The hydrodynamic coupling between two rotors increases monotonically with h; for an isolated pair, this manifests in more robust synchronization at larger heights.For a chain of rotors, increasing h has the reverse effect, disrupting the stable MW with wobbles and punctuating it with periodic phase defects (Fig. 2).The hydrodynamic coupling between every pair of rotors in the chain grows as h is increased.For just two rotors, Eq. (1) shows that equivalent changes to the hydrodynamic coupling can be achieved through modification of the mean interparticle separation l.For the chain of 6 rotors, in which longer range hydrodynamic interactions also occur, changes to h and l are no longer equivalent.The peculiar dynamics observed arise from a change in the relative contributions of interactions with different neighbors.The no-slip wall has the effect of screening the hydrodynamic interactions in a way that qualitatively changes as a function of β = 2h/ .This is an important determinant of MW stability, as observed also in simulations of colloidal "rowers" [16].Figure 2c shows the magnitude of the coupling of a given oscillator with its jth nearest neighbor, estimated with Eq. ( 1), normalized by the total interaction strength with the first 5 neighbors.Although all pairwise couplings grow monotonically with h, the relative magnitude of the nearest neighbor interactions actually diminishes.Conversely, the relative importance of all others increases with h.Hydrodynamic disturbances parallel to the wall decay as u ∼ r −j where j = 1 and 3 represent the far (β 1) and near (β 1) asymptotic limits [16].For the end rotor the magnitude of the coupling with the nth nearest neighbor, normalized by the total coupling strength is S(n) = n −j / 5 i=1 i −j .For β 1 the interactions are dominated by nearest neighbor, with S(1) = 0.84, while for β 1, S(1) = 0.44 (see black curve in Fig. 2c).We test the hypothesis that the breakdown of the traveling wave is due to long-range hydrodynamic interactions through simulations in which interactions are truncated at nearest neighbors, and find the abundance of defects is significantly reduced.Importantly, the dynamics are nearly insensitive to h, with a maximum relative variation in end-to-end drift speed of just 3% between h = 5 µm and 1000 µm (Fig. 4) [28].
Additional numerical simulations permit the wider ex-ploration of parameter space.Figure 4a shows the average end-to-end phase drift per beat as a function of λ and h.We also compute the complex order parameter [16,31].Using the average values Ā and Ψ for t > 200 s [Fig.4b,c], and following the experimental path (dotted line), we see the stable traveling wave at small h 3 developing and, as h increases, defects and wobbles along the whole chain first 4 and then localized to the initial half of the chain 2 , with the remaining three oscillators constantly phase-locked.At values of λ smaller than the experimental one, however, we observe richer dynamics.For λ 2.5 pN/µm, the pattern morphs continuously between different types of complete synchronization, going from a MW 3 to a chevron-like pattern 1 .These transitions happen without the emergence of defects [11].For 2.5 pN/µm < λ < 3 pN/µm the system shows reentrant behavior with defects only at intermediate heights, separating a MW region from a chevron-like region.The order parameter angle | Ψ| (Fig. 4c) identifies clearly the stable MW (yellow/orange) and chevron (dark blue) regions of parameter space.For a fixed h 50 µm, increasing λ results in a monotonic decrease in Ā owing to the reduced rotor compliance.Conversely, the end-to-end phase drift exhibits a strong peak around λ = 4.5 pN/µm, where the rotors slip approximately one beat in every five, despite an intrinsic frequency difference of just 5%.These nontrivial dynamics emerge due to the combination of phase slips induced by long-range interactions, and rapid healing of phase defects through orbit compliance.The complete absence of these features from the simulations with nearest neighbor coupling alone (Fig. 4d) highlights the role played by competition between interactions at different ranges.Changing h is then a simple and accessible way to mod-ulate their relative strength (see Fig. 2).
Large arrays of cilia are synonymous with no-slip boundaries, and in many cases, the spacing between these organelles is comparable to their length [12], so that effectively h/ ∼ 1 (see Fig. 4a).Our results suggest that flagella of Volvox may then be balancing the need to extend out into the fluid enough to generate a vigorous thrust, with the screening of long-rage hydrodynamic interactions necessary to stabilize MWs on the colony surface.As a result, ensembles of flagella in Volvox [11] (but see also numerical simulations [20]) may operate in a regime naturally prone to the emergence of metachronal phase defects, which are indeed observed experimentally [12].
We are grateful to T.J. Pedley and D. Bartolo for useful discussions.This work was supported by a Human Frontier Science Program Cross-Disciplinary Fellowship (DRB), a Wellcome Trust Senior Investigator Award (REG), and the EU ERC CoG Hydrosync (PC).

SUPPLEMENTARY MATERIAL Single rotor force calibration
For a single bead of radius a in a viscous fluid, situated at a distance h from an infinite no-slip boundary, the external applied force F is related to its velocity v according to F = ζ • v, where ζ is the anisotropic drag matrix, given by [14] The coefficient ζ 0 = 6πµa is the drag on the sphere in an unbounded fluid of viscosity µ (equivalent to setting h → ∞).We are interested only in trajectories which are parallel to the no-slip wall (v • e z = 0).For a constant applied driving force F dr , the sphere's speed v = |v| is given by implying a monotonic increase of the sphere's speed with h for a given F dr .Each set of experiments involves studying the colloidal oscillators at a number of different heights h.For each set, the center of the trajectory, its radius and the driving and radial forces are calibrated, for each individually loaded rotor, at the height of h = 22 µm.These are then checked for independence on h. Figure S1 shows, after a full calibration, the speed of an individually loaded colloidal oscillator at 6 different heights together with the prediction from Eq. (S2) using a constant driving force.The two agree well for F dr = 2.23 pN.

Hydrodynamic interaction of two rotors at an arbitrary distance from a no-slip plane
The fluid disturbance produced by the motion of a sphere parallel to a no-slip wall depends on its height above the planar boundary.For two such spheres situated in the fluid, it is important to calculate the strength of the hydrodynamic interactions between them, and the subsequent effects on their dynamics.We consider two spheres of radius a driven around circular orbits of radius R 0 which are parallel to a no-slip wall.The orbit's centers are located at positions (x, y, z) = (0, 0, h) and ( , 0, h) respectively.The plane z = 0 represents the noslip boundary, with the semi-infinite domain z > 0 filled with fluid of viscosity µ. (ê φ1 , êR1 ) and (ê φ2 , êR2 ) are the unit vectors of the local cylindrical frame of reference of each single rotor.This reference frame is centered at the center of the orbit.The displacement of each sphere from the center of its trajectory will be expressed in its cylindrical frame of reference as (R 1 , φ 1 ) and (R 2 , φ 2 ) respectively.Each rotor is subject to a constant tangential driving force F i = F i êRi , and also to a radial spring force with stiffness λ, which suppresses excursions from the equilibrium radius R 0 .The spring stiffness is assumed to be large enough that the radial degree of freedom is slaved to the angular degree of freedom.That is, knowing (φ 1 , φ 2 ), we know the instantaneous value of (R 1 , R 2 ).It will be our goal to derive the equations of motion of the spheres, without making any assumptions about the relative magnitudes of and h.
We assume that sphere radius and trajectory radius are both small compared to other length scales (a, R 0 h, ).Correspondingly, the two orbits are sufficiently far from each other that we can neglect the variation in the relative separation between the spheres as they move.The separation vector will always be taken to be êx .We will write the relation between the sphere's angular velocity ω i and the tangential driving force F i , as F i = ζ 0 ζ w R 0 ω i , where ζ 0 = 6πµa is the bulk drag coefficient, and ζ w is the correction due to the presence of the wall.For sufficiently small a/h this correction can be written as ζ w = 1+(9a/16h)+O((a/h) 3 ).The first thing required is the generic expression of the 'Blakelet', i.e. the Stokeslet on top of a bounding wall.This is given by [32]: Here r = ( , 0, 0), r = |r|; R = ( , 0, 2h), R = |R|; α ∈ {1, 2}.In our case, the point force will be in the (x, y) plane.Furthermore, we are only interested in the component of the velocity in the same plane, since the trajectories are constrained to lie at z = h.The background flow is due to the motion of the other sphere, which is: In the limit as h → 0, this reduces to Eq. ( 16) in [13].
Being interested only in the || component, and because Calling, as in [13], s = F/8πµ we get where r = r/r = êx , called n 21 in [13]; β = 2h/ , and . (S8) Notice that, due to the nearby wall, the strength of the Stokeslet s is written in terms of the sphere's velocity as: The derivation of the equations of motion follows the same procedure outlined in Appendix A of [13].The only things we need to calculate are: The rest of the calculation can be carried out in exactly the same way as in [13] and the final result is where now ρ = 3aζ w /8 , α = ωζ 0 ζ w /λ, and For example, this means that the phase difference χ = φ 2 − φ 1 and phase sum Φ = φ 1 + φ 2 evolve according to χ To the first order in the small quantities ∆ω = ω 2 − ω 1 and ρ, and averaging over a "natural" timescale of the fast variable Φ we get These functions can be rewritten as where D(χ) = D 0 sin χ and S(χ) = S 0 cos χ.As in [13], time has been rescaled according to the mean angular speed ω, and both ω i are measured in units of ω.Rewriting Eq. (S18) in dimensional units yields χ = ∆ω − 3a 4 (S24) For each rotor (now indexed by i ∈ {0, 1}), the intrinsic angular frequency is given by ω i = F i /(ζ 0 ζ w R 0 ) and so the above equations reads: (S25) Since a detuning factor D is included so that the driving force is F i = F dr D i−1/2 , the above equation can be written as (S26) The threshold value of D beyond which the rotors' phase difference will drift can be calculated explicitly.

Varying chain length
In order to assess the generality of the results presented in the main text, we used numerical simulations to explore the effect of changing the number of rotors, N , present in the linear array (see Fig. 1a).Figure S2 shows the average phase drift (measured in beats per beat) with respect to the first rotor, along chains of different length, N ∈ {2, 15}.Each chain has a fixed detuning of 5% between the end rotors.For each height h = 10 µm and h = 100 µm, simulations were conducted with full hydrodynamic coupling and nearest neighbor coupling only.
The results for N = 6 are representative of the dynamics across a range of chain lengths.For chains in which rotors are coupled through nearest neighbor interactions, the rotors tend to phase-lock in clusters of 2-5 rotors.As discussed in the main text, the nearest neighbor results are fairly insensitive to changes in h, shown here by the similarly between the results of Fig. S2b and d.In stark contrast, the chains in which rotors are fully coupled to one another through hydrodynamic interactions exhibit qualitatively different behavior at different heights.

Truncation of hydrodynamic interactions
Figure S3 shows the results of deterministic numerical simulations, with hydrodynamic interactions truncated to be nearest neighbor in nature (see also Fig. 4d).The dynamics are almost completely insensitive to changes in h, across several orders of magnitude.

FIG. 2 .
FIG. 2. (color online).Results for the linear array of driven colloidal oscillators.(a) Kymographs showing sin φi at three heights above the wall.With increasing h, the traveling wave becomes frustrated, with the introduction of wobbles (arrows) and phase defects (circles).(b) Numerical results from model.(c) Fraction of total coupling corresponding to interacting with different neighbors, as a function of h.Shaded red region represents experimental parameter regime.

FIG. 3 .
FIG. 3. (color online).Experimental phase dynamics.(a,b) Phase difference relative to the first rotor, φi − φ0, at h = 11.7 µm.(c) Wobbles are characterized by their magnitude W (radians) and timescale τ / T (normalized by rotor period), shown as a function of h in panels (d) and (e).(f) Probability that a propagating wobble ends at rotor i, resulting in a slip.

FIG. 4 .
FIG. 4. (color online).(a) Average phase drift per beat between end oscillators (measured in beats) as a function of height above the wall and radial spring stiffness.Shown also are four representative kymographs.(b) Time-averaged amplitude Ā and (c) angle | Ψ| of the complex order parameter Z = Ae iΨ .The axes are the same as in a.(d) Same as a but with hydrodynamic interactions truncated to nearest neighbor.Parameters used include a = 1.74 µm, = 9.19 µm, R = 1.59 µm, viscosity µ = 6 mPas.Simulations correspond to 0 ≤ t ≤ 2000 s.The dashed white line shows the value of λ corresponding to Fig. 2.

FIG. S1 .
FIG. S1.Calibration of an individual colloidal oscillator, moving under the influence of a harmonic potential in an optical tweezer.Experimental results (dots) are shown alongside the prediction of Eq. (S2), with which the driving force of F dr = 2.23 pN can be extracted.
FIG. S2. (color online).Average phase drift with respect to the first rotor, for chains of different lengths N ∈ {2, 15}, at two different heights h ∈ {10 µm, 100 µm}, and subject to either full hydrodynamic coupling or nearest neighbor interactions only.The end-to-end detuning is fixed at 5% in each case, the radial spring stiffness is λ = 4.5 pN/µm, and all other parameters are as in Fig.4.
FIG. S3. (color online).Kymographs showing the phase sin φi along the linear chain of model rotors, coupled hydrodynamically through the Blake tensor, but with interactions artificially restricted to be nearest neighbor.The radial spring stiffness is λ = 4.5 pN/µm and all other parameters are as in Fig. 4.