Role of spatial heterogeneity on the collective dynamics of cilia beating in a minimal 1D model
RRole of spatial heterogeneity on the collective dynamics of cilia beating in a minimal1D model
Supravat Dey, ∗ Gladys Massiera, † and Estelle Pitard ‡ L2C, Univ Montpellier, CNRS, Montpellier, France. § Cilia are elastic hairlike protuberances of the cell membrane found in various unicellular organismsand in several tissues of most living organisms. In some tissues such as the airway tissues of the lung,the coordinated beating of cilia induce a fluid flow of crucial importance as it allows the continuouscleaning of our bronchia, known as mucociliary clearance. While most of the models addressingthe question of collective dynamics and metachronal wave consider homogeneous carpets of cilia,experimental observations rather show that cilia clusters are heterogeneously distributed over thetissue surface. The purpose of this paper is to investigate the role of spatial heterogeneity on thecoherent beating of cilia using a very simple one dimensional model for cilia known as the rowermodel. We systematically study systems consisting of a few rowers to hundreds of rowers and weinvestigate the conditions for the emergence of collective beating. When considering a small numberof rowers, a phase drift occurs, hence a bifurcation in beating frequency is observed as the distancebetween rowers clusters is changed. In the case of many rowers, a distribution of frequencies isobserved. We found in particular the pattern of the patchy structure that shows the best robustnessin collective beating behavior, as the density of cilia is varied over a wide range.
PACS numbers: 87.16.Qp,05.45.Xt,47.63.-b
I. INTRODUCTION
Cilia are elastic hairlike protuberances of the cell mem-brane found in various unicellular organisms and in sev-eral tissues of most living organisms. As a propulsor,a cilium is periodically beating in a succession of powerand recovery strokes, propelled by its internal molecularmotors [1]. Propulsion acts either on the microorgan-ism itself such as for
Paramecium or Volvox algae, orto induce the flow of the surrounding fluid, as this thecase for airway, brain or oviduct tissues [2–4]. In theairway tissues of the lung, this fluid flow induced by thecoordinated beating of cilia is of crucial importance asit allows the biological function of cilia to help to expelthe mucus and the impurities out of the airways, knownas mucociliary clearance [4].Cilia are often observed as scattered clumps and beatin coordinated manner in the form of metachronal waveskeeping a constant phase difference with adjacent cilia[5–7]. This type of large scale coordinated beating pat-tern is of great importance for efficient propulsion [8–12].It is now well established that hydrodynamic interactionsplay a crucial role [7, 10, 13, 14] for the emergence of suchlarge scale metachronal waves.While most of the models addressing the question ofcollective dynamics and metachronal wave consider ho-mogeneous carpets of cilia, in real samples, cilia formpatches. In cultured and in vivo airways epithelium tis-sues, cilia patches are heterogeneously distributed overthe surface [15, 16]. For example, observations per-formed on human bronchial epithelium cultures showthat cilia are distributed in clusters containing ∼ −
200 cilia [17], and that these clusters are separated ina random manner. The mucociliary dysfunction due to impaired coordinated beating of cilia is poorly under-stood and could be related to the spatial heterogeneityof healthy cilia distribution. The main focus of this workis to study the collective behavior of a system of manyhydrodynamically coupled beating cilia with heteroge-neous spatial configurations.Several theoretical models with various levels of com-plexity investigate the relation between the hydrody-namic coupling and the metachronal synchronization[14, 18]. We will concentrate our study on a minimalmodel, where the beating pattern is simple and is com-posed of a few degrees of freedoms allowing us to un-derstand the sole role of hydrodynamical coupling in thebeating synchronization of an array of cilia. There aretwo classes of well studied minimal models. In one classof models, a cilium is described as a rower: a spher-ical bead oscillates between two distinct states where,in each state, the bead moves in a specific driving po-tential, mimicking the power and recovery strokes of areal cilium [14, 19–25]. In another class models, cilia(known as rotors) are considered as spherical beads or-biting on rigid or flexible two dimensional trajectoriesunder a driving torque [13, 18, 26–28]. In more complexmodels, cilia are modelled as actively driven semiflexiblefilaments with more realistic beating pattern [10, 29, 30].In such models, the hydrodynamic coupling between ciliacan lead to metachronal waves in systems of cilia in one-and two-dimensional lattices. In this paper, to under-stand the role of spatial heterogeneity, we consider theframework of the rower model.In this work, we will present new results on how het-erogeneity of cilia position influences the stability androbustness of coherent beating of cilia. We focus on oneof the simplest models for cilia beating, the rower model a r X i v : . [ c ond - m a t . s o f t ] D ec [19] in one dimension. We carefully study the effect ofspatial heterogeneity in systems with few rowers as wellas in systems with a large number of rowers with variouskinds of heterogeneities. First, we study the 3-rowerscase in great detail as it is the minimal way to introduceposition irregularities. We find phase drifting and bifur-cation of the frequencies of beating when the distancebetween the second and third rower is varied, keeping thedistance between the first and second rowers fixed. Thephenomenon of phase drifting is also present in a systemwith more than 3 rowers, and we study different spatialconfigurations with a finite number of rowers. This al-lows us to identify a crossover distance, above which theseparation between consecutive clusters is large enoughto decouple their dynamical behavior. In the case of alarge number of coupled rowers, we consider several typesof spatial heterogeneity. We find that when the cilia arespatially distributed on randomly clustered configura-tions, their dynamical behaviour is characterized by arobust average common beating frequency, that dependsonly weakly on the density of cilia. This seems to corre-spond to preliminary results on human epithelial tissues[31]. On the other hand, the collective behaviour ob-served for other types of spatial heterogeneities is strik-ingly different. II. THE ROWER MODEL
The rower model for cilia, proposed by Cosentino et al. in 2003 [19], experimentally realized by driven colloidsin viscous fluids [22], remains one of the basic model tostudy hydrodynamic synchronization [14]. In this model,the complex structure of a cilium is coarse-grained as aspherical bead, and the periodic beating pattern is de-scribed as an oscillating linear motion of a bead in aviscous fluid. In order to create a sustained oscillat-ing motion of the bead, two different driving potentialsare used, and a mechanism for geometrical switchingbetween these potentials is employed. In Fig. 1 (left),we show a schematic diagram of the geometric switchand the harmonic potentials. The bead motion undereach driving potential corresponds to a specific state ofthe bead σ = ± ± y direc-tion). If the bead reaches a particular limiting ampli-tude y = ± s , it switches to the other driving potentialand consequently reverses the direction of its motion.Hence, it successfully creates a sustained oscillating mo-tion. These two states of the bead mimic the power andrecovery stroke of beating of cilia.We consider an array of N rowers on a one dimen-sional lattice of size L (see Fig. 1 (right)). The spacingbetween two consecutive lattice sites is d = 1. If thereis no heterogeneity, all the lattice sites are occupied andhence, L = N . In the case where rowers are placed het- erogeneously on the lattice, one allows for empty latticesites and L > N . Let us call the positions of the rowers x , x , .., x n , .., x N where x < x , .., < x n , .., < x N . Twodynamical variables σ n (state) and y n (displacement) de-scribe the motion of the n th rower.A bead switches between the two states if the dis-placement y reaches the maximum amplitude ± s , i.e.if y n ( t ) = ± s then σ n ( t ) = − σ n ( t ) and σ n ( t ) = σ n ( t )otherwise. We choose the two driving potentials to beharmonic V ( y n , σ n ) (Fig. 1 (left)). The external drivingforce f n for the n th rower is given by: f n = − ∂V ( y n , σ n ) ∂y n = − ( ky n − σ n ) . (1)In our study, the stiffness constant of the potential k isassumed to the same for all the rowers. This externaldriving force on the bead is a simple approximation ofthe complex internal active force of a real cilium.Cilia motion corresponds to a low Reynolds numberregime and we are interested by the far-field hydrody-namic regime: the size and displacements of the beadsare small compared to the lattice spacing. In this limit,the hydrodynamic coupling between the rowers is givenby Oseen tensor [32, 33]. The velocity at any instantof time v mn acting on the n th rower induced by activeoscillation of rower m is given by: v mn = O ( m, n ) f m . (2)The Oseen coupling between any two rowers m and n is O ( m, n ) = 1 / (8 πηd mn ), where d mn is the distance be-tween the sites ( d nm = | x n − x m | ), and η is the viscosityof the fluid medium. Here, we note that as our focusin this paper is on the coherence of cilia beating (noton the flow of surrounding fluid) and therefore the driv-ing forces for both states are chosen to be symmetric forsimplicity (as in [22, 24, 25]). However, one can makethe motion of the two states asymmetric by consideringdifferent driving forces [34] or different drag coefficientsfor both states [19].Hence, in the overdamped limit, the dynamical equa-tions for this system are given by,˙ y n = 16 πηa f n + N (cid:88) m (cid:54) = n v mn = 16 πηa f n + N (cid:88) m (cid:54) = n f m πηd mn . (3)Here, the factor 6 πηa is the viscous drag coefficient ofa spherical bead with radius a . For computation pur-poses, we choose 6 πηa = 1, hence 1 / (8 πη ) = 3 a/ ≡ α .Replacing a and η , Eq. (3) can be rewritten as,˙ y n = f n + α N (cid:88) m (cid:54) = n f m d mn . (4) y−s s V(y) σ = σ =− σ=1σ=−1 x ys −s σ=1σ=−1 FIG. 1. A schematic diagram of the model.
Left
The harmonic potentials correspond to the two different states ( σ = ± y = ± s . The driving potentials and the switchingmechanism form a simple description of real beating. Middle
A schematic representation of simplified beating pattern of acilium.
Right N rowers are placed on a regular one dimensional lattice. The positive y direction is σ = 1 and negative y direction is σ = −
1. We create heterogeneity in rowers’ position by keeping empty lattice sites in-between occupied sites.These are shown in Fig. 11.
Note that in a realistic physical situation, the couplingstrength α is always a positive quantity. In the absenceof any hydrodynamic coupling, (i.e. α = 0), it is easyto calculate the natural frequency of the rowers ω . Itdepends on the force constant k and the value of limitingdisplacement s ; its analytical expression is given by [19,22]: ω = 2 πT = πk ln (cid:104) ks − ks (cid:105) . (5)As mentioned before, this model has been studied forsystems of several rowers on regular lattices in one andtwo dimensions [19, 34]. The collective beating of therowers lead to metachronal waves as a result of anti-phase synchronization of neighbours for k >
0. When k < α < α , k and s already used in the original paper[19]. The force constant for the harmonic driving po-tential is k = 1 .
0. A single rower oscillates between − s to + s with s = 0 .
8. The value of α is taken to be 0 . N rowers, we use theEuler method with an integration step h between 10 − and 10 − , depending on the situation. We consider herethe deterministic case, hence simulations are run with-out thermal noise. The case of adding thermal noise willbe considered elsewhere [36]. In the following, we try to understand the role of spa-tial heterogeneity on the collective behavior of cilia beat-ing. First, we study a 3-rower system which is the min-imal sets to study heterogeneity effect. We consider 3rowers on a one dimensional lattice such that the two firstrowers occupy two consecutive lattice sites whereas theposition of the third rower increases (see Fig. 4). Thenwe study the effect of heterogeneity in systems composedof 4, 5, and 10 rowers, where we find that the resultingdynamics obeys a general scenario (section III). Finally,we study systems of a large number of rowers with threetypes of heterogeneities (section IV)corresponding to dif-ferent degrees of randomness — (i) the regular clusteredcase, (ii) random case, and (iii) the random clusteredcase. How to characterize collective behavior?
A collection of oscillators can display different emerg-ing features through coherence in phases and frequencies[37, 38]. Most common phenomena of phase coherenceare synchronization and phase-locking. In order to char-acterize the phase coherence, we define a phase variable φ n for n th rower using the following prescription (givenby Stark et al. [34]): φ n = 2 πm n + π σ n y n s , if y n > σ n = 1 , = 2 πm n + π σ n y n s + 2 π, if y n < σ n = 1 , = 2 πm n + π σ n y n s + π, if σ n = − . (6)The phase number m n ( ∈ N ) is increased by 1 afterrower n completes a full cycle. The piecewise linear re-lation between phase and displacement ensures that one FIG. 2. The periodic beating of a single rower is plottedas function of time t . The red curve is the displacement y of the rower, and the blue curve is its phase representationaccording to Eq. (6). full cycle in y is equivalent to a rotation of 2 π in φ . InFig. 2, we plot the time evolution of displacement y andcorresponding phase φ for a single rower. In the case ofperfect synchronization, all oscillators oscillate in phase;it is defined as φ j = φ j +1 for all j . For phase-lockingsystems, neighbouring oscillators maintain a constantnonzero phase difference δ , φ j = φ j +1 + δ . A nonzero δ leads to the formation of traveling wave in a systemof many oscillators. This wave is known as metachronalwave in the literature [6]. The degree of phase coherencecan be measured by a complex order parameter Z , whichis defined as [34], Z = A e i Φ = 1 N − N − (cid:88) n =1 e i ∆ φ n , (7)where ∆ φ n = φ n +1 − φ n and N is the total numberof oscillators. The system is maximally coherent when A = 1, and has no coherence for A = 0 (see AppendixI). When ∆ φ n = δ for all n , A = 1. This is true forany constant δ including zero i.e, for both phase lockedand in phase synchronization solutions. Therefore, forboth perfect metachronal wave and fully synchronizedstates the value of the order parameter A is 1. In orderto distinguish the synchronized states from metachronalwaves, one should also measure Φ, which gives the aver-age angle of phase difference. For a synchronized stateΦ = 0, and for a metachronal wave, Φ = δ (cid:54) = 0.Beside phase coherence, another observed emergingbehavior in a system of coupled oscillators is frequencylocking. In this case, the phases of oscillators can be dif-ferent but they oscillate with a common frequency whichis different from their natural frequency (the frequencyin the absence of coupling) [37, 38]. We define the aver-age frequency ω i of a rower i in the following way [38], ω i = lim t →∞ t (cid:90) t + tt ˙ φ i dt = lim t →∞ φ i ( t + t ) − φ i ( t ) t . (8) Here, t is sufficiently large so that the system hasreached a dynamical steady state at that time. In gen-eral, the frequencies of the oscillators are distributed overa distribution P ( ω ). For a perfectly frequency lockedsystem, the probability distribution of frequencies P ( ω )is a δ -function. III. RESULTS FOR FINITE-SIZE SYSTEMS
In this section, we systematically study arrays con-sisting of a few rowers (from N = 2) to many row-ers ( N = 100). We compute the order parameter A (Eq. (7)), and the distribution of beating frequencies ω i (Eq. (8)) to characterize the collective behavior of thebeating. Our objective is to understand the role of thespatial heterogeneity of cilia on the stability and robust-ness of the synchronization. This heterogeneity is intro-duced as soon as the number of rowers is 3.
1. 2-rowers case
The dynamical behavior of a system two rowers arewell understood. It is known that two rowers oscillate inopposite phase [19, 22, 34], and the collective frequencyof the oscillation depends on the separation and hydro-dynamic strength [22]. Here, we revisit the two-rowercase mainly to study collective frequency for differentseparating distances d . As we will see in the next sec-tions, this study will be useful to understand a systemwith more rowers.In the inset of Fig. 3(a), the phase difference∆ φ ( t ) = φ ( t ) − φ ( t ) is plotted as a function of time t for α = 0 .
1. At t = 0, we start the simulation withan arbitrary phase difference. We observe that within asmall transient time the 2 rowers reach to a perfect anti-phase synchronization state (∆ φ ( t ) = π ). Let τ tr bethe average time required to reach a perfect anti-phasesynchronized state. In our simulation, we compute τ tr by the time t at which ∆ φ ( t ) reaches the value 0 . π and then averaging the data over many random initialconfigurations( ∼ τ tr linearly in-creases with d when α is constant (or with 1 /α when d is constant) (see Fig. 3(a)). This is consistent withthe fact that the larger the interaction ( α/d ), the fasterthe rowers reach the synchronization state.What is the collective frequency of the anti-phase syn-chronization for two rowers? One can derive the expres-sion for the collective frequency ω coll, ( d ), by assumingthe anti-phase solution in the dynamical equation Eq. (4)[22]. The collective frequency ω coll, ( d ) is a functionof d and α , and is given by: ω coll, ( d ) = ω , = ω (1 − | α | /d ) , (9)where ω is the natural frequency of the rower when theinteraction is absent (see Eq. (5)). This matches exactlywith the numerical findings, as shown in Fig. 3(b). (a) τ t r d π ∆ φ ( t ) t (b) ω d ω (1+ α /d ) ω (1- α /d ) FIG. 3. (a) The transient time ( τ tr ) to reach the anti-synchronized state from an arbitrary initial condition is plot-ted against separation d . τ tr increases linearly with d . In-set: The phase difference between the rowers, ∆ φ , is plottedagainst time t . At t = 0, the rowers start from an arbitraryinitial condition, and after τ tr , ∆ φ converges to π , i.e. therowers beat in exact anti-phase synchronization. (b) Collec-tive frequency ω coll, ( d ) as a function of d for differentkinds of initial conditions. The data with solid squares areobtained when the two rowers start from a initial conditionwhich is different from an exact in-phase configuration whilethe data with empty squares are obtained with in exact in-phase initial condition. Their frequencies converge to ω forlarge d . A special case arises for a specific initial condition,which is specific to the deterministic (zero noise) system.At t = 0, if the rowers are in exact in-phase ( y = y and σ = σ ) there will not be any anti-phase synchroniza-tion, and the rowers will continue to beat in phase for-ever. The analytical expression for ω coll, of the beatingis given by: ω coll, ( d ) = ω , = ω (1 + | α | /d ) . (10)In Fig. 3(b), we plot the frequency as a function of d when the rowers start from an exact in-phase initial con-dition (empty squares): the data points exactly matchthe analytical expression given by Eq. (10). However,this solution is very unstable. If the simulation startswith an initial condition which is slightly different fromthose specific initial conditions, one reaches the usualsolution discussed above.
2. 3-rowers case
A 3-rowers system is the minimal setup to study thespatial heterogeneity of cilia. It shows an interesting phe-nomenon called phase drifting in which the third roweroscillates with a different phase than the first two rowers,typical of a dynamical bifurcation. This phenomenon is familiar in other coupled non-linear systems of oscilla-tors when the coupling strength or the noise strengthis varied (Adler systems [37, 39–41]), and has also beenreported for rotors near a wall as the distance from thewall is increased [42].We consider the 3-rowers system in one dimension,schematically shown in Fig. 4, where the first two rowersoccupy two consecutive lattice sites and the third roweris placed on a site which is a distance d apart from thesecond rower. The lattice constant is d = d = 1. Weconsider different cases for different values of d .
21 3 d FIG. 4. The minimal set up to study the spatial inhomo-geneity of cilia position. The first two rowers occupy the twoconsecutive lattice sites d = 1. The third is placed d distance apart from the second rower. For a given value of d , we solve the dynamical equa-tions for 3 rowers numerically. We find that, in the dy-namical steady state, the first two rowers always oscil-late in anti-phase with each other and the behavior ofthe third depends on the value of d . In Fig. 5(a), weplot the phase difference ∆ φ = φ − φ as a function oftime t for various d . We see that, for d < d c = 2, thephase difference ∆ φ is constant and for d ≥ d c , ∆ φ grows and is modulated in time. The latter implies thatthe phase locking between the second and third is lostdue to the appearance of phase drifting [37, 39, 40].As a result, for d < d c all 3 rowers oscillate with thesame frequency and for d ≥ d c , rower 1 and 2 oscillatetogether with the same frequency but rower 3 oscillateswith a different frequency. In Fig. 5(b), we plot thefrequency of the 3 rowers ω , ω , and ω as a function d . For d = 1, all the frequencies of all rowers are thesame and for d ≥ ω = ω (cid:54) = ω .Hence, for the set of parameters used in our simula-tion ( k = 1, s = 0 . α = 0 . d c is2. We simulated the 3-rowers system for several valuesof hydrodynamical coupling α = 0 . , .
05, and 0 .
01 withthe same value of the driving parameters k and s andfind d c = 2 in all cases. In this range of investigatedparameters, d c is independent of α . By studying severalcases with different sets of driving forces of the rowersfor a given α , we found that the value of d c depends onthe force constant k , and on the amplitude of the oscil-lation s i.e. d c = d c ( k, s ) (see Appendix II). This resultis very interesting. It means that hydrodynamic inter-action strength does not have any role on determining d c . The value of d c is therefore completely determinedby the internal activity of the cilia and not by the fluidviscosity. However, this result is weakened for N > φ - φ t (a) d = 1 2 3 10 (b) ω ω ω ω coll,2 + α /d ω - α /d FIG. 5. (a) The phase difference between the rower 2 and 3,∆ φ versus time t for different d . ∆ φ remains constantfor d < d c . At d c , a bifurcation occurs and ∆ φ grows intime for d ≥ d c . Here, the value of d c is 2. (b) Frequenciesof the three rowers as a function of d . Bifurcation occursat d = d c = 2. For very large d , rower 3 becomes al-most independent of the other rowers and oscillates with itsnatural frequency ω while the first 2 rowers oscillate withthe collective frequency of a 2-rowers system ω coll, . Thedata here is obtained from a single configuration. Indeed, wechecked that the steady state does not depend on the initialconfiguration if the simulations does not start with a peculiarcondition (such as: all rowers are in phase or all are in antiphase at t=0). τ d r i ft d (a) -5 -4 -3 -2 -1 α τ d r i ft - c on s t d (b) α=0.1α=0.05α=0.01 FIG. 6. (a) The log-log plot of ατ drift − const for α = 0 . , . .
01. (b) This scaling makes all the curves fall into asingle curve f ( d ). The function f ( d ) s well fitted by 1 /d .To compute the power laws more accurately, we choose verysmall step size for the Euler integration method, h = 10 − . For N = 4, we observe d c depends on α as well as initialconfiguration of rowers (see Appendix II).The values of ω , ω , and ω depends on d . Aswe have noted above, ω = ω for all d , but ω be-haves differently for d ≥ d c . For very large d , thehydrodynamic coupling of the 3rd rower due to the firsttwo can be neglected, hence the frequencies of the row-ers become independent of d . In this limit, the 3rdrower beats with the natural frequency ω and the firsttwo rowers oscillate with the collective frequency of tworowers ω coll, (see section III.1). On the contrary, forsmall d , the dependence of the rower frequencies on d is non-trivial. We use the following fits: ω ( d ) = ω ( d ) (cid:39) ω coll, + α/d , and ω ( d ) (cid:39) ω − α/d .In Fig. 5(b), we plot these functions (solid for ω , anddotted line for ω ). The two lines match nicely with nu-merical data (points). Note that ω saturates rapidly to ω ( ω − ω (cid:39) α/d ), while the decay of ω to ω coll, isslow ( ω − ω coll, (cid:39) α/d ). The expression of ω ( d )and ω ( d ) can be understood by the following way.Consider that the 3 rowers oscillate with the same fre-quency, the first two being exact anti-phase, and thirdone being exact anti-phase with the 2nd one. Althoughthis is obviously a very crude assumption, one can useEq. 4 and solve it for the 1st rower: ω = ω (1 − α + α/ ( d + 1)) (cid:39) ω coll, + ω α/d . (11)Solving the equation for the 3rd rower leads to: ω = ω (1 − α/d + α/ ( d + 1)) (cid:39) ω − ω α/d . (12)The growth of ∆ φ is not uniform and follows a peri-odic pattern in time (Fig. 5(a)). For d = 2, the phasedifference rather shows sharp 2 π jumps followed by al-most flat regimes. With the increase of d , the growthrate of the phase difference becomes faster and steps dis-appear. In order to characterize the nature of the phasedrift, we define the time scale τ drift ( α, d ) as the timeneeded for the phase difference ∆ φ to increase by 2 π (i.e. the time period of ∆ φ ), and we compute it fora given internal activity of rowers (i.e for fixed k and s ). We find that α τ drift ( α, d ) = f ( d ) + const , wherethe function f does not depend on α . In Fig. 6, weplot ατ drift − const in log-log scale and we show thatthe curves for different α collapse into a single curve.The scaling function f ( d ) decays algebraically to zeroas 1 /d . This collapse and the 1 /d decay can be ra-tionalized using the following simple arguments and ap-proximations. From the definition of the drift time, it iseasy to see that τ drift ( d ) ∼ / ( ω ( d ) − ω ( d )) (cid:39) / ( α d ) + const/α , using the fitting expressions for ω ( d ) and ω ( d ).
3. 4-rowers case
The phenomena of phase drifting and bifurcation isalso observed for a system of 4 or more rowers. Here, westudy the case of 4 rowers. We divide the 4 rowers intotwo subgroups in order to study different configurations.For a given internal structure (i.e. for given k and s ), oneobtains a variety of behavior, as the value of d c dependson the number of rowers in each groupCase I — Asymmetric case: The first group of rowershas 3 consecutive rowers with equal gap, d = d = 1.The second group consists only of the 4th rower and thedistance separating the two subgroups is d . In Fig. 7(b), we plot the frequencies for rowers as a function of d . A bifurcation similar to the one obtained for threerowers is observed for d ≥ d c . The critical distance forthis case is different from the 3-rowers system. Unlike for N = 3, the value of d c depends α and initial configura-tion of rowers (see Appendix II). For very large d , the3-rowers group behaves independently of the 4th rowerwhich oscillates at ω , the natural single rower frequency.Case II — Symmetric case: The first group has 2 con-secutive rowers with gap, d = 1. The second groupconsists of the 3rd and 4th rower with gap d =1. Here,all rowers oscillate with the same frequency and as thedistance between two groups d is very large the fre-quency asymptotically saturates to the frequency of 2-rowers system. Here the two clusters independentlybeat with the same frequency characteristic of their size(Fig. 7 (a)). (a) d ω ω ω ω ω coll,2 d (b) d ω ω ω ω ω coll,3 ω FIG. 7. Angular frequencies of 4 rowers for different gapsbetween clusters; α = 0 .
4. 10-rowers case
We divide the 10 rowers set into two subgroups andstudied the 5 different following cases.Symmetric case (5-5): each subgroup consists of 5 row-ers equally spaced, but the distance d between the twogroups can vary. We observe no bifurcation in this sym-metric case (see Fig. 8 (a)).There are 4 asymmetric cases: 6-4, 7-3, 8-2 and 9-1groups. In all these cases, a bifurcation occurs. In Fig. 8(b), (c), (d), (e), results for the collective frequencies ofthe 2 clusters are shown. The critical distances afterwhich bifurcation occurs are different in each case.From the above study of finite size systems, we learntwo instructive features. (i) As expected from symme-try, there is no bifurcation if the sizes of the 2 clustersare equal, instead all rowers beat at the same frequencywhatever the distance between the 2 clusters. (ii) Thedistance where the bifurcation occurs is maximal for non-trivial cluster sizes and is minimal for strikingly differentsizes. (a) ω ω ω coll,5 (b) ω ω ω coll,6 ω coll,4 (c) ω ω ω coll,7 ω coll,3 (d) ω ω ω coll,8 ω coll,2 (e) d ω ω ω ω coll,9 FIG. 8. From (a) to (e), the sizes of the 2 clusters arerespectively 5-5/ 6-4/ 7-3/ 8-2/ 9-1. When the size of thetwo clusters is comparable (5-5 and 6-4) the solution of thedynamical equations depends on initial conditions. The dataare then averaged over 1000 initial conditions.
5. Link between the order parameter A and therowers phases φ As seen before, local observables such as phase dif-ferences between rowers ∆ φ ij or individual frequencies ω i are necessary to describe in detail the good or poorcoherence state of an assembly of rowers. If one nowlooks at the global observable A , how will A be relatedto the dynamical state of the system? The relation be-tween ∆ φ ij between consecutive rowers and A is givenby Eq. (7). However, by looking only at the value of A , it is not possible to infer the dynamical state of thesystem (i.e ∆ φ ij ). The study of the 3+1-rowers case isinstructive in this respect and is reported here. As seenin Fig. 7 (b), d c = 3 in this case. We show in Fig. 9 that3 distinct dynamical phases can be identified by compar-ing the time averaged value of A with the rowers phase (a) < A > d -4-2 0 2 4 15 30 45 60 75 (b) t ∆φ ∆φ ∆φ -4-2 0 2 4 15 30 45 60 75 (c) t ∆φ ∆φ ∆φ -4-2 0 2 4 15 30 45 60 75 (d) t ∆φ ∆φ ∆φ FIG. 9. (a) The time averaged order parameter (cid:104) A (cid:105) , av-eraged over a large time window after the system reachesdynamical steady state for an arbitrary initial condition forthe 3+1-rowers case, is plotted against separation d . Thetime evolution of phase differences between consecutive row-ers for d = 1, d = 2, and d = 50 are plotted in (b), (c),and (d) respectively. differences. If d = 1, the group of rowers is compactand oscillate in almost anti phase, hence (cid:104) A (cid:105) is maximaland close to 1. If d = 2, the system stands right beforethe bifurcation ( d c = 3), and the rowers state is dynam-ically disordered as one can see on Fig. 9 (c). If d ismuch larger than d c (and this is illustrated by lookingat d = 50, the two clusters of rowers oscillate indepen-dently but the 3 rowers in the first group are coordinatedin almost anti phase, leading to A reaching a saturationvalue after the dip at d = 2.
6. Collective beating frequency of rowers on aregular array
For a regular array of N equally spaced rowers, we ob-serve that all rowers oscillate at same frequency (excepta few rowers close to the boundaries for very large N ,where N is the number of total rowers) showing phasecoherence both for α < α > ω coll dependson N . We observe that for α < ω coll strongly dependson N (see Appendix IV). For α >
0, however ω coll stabi-lizes to a constant value. This confirms that metachronalwaves in mucociliary systems can be stable even in large samples. ω c o ll Nsimulationtheory
FIG. 10. The collective frequency is plotted as a function oftotal number of rowers N . The rowers are placed on a regularlattice with lattice constant d = 1. The data are averagedover 1000 initial configurations. It reaches a stable frequencyfor N > N c ( ∼ The value of the collective beating frequency ω coll ofthe rowers oscillates for small N (total number of row-ers), and it eventually converges to a constant value forlarge N . We can understand this behavior from a verysimple picture. Consider the following dynamical steadystate for our oscillating rowers in which a rower beats inperfect anti-phase with its nearest neighbours. This im-plies y ( t ) = − y ( t ) = y ( t ) = − y ( t ) = .., = y N − ( t ) = − y N ( t ). In this particular case, we can easily solve thedynamics (Eq. (4)) and calculate ω coll . For N numbersof total rowers, the collective frequency is given by, ω coll,N = ω [1 − α (1 −
12 + 13 + .. + ( − N N − . (13)The coefficient of α is the alternate harmonic (conver-gent) series. In Fig. 10, we plot both the collective fre-quency computed from theory (Eq. (13)), and the oneobtained from simulation. For N = 2, the values ofcollective frequency computed from these two differentmethods match exactly; the assumption of perfect anti-phase synchronization is indeed true for N = 2 (seenpreviously, Fig. 3(a)). However, for N >
2, this is notexactly the case because the assumption of perfect anti-phase synchronization is no more valid. We believe thatthe influence of boundaries is indeed important even forlarge N , as seen on Fig. 10.We have simulated systems with small number of row-ers for another value of α (= 0 .
2) and have observed thatthe qualitative behavior of our results is independent ofcoupling strength.
IV. COLLECTIVE BEHAVIOR OF MANYSPATIALLY HETEROGENEOUS ROWERS
In the previous section, we have studied how spatialheterogeneity of rowers affect the coherent beating whenthe number of rowers is small. We observed that thespatial arrangement of rowers can lead to phase drifting,phase incoherence, and bifurcation in frequency whenthe separation between two clusters is greater than acritical distance. We now turn to asking the questionof how spatial heterogeneities in the position of rowersaffect the coherent beating in a system of hundreds ofrowers.On the one dimensional lattice, we will consider threetypes of heterogeneities — (i) regular clustered config-urations, (ii) random configurations, and (iii) randomclustered configurations, which is an intermediate situa-tion between (i) and (ii).We generate these types of heterogeneities for differentvalues of the density of rowers ρ , and study the dynam-ical properties of the corresponding systems.(i) In the case of regular clustered heterogeneity, clus-ters of fixed number of rowers are placed on a lattice,separated by regular gaps. A typical lattice configura-tion of rowers is shown in Fig. 11(i). We call cn thenumber of rowers in a cluster, and gl the gap length be-tween two consecutive clusters. We refer to such clustersas “cn-gl”. One way to generate configurations of differ-ent densities for a given cluster length cn is to vary thegap length gl between two clusters. The density is givenby, ρ = cncn + gl − . (14)(ii) For random heterogeneity (see Fig. 11(ii)), N row-ers are placed randomly on a one dimensional lattice ofsize L . The density of rowers is given by ρ = N/L . For afixed N , we generate configurations with different valuesof ρ .(iii) In the case of random clustered heterogeneity,we introduce randomness in the sizes of the clusters aswell as in the gap lengths between two clusters. Here,cluster lengths cn are chosen randomly between cn min and cn max using uniform random distribution. The gaplengths gl are chosen randomly between gl min and gl max using uniform random distribution. We refer to this typeof heterogeneity as [ cn min , cn max ] − [ gl min , gl max ]. Atypical lattice configuration is shown in Fig. 11(iii). Inorder to generate configurations with different densitieskeeping fixed cn min , cn max , and gl min , we vary gl max .The average density can be computed numerically, aver-aged over a large number of realizations: ρ = (cid:104) cncn + gl − (cid:105) . In our computation, we will use 2- gl and 3- gl regularclustered heterogeneities. For random clustered configu- FIG. 11. Different types of spatially heterogeneous row-ers configurations (rowers are the arrows). (i) Regular clus-tered: clusters of fixed number of rowers are placed on thelattice regularly spaced. In this diagram, the number of row-ers in a cluster is cn = 2, and the length of the gap be-tween two clusters of rowers gl = 3. We refer to this latticeas “2-3”. (ii) Random: rowers are randomly placed on thelattice. (iii) Random clustered: the number of rowers in acluster is not fixed, and is chosen uniformly in the interval[ cn min , cn max ]. The gap lengths between two consecutiveclusters are also randomly taken from [ gl min , gl max ]. In thisdiagram, the minimum and maximum length of clusters ofrowers are cn min = 2 and cn max = 4. The minimum andmaximum length of a gap between two clusters are gl min = 2and gl max = 5. This structure is referred to as “[2,4]-[2,5]”. rations, we will use the following heterogeneities: [2,4]-[2, gl max ] and [3,5]-[2, gl max ].We measure the distribution of beating frequency P ( ω ) as a function of ρ in the different types of hetero-geneities discussed above. We use N = 100 and α = 0 . P ( ω ), we consider several spatialstructures and initial conditions of rowers. In the caseof random and regular clustered, we take 100 differentspatial configurations of rowers. For regular clusteredheterogeneity, as the spatial configuration of rowers ona lattice is fixed, we consider 100 different initial config-urations of the rowers.
1. Case ( i ) : Regular clustered We consider “2-gl” regular clustered configurations. InFig. 12, we plot the spatial profile of the average beat-ing frequency. We observe that for a given ρ all rowersoscillate with a single collective frequency. The value ofthis collective frequency decreases as ρ is decreased (byincreasing gl ). In the large gl limit, the inter clusters in-teraction can be neglected. As a consequence, we see inFig. 12 that the frequency of the rowers for gl = 20 is thesame as the collective frequency of a cluster of two row-ers separated by a distance d = 1, ω coll, (section IV.1).0 ω rowers’s no2-202-52-22-1 ω coll,2 FIG. 12. (i) Regular clustered case. Spatial profiles of therowers frequency for “2-gl” regular clustered heterogeneityfor different gap lengths gl = 1 , , gl values are given by ρ = 1, 0.67, 0.33,and 0.1 respectively. For a given gap length, all rowers in alattice oscillate approximately in a single collective frequency.The value of the collective frequencies depend on the valueof gl . For gl = 20, the inter-clusters hydrodynamic couplingcan be neglected, and in this limit, rowers oscillate with thecollective frequency of a cluster of two rowers, ω coll, . However, for “3-gl” regular clustered configurations, ω increases as ρ is decreased, and reaches ω coll, at small ρ (see Fig. 15(a)).
2. Case ( ii ) : Random ω rowers’ no (a) P ( ω ) ω (b) ρ =0.1 ρ =0.4 ρ =0.75 ρ =1.0 FIG. 13. (ii) Random case. Frequency plots for randomlypositioned rowers with density ρ = 0 . , . , .
75, and for aregular lattice with ρ = 1. (a) The frequency of the rowers, ω is plotted as a function of rower’s number. For ρ = 1,all the rowers oscillates in a single collective frequency. For ρ <
1, we see the whole system breaks into different parts,and these different local groups (synchronized island) oscil-late with different frequencies. The number of such groupsincreases with decreasing ρ . (b) The probability distribu-tions P ( ω ) are plotted for the same values of ρ (as in (a)).For ρ = 1, it is a δ -function (data is scaled by a arbitrarynumber for best visualization). As ρ is decreased, the widthof the distribution increases, while the average ω increases.For ρ = 0 . , . In Fig. 13 (a), we have plotted the steady state fre-quency profile as a function of rowers’ number for dif-ferent densities in the random case. For a given ρ , weconsider a configuration of randomly placed rowers on alattice and evolve this system with a random initial con-dition of rowers. For ρ = 1 (perfect regular array), all therowers oscillates with the same frequency. This resultsin a δ -function for P ( ω ) (Fig. 13 (b)). For ρ <
1, variousclusters of rowers start to oscillate at different frequen-cies, leading to a finite width in P ( ω ) plot (Fig. 13 (b).Let us call synchronized island a group of consecutiverowers (connected by next occupied lattice sites) beatingwith the same frequency. We find that a synchronized is-land consists of a few clusters ( ∼ ρ (except for very small density where mostrowers beat at their natural frequency) and so does thewidth of the distribution P ( ω ). From our study on fi-nite numbers of rowers, we know that if the separationbetween two groups is less than a critical distance theyoscillate with the same frequencies. The value of the fre-quencies depends on the separation between two groupsand number of rowers in each group. This leads to anon-trivial distribution of the beating frequencies.For small densities, the frequency spectrum will bedominated by small clusters characteristic frequencies.Indeed, for ρ = 0 .
1, the position of the distant peak at ω (cid:39) ω is due to isolated single rowers, while the otherpeaks at small ω are due to a collective effect of therowers. For larger densities, the frequency spectrum isgetting even more complex and broad.
3. Case ( iii ) : Random clustered We consider [2,4]-[2, gl max ] random clustered hetero-geneous configurations. The frequency profile and P ( ω )are plotted for gl max =2, 4, 6, 8 and 24 in Fig. 14. As inthe case of random heterogeneity, we observe that vari-ous clusters of rowers oscillate with different frequenciesand the number of such clusters increases with gl max .This leads to finite widths in P ( ω ) (see Fig. 14 (b)). Aremarkable feature is that the mean of the distributionis almost independent of the density ρ ( gl max ). For [2,4]-[2, gl max ] configurations such that gl max >>
1, the av-erage gap between two consecutive clusters is large andthe effect of hydrodynamic couplings between them willbe negligible in many cases. Indeed the peaks in P ( ω )correspond to the collective frequencies of those clusters.The study with [3,5]-[2, gl max ] configurations also leadsto the same qualitative results (Fig. 15). We also findthat a synchronized island consists of several numbersof clusters. This number increases as a function of thedensity of rowers, and for large densities, this numberis very high compared to the random heterogeneity case1(see Appendix III). ω rowers’ no (a) P ( ω ) ω (b) [2,4]-[2,24][2,4]-[2,8][2,4]-[2,4][2,4]-[2,2] FIG. 14. (iii) Random clustered case. Frequency plots for“[2,4]-[2, gl max ]” lattices for gl max = 2 , , , and 24. The den-sities corresponding to these gl max are given by ρ = 0.75, 0.6,0.43, and 0.2 respectively. (a) The frequency of the rowers, ω is plotted as a function of rower’s number for a given spatialstructure of rowers’ position. The whole system breaks intosynchronized islands, which oscillate with different frequen-cies. The number of such groups increases with increasing gl max (i.e. with decreasing ρ ). (b) The probability distribu-tions P ( ω ) are plotted for the same values of gl max as in (a).Remarkably, the mean and the width do not strongly dependon gl max .
4. Comparing different structures with same density
In the last sections, we have considered each type ofheterogeneity separately and have investigated the na-ture of P ( ω ) for various densities. As experimental sam-ples can present a variety of spatial heterogeneities, wenow compare the results for different heterogeneities fora given density. Moreover, many studies only report ex-perimental values of ρ and not the precise spatial distri-bution of cilia. Hence comparing theoretically the dy-namical behavior of different configurations of cilia withfixed density will provide different types of scenarios thatwe hope we can compare to experimental observationsand check whether this is consistent with the observedspatial arrangement of cilia.
1. Frequency spectrum for different structures
We compare the mean frequency (cid:104) ω (cid:105) and the coeffi-cient of variation C v = (cid:112) (cid:104) ω (cid:105) − (cid:104) ω (cid:105) / (cid:104) ω (cid:105) of the distri-butions P ( ω ). In Fig. 15(a), we plot (cid:104) ω (cid:105) as a function ofdensity for different spatial heterogeneities of rowers. Forrandom heterogeneity, (cid:104) ω (cid:105) decays monotonically withdensity ρ . In the limit of ρ →
1, all structures look likeregular lattices and consequently all reach the collectivefrequency for many rowers. In the other limit ρ → (cid:104) ω (cid:105) converge to specific values depending on structures. Inthe case of random structures, (cid:104) ω ( ρ → (cid:105) = ω . In thecase of regular clusters, (cid:104) ω ( ρ → (cid:105) = ω coll, for “2-x” and (cid:104) ω ( ρ → (cid:105) = ω coll, for “3-x”. In the case of ran-dom clusters, (cid:104) ω ( ρ → (cid:105) lies between ω coll, and ω coll, for “[2,4]-[2,x]” and between ω coll, and ω coll, for “[3,5]-[2,x]”. Note that (cid:104) ω (cid:105) is almost independent of density inthe random clusters case, in contrast with the results forother heterogeneity types. The standard deviation C v ,Fig. 15(b), shows that the frequency distribution is thebroader for random heterogeneity.The average frequency ω appears to be a robust quan-tity for randomly clustered configurations of cilia, inde-pendently of the surface coverage ρ . This most probablycorresponds to the type of surface coverage observed inin vivo samples. Consequently, this robustness impliesthat even at low densities ( ρ (cid:39) . (a) < ω > ρ random2-gl3-gl[2,4]-[2,gl max ][3,5]-[2,gl max ] C v ρ (b) FIG. 15. Mean frequency and its fluctuation as a function ofdensity ρ for α = 0 . N = 100. For random and randomclustered lattice, the data are averaged over 100 spatial con-figurations. For regular clustered lattice, the data are aver-aged over 100 initial configuration of rowers. For these plotswe consider “2-x” and “3-x” regular clustered lattice, and“[2,4]-[2,x]’ and “[3,5]-[2,x]” random clustered lattice. (a) (cid:104) ω (cid:105) against ρ (with standard errors, errorbars are often smallerthan point sizes), and (b) normalized variance C v versus ρ .
2. Order parameter for different structures
The results in the previous sections show that hetero-geneity in the rowers positions lead to partial loss of fre-quency locking. Hence, all rowers do not oscillate with acommon global frequency. But, some rowers do oscillatewith a common frequency locally. This common localfrequency is different in the different parts of the lat-tice and depends on the local environment of the rowersin a non-trivial manner. The local frequency may alsodepend on the initial condition of the rowers.Hence, it is complementary to investigate whetherthere is any phase coherence among the rowers even inthis partially frequency locked system. This might give2an idea about the local phase coherence of the rowersthat oscillate with the same frequency. In order to esti-mate phase coherence, we measure the order parameter A (defined in Eq. (7)) and compute its average value, (cid:104) A (cid:105) for different values of ρ and different lattice struc-tures, in the same spirit as what we did in the 3+1-rowerscase in Section II.5. (cid:104) A (cid:105) is computed after steady statetime average and ensemble average (average over severalspatial configurations and/or several initial conditions ofrowers).For ρ = 1 the order parameter A (cid:39) .
9, which is closebut not equal to 1. The definition of A (Eq. (7)) showsthat its value depends on the distribution of phase differ-ence of neighbours ∆ φ j = φ j +1 − φ j . The phase lockedangle is not perfectly π , rather we see a finite width ofthe distributions of ∆ φ j which is (almost) independentof system sizes. This leads a l value slightly less than 1(see Appendix I).Not surprisingly, the regular clustered configurationsconserve the phase coherence, except at very low den-sities. However, it is remarkable that the random clus-tered case is more coherent at all densities than the purerandom case. In this case, which is likely to be the ex-perimentally relevant case, spatially separated clustersare internally coherent, and all rowers beat with almostthe same average frequency whatever its position andthe density of the sample. < A > ρ randomregular clusteredrandom clustered FIG. 16. Order parameter < A > as a function of density ρ for different structures for N = 100. For random and randomclustered lattice, the data are averaged over 100 spatial con-figurations. For regular clustered lattice, the data are aver-aged over 100 initial configuration of rowers. For these plotswe consider “2-x” regular clustered lattice and “[2,4]-[2,x]”random clustered lattice. V. CONCLUSION
While ciliated epithelium are most often modeled aslarge homogeneous carpets, they are observed experi- mentally to be inhomogenous and often rather sparse.However, cilia on such surfaces show coordination to alarge extent, and thanks to this coordination fulfill theirbiological role. We have addressed this issue here, in asimplified model of a 1D heterogeneous array of rowers.Hydrodynamic coupling is long-range and leads to col-lective coordination of the rowers phases in the homoge-neous array previously studied. When spatial hetero-geneity is introduced, gaps where rowers are absent leadto unperfected or lack of coordination. This is true inthe simplest case of a very small number of rowers: if thegap length between two small clusters of rowers exceedsa critical distance, phase drifting is observed and ulti-mately leads to decoupling of rower clusters which are in-ternally coordinated. Interestingly this critical distancedoes not depend on the prefactor in the hydrodynamicinteraction, which is proportional to the viscosity; it onlydepends on the internal parameters of the rower motility(but is likely to depend on the range of the interactions).Looking at the coordination of two clusters of rowersas the spatial gap between them increases paves the wayto understanding more complex arrays of rowers. Wefind that when the two clusters contain the same num-ber of rowers, all the rowers beat at the same frequencywhatever the gap. If they do not contain the same num-ber of rowers, each cluster bifurcates to its own intrinsiccharacteristic frequency after the gap length exceeds acertain value. The distance after which the clusters aretotally decorrelated (their frequency reach the one theywould have in the absence of the other cluster) varies ina non trivial way with the size of the patches.Hence the study of more complex arrays of rowers con-taining a large number of clusters of various sizes revealsa distribution of frequencies based on the physics de-scribed above. Not surprisingly, a regular configurationof clusters of similar size will show a delta distributionof frequencies. In contrast, a totally random configu-ration will produce a wider distribution of frequenciesif the density of rowers is small, while the average fre-quency will decrease. In between those two cases, a ran-domly clustered configuration (consistent with experi-mental observations) has a distribution of frequenciesthat depend barely on the density of rowers, while theaverage beating frequency is close to a constant when thedensity is varied. Consequently a realistic carpet of ciliawith density typically ρ = 0 . A which quantifiesthe coordination of phase differences between neighbor-ing rowers (existence of a metachronal wave) are com-plementary in providing a description of the dynamicalstates of the rowers.While an experimental situation is obviously subjectto thermal noise, our study and many others [19, 25, 34]are deterministic. In this approximations, the model3presents a few complications as the dynamical states maymarginally depend on initial conditions, as discussed inthe text. This complication is removed as one would in-clude noise, and we believe the results stay qualitativelythe same [24]. Besides the thermal noise, there is an ac-tive source of noise which is intrinsic to the system. Inthe case of the flagellar beating for Chlamydomonas andsperm cells, it has been demonstrated that the later hasa dominating contribution in their synchronized beating[43, 44]. It would be interesting to study the role of spa-tial heterogeneity in the presence of such intrinsic noisein future.Moreover, other studies have considered a variationof hydrodynamical coupling, by taking into account thepresence of a wall on which the rowers are attached (how-ever without including noise [34, 42] ). This variation in-troduces remarkable changes in the observed dynamicalstates of the cilia, like the coexistence of phase-lockedand desynchronized clusters, known as chimera states.This is likely to happen in homogeneous arrays as stud-ied in the above cited papers, and is likely to happentoo in the heterogeneous case, though it has not beenstudied so far.Along the same line, it will be interesting to study2D carpets of rowers, as well as a realistic experimentalspatially-resolved sample. In this case, we believe thatthe orientation of beating of the rowers will also be re-sponsible for the spatial correlations of the dynamicalstates in 2D [36]. Acknowledgments - This work has been carried outthanks to the support of the LabEx NUMEV project(n o ANR-10-LABX-20) and of the MUCOCIL project(n o ANR-13-BSV5-0015) both funded by the Investisse-ments d’Avenir French Government program, managedby the French National Research Agency (ANR). Wethank our collaborators (D. Donnarumma, M. Jory, A.Bourdin, I. Vachier, A. Fort-Petit, D. Gras, P. Chanez,A. Viallat, K. Khelloufi) for fruitful exchange on theirknowledge of bronchial epithelium. We also wish tothank P. Cicuta, L. M. Cosentino, and B. Friedrich foruseful discussions on the model and results. ∗ [email protected] † [email protected] ‡ [email protected] § Present address of SD:
School of Physics and Astronomy,Rochester Institute of Technology, Rochester, New York14623, USA[1] J. Gray,
Ciliary Movements (Cambridge Universiy Press,1928).[2] I. Iba˜nez Tallon, N. Heintz, and H. Omran, HumanMolecular Genetics , R27 (2003).[3] K. Sawamoto et al. , Science , 629 (2006).[4] D. Smith, E. Gaffney, and J. Blake, Respiratory physi- ology and Nurobiology (2008).[5] C. Brennen and H. Winet, Annual Review of Fluid Me-chanics , 339 (1977).[6] M. Sanderson and M. Sleigh, J. Cell Sci. , 331 (1981).[7] D. R. Brumley, M. Polin, T. J. Pedley, and R. E. Gold-stein, Phys. Rev. Lett. , 268102 (2012).[8] S. Gueron and K. Levit-Gurevich, Proc. Natl. Acad. Sci.USA , 12240 (1999).[9] N. Osterman and A. Vilfan, Proc. Natl. Acad. Sci. USA , 15727 (2011).[10] J. Elgeti and G. Gompper, Proc. Natl. Acad. Sci. USA , 4470 (2013).[11] S. Michelin and E. Lauga, Physics of Fluids , 111901(2010).[12] K. Y. Wan and R. E. Goldstein, Proc. Natl. Acad. Sci.USA , E2784 (2016).[13] R. Golestanian, J. M. Yeomans, and N. Uchida, SoftMatter , 3074 (2011).[14] N. Bruot and P. Cicuta, Annual Review of CondensedMatter Physics , 323 (2016).[15] M. G. Gabridge, C. C. Agee, and A. M. Cameron, TheJournal Of Infectious Diseases , 9 (1997).[16] H. R, B. A. aun Knauf S, K. FJ, and B. M, J MedPrimatol (2014).[17] J. Blake, Journal of Biomechanics , 179 (1975).[18] R. Golestanian, J. M. Yeomans, and N. Uchida, SoftMatter , 3074 (2011).[19] M. C. Lagomarsino, P. Jona, and B. Bassetti, PhysicalReview E , 021908 (2003).[20] M. Leoni, J. Kotar, B. Bassetti, P. Cicuta, and M. C.Lagomarsino, Soft Matter , 472 (2009).[21] M. Leoni, B. Bassetti, J. Kotar, P. Cicuta, and M. C.Lagomarsino, Physical Review E , 036304 (2010).[22] J. Kotar, M. Leoni, B. Bassetti, M. C. Lagomarsino,and P. Cicuta, Proceedings of the National Academy ofSciences , 7669 (2010).[23] N. Bruot, L. Damet, J. Kotar, P. Cicuta, and M. C.Lagomarsino, Physical review letters , 094101(2011).[24] L. Damet, G. M. Cicuta, J. Kotar, M. C. Lagomarsino,and P. Cicuta, Soft Matter , 8672 (2012).[25] G. M. Cicuta, E. Onofri, M. C. Lagomarsino, and P. Ci-cuta, Physical Review E , 016203 (2012).[26] T. Niedermayer, B. Eckhardt, and P. Lenz, Chaos:An Interdisciplinary Journal of Nonlinear Science ,037128 (2008).[27] N. Uchida and R. Golestanian, Phys. Rev. Lett. ,178103 (2010).[28] N. Uchida and R. Golestanian, Phys. Rev. Lett. ,058104 (2011).[29] S. Gueron, K. Levit-Gurevich, N. Liron, and J. Blum,Proc. Natl. Acad. Sci. USA , 6001 (1997).[30] B. Guirao and J. F. Joanny, Biophys. J. , 1900 (2007).[31] M. K. KHELLOUFI, Physique de la dynamique mu-cociliaire / Dispositif d´etude de la migration cellulaire3D Application `a lasthme et `a la BPCO , PhD thesis,Aix Marseille Universit´e, 2015.[32] J. Happel and H. Brenner,
Low Reynolds number hydro-dynamics (Martinus Nijhoff, Boston, 1983).[33] J. Dhont,
An Introduction to Dynamics of Colloids (El-sevier, 1996).[34] C. Wollin and H. Stark, Eur. Phys. J. E , 42 (2011).[35] Y. Roy, V. Sivathanu, and S. K. Das, Comp. in Bio. and Med. , 1758 (2013).[36] S. Dey, G. Massiera, and E. Pitard, in preparation .[37] A. Pikovsky, M. Rosenblum, and J. Kurths, Syn-chronization: a universal concept in nonlinear sciences (Cambridge university press, 2003).[38] J. A. Acebr´on, L. L. Bonilla, C. J. P´erez Vicente, F. Ri-tort, and R. Spigler, Rev. Mod. Phys. , 137 (2005).[39] R. Adler, Proc. IRE (1946).[40] R. Stratonovich, Topics in the Theory of Random Noise (Gordon and Breach, New York, 1963).[41] B. M. Friedrich, The European Physical Journal SpecialTopics , 2353–2368 (2016).[42] D. R. Brumley et al. , Phys. Rev. Fluids , 081201 (2016).[43] R. E. Goldstein, M. Polin, and I. Tuval, Phys. Rev. Lett. , 168103 (2009).[44] R. Ma, G. S. Klindt, I. .H. Riedel-Kruse, F. J¨ulicher, andB. M. Friedrich, Phys. Rev. Lett. , 048101 (2014). Appendix I. The order parameter A.
As we have discussed in the main text, the degree ofsynchronization in the system of phase oscillators can bemeasured by the following complex order parameter, Z = A e iΦ = 1 / (N − N − (cid:88) j=1 e i∆ φ j , (A1)where ∆ φ j = φ j +1 − φ j are phase difference of nearestneighbours and N is the number of the oscillators. Themagnitude A describes the phase-coherence, and polarangle Φ indicates average phase-difference of neighbours.The system show a maximal coherence when A = 1, andno coherence for A = 0. In the large N limit, we canwrite the above summation by integration, A e iΦ = (cid:90) π d(∆ φ ) P(∆ φ ) e i∆ φ , (A2)where P (∆ φ ) is the probability distribution of phase dif-ference of nearest neighbours.Let us discuss three cases in detail for different types ofphase coherence. Perfect phase locking :
In this case,the neighbouring pairs are phase locked to a constantangle δ i.e, ∆ φ j = δ . Therefore, P (∆ φ ) = p ( δ − ∆ φ ).From Eq. (A2) we get A = 1, and we also see that A isindependent of δ . No synchronization :
The proba-bility distribution of ∆ φ is totally random. In this case, P (∆ φ ) = 1 / (2 π ) which leads A = 0. Partial phaselocking :
The width (standard deviation) of the distri-bution of ∆ φ is neither zero (as in delta function) nor itpossess a maximal width (as in uniform distribution). Itrather has a finite width in between two extreme cases.Let us assume such a distribution by a normal distribu-tion with standard deviation σ , P (∆ φ ) = 1 √ πσ e − φ / σ . (A3) (k, s) α d c for N = 3 d c for N=4(1.0, 0.8) 0.01 2 (3-5) 3.51 ± ± ± ± ± ± ± ± ± ± ± ± ± ± d c for five different setsof internal parameters of cilia ( k, s ) for N = 3 (column 3)and 4 (asymmetric case, column 4). For each set of ( k, s ),we consider three different values of hydrodynamic couplingstrength α = 0 .
1, 0 .
05, and 0 .
01. For the computation ofeach d c , 1000 initial configurations are used. The integrationstep h = 10 − . As the value of the lattice constant in ourcomputation is 1, the value of d c is an integer. For N = 4, d c is sensitive to initial configuration of rowers and its value liesbetween ( d max − d min ), as shown in column 4. The averagevalues with error are also presented. For the above distribution, we can derive the expressionfor A using Eq. (A2). After calculating the Gaussianintegral we get, A = e − / σ . (A4)Please note that A → σ →
0, and A → σ → ∞ . Appendix II. Critical distance for variousparameters
In Table. A1, we present the value of critical distancebetween two clusters at which phase drift behavior ap-pears for five different sets of internal parameters of cilia( k, s ) for N = 3 and 4 (asymmetric case). For N = 3,we observe that the value of d c is depends on k and s but, independent of α and initial conditions. However,for N = 4 the value of d c is quite sensitive to α andinitial conditions. Appendix III. Synchronized island
Let us recall the definition of a synchronized island: agroup of consecutive rowers (connected by next occupiedlattice sites) beating in a common frequency. We havemeasured this quantity for random and random clustered5heterogeneity for the same data provided in Fig. 13 andFig. 14 in the main text. In Fig. 17, we plot the averagenumber synchronized islands ( (cid:104) N cluster (cid:105) ) and the ratioof the average number of clusters to the average num-ber synchronized islands ( (cid:104) N cluster (cid:105) / (cid:104) N island (cid:105) ). We ob-serve that for all heterogeneity (cid:104) N island (cid:105) decreases withincreasing ρ and the value of (cid:104) N island (cid:105) is large for ran-dom heterogeneity compared to random clustered het-erogeneity, suggesting that the synchronization is morevulnerable to random heterogeneity. From Fig. 17 (b),it is clear that a synchronized island consists of severalclusters. This number increases with decreasing ρ . (a) < N i s l a nd > ρ random[2,4]-[2,gl max ][3,5]-[2,gl max ] (b) < N c l u s t e r > / < N i s l a nd > ρ random[2,4]-[2,gl max ][3,5]-[2,gl max ] FIG. 17. Plots of synchronized island for different het-erogeneities. (a) (cid:104) N cluster (cid:105) is plotted against ρ . (b) (cid:104) N cluster (cid:105) / (cid:104) N island (cid:105) against ρ . Appendix IV. The attractive case α < . In the main text, we have discussed the effect of spa-tial heterogeneity in the case when two rowers oscil-late in anti-phase, and many rowers oscillate collectivelyas metachronal waves. Here, we discuss the effect ofspatial heterogeneity, in the case when rowers show in-phase synchronization. While for a system of rowers,the anti-phase solution can be achieved using a realistic ∆ φ ( t ) t 0 FIG. 18. In phase synchronization for a two rower system.The phase difference ∆ φ = φ − φ is plotted as a functionof time t . At t = 0 the rowers start from an arbitrary initialcondition. The phase difference vanishes as t increases. (a) d ω ω ω ω collective (2) ω -5 -4 -3 -2 -1 (b) | α | τ d r i ft - c on s t d α=−0.2α=−0.02α=−0.002 ) FIG. 19. (a) Angular frequencies of three rowers as a func-tion of d is plotted for α = − .
02. Bifurcation occurs for d ≥
4. For very large d , rower 3 become almost inde-pendent of others and oscillates with its natural frequency ω while first two rowers oscillates with the collective fre-quency of two rower system. (b) The log-log plot of τ drift for α = − . , − .
02 and − . | α | . This scaling makes all the curves fall into asingle curve g ( d ). The function g ( d ) shows 1 /d decay. ω c o ll e c t i v e N α =-0.02 FIG. 20. The collective frequency is plotted as a function oftotal number of rowers N . The rowers are placed on a regularlattice with lattice constant d = 1. The data are averagedover 1000 initial configurations. It decays with system size N . hydrodynamic interaction ( α >
0) (as discussed in themain text), the in-phase synchronization can be realizedthrough a “non-realistic” hydrodynamic coupling withnegative α [A19] or using a negative force constant k ofharmonic driving force of rowers [A34]. We investigatethe case of in-phase synchronization using the same dy-namical equation Eq. (4) but with negative α . Here, wepresent the results for α = − .
02 and the same valuesof the force constant k and amplitude s as mentioned inthe main text are used. Few rowers
In Fig. 18 we show that two rowers, initially havingarbitrary phases, oscillate in the same phase after a tran-6 < ω > ρ (a) random2-gl3-gl[2,4]-[2,gl max ][3,5]-[2,gl max ] C v ρ (b) FIG. 21. Mean frequency and its fluctuation as a function ofdensity ρ for N = 100. For random and random clusteredheterogeneity, the data are averaged over 100 spatial configu-rations. For regular clustered heterogeneity, the data are av-eraged over 100 initial configuration of rowers. For these plotswe consider “2-x” and “3-x” regular clustered, and “[2,4]-[2,x]’ and “[3,5]-[2,x]” random clustered heterogeneity. (a) (cid:104) ω (cid:105) against ρ , and (b) the coefficient of variation C v against ρ . sient period. For a three rower system, the phenomenaof phase drifting has seen as the distance between 2ndand 3rd rowers d is varied. In Fig. 19 (a), we plot thefrequency of all the 3 rowers. Note that a bifurcationoccurred for d ≥ d c (= 4). In Fig. 19 (b), we showthe drift time τ drift for different values of negative α .As in the case of positive α , the similar scaling of τ drift holds for negative α : | α | τ drift ( α, d ) = g ( d ) + const .However, the decay of scaling function g ( d ) is differ-ent from the scaling function f ( d ) for positive α (seemain text). Here, the scaling function g ( d ) decays as g ( d ) ∼ /d which is faster than the decay of f ( d ). Many rowers
We first consider N rowers on a regular lattice withlattice constant d = 1. In this case, all the rowers oscil-late with a single common frequency ω coll,N . The phaseof a rower and its adjacent neighbors are in-phase syn-chronized. In Fig. 20, we plot collective frequency asa function of system size N . We observe that ω coll,N decays with N .Next, we present the results of three different typesof heterogeneous lattices (regular clustered, random andrandom clustered) for α <
0. We find that the spatialinhomogeneity in rowers’ position leads to fluctuation inthe beating frequency of the rowers and reduces the or-der in phase coherence. In Fig. 21(a), we plot averagefrequency of the rowers. We observe that (cid:104) ω (cid:105) decreasesas a function of the density ρ for all the types of het-erogeneous lattices. In Fig. 21(b) we plot the coefficientof variation C v which is a normalize fluctuation of fre-quency. We observe that for a fixed ρ , C v depends onthe type of heterogeneity. For regular clustered C v is < A > ρ randomregular clusteredrandom clustered FIG. 22. Order parameters as a function of density ρ fordifferent structures for N = 100. For random and randomclustered lattice, the data are averaged over 100 spatial con-figurations. For regular heterogeneity, the data are averagedover 100 initial configuration of rowers. For these plots weconsider “2-x” regular lattice and “[2,4]-[2,x]” random clus-tered lattice. For regular lattice, the phases of the rowerscoherent for all the density. For random clustered and ran-dom lattice, the synchronization of phases gets destroyed atlower densities. less and for random heterogeneity C v is high, while forrandom clustered it is intermediate.The phase coherence can be measured by the orderparameter A (defined in the main text). The order pa-rameter A is plotted as function of density ρρ