Elastohydrodynamic Synchronization of Adjacent Beating Flagella

It is now well established that nearby beating pairs of eukaryotic flagella or cilia typically synchronize in phase. A substantial body of evidence supports the hypothesis that hydrodynamic coupling between the active filaments, combined with waveform compliance, provides a robust mechanism for synchrony. This elastohydrodynamic mechanism has been incorporated into `bead-spring' models in which flagella are represented by microspheres tethered by radial springs as internal forces drive them about orbits. While these low-dimensional models reproduce the phenomenon of synchrony, their parameters are not readily relatable to those of flagella. More realistic models which reflect the elasticity of the axonemes and active force generation take the form of fourth-order nonlinear PDEs. While computational studies have shown synchrony, the effects of hydrodynamic coupling between nearby filaments governed by such models have been theoretically examined only in the regime of interflagellar distances $d$ large compared to flagellar length $L$. Yet, in many biological situations $d/L \ll 1$. Here, we first present an asymptotic analysis of the hydrodynamic coupling between two filaments in the regime $d/L \ll 1$, and find that the form of the coupling is independent of the details of the internal forces that govern the motion of the filaments. The analysis is like the localized induction approximation for vortex filament motion, extended to the case of mutual induction. To understand how the coupling mechanism leads to synchrony of extended objects, we introduce a heuristic model of flagellar beating, a single fourth-order nonlinear PDE whose form is derived from symmetry considerations, the physics of elasticity, and the overdamped nature of the dynamics. Analytical and numerical studies of this model illustrate how synchrony between two filaments is achieved through the asymptotic coupling.


I. INTRODUCTION
In nearly all of the contexts in biology in which groups of cilia or flagella are found they exhibit some form of synchronized behavior.At the level of unicellular organisms this often takes the form of precise phase synchrony as in the breaststroke beating of biflagellated green algae [1], but it has also been known since the work of Rothschild [2] that swimming sperm cells can synchronize the beating of their tails when they are in close proximity.In multicellular organisms such as Paramecium [3] and Volvox [4,5], and in the respiratory and reproductive systems of higher animals one often observes metachronal waves, which are long-wavelength modulations in the beating of ciliary carpets.There are three primary dynamical behaviors of ciliary groups, two involving beating waveforms that have a 'power stroke' in which the filament pivots as a nearly straight rod, followed by a 'recovery stroke' in which it is strongly curved, and a third in which the flagella beating is undulatory.In the first two cases it is useful to categorize the different geometries on the basis of the orientation of the power strokes of adjacent flagella.If we follow reference points along each flagellum -say, each center of mass -then they will move with either parallel (cilia) or anti-parallel (biflagellate) angular velocities (Figs.1a & b).In the undulatory case (Fig. 1c), found in mutants of Chlamydomonas and during the 'photoshock response' [6], nearby flagella beat parallel to each other.
Based on the observations of Rothschild [2] on synchronized swimming of nearby sperm cells, Taylor [7] investigated the possibility that hydrodynamic interactions could lead to synchrony.His 'waving sheet model' considered two infinite parallel sheets each in the shape of a prescribed unidirectional sinusoidal traveling wave.Examining the rate of viscous dissipation as a function of the phase shift between the two waves, he found that the synchronized state had the least dissipation.While highly plausible as an explanation of synchronization, this model does not include a dynamical mechanism by which the synchronized state is achieved from arbitrary initial conditions.Subsequent work [8] has shown that adding waveform flexibility to the model yields a true dynamical evolution toward synchrony, and this has been confirmed by experiment [9].The recognition that hydrodynamic interactions alone are insufficient to generate dynamical evolution toward synchrony, and that some form of generalized flexibility is necessary had already been seen in the study of rotating helices as a model for bacterial flagella [10].This notion of 'orbital compliance' was subsequently incorporated into several variants of bead-spring models [11][12][13] of ciliary dynamics in which each beating filament is replaced by a moving microsphere which is driven around an orbit by internal forces and allowed to deviate by a radial spring.Under the assumption that radial motions evolve rapidly relative to azimuthal ones, these models generically yield a nonlinear ODE for the phase difference between the oscillators that takes the form of the Adler equation [14].
While these models lead to a microscopic interpretation of the generic Adler equation, they are most appropriate for the situation in which the distance d between the flagella is large compared to their length L, where a far-field description in terms of Stokeslets is valid [15].But many of the most interesting situations, such as the parallel and undulating geometries in Figs.1b&c, are in precisely the opposite limit, d/L 1, while still in the regime a d, where a is the filament radius.In this limit it is clear that a proper description of the entire filament is necessary, because there are very strong near-field interactions all along them, and therefore a representation by a single point force would not be realistic.While computations incorporating microscopic models of flagella embedded in a viscous fluid show that synchronization does indeed occur through hydrodynamic interactions in this regime [16,17], it was only in subsequent work that the formally exact nonlocal description of hydrodynamic interactions in multi-filament systems was presented [18].Taking advantage of the separation of scales a d L, we present in Section II an asymptotic derivation of the leading-order hydrodynamic coupling between two filaments.In particular, we find that the relevant small coupling parameter is = ln(L/d) ln(L/a) .The analysis leading to this result is reminiscent of the 'localized induction approximation' in vortex filament dynamics [19].
At present, there is no single generally-accepted microscopic model for eukaryotic flagellar beating, although recent studies have begun to address the relative merits of several promising candidate models [20][21][22], which typically consist of a pair of coupled equations, one for the filament displacement, incorporating filament bending elasticity and viscous drag, and the other for the active bending forces associated with molecular motors.To illustrate how the hydrodynamic coupling and waveform compliance lead to synchronization, in Section III we introduce and analyze a heuristic single PDE, of the form where N is a nonlinear operator, which displays selfsustained finite-amplitude wavelike solutions.Section IV considers, both analytically and numerically, the dynamics of two filaments of the type introduced in Section III interacting through the coupling derived in Section II.The concluding Section V discusses possible applications of the model.

II. ASYMPTOTICS
We consider two slender filaments of length L undergoing some waving motion, as shown in Fig. 2. Their mean separation is d and we assume that their waving amplitude is at most of order d.We seek to derive, in the linear regime of small amplitude displacements, the forces resulting from hydrodynamic interactions in the asymptotic limit d L.
Let h 1 and h 2 denote the vertical displacements of the filaments from their mean positions.We present the derivation of the force on one filament only, say filament 1, as the dynamics of the other can be deduced by symmetry.In the linear regime, it is only necessary to consider the balance of forces in the vertical direction.Using the framework of resistive-force theory (RFT) [23], the vertical component F 1 of the hydrodynamic force per unit length acting at the point (x 1 , h 1 ) of filament 1 is where ζ ⊥ is the drag coefficient for motion normal to the filament.We now proceed to calculate the flow u induced by filament 2. In the linear regime, this flow arises from a superposition of hydrodynamic point forces acting in the y direction along filament 2. Neglecting end effects, the flow u 2→1 y can thus be written as the following integral of Stokeslets: where µ is the dynamic viscosity, 1 the identity, and r the vector that points from (x 2 , d + h 2 ) on filament 2 towards (x 1 , h 1 ) on filament 1 (Fig. 1).In the linear regime, the force density is simply given by that from RFT as where the positive sign indicates that the force f 2 acts on the fluid.We thus have into Eq.( 4) and linearizing for small h 1 and h 2 , we obtain In order to compute the asymptotic value of Eq. ( 5) in the limit d L we introduce the dimensionless lengths , and obtain where for convenience we have dropped the prime on the integration variable, while all other quantities in (6) By construction |x 1 − x 2 | is at least δ 1 and thus it is possible to neglect 1 in the nonlocal integral to obtain Provided that x 1 (1 − x 1 ) δ 2 , then in the limit δ → 0, the integral in Eq. ( 8) diverges logarithmically as Similarly, and with the change of variables ∆ = x 2 − x 1 , the local part of the integral can be written as where h 2 is now understood to be a function of ∆.In the limit δ → 0, the term ∂h 2 /∂t takes its value near ∆ = 0, and the remaining integral can be computed exactly, yielding Because 1 δ the result to order 1 reads Remarkably, the ln δ divergence from Eq. ( 9) exactly cancels out the one from Eq. ( 12) producing a result which is independent of the particular choice made for the cutoff δ.Thus, the final expression for I reads Thus, the vertical component of the hydrodynamic force on filament 1 is where we have used ζ ⊥ ∼ 4πµ/ ln(L/a), with a the radius of the flagella, and As in nearly all applications of slender body hydrodynamics, the parameter is asymptotically small only in the unphysical case when the aspect ratio is exponentially large.However, it is well-known that use of the leading order approximation of slender body hydrodynamics for larger values of is robust [24,25].
Finally, it is important to note that even in the case when the filaments are in phase and h 1 (x 1 , t) = h 2 (x 2 , t) = h(x, t) everywhere, and where it is not appropriate to evaluate Eq. ( 16) at close contact between the two flagellar filaments (d = 2a).because the induced flow was computed as a superposition of Stokeslets only.This is a good approximation only when all other singularities present have decayed away, in particular the (potential) 1/r 3 source dipole which arises from the finite-size of the flagella.Thus, the approximation requires that d a from every point of flagella 1 to every point in flagella 2 (and vice-versa).In other words, the result ( 14) is valid only within the limits a d L. For example, for eukaryotic flagella with L ∼ 50 µm and a ∼ 0.1 µm, then the analysis is valid when the flagella are separated by a few microns, in which case decreases from 0.5 for d ∼ 2 µm to 0.25 for d ∼ 10 µm.
The results above imply that if each of two nearby filaments is governed by an equation of the form where {c i } are the parameters that differentiate the flagella, then the coupled pair evolves according to As we have only computed the hydrodynamic interaction to order , it is appropriate to consider the leading-order form of ( 17) as ∂h i /∂t = N ci (h i ) + N cj (h j ) for i, j = 1, 2 and i = j.

III. PHENOMENOLOGICAL MODEL OF A SINGLE FILAMENT
A. Background and model Here, to represent the situation in which a flagellum is attached to an organism's surface, we focus on the case of a finite beating filament, say 1, pinned at its left end to a fixed support, with a free right end.As with all models for systems of this type, and the analysis in the previous section, we focus on low Reynolds number dynamics.The structure of the most general equation of motion for a filament arising from balancing its tangential and normal forces and bending moments is well known [20,21,26,27].Under the further assumption of linear filament elasticity and resistive force theory, and assuming that the filament deviates only slightly from straight, the linearized equation of motion for the tangent angle ψ(s, t) as a function of the arc length s and time t takes the form ζ ⊥ ψ t = af ss − EIψ ssss , where E is the Young's modulus and I the moment of inertia, per unit density, of the filament cross section about the axis of rotation, and f is the active bending moment.Recognizing that within this approximation ψ ∂ x h 1 we obtain where A = EI is the bending modulus of the filament.The distinction between different models of active bending is to be found in the particular form of f , which may be coupled back to the geometry (e.g.ψ and its derivatives) through an auxiliary equation of motion.
Recent work [21] has studied the linearized dynamics of the unstable modes that arise in three models of the form (18), known as sliding-control [28], curvature-control [27], and the geometric clutch [29].The sliding control model, whose equation of motion does not explicitly break left-right symmetry, was shown not to exhibit base-to-tip propagating modes of the kind seen in experiment.In contrast, both the curvature-control and geometric clutch models, which do display modes with the qualitatively correct behavior, have a symmetry-breaking term ∂ xxx h 1 .
Furthermore, as the dynamics is translational invariant in y, there can be no terms explicitly dependent on h 1 itself.
Motivated by these results, and with an eye toward the simplest PDE that will encode a characteristic wavelength, amplitude, and frequency of flagellar beating, we propose a form in which the left-right symmetry is broken with the derivative of the lowest order possible, where G is a nonlinear amplitude-stabilizing function of κ ∂ xx h 1 , the filament curvature.The presence of this advective term does not reflect any explicit fluid flow, rather it encodes processes internal to the filament that break symmetry.If we assume that the filament has h 1 → −h 1 symmetry, then G will be an odd function of its argument.
Expressing G as a Taylor expansion up to cubic order we arrive at the model of interest, henceforth called the advective flagella (AF) model where A, B, C and D are heuristic parameters of the model.We can now introduce dimensionless variables ξ = αx, τ = t and h = Hh 1 , where α, and H are constants.Direct substitution of these new variables into (20), with the , and the dimensionless PDE where 0 ≤ ξ ≤ Λ, with Λ = αL.The length c is the well-known elastohydrodynamic penetration length that arises in the study of actuated elastic filaments [30,31], and thus Λ is the so-called Sperm number.Of the many possible boundary conditions we adopt the simplest: hinged at the left end (h(0, τ ) = h ξξ (0, τ ) = 0) and free at the right end (h ξξ (Λ, τ ) = h ξξξ (Λ, τ ) = 0).In this rescaled form, the dynamics of a single filament is completely specified by c and Λ.Note that the dynamics ( 21) does not enforce filament length conservation beyond linear order, but this should not introduce any qualitative changes in the results below.
The linear operator in (21) is that of the advective version [32] of the Swift-Hohenberg model [33] without the term linear in h, and is identical to that in the phenomenological model of dendrite dynamics introduced by Langer and Müller-Krumbhaar [34] and studied by Fabbiane, et al. in the context of control theory [35].Heuristically, we recognize that the second and fourth derivatives naturally select a most unstable length scale for a linear instability of the state h = 0, the advective term leads to rightward wave propagation, and the nonlinearity leads to amplitude saturation.The intuitive understanding that these are the minimum necessary, but also sufficient, ingredients for a model of the beating of a single flagellum, is proven to be correct by the numerical and analytic work presented in the following subsections.Moreover, in Section IV, we will also show that the dual features of linear terms −2h ξξ − h ξξξξ which lead to a band of unstable modes, and amplitude saturation through nonlinearity will allow for waveform compliance when two filaments are hydrodynamically coupled.

B. Numerical studies of a single filament
In this section we present the dynamical behavior of the AF model obtained through numerical computations.
The model (21) was solved numerically using an implicit finite-difference scheme described in detail by Tornberg and Shelley [18], including one-sided stencils for derivatives at the filament endpoints.We shall see that for any choice of c and Λ the model selects a wavelength λ (wavenumber k = 2π/λ), frequency ω and maximum amplitude A. Of particular interest are the limiting regimes Λ/λ ∼ 1 and Λ/λ 1.The former is the regime seen in many experiments (Fig. 1), while the latter corresponds to a semi-infinite system whose behavior far from the pinning wall approximates a traveling wave.
Following a short transient, we find that the filament evolves toward a periodic waveform of nonuniform, finite amplitude.Figure 3a shows overlays of h(ξ, τ ) at various times during a full cycle for Λ = 15 and three values of c.
We observe that for this length, and c = 1 (top waveform) the oscillation amplitude has not reached saturation at the free end.Thus, even larger values of c (middle and bottom) naturally advect the waveform even faster to the right, resulting in a smaller maximum amplitude at the free end.As seen by comparing the two examples in Fig. 3b with the top in Fig. 3a, where all three filaments have c = 1 but different length, for Λ = 10 the amplitude envelope is clearly linear in ξ, for S = 15 there is the beginning of a rollover, and for Λ = 30 there is clear saturation.For Λ/λ ∼ 1 the maximum amplitude of the waveform is reached at the free end of the filament, as seen in the natural biological waveforms and also in the models described above [21].In general, for a given value of c, as Λ grows relative to λ, the filament displays two distinct regions: a transition region adjacent to the wall where the amplitude grows along ξ up to a characteristic L, and a region beyond where the amplitude saturates and the oscillations approximate a traveling wave.The 'healing length' L increases with c.
In the regime Λ/λ 1, and far from the wall, the amplitude A, wave vector k, and frequency ω of the approximate traveling wave are determined by c only.The solid symbols in Figure 4 indicate those numerical results.To put these results in context, note that a linear stability analysis of the operator −2h ξξ − h 4ξ leads to a growth rate which is maximized at k * = 1 where, without loss of generality, we consider only the positive branch of solutions.We see from Fig. 4a that the value of k selected by the system is always less than 1, but k → 1 as c increases, and also, from Fig. 4b, that the selected frequency is consistent with the relation ω ck.Finally, the saturated amplitude A of the waveform far from the wall is a strongly decreasing function of k (Fig. 4c).
C. Approximate analytical solution

Amplitude saturation
Far away from the wall the solution of ( 21) is well approximated by a traveling wave with a time dependent amplitude of the form Substituting ( 22) into (21), rewriting the cubic cosine term with the multiple angle formula, ignoring terms involving O[3(kξ − ωτ )], and equating the coefficients of the sine and cosine to zero yields an evolution equation for A(τ ) and a relationship between the wavenumber k and the angular frequency ω, namely The evolution equation for the amplitude A(τ ) is of the Bernoulli type and can be easily solved after making the change of variables z = A −2 , with solution of the form z = Be βτ + γ, where B, β, and γ are constants, corresponding to a time dependent amplitude given by Provided 0 < k < √ 2, A has a finite value as τ → ∞, namely A ∞ (k) is shown as the solid line in Fig. 4c, and it is in excellent agreement with the numerical data for A expressed now as a function of the numerically selected k.We note that the steady-state solution can be found exactly in terms of elliptic integrals by transforming into a moving coordinate system with X = ξ − cτ , h(ξ, τ ) → F (X) and solving However, the slight increase in accuracy that this approach produces comes at the expense of a lack of clarity compared to the one-mode approximation (22).

Wavenumber and frequency selection
While the traveling wave solution is a good approximation far away from the wall, it only provides an asymptotic relationship between the wavenumber and the frequency but does not yield any information about the possible numerical values of k and ω.In order to find these numerical values the linearized evolution equation may be used.Besides the trivial solution h(ξ, τ ) = h 0 where h 0 is a constant, (26) has a solution of the form h(ξ, τ ) = exp(Ωτ + Kξ), where Ω = σ + iω and K = p + ik.In the present case, the only solutions of interest are those with a positive temporal growth rate σ > 0 because any other solution will decay to h(ξ, τ ) = h 0 .Moreover, we expect the solution corresponding to the maximum growth rate to be the one that will control the evolution of the shape h(ξ, τ ).Hence, to find the extremal values of K and Ω we replace h(ξ, τ ) into ( 26) and then, from the characteristic equation, calculate the derivative of Ω respect of K: After separating real and imaginary parts they become: where we have already set the real part of the derivative of Ω equal to zero and note that the complex part of the derivative is identical to calculating ∂ω /∂p.In general, ∂ω /∂p is also set to zero to find the boundary between convective and absolute instability [36,37].Here this second constraint has been relaxed so that it is possible to find the curve ω (k) along which the system is most unstable.In fact, rewriting equation (28c) as −12pk 2 = −c − 4p − 4p 3 and replacing this expression into (28b) yields ∂ω /∂p = −8k 3 .Equating the two expressions for ∂ω /∂p then gives Replacing ( 29) back into equations (28c-a) gives the expressions for the velocity c, the angular frequency ω = −ω and the growth rate σ as functions of the wave number k: From these results it becomes clear that the only solutions which satisfy the hypothesis for h(x, t) are those with values k ≥ 1/ √ 3 because for values of k < 1/ √ 3 both, c and ω are complex.In particular, the value k = 1/ √ 3 corresponds to c = ω = 0.As seen in Fig. 4b the numerical counterpart to this boundary is k = 1/2.While this represents a discrepancy of approximately 15 percent between the analytic prediction and numerical results, it can be seen in the figure that the basic trend in the selected frequency is captured.Figure 4a shows the comparison between the numerical k(c) and the analytical result above, and in this case larger the value of c the closer the two become.Because the analytic expression (25) for the maximum amplitude as a function of k was obtained using the full nonlinear equation, this prediction shows almost perfect agreement with the numerical results (Fig. 4c).
Another value of interest is the critical k * with zero growth rate, and which can be found by solving the bi-quadratic equation 4k 4 − 4k 2 + 1 3 = 0 which has only one real root above the threshold for k: All solutions in the interval 1/ where σ(k * ) = 0 and σ(1/ √ 2) = 2/3; while all solutions for which k > k * have a negative growth rate and, hence, relax back to the solution h = h 0 .The critical angular frequency and velocity that correspond to σ(k * ) = 0 are The numerical value of the velocity for which the oscillatory solutions "disappear" is c * 3.2, and the analytic prediction (32b) is well within the acceptable limits of agreement for the approximation.Moreover, the analytic prediction for the functional form of the ratio ω(k)/c(k) is qualitatively correct and given by which as expected, tends to ω/c → k as k grows.

Frequency-amplitude relation
Before proceeding to numerical studies of coupled filaments we are in a position to see how synchronization of two nearby filaments may occur.First, we note that when two filaments beat in synchrony the fluid gap between them is nearly constant, whereas when they are out of synchrony the gap varies with position.That variation produces fluid forces that will deform the filaments such that their local amplitude and frequency will be altered.As seen in bead-spring models [11], for example, a stable synchronous state may occur when the frequency is a decreasing function of amplitude.In the case of the approximate traveling-wave states we have discussed in this section, the results above can be used directly to calculate ω(A), as shown in Fig. 5a, To see how the synchronization mechanism operates within the present model, and anticipating the numerical results presented below, consider the configuration of two filaments shown in Fig. 5b, each traveling to the right, where the black arrows indicate the local direction of motion of the filaments at two distinct coordinates, and the blue arrows indicate the direction of the fluid flow induced by filament 1 (h 1 (x, t)) on 2. At the point x a the fluid flow acting on filament 2 will push it further down, whereas at point x b that flow will pull it upwards.The net effect is that the local wave amplitude of filament 2 will be decreased, and by the relationship in Fig. 5a its frequency will increase, moving it faster to the right and hence catching up with filament 1.Similar considerations show that the effect of filament 2 on 1 is to increase its amplitude, hence to decrease its frequency.Thus, the stable state is the in-phase synchronized one.This elastohydrodynamic mechanism is the continuum analog of that which operates in bead-spring models.
IV. TWO COUPLED FILAMENTS

A. Symmetry and stability considerations
In this section we discuss the coupled dynamics of two nearby filaments using the model of Section III with the coupling derived in Section II.When two coupled filaments have the same intrinsic speeds c 1 = c 2 = c then they obey where, following Eq.21, the nonlinear operator N c ( ĥi ) is By direct substitution, it is straightforward to show that Eqs.35, with the operator 36, have two exact solutions S 1 : ( ĥ1 (ξ, τ ), d + ĥ1 (ξ, τ )) (sinuous solution), and S 2 : ( ĥ1 (ξ, τ ), d − ĥ1 (ξ, τ )) (varicose solution), where d is a dimensionless constant, provided ĥ1 (ξ, τ ) satisfies the nonlinear, autonomous equation which is the rescaled PDE that governs a single isolated filament.Note that here the plus (minus) sign corresponds to the sinuous (varicose) configuration.The effective times τ = τ /(1 ∓ ) correspond to faster motion than an isolated filament in the sinuous case and slower motion in the varicose arrangement.This can be understood as a result of decreased and increased viscous dissipation in the two cases, respectively, consistent with the results of the waving-sheet model.
To study the stability of these solutions a small perturbation is introduced.Because our focus is on the difference between the filament positions, it is sufficient to examine perturbed solutions of the form ĥ1 (ξ, τ ), d ± ĥ1 (ξ, τ ) + η(ξ, τ ) where |η(ξ, τ )/ ĥ1 (ξ, τ )| 1.Then, the linearized equation of motion for η is Note that here the upper (lower) sign corresponds to the sinuous (varicose) configuration.The dynamics (38) is close in form to the linearized dynamics of a single filament, with one crucial difference: and in the original operator the coefficient of the second spatial derivative was negative (namely −2) which corresponds to the "anti"-diffusion that produces the instability, whereas now, the base solution parametrically forces the perturbation through the coefficient of the second derivative.Thus, depending on the characteristics of the solution ĥ1 , and only in the sinuous case, it is possible to find regions with a positive effective diffusion coefficient.
As pointed out in the previous section, the asymptotic solution in time of (37) far away from the origin is a traveling as the average coefficient of the diffusive term in Eq. 38.This coefficient is positive (hence stabilizing) for the sinuous configuration and negative (destabilizing) for the varicose one.In the presence of the always stabilizing influence of the fourth derivative, this implies that the only linearly stable configuration is the sinuous one.
When the speeds of the filaments are slightly different, c 1 = c − (∆c/2) and c 2 = c + (∆c/2) with ∆c/c 1, a similar analysis to the one described above makes it possible to find an approximate solution that corresponds to a quasi-sinuous configuration, and for which the first order solution is the same as when ∆c = 0.The initial system that governs the motion is and the proposed solution is S : ( ĥ1 = h − η, ĥ2 = h + η) where η/ h 1. Direct substitution and some simple algebra yield the evolution equation for h, ( Keeping only linear terms in η this reduces to the original autonomous one for a single filament, but with coefficient 1 − in front of the temporal derivative.The associated linearized equation for η is similar to (38), but also has a forcing Replacing into (42) the traveling wave solution (22) we obtain Notice that the diffusion coefficient and the forcing are out of phase.When | cos(kξ − ωτ )| is large, the diffusion coefficient is positive and η decreases making ĥ1 ĥ2 ; these are the regions where the forcing is very small and tends to leave the system unperturbed.In the regions of space and time where the cosine is very small and the diffusion coefficient takes a negative value, the forcing is strong and in opposition to the growth of η, and even though the filaments can never be in absolute synchrony, as when the two velocities are the same, the forcing keeps the asynchrony to a minimum.

B. Numerical results
Numerical studies of the coupled dynamics (35) of two filaments with the same speed parameter c show robust synchronization for a broad range of initial conditions.Figure 6 shows the evolution toward synchrony of a pair of A convenient quantity to characterize the degree of synchrony of two filaments is the L 2 norm of the difference of their displacements, D = ĥ1 (ξ, τ ) − ĥ2 (ξ, τ ) 2 .Figure 7 shows the temporal behavior of D displacements, for three  different values of the coupling constant .While the rate of synchronization increases with , and the details of the decay of the norm differ, synchrony occurs in all cases.Further analysis shows that, when averaged over the fast oscillation, the approach to synchrony is exponential, as would be expected from the fact that the linearized dynamics close to synchrony is first order in time.In the physical regime, with ∼ 0.2-0.5, we see in Fig. 7 that synchrony occurs in a matter of a few beat cycles.This rapid synchronization is often seen in biological systems, including Chlamydomonas flagella subjected to hydrodynamic perturbations [38] and during the photoshock response [39].Finally, we discuss the case with slightly different values of the speed parameter c.We expect that the two filaments will frequency lock but display a finite phase shift.This is borne out in the numerical studies, as shown in Fig. 8 for c 1 = 0.9 and c 2 = 1.0,where the upper filament (2) leads the lower one (1), while the two display identical frequencies.8. Synchronization with unequal advective coefficients.Overlaid waveforms when frequency locking has been achieved, for c1 = 0.9 and c2 = 1.0 and Λ = 10.The faster flagellum (upper) is phase-shifted forward with respect to the slower one.

V. DISCUSSION
In this work we have presented two main results.First, under the assumptions of coplanarity and small-amplitude deformations, we have derived the leading order coupling term that describes the hydrodynamic interaction between two nearby slender bodies whose separation d lies within the asymptotic limit a d L. This is a very general result formally expressed simply in terms of the velocities of each filament, whatever their microscopic origins.Second, we have applied this result to a model of flagellar beating that has the minimum necessary features to capture the essence of the system: self-sustained oscillations, broken left-right symmetry, bending elasticity, and waveform amplitude saturation.Analytical and numerical studies of the model for the case of single filaments illustrate the mechanisms of wavelength, frequency, and amplitude selection, while those for coupled pairs display the basic elastohydrodynamic mechanism of synchrony.
We emphasize that because of the extremely general form of the inter-filament coupling term, it can immediately be used to extend any particular model of the beating of a single filament to the coupled dynamics of two or more.
A worthwhile extension would be the generalization of this result to include the presence of a no-slip wall at which filaments are anchored, with the goal of understanding metachronal waves.Finally, we expect the approach to flagellar dynamics outlined here, based on long-wavelength expansions and symmetry considerations, and reduced to a single autonomous PDE, will prove useful in the analysis of experimental waveforms and dynamics.

FIG. 2 .
FIG. 2. Notation for the calculation.Two filaments, each of length L, are separated by a mean distance d.The displacements from the mean of two arbitrary points on the curves, separated by the vector r, are h1(x1, t) and h2(x2, t).

FIG. 3 .
FIG. 3. Numerical results for the AF model of a single filament.(a) Overlaid waveforms during a single beating period for filament length Λ = 15 and three values of c.(b) As in (a) but for fixed c = 1 and Λ = 10 and 30.

FIG. 4 .
FIG. 4. Waveform selection in the AF model.Numerical (red circles) and analytical (blue lines) results for long filaments: (a) selected wave vector k vs. advective parameter c, (b) selected frequency ω vs. k, and (c) saturated amplitude A vs. k.

FIG. 5 .
FIG. 5. Mechanism of synchronization.(a) Amplitude dependence of traveling-wave frequency ω.(b) Schematic of two nearby filaments, with black arrows indicating the direction of motion of points on the filaments, and blue arrows showing the direction of fluid motion induced by filament 1.

filaments with c = 1 ,
Λ = 10, and = 0.1, computed over a total time of T = 60.Panel (a-c) show the two filaments at early (2 ≤ τ ≤ 10.7), intermediate (16.3 ≤ τ ≤ 25.6), and late (48.7 ≤ τ ≤ 57.3) times and at four equally spaced time intervals within each corresponding cycle.From the clearly disordered pattern in (a) the filaments evolve to a fully-synchronized state in (c).

6 .
Numerical synchronization dynamics for two filaments with Λ = 10 and c = 1.Panels (a)-(c) show overlaid waveforms at four points within a single oscillation cycle at early (a), middle (b) and late (c) times.

7 .
Approach to synchrony.The L2 norm of the difference in waveform between two filaments as a function of time, for Λ = 10, c = 1, and indicated values of . FIG.
are still dimensional.This integral has two contributions: a local integral, I L , where x 2 is close to x 1 and a nonlocal one, I N L where |x 2 − x 1 | 1 .To evaluate I N L , we choose an intermediate length scale δ satisfying 1