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 that generate a variety of collective motions, which can be quantified with high-speed imaging in microfluidic environments [4-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-13], akin to the “Mexican wave” propagating around a packed stadium. Many surrogate models for flagellar dynamics have been proposed [13-24], typically with a set geometry that fixes the range and coupling between oscillators.

Here we relax this condition and study a linear array of colloidal oscillators [25] driven in circular trajectories at a controllable height above a no-slip wall. Originally introduced as a mathematically convenient minimal model for synchronization at low Reynolds numbers [15], colloidal rotors have been experimentally shown to reproduce the time-dependent flow field associated with a beating flagellum down to distances comparable to its size (~10 μm) [9,26]. When generalized to include waveform flexibility [14,27,28] they are also capable of capturing interflagellar synchronization in bulk [9]. The system of colloidal rotors studied here can be modified continuously from being primarily coupled 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 [18], these dynamics are reminiscent of chimera states, in which oscillators split into phase-locked and desynchronized clusters [29,30].

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 [31,32]. The OTs describe a planar array of circular trajectories [Fig. 1(a)] 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) ⩽ 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 another common situation where the ciliary beating plane is perpendicular to an organism’s surface.

Fig. 1

Experimental setup and results. (a) Microspheres of radius a = 1.74 μm, situated at a distance h above a no-slip boundary are driven by time-sharing optical tweezers in circular trajectories of radius R = 1.59 μm and center-to-center separation ℓ = 9.19 μm. (b) Average phase drift for a rotor pair vs detuning D for h = 4.2 (), 6.7 (), 11.7 (), 16.7 (), 31.7 (), 51.7 () μm. (c) Phase diagram showing experimental regions of synchrony (blue) and drift (red), the boundary from hydrodynamic simulations (dashed), and theory from Eq. (2) (solid).

The oscillators are imaged using 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 rotor positions are measured using an algorithm that correlates the image intensity I (x, y) with a rotationally symmetric kernel image K(x, y) constructed from a real colloid. By fitting the two-dimensional (2D) cross-correlation function C(x0,y0) = ∑( I (x,y) × K(x − x0, y − y0) with a 2D parabola and maximizing this function, the rotor positions (x0, y0) are extracted with subpixel resolution and used to track their phases {ϕ(t)} over time [Fig. 1(a)].

The rapid feedback loop between colloid and trap positions facilitates the arbitrary placement of the OTs with respect to the colloidal particles. The trap positions are maintained at a constant radius R and fixed angular distance ahead of the colloids. Consequently, 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 = FdrD leading to rotation. The period of rotation is therefore not fixed, and it is this degree of freedom that permits synchronization of interacting particles. The choice of λ reflects estimates of the bending rigidity of flagella, κ ~ 4 × 10−22 N m2, and their length, L. From λ = κ/L3 [14], and for typical values of L, values of should represent typical flagella. Unavoidable delays in the OT’s feedback response introduce a mismatch between the parameters used in experiments and simulations, which is corrected by increasing the simulation value of λ by a constant factor γ relative to the experimental one. The previously reported value of γ = 2.21 [32] is adopted throughout this paper, which results in quantitative agreement of simulations with the present experiments [see Fig. 1(c)].

Isolated oscillators rotate with a height-dependent angular velocity ω = F0ζ, where ζ0 = 6π μa is the sphere’s bulk drag coefficient, and accounts for the presence of the wall [33]. Experimental results are compared with deterministic hydrodynamic simulations in which colloids are treated as point-like particles above a no-slip boundary, and therefore coupled through the so-called “Blake tensor” [11,34]. Before each experiment we calibrate Fdr ≃ 2.23 pN (see Supplemental Material [35]; 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.

Consider first two rotors separated by a distance ℓ. For rotors with instantaneous positions {} and velocities {}, the hydrodynamic drag on the ith rotor is given by where is the net external force acting on the jth sphere and G(,) 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 [13,14]. The timescale for synchronization is proportional to the spring constant λ [see Eq. (1)] and also depends on the strength of hydrodynamic interactions between rotors, which is 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. 1(b) 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. 1(c), where the experimental synchronization boundary is based on a threshold of five slips in the whole experiment The results of individual experiments are classified based on this threshold, and are represented as either red (drift) or blue (synchronized) points in Fig. 1(c). 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. 1(c)] up to small shift in D, which could come from the finite value of a/h 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 [14,35]. As phase-locking is slow compared to the rotation period, we average over this fast timescale and find where From Eq. (1), the average phase drift for non-phase-locked states reads Given the functional form of the frequency detuning, F = FdrD, Eq. (2) can be solved explicitly to yield the critical detuning D*(h) [solid line in Fig. 1(c)]. The theoretical and numerical solutions for the boundary in Fig. 1(c) slightly under- and overestimate the data, respectively, owing to neglect of temporal variations in the interparticle 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 six rotors, with the ith rotor centered at = (il,0,h). This is the longest controllable chain with our active-feedback-based OTs. Linear arrays of colloidal oscillators have been shown to capture the dynamics of two-dimensional arrays [13], so this simplified geometry will be the focus here. The dynamics are studied experimentally as a function of h, but numerical simulations allow wider exploration of parameters, including changes in the radial stiffness λ, which governs the coupling strength [9,11,13,14,32] as in Eq. (1). In both experiments and simulations we introduce a mild frequency bias D = 1.01, typical also of Volvox colonies [13], 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 2(a) 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. 2(b) 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.

Fig. 2

Results for the linear array of driven colloidal oscillators, shown schematically in grey (not to scale). (a) Kymographs showing sin ϕ 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. The shaded red region represents the experimental parameter regime.

The phase dynamics of wobbles and defects are shown in Fig. 3(a) for h = 11.7(4) μm. The first 25 seconds of the time series show fluctuations in ϕ − ϕ0 (wobbles), even while the system is frequency locked. Fluctuations start at the first oscillator pair and travel unidirectionally along the chain [Fig. 3(b)] with a preserved, soliton-like signature [36]. Occasionally, they terminate within the chain with a slip [Fig. 3(a)]; these are the phase defects observed in kymographs. Both wobbles and defects are characterized by initial excursions of amplitude W and recurrence time τ [Fig. 3(c)], which depend on h [Figs. 3(d) and 3(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. 3(d)], mirroring the increased probability that a wobble will terminate in a slip within the chain, causing a defect [Fig. 3(f)]. 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 [35].

Fig. 3

Experimental phase dynamics. (a),(b) Phase difference relative to the first rotor, ϕ − ϕ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.

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 six 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” [18]. Figure 2(c) shows the magnitude of the coupling of a given oscillator with its nth nearest neighbor, estimated with Eq. (1), normalized by the total interaction strength with the first five 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− where j = 1 and 3 represent the far (β ≫ 1) and near (β ≪ 1) asymptotic limits [18]. For the end rotor the magnitude of the coupling with the nth nearest neighbor, normalized by the total coupling strength, is For β ≪ 1 the interactions are dominated by nearest neighbor, with S(1) = 0.84, while for β ≫ 1, S(1) = 0.44 [see the black curve in Fig. 2(c)]. 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) [35].

Fig. 4

(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. 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, and viscosity μ = 6 mPas. Simulations correspond to 0 ≤ t ≤ 2000 s. The dashed white line shows the value of λ corresponding to Fig. 2. Note the different color scales used throughout.

Additional numerical simulations permit the wider exploration of parameter space. The range of λ values studied here corresponds to the estimated values λ = κ/L3 based upon ciliary lengths of L ~ 4−10 μm [14]. Figure 4(a) shows the average end-to-end phase drift per beat as a function of λ and h, and enables analysis of many numerical simulations without looking at the individual kymographs. The area of solid blue corresponds to specific parameter combinations for which complete phase-locking occurs. However, from the drift alone, one cannot distinguish between a linear traveling wave ③ and a chevron phase profile ①. For this we compute the complex order parameter where χ = ϕ − ϕ [18,37]. Note the use of the pairwise phase differences, not the individual rotor phases. Looking at the mean value of [see Fig. 4(c)], the region of parameter space corresponding to complete phase-locking can be decomposed into chevron and MW regions. Using the average values Ā and for t > 200 s [Figs. 4(b) and 4(c)], we see that as h is increased at the experimental value of λ (white dotted line), the stable traveling wave at small h shifts to a profile with defects and wobbles, initially along the whole chain, and then localized to one-half of the chain with the remaining three oscillators constantly phase-locked.

At values of λ smaller than the experimental one, however, we observe qualitatively different dynamics. For λ ≲ 2.5 pN/μm, the pattern morphs continuously between different types of complete synchronization as h is increased, going from a MW ③ to a chevron-like pattern ①. These transitions happen without the emergence of defects [11]. For 2.5 < λ < 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. 4(c)] 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. 4(d)] highlights the role played by competition between interactions at different ranges. Changing h is then a simple and accessible way to modulate 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 [13], so that effectively h/ℓ ~ 1 [see Fig. 4(a)]. 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 [22]) may operate in a regime naturally prone to the emergence of metachronal phase defects, which are indeed observed experimentally [13].