# Conditions for metachronal coordination in arrays of model cilia

^{a}Department of Living Matter Physics, Max Planck Institute for Dynamics and Self-Organization, Göttingen 37077, Germany;^{b}CAS Key Laboratory for Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China;^{c}School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100190, China;^{d}Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Oxford OX1 3PU, United Kingdom;^{e}School of Mathematics, University of Bristol, Bristol BS8 1TW, United Kingdom;^{f}Department of Physics, Tohoku University, Sendai 980-8578, Japan

See allHide authors and affiliations

Edited by Charles S. Peskin, New York University, New York, NY, and approved July 7, 2021 (received for review February 15, 2021)

## Significance

Motile cilia can coordinate with each other to beat in the form of a metachronal wave, which can facilitate the self-propulsion of microorganisms such as *Paramecium* and can also be used for fluid transport such as mucus removal in trachea. How can we predict the collective behavior of arrays of many cilia coordinated by hydrodynamic interactions, and in particular, the properties of the emerging metachronal waves, from the single-cilium characteristics? We address this question using a bottom-up coarse-graining approach and present results that contribute to understanding how the dynamical self-organization of ciliary arrays can be controlled, which can have significant biological, medical, and engineering implications.

## Abstract

On surfaces with many motile cilia, beats of the individual cilia coordinate to form metachronal waves. We present a theoretical framework that connects the dynamics of an individual cilium to the collective dynamics of a ciliary carpet via systematic coarse graining. We uncover the criteria that control the selection of frequency and wave vector of stable metachronal waves of the cilia and examine how they depend on the geometric and dynamical characteristics of a single cilium, as well as the geometric properties of the array. We perform agent-based numerical simulations of arrays of cilia with hydrodynamic interactions and find quantitative agreement with the predictions of the analytical framework. Our work sheds light on the question of how the collective properties of beating cilia can be determined using information about the individual units and, as such, exemplifies a bottom-up study of a rich active matter system.

Motile cilia are hair-like organelles that beat with a whip-like stroke that breaks time-reversal symmetry to create fluid flow or propel swimming microorganisms under low Reynolds number conditions (1⇓–3). The beat is actuated by many dynein motors, which generate forces between microtubules that cause the cilium to bend in a robust cyclic manner with moderate fluctuations (4, 5). On surfaces with many cilia, the actuating organelles can coordinate with each other and collectively beat in the form of metachronal waves, where neighboring cilia beat sequentially (i.e., with a phase lag) rather than synchronously (6). The flows created from this coordinated beating are important for breaking symmetry in embryonic development (7, 8), creation of complex and dynamic flow patterns for the cerebrospinal fluid in the brain (9, 10), and providing access to nutrients (11). In microorganisms such as *Paramecium* and *Volvox*, the metachronal beating of cilia provides propulsion strategies in viscous environments (12, 13). It has been shown that depending on the parameters, beating ciliary carpets can exhibit globally ordered and turbulent flow patterns (14), which can be stable even with a moderate amount of quenched disorder (15), and that metachronal coordination optimizes the efficiency of fluid pumping (16, 17). Natural cilia have inspired various designs of artificial cilia (18⇓⇓⇓–22), which may be used for pumping fluid (23, 24) and mixing (25), or fabrication of microswimmers (26).

Hydrodynamic interactions have been shown to play a key role in coordinated beating of cilia (27, 28) and mediating cell polarity control (29). To achieve synchronization between two cilia via hydrodynamic interactions, it is necessary to break the permutation symmetry between them [e.g., by exploiting the dependence of the drag coefficient on the distance from a surface (30), flexibility of the anchoring of the cilia (31), nonuniform beat patterns (32, 33), or any combination of these effects (34)]. In addition to the hydrodynamic interactions, the basal coupling between cilia can also facilitate the coordination (35⇓–37).

How can we predict the collective behavior of arrays of many cilia coordinated by hydrodynamic interactions, and in particular, the properties of the emerging metachronal waves, from the single-cilium characteristics? Extensive numerical simulations using explicitly resolved beating filaments (16, 17, 27, 38⇓–40) and simplified spherical rotors (13, 14, 41, 42) have demonstrated that metachronal coordination emerges from hydrodynamic interactions. However, insight into this complex many-body dynamical system at the level that has been achieved in studies of two cilia is still lacking. Here, we propose a theoretical framework for understanding the physical conditions for coordination of many independently beating cilia, which are arranged on a substrate in the form of a two-dimensional (2D) array immersed in a three-dimensional (3D) fluid. We uncover the physical conditions for the emergence of stable metachronal waves and predict the properties of the wave in terms of single-cilium geometric and dynamic characteristics.

## Results

We use a simplified model of a cilium as a force monopole moving along a circular trajectory of radius a above a substrate, represented by a sphere of radius b driven by a force *A*. To theoretically study metachronal coordination, we consider such model cilia on a lattice of spacing ℓ in the x–y plane as shown in Fig. 1*B*. We examine the role of the geometric parameters in determining the collective mode of coordination. We parametrize the orientation of the cilia by the angle θ that the plane of the circular trajectory makes with the *A*). The position of the sphere representing the force monopole along the trajectory, which we parametrize by the polar angle *i*th sphere, is

### Dynamical Equations.

Each cilium is driven independently by a tangential force acting on the bead. The magnitude of the force,

The intrinsic angular speed of each cilium given as**3**, we introduce a coordinate transformation

Using the translational invariance along the substrate, we can express the Blake tensor in the 2D Fourier space *SI Appendix* has the details of the derivation) and recast Eq. **3** in terms of the new coordinate. We then use a separation of timescale between the mean phase and the phase difference to simplify the dynamics (*SI Appendix* has some details of the calculations, and Movies S1–S12 show how a metachronal wave emerges in our simulation). By changing the notation from *SI Appendix* has some details of the calculations, and Movies S1–S12 show how a metachronal wave emerges in our simulation)*Materials and Methods*). Cilium *SI Appendix*, Fig. S1). Subtracting the velocity of

The compact form of Eq. **9** allows us to systematically investigate the conditions under which the array of cilia can admit stable metachronal wave solutions and what determines the direction of propagation and the wavelength of the wave. As we shall see below,

### Dispersion Relation.

Let us now consider a situation where the cilia beat in coordination and generate a metachronal wave of frequency ω and wave vector k. We describe the traveling wave as **9** at the zeroth order and making use of the identity

The dispersion relation is plotted in Fig. 2 for various choices of cilia orientation. Note that the resulting frequencies for all modes are somewhat larger than the single-cilium frequency *A*, namely

### Stability Criterion.

Satisfying the dispersion relation provides a necessary condition for frequencies and wave vectors to represent metachronal waves. However, it does not guarantee that the wave is a stable solution to Eq. **9**. To check the stability of a solution, we can expand Eq. **9** in terms of

In Fig. 2, the wave vectors corresponding to linearly stable solutions of Eq. **9** are shown as (blue) dots. The explicit expression for S allows us to make predictions about the necessary criteria for the stability of the waves. For example, when *A*–*C*). Increasing the angle χ, which amounts to tilting the ciliary beating orbit away from the z axis, further accentuates this feature while allowing for the direction of propagation to deviate from the direction of beating, leading to the formation of dexioplectic or laeoplectic metachronism (Fig. 2 *D* and *E*).

To examine the role of the underlying lattice structure, we consider a triangular lattice, which is characterized by reciprocal lattice vectors *F* shows for the case of *Materials and Methods* and *SI Appendix*, Fig. S1). We thus find that tuning the orientation of the cilia trajectories and controlling the positioning of the cilia in the array provide the possibility to generate metachronal waves with desired wavelengths and directions of propagation.

We note that in our current formulation, the stability criterion is degenerate with respect to the direction of propagation (i.e., **8** guarantees that it will also describe the stability of the modes in terms of the original ϕ coordinate.

### Agent-Based Simulation.

To support the validity of the above analytical description, which is analyzed within the framework of linear stability analysis, we perform numerical simulations based on the governing dynamical equations of the cilia (Eq. **3**). The examples of the time evolution are presented in Fig. 3 for an *A*–*E*) with different tilting angles θ and χ, as well as a triangular lattice with *F*). The cilia are initiated with the same phase

As can be seen in Fig. 3 and Movies S1–S12 (*SI Appendix* has some details of the calculations, and Movies S1–S12 show how a metachronal wave emerges in our simulation), in all cases the cilia coordinate with each other and form a metachronal wave after a transient period. For example, in the case of the cilia on a square lattice with the tilting angles of the trajectory as *A*), the cilia beat in the form of the metachronal wave with the wave vector

These simulations do not use the separation of timescales that was used in the theoretical analysis. The emergence of metachronal waves is over

## Discussion

We have constructed a theoretical framework to study metachronal waves in ciliary arrays, where each cilium is driven independently with the same beat pattern and interacts with the others via hydrodynamic interactions for arbitrary geometric configurations. We calculate the dispersion relation of the system, relating the propagation frequency and the wave vector of the metachronal wave, and observe a relatively small frequency range, with higher frequencies at longer wavelengths. We have found that stable waves correspond to finite domains of wave vector, which are selected with relatively well-defined orientation of propagation that is determined by the geometric characteristics of the ciliary beating pattern and the lattice structure. Our results allow us to predict the role of the different harmonics in the moment decomposition of the beat pattern and the friction, which in turn, can be used to make predictions about control of metachronal waves using external cues, as has been demonstrated in the case of phototaxis of *Chlamydomonas* (45).

Although quantitative studies of the beating patterns of cilia in vivo are still lacking, a precise 3D tracking of the motion of a single cilium was recently performed using murine tracheal cilia attached onto a glass surface (46). In a 100 μM adenosine triphosphate (ATP) solution, the demembranated and reactivated cilium had the beating frequency

## Materials and Methods

### Synchronization of Two Cilia.

Consider two cilia rotating in

In the case of **22**,

## Data Availability

All study data are included in the article and/or supporting information.

## Acknowledgments

We thank Andrej Vilfan and Masao Doi for fruitful discussions. This work has received support from the Max Planck School Matter to Life and the MaxSynBio Consortium, which are jointly funded by the Federal Ministry of Education and Research of Germany, and the Max Planck Society. F.M. received partial support from the Alexander von Humboldt Foundation, Strategic Priority Research Program of Chinese Academy of Sciences Grant XDA17010504, and National Natural Science Foundation of China Grant 12047503. R.R.B. acknowledges a doctoral scholarship from the Engineering and Physical Sciences Research Council (EPSRC) and a University of Bristol Vice-Chancellor’s Fellowship.

## Footnotes

↵

^{1}F.M. and R.R.B. contributed equally to this work.- ↵
^{2}To whom correspondence may be addressed. Email: ramin.golestanian{at}ds.mpg.de.

Author contributions: F.M., R.R.B., N.U., and R.G. designed research; F.M., R.R.B., N.U., and R.G. performed research; F.M., R.R.B., and R.G. analyzed data; and F.M., R.R.B., N.U., and R.G. wrote the paper.

The authors declare no competing interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at https://www.pnas.org/lookup/suppl/doi:10.1073/pnas.2102828118/-/DCSupplemental.

- Copyright © 2021 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

## References

- ↵
- J. Gray

- ↵
- ↵
- ↵
- ↵
- ↵
- E. W. Knight-Jones

- ↵
- ↵
- A. Takamatsu,
- K. Shinohara,
- T. Ishikawa,
- H. Hamada

- ↵
- R. Faubel,
- C. Westendorf,
- E. Bodenschatz,
- G. Eichele

- ↵
- N. Pellicciotta et al

- ↵
- M. B. Short et al

- ↵
- S. L. Tamm,
- T. M. Sonneborn,
- R. V. Dippell

*Paramecium*. J. Cell Biol. 64, 98–112 (1975). - ↵
- ↵
- ↵
- N. Uchida,
- R. Golestanian

- ↵
- N. Osterman,
- A. Vilfan

- ↵
- J. Elgeti,
- G. Gompper

- ↵
- ↵
- M. Vilfan et al

- ↵
- ↵
- T. Sanchez,
- D. Welch,
- D. Nicastro,
- Z. Dogic

- ↵
- F. Meng,
- D. Matsunaga,
- J. M. Yeomans,
- R. Golestanian

- ↵
- ↵
- ↵
- D. Matsunaga et al

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- N. Uchida,
- R. Golestanian

- ↵
- A. Maestro et al

- ↵
- N. Narematsu,
- R. Quek,
- K. H. Chiam,
- Y. Iwadate

- ↵
- K. Y. Wan,
- R. E. Goldstein

- ↵
- ↵
- S. Gueron,
- K. Levit-Gurevich,
- N. Liron,
- J. J. Blum

- ↵
- ↵
- Y. Ding,
- J. C. Nawroth,
- M. J. McFall-Ngai,
- E. Kanso

- ↵
- C. Wollin,
- H. Stark

- ↵
- A. Ghorbani,
- A. Najafi

- ↵
- ↵
- ↵
- ↵
- T. A. Katoh et al

- ↵
- K. Ikegami,
- S. Sato,
- K. Nakamura,
- L. E. Ostrowski,
- M. Setou

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Applied Physical Sciences

- Biological Sciences
- Biophysics and Computational Biology