Hydrodynamic synchronization of spontaneously beating filaments

TL;DR

This study models how hydrodynamic interactions and biochemical feedback lead to synchronization of beating filaments like sperm and cilia, explaining experimental observations and the roles of waveform symmetry.

Contribution

It introduces a geometric feedback model that captures elastohydrodynamic phase synchronization across various filament waveforms, including asymmetric and symmetric patterns.

Findings

Both in-phase and anti-phase synchronization can occur in asymmetric beats.

Symmetric waveforms tend to synchronize in-phase.

Biochemical noise can cause phase slips in synchronization.

Abstract

Using a geometric feedback model of the flagellar axoneme accounting for dynein motor kinetics, we study elastohydrodynamic phase synchronization in a pair of spontaneously beating filaments with waveforms ranging from sperm to cilia and Chlamydomonas. Our computations reveal that both in-phase and anti-phase synchrony can emerge for asymmetric beats while symmetric waveforms go in-phase, and elucidate the mechanism for phase slips due to biochemical noise. Model predictions agree with recent experiments and illuminate the crucial roles of hydrodynamics and mechanochemical feedback in synchronization.

Equations

$$f_m(s,t)=\rho (n_{+}{F}_{+}+n_{-}{F}_{-})-K\Delta ,$$

$$\displaystyle\begin{split}\sigma_{ss}-\Big{(}1+\frac{c_{\parallel}}{c_{\perp}}\Big{)}N_{s}\phi_{s}-N\phi_{ss}-\frac{c_{\parallel}}{c_{\perp}}\sigma\phi_{s}^{2}&\\ =c_{\parallel}\big{(}\phi_{s}u^{d}_{n}-&\,\partial_{s}u^{d}_{t}\big{)},\end{split}$$

$$\displaystyle\begin{split}N_{ss}-\frac{c_{\perp}}{c_{\parallel}}N\phi_{s}^{2}+\sigma\phi_{ss}+\Big{(}1+\frac{c_{\perp}}{c_{\parallel}}\Big{)}\sigma_{s}\phi_{s}&\\ =c_{\perp}\big{(}\mathrm{Sp}^{4}\phi_{t}-u^{d}_{t}\phi_{s}-&\,\partial_{s}u^{d}_{n}\big{)},\end{split}$$

$${\varphi }_{ss}+{\mu }_a{f}_m+N=0,$$

$${\mathbf{u}}^d(s_{\alpha })=\mathbf{K}[{\mathbf{f}}_e^{\alpha }](s_{\alpha })+{\int }_0^1\mathbf{G}(s_{\alpha };s_{\beta })\cdot {\mathbf{f}}_e^{\beta }(s_{\beta })d{s}_{\beta },\text{vspace}-0.cm$$

$$\begin{array}{rl}{\partial }_t{n}_{\pm }=& \eta (1-n_{\pm })\\ & -(1-\eta )n_{\pm }\exp \left[f^{\ast }(1\mp \zeta {\Delta }_t)\right]+\xi (s,t),\end{array}$$

$$\dot{\delta }=ϵ\sin \delta +\alpha \sin 2\delta ,\text{vspace}-0.1cm$$

Full text

Studies on flagellar synchronization date back to observations by Rothschild Rothschild (1949) on nearby swimming sperms and subsequent theoretical work by Taylor Taylor (1951), who proved that dissipation for two swimming sheets is minimized for an in-phase configuration. While biology is often not driven by dissipation principles, it has long been hypothesized that hydrodynamic interactions play a central role in synchronization Elfring and Lauga (2009) and in collective behaviors such as metachronal waves in ciliary arrays Golestanian et al. (2011). Over the last two decades, experiments Rüffer and Nultsch (1985); Polin et al. (2009); Goldstein et al. (2009); Leptos et al. (2013); Wan et al. (2014); Brumley et al. (2014) using micropipette-held Chlamydomonas have revealed that elastohydrodynamic interactions may indeed be at play in causing its two flagella to synchronize their breaststrokes, with periods of asynchrony thought to arise due to biochemical noise Brumley et al. (2014); Pikovsky et al. (2003).

Theoretical progress in understanding synchronization is complicated by the intricate internal structure and actuation of the flagellum core, or axoneme. In presence of ATP, thousands of dynein molecular motors act in concert to bend the structure and drive spontaneous beats. Much work has gone into developing minimal models that neglect this biological complexity and coarse-grain flagella as microspheres driven on compliant or tilted orbits Vilfan and Jülicher (2006); Niedermayer et al. (2008); Uchida and Golestanian (2011); Maestro et al. (2018). More detailed numerical models have relied on pre-imposed internal or external actuations to analyze metachronal waves Gueron et al. (1997); Elgeti and Gompper (2013) or the bistability of elastic filaments Guo et al. (2018), yet these descriptions poorly capture experimental waveforms Not . Only recently have there been attempts to study the role of hydrodynamics in simplified models of active microfilaments Goldstein et al. (2016).

While the detailed process leading to spontaneous flagellar oscillations remains controversial, several mechanisms have been proposed ranging from flutter-like instabilities Bayly and Dutcher (2016); De Canio et al. (2017); Ling et al. (2018) to dynamic internal tension Han and Peskin (2018) and geometric control of dynein kinetics Lindemann (1994); Brokaw (1971); Bayly and Wilson (2014, 2015); Sartori et al. (2016); Oriola et al. (2017); Sartori (2019).  The role of dynein motors in driving oscillations of microtubules has also been established in purified in vitro systems Vilfan et al. (2019). Building on previous sliding and curvature control models Riedel-Kruse et al. (2007); Oriola et al. (2017); Sartori et al. (2016); Sartori (2019), we recently developed Chakrabarti and Saintillan (2019) a microscopic description for an active elastic flagellum accounting for internal dynein motor kinetics, which produces spontaneous oscillations following saturation of a Hopf bifurcation and can generate a variety of beating patterns observed in nature. In this Letter, this microscopic model is employed to analyze the temporal dynamics and synchronization of a pair of spontaneously beating filaments. Our results explain synchronization in various situations with trends consistent with experiments Brumley et al. (2014), and underscore the crucial roles of hydrodynamic interactions, mechanochemical feedback, and biochemical noise.

*Model formulation.—*The flagellar axoneme is a 3D structure with circular cross-section composed of 9+2 pairs of microtubules arranged in a cyclic fashion Alberts et al. (2013). Following past models Machin (1958); Brokaw (1965, 1971, 1979); Oriola et al. (2017); Chakrabarti and Saintillan (2019); Camalet and Jülicher (2000), we idealize this structure as a planar projection with diameter $a$ and length $L$, where microtubules are represented by two polar filaments ${\bf x}_{\pm}$ clamped at the base [Fig. 1(a)]. We seek an evolution equation for the centerline ${\bf x}(s,t)$ parametrized by arclength $s\in[0,L]$. Internal actuation arises from dynein motors that extend from each filament and stochastically bind with the opposite one. In presence of ATP, these motors move along the filaments and generate internal shear forces resulting in an arclength mismatch known as the sliding displacement Camalet and Jülicher (2000): $\smash{\Delta(s,t)=\int_{0}^{s}\left(|\partial_{s^{\prime}}{\bf x}_{-}|-|\partial_{s^{\prime}}{\bf x}_{+}|\right)ds^{\prime}=a\phi(s,t)}$, where ${\phi}(s,t)$ is the tangent angle and $(\cdot)_{s}$ denotes arclength derivative. This sliding is resisted by internal protein linkers, or nexin links, modeled as linear springs of stiffness $K$. Both dynein motors and nexin links result in equal and opposite force densities along the filaments:

[TABLE]

where $\rho$ is the mean motor density, $n_{\pm}$ are the fractions of motors in the bound state, and $F_{\pm}$ are the associated loads. These sliding forces generate internal moments $\smash{M(s,t)=B\phi_{s}-a\int_{s}^{L}f_{m}(s^{\prime},t)ds^{\prime}}$, where the first term captures the passive elastic response of the structure modeled as an inextensible Euler-Bernoulli elastica with bending rigidity $B$ du Roure et al. (2019). Dynein motor populations evolve as $\partial_{t}n_{\pm}=\pi_{\pm}-\epsilon_{\pm}$ where $\pi_{\pm}$ and $\epsilon_{\pm}$ are the attachment and detachment rates, respectively. The attachment rate is proportional to the fraction of unbound motors: $\pi_{\pm}=\pi_{0}\left(1-n_{\pm}\right)$. The detachment rate depends linearly on the fraction of bound motors and exponentially on the carried load Svoboda and Block (1994); Müller et al. (2008): $\epsilon_{\pm}=\epsilon_{0}n_{\pm}\exp(\pm F_{\pm}/f_{c})$, where $f_{c}$ is a critical load above which rapid unbinding occurs. To complete the model with the appropriate geometric feedback, we specify a force–velocity relation for the dyneins Oriola et al. (2017); Shelley (2016): we assume that the motors have a velocity $v_{0}$ at zero load that decreases linearly with sliding velocity $\Delta_{t}\equiv a\phi_{t}$ and are able to carry a load $f_{0}$ when stalled, yielding the expression: $F_{\pm}=\pm f_{0}(1\mp\Delta_{t}/v_{0})$.

Motion of the centerline is governed by the force balance for an elastic rod in viscous flow Antman (1995): $\partial_{s}\mathbf{F}_{e}+{\bf f}_{v}=0$, where $\mathbf{F}_{e}(s,t)=\sigma\boldsymbol{\mathrm{\hat{t}}}+N\boldsymbol{\mathrm{\hat{n}}}$ is the elastic force with tension $\sigma$ and normal force $N$, and ${\bf f}_{v}(s,t)$ is the viscous force density captured by nonlocal slender body theory Keller and Rubinow (1976); Tornberg and Shelley (2004). This is accompanied by a moment balance in the $z$ direction: $M_{s}+N=0$. Scaling lengths by $L$, sliding displacement by $a$, time by the correlation timescale $\tau_{0}=1/(\pi_{0}+\epsilon_{0})$, elastic forces by $B/L^{2}$, and motor loads by $\rho f_{0}$ produces four dimensionless groups, of which two are of primary interest: (i) the sperm number $\smash{\mathrm{Sp}=L\left(8\pi\nu/B\tau_{0}\right)^{1/4}}$, where $\nu$ is the viscosity, compares the relaxation time of a bending mode to the motor correlation time; (ii) the activity number $\smash{\mu_{a}=a\rho f_{0}L^{2}/B}$ compares motor-induced sliding forces to characteristic elastic forces. The two other dimensionless groups are: $\smash{\mu=Ka^{2}L^{2}/B}$ and $\zeta=a/(v_{0}\tau_{0})$ Chakrabarti and Saintillan (2019). With these scalings, the dimensionless equations for $\phi(s,t)$ read

[TABLE]

[TABLE]

where the first two equations are force balances in the tangential and normal directions while the third is the moment balance. The dimensionless tangential and normal drag coefficients $c_{\parallel,\perp}$ derive from local terms in slender body theory Chakrabarti and Saintillan (2019) and satisfy $c_{\perp}/c_{\parallel}\to 2$ for infinitely slender filaments. Hydrodynamic interactions are captured by the disturbance velocity $\mathbf{u}^{d}$, with projections $u_{t}^{d}$ and $u_{n}^{d}$ in the tangential and normal directions, respectively. Given two filaments indexed by $\{\alpha,\beta\}$, the flow is obtained as

[TABLE]

where ${\bf f}_{e}\equiv\partial_{s}\mathbf{F}_{e}$ is the elastic force density. The first term in Eq. (5) is the finite-part integral of slender body theory Keller and Rubinow (1976); Tornberg and Shelley (2004); Chakrabarti and Saintillan (2019) and captures hydrodynamic interactions within a filament. The second term accounts for the flow induced by the other filament, with the Green’s function $\mathbf{G}(s_{\alpha};s_{\beta})$ given by the Oseen tensor. These equations are supplemented by clamped boundary conditions at $s=0$, and moment- and force-free conditions at $s=1$ Chakrabarti and Saintillan (2019). In dimensionless form, the evolution equation for the bound motor populations reads

[TABLE]

where $\eta=\pi_{0}/(\epsilon_{0}+\pi_{0})$ is the fraction of time spent by motors in the bound state and ${f}^{*}=f_{0}/f_{c}$ is the ratio of the stall load to the characteristic unbinding force. The last term accounts for biochemical noise with $\langle\xi(s,t)\rangle=0$ and $\langle\xi(s,t)\xi(s^{\prime},t^{\prime})\rangle=2\Lambda\delta(s-s^{\prime})\delta(t-t^{\prime})$, where $\Lambda$ is an effective temperature. These governing equations are solved numerically as outlined in Chakrabarti and Saintillan (2019).

Spontaneous oscillations.—We first describe the dynamics of isolated filaments, with model parameters estimated from experiments Riedel-Kruse et al. (2007); Oriola et al. (2017); Chakrabarti and Saintillan (2019); see Supplemental Information Not for details. With a choice of $L\sim 50\,\mu$m for human sperm, $B\sim 0.9-1.7\times 10^{-21}\,$N m2, $f_{0}\sim 1-5\,$pN, $\tau_{0}\sim 50\,$ms and $\rho\sim 10^{3}\,\mu$m-1, we estimate $\mathrm{Sp}\sim 8-20$ and $\mu_{a}\sim 2-10\times 10^{3}$ and explore beating patterns in this range. For a given sperm number, a Hopf bifurcation occurs beyond a critical activity level $\mu_{a}^{c}$ and gives rise to spontaneous traveling waves Not ; Chakrabarti and Saintillan (2019). Close to the bifurcation, the waves propagate from the free end towards the base as previously seen in other simulations of sliding control models Oriola et al. (2017); Riedel-Kruse et al. (2007) but in disagreement with typical sperm beating patterns. However, far from the bifurcation, nonlinearities give rise to a reversal in the direction of propagation Chakrabarti and Saintillan (2019), with sperm-like waveforms shown in Fig. 1(b) that resemble experiments Brokaw (1965) and have beating frequencies $f\sim 10-15\,$Hz. In the following discussion, we focus on this anterograde propagation regime as it is biologically most relevant.

Asymmetric beating patterns more typical of cilia can be captured by setting different attachment and detachment rates for the motor populations on $\mathbf{x}_{\pm}$ Chakrabarti and Saintillan (2019). This bias in the kinetics allows the flagellum to bend in one direction preferentially, resulting in asymmetric power and recovery strokes as shown in Fig. 1(c). The flagella of wildtype Chlamydomonas also have a static mode of deformation Sartori et al. (2016) that we account for using a spontaneous shape $\phi^{0}(s)$. To better approximate their asymmetric breaststrokes [Fig. 1(d)], a curvature control mechanism is introduced along with the biased kinetics Chakrabarti and Saintillan (2019) that uses a generalized Bell’s law for the dynein detachment rate: $\epsilon_{\pm}=\epsilon_{0}n_{\pm}\exp\left(F_{\pm}/f_{c}\pm(\kappa(s)-\phi^{0}_{s}(s))/\kappa_{c}\right)$, where $\kappa(s)$ is the curvature and $\kappa_{c}$ is the threshold value for rapid dissociation. The subtraction of the zero mode in the curvature control follows Sartori et al. Sartori et al. (2016), who suggested that motor forces respond to derivatives of curvature rather than curvature itself. Accounting for the short length $L\sim 6-15\mu$m Sartori et al. (2016) of cilia and Chlamydomonas flagella with $B\sim 0.5-5\times 10^{-22}\,$N m2, we estimate $\mathrm{Sp}\sim 2-3$ and measure spontaneous frequencies of $f\sim 10-20$ Hz.

*Pair synchronization.—*We first focus on the synchronization of pairs of sperms placed side by side as shown in Fig. 1(e). We initialize the simulation in absence of inter-filament hydrodynamic interactions (HI) by letting spontaneous oscillations reach steady state after saturation of dynein kinetics. The initial configuration is chosen such that the filaments are almost in antiphase (AP) [top panel of Fig. 1(e)]. We then switch on HI and, after several periods, the sperms go in-phase (IP) and remain phase-locked thereafter [bottom panel of Fig. 1(e); see movies in Not ]. The key role of hydrodynamics in this process is best illustrated by Fig. 2, showing the evolution of the bound motor populations $\smash{n_{+}}$ at $s=1/4$ on both filaments (the behavior is identical for $n_{-}$ and at other locations). Before HI are switched on, motor populations are uncoupled and undergo periodic oscillations in antiphase with cusp-shaped waveforms typical of motors far from equilibrium Jülicher and Prost (1995) and only a small fraction of bound motors at any given time. Once HI start acting, both the phase and amplitude of the motor populations change. This is attributed to elastic deformations of the filaments in their induced flow fields, which feed back to the kinetics through the change in sliding displacement and velocity. As seen in Fig. 2, the two motor populations rapidly go in phase with a marginally increased amplitude, resulting in spontaneous IP synchronization of the beating patterns. The cartoon in Fig. 2 highlights this cyclic process fundamental to elastohydrodynamic synchronization, by which HI affect beating patterns via geometry-dependent motor kinetics. This feedback is most dramatic when the filaments are closeby and sufficiently flexible.

A similar mechanism is at play for asymmetric ciliary beats in Fig. 1(f, g). When the power strokes of the two cilia indicated by red arrows point in the same direction, an IP beat emerges with net unidirectional pumping of the fluid [Fig. 1(f)]. When the power strokes are in opposite directions, our model leads to AP synchronization with beating patterns resembling a ‘freestyle’ swimming gait as shown in Fig. 1(g). Similar AP patterns are obtained for Chlamydomonas beats. These observations hint at the hypothesis Wan and Goldstein (2016) that the IP breaststrokes seen in wildtype cells result from elastic basal couplings between the two flagellar axonemes rather than from HI alone. Indeed, experiments with vfl mutants that are deficient in these filamentary connections Wan and Goldstein (2016) or with Volvox cells held in separate micropipettes Brumley et al. (2014) have shown AP synchronization for power strokes with opposite orientations, consistent with our model findings. Note that in the case of swimming or even weakly clamped cells flagellar synchronization can also happen through a rocking motion of the cell body independent of HI or in absence of basal coupling Friedrich and Jülicher (2012); Polotzek and Friedrich (2013); Geyer et al. (2013); Bennett and Golestanian (2013a, b). The relative importance of these mechanisms remains to be explored for the various asymmetric waveforms Bennett and Golestanian (2013b) arising in our model.

For a more quantitative analysis of synchronization, we introduce a definition of the phase $\psi$ of a waveform. To this end, we perform the Hilbert transform of the continuous periodic time series $\beta(t)=\phi(1/2,t)$, providing the analytic continuation $\smash{\zeta(t)=\beta(t)+\mathrm{i}\widehat{\beta}(t)}$ where \smash{\widehat{\beta}(t)\equiv(1/\pi)\mathchoice{{\vbox{\hbox{\textstyle-}}\kern-4.86108pt}}{{\vbox{\hbox{\scriptstyle-}}\kern-3.25pt}}{{\vbox{\hbox{\scriptscriptstyle-}}\kern-2.29166pt}}{{\vbox{\hbox{\scriptscriptstyle-}}\kern-1.875pt}}\!\int_{-\infty}^{\infty}\beta(\tau)/(t-\tau)\,d\tau}. The phase of the waveform is then calculcated as ${{\psi}(t)=\arctan[\widehat{\beta}(t)/\beta(t)]}$, and we use an appropriate geometric gauge to define a true phase that grows monotonically with time Kralemann et al. (2008). The phase difference $\delta(t)=\psi_{1}-\psi_{2}$ for two nearby sperms going from AP to IP is shown in Fig. 3(a) and decays to zero over the course of several periods. In spite of the complexity of the governing equations in presence of HI, the phase difference is well described by a simple low-dimensional Adler equation as in past experiments with Chlamydomonas Goldstein et al. (2009) and in minimal rotor models Niedermayer et al. (2008); Uchida and Golestanian (2011); Maestro et al. (2018); Not . Here, we seek a two-parameter equation of the form

[TABLE]

where constants $\epsilon,\alpha$ are estimated numerically. A solution to this equation follows the numerical data very well in Fig. 3(a). In all our computations, we find that $|\epsilon|\gg|\alpha|$ and thus define $|\epsilon|$ as the effective coupling strength in accordance with minimal models of rotors Not . When plotted as a function of interflagellar distance $d$ in Fig. 3(b), $|\epsilon|$ shows a far-field algebraic decay of $1/d$ over the limited range of accessible values, which is a signature of the dominant Stokeslet HI and can be rationalized from a simple rotor model Not . A slower decay is seen at short separations, where complex near-field interactions take place. Stronger coupling arises for symmetric spermlike beats than for ciliary beats, primarily due to the longer lengths of sperm flagella. For cilia, we also find that $|\epsilon|_{\mathrm{IP}}>|\epsilon|_{\mathrm{AP}}$ in agreement with experiments Brumley et al. (2014), which can be attributed to the fact that filaments spend more time close to one another during IP beats and thus interact more strongly.

Intrinsic to the kinetics of molecular motors is biochemical noise, which alters the precise notion of synchronization. To probe its effects, we study the long-time statistics of the phase difference in presence of noise for spermlike waveforms in Fig. 4(a). Fluctuations follow a Gaussian distribution centered around the mean IP configuration of $\delta=0$, with a variance scaling linearly with separation distance $d$. This is a consequence of the $1/d$ decay of the coupling strength $|\epsilon|$ and is further corroborated by the collapse of the distributions under the rescaling $\delta\to\delta/\sqrt{d}$ in Fig. 4(b) Brumley et al. (2014). We model the noisy phase dynamics by a stochastic Adler equation $\dot{\delta}=\epsilon\sin\delta+\chi(t)$ with $\langle\chi(t)\rangle=0$ and $\langle\chi(t)\chi(t^{\prime})\rangle=2D\delta(t-t^{\prime})$, where $D$ is the phase diffusivity with units of s*-1*. Associated with the Adler equation is a Fokker-Planck description for the probability distribution $P(\delta)$ of the phase difference, with steady-state solution given by $P(\delta)=\exp(-\epsilon\cos\delta/D)/2\pi I_{0}(|\epsilon|/D)$, where $I_{0}$ is the modified Bessel function of order zero and where we estimate $D$ numerically Chakrabarti and Saintillan (2019); Stratonovich (1967). The interaction potential $U(\delta)=-\ln P(\delta)$, which is $2\pi$-periodic, is shown in Fig. 4(d) for increasing noise levels. When noise is weak, the filaments remain phase-locked and fluctuate around the IP configuration, which translates into a deep potential well at $\delta=0$. With increasing noise, the potential well flattens as deviations from perfect AP synchrony become more frequent and intense. Occasionally, accumulated noise allows the filaments to gather a complete phase of $2\pi$, causing them to ‘slip’ towards $\delta\pm 2\pi$. These slips are visible in the phase trajectories of Fig. 4(c) and can be interpreted as thermally assisted hops between neighboring wells in the flattened periodic potential. In absence of frequency mismatch, slips are equally probable in $\pm 2\pi$, and the stochastic Adler model predicts a frequency of $G=D|I_{0}(|\epsilon|/D)|^{-2}/4\pi^{2}$ Stratonovich (1967).  Using the computed value of $D$, this prediction indeed provides a quantitative estimate of the mean frequency of slips in full simulations; see Supplemental Information Not .

Concluding remarks.— We have used an idealized planar model of the flagellar axoneme that captures the essential physics of internal dynein activity and produces spontaneous oscillations similar to those seen in nature Not to elucidate elastohydrodynamic synchronization of nearby flagella and cilia. Our simulations underscore the essential roles of hydrodynamic interactions and associated mechanochemical feedback in enabling synchronization. Our model predictions for various beating patterns and orientations all agree with experiments and give credence to a combination of sliding and curvature control mechanisms for the generation of spontaneous beats. We were also able to reproduce experimentally observed phase slips induced by biochemical noise. Future studies with our model will probe the role of elastic basal couplings Wan and Goldstein (2016), swimming cells that are free to adjust phase by sliding past one another Yang et al. (2008) or by rotational motion of their body Geyer et al. (2013); Bennett and Golestanian (2013b), and emergent dynamics in large-scale ciliary arrays.