Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Synchronization in carpets of model cilia: numerical simulations and continuum theory
(USC Thesis Other)
Synchronization in carpets of model cilia: numerical simulations and continuum theory
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
SYNCHR ONIZA TION IN CARPETS OF MODEL CILIA:
NUMERICAL SIMULA TIONS AND CONTINUUM THEOR Y
b y
An up V Kanale
A Dissertation Presen ted to the
F A CUL TY OF THE USC GRADUA TE SCHOOL
UNIVERSITY OF SOUTHERN CALIF ORNIA
in P artial F ulfillmen t of the
Requiremen ts for the Degree
DOCTOR OF PHILOSOPHY
(Mec hanical Engineering)
Decem b er 2021
Cop yrigh t © 2021 An up V Kanale
T o A pp a .
F or giving me the gr e atest gift of al l— curiosity.
ii
A c kno wledgmen ts
प्रथमवयिस दत्त ं तोयमल् प ं स्मरन् तः
िशरिस िनिहतभारा नािरक े ला नराणाम ् ।
सिललमम ृ तकल् प ं दद्य ु राजीवनान् त ं न िह क ृ तम ु पकार ं साधवो िवस्मरिन् त ।।
This PhD has b een a steep learning curv e for me, p ersonally and professionally , and
w ould not ha v e b een p ossible without help from man y w onderful p eople in m y life. I w ould
lik e to tak e this opp ortunit y to express m y gratitude to them.
Firstly , I w ould lik e to thank m y advisor Prof. Ev a Kanso for her guidance and patience
throughout the course of this PhD. I ha v e learn t a great deal from her not only ab out pursuing
researc h problems, but also comm unicating scien tific researc h in the most effectiv e manner.
Her high standard of preparing and deliv ering presen tations, and her ey e for visuals are
things I will alw a ys lo ok up to, and will hold me in go o d stead for the rest of m y professional
life. I could not ha v e undertak en this w ork without her supp ort.
I thank our researc h collab orator Dr. Sebastian F ürthauer for b eing patien t with me
as I w as learning new to ols and tec hniques outside of m y comfort zone. Sebastian’s depth
of kno wledge as a ph ysicist, his c heerfully curious p ersonalit y , and encouraging w ords– “Ah
iii
y ou need to mak e as man y mistak es as y ou can as fast as y ou can. That’s ho w y ou learn the
most!”– will con tin ue to b e a source of inspiration for me.
I w ould lik e to thank Divy a, Appa, Amma and Vibha for putting up with me in m y diffi-
cult times. Their constan t, un w a v ering supp ort and encouragemen t help ed me sail through
some truly difficult times in m y PhD. I’m esp ecially thankful to Divy a for nev er losing faith
in m y abilities.
Thanks m ust also go to m y committee mem b ers– Profs. Luhar, Newton, Nakano and
Bermejo-Moreno for excellen t feedbac k on m y prop osal and constan t encouragemen t. Be-
sides, the non-linear dynamics class with Prof. Newton and the computational ph ysics class
with Prof. Nakano ha v e b een the t w o most enjo y able and illuminating courses I ha v e tak en
in m y time at USC. Mitul has b een a v ery imp ortan t figure in m y life as a grad studen t– first
as m y MS researc h advisor, then as the instructor for the class I graded for (whic h help ed me
pa y m y ren t!) and most imp ortan tly , as an approac hable, wise presence I ha v e alw a ys felt
comfortable reac hing out to. I also thank Prof. Bermejo-Moreno for his suggestions during
m y qualifying exam.
I’m thankful to F eng, not only for helping me dev elop the idea of Kuramoto ellipse, but
also for man y con v ersations ranging from science, to philosoph y , to gym w ork out sessions!
With his immense depth and breadth of kno wledge, I lo ok forw ard to seeing exciting researc h
coming from him in the near future. I thank Hanliang for b eing m y first researc h men tor,
and helping me get started on the pro ject.
I thank all the other mem b ers in the Kanso lab that I ha v e crossed paths with. My
in teractions with Brendan ha v e alw a ys b een uplifting and ha v e had an impact on m y view
of researc h philosoph y . Sina has b een a great labmate and friend throughout these 5 y ears.
iv
I thank Y usheng, Basile, Yi, Lionel, Chenc hen, Jingyi and Haotian for main taining a w arm
and w elcome lab atmosphere.
I w ould lik e to thank Girish anna, Reshma attige and the kids. I ha v e nev er felt homesic k
in m y 7 y ears in the US, and it is solely due to their w elcoming presence nearb y . I w ould
also lik e to thank m y in-la ws for their constan t supp ort, encouragemen t and patience.
I w ould lik e to thank Christoph, Bo om, V amsi and Shan tan u for b eing fan tastic p eers
at AME. I will miss their compan y and con v ersations with them ab out researc h, academia,
science, and life in general.
Graduate life can get quite isolat ing, but I consider m yself fortunate to ha v e had the
compan y of man y great friends. I thank Daka, Na v een and Krishna for man y in tellectually
stim ulating con v ersations and trading advise ab out almost ev erything under the sun. I thank
Bhan u, Sra vy a and Akanksha for the relaxing w eek end get-togethers. I w ould lik e to thank
Harshitha, Nit y a, Bala, Aniruddha, Dino, Gautam, Jitin, Praneet, Shalini, P o o ja, Divy a
Ramesh, Shashank for their supp ort.
Somewhat uncon v en tionally , I w ould also lik e to thank p eople who I ha v e nev er met,
but ha v e had a great impact on m y e n th usiasm for science. I thank Profs. Stev en Strogatz,
Man u Prakash and Ra y Goldstein for b eing inspirational figures to lo ok up to as a y oung
researc her, and finding and pursuing suc h exciting researc h questions. Finally , I w ould lik e
to thank Gran t Sanderson of 3blue1bro wn for setting the highest standard of comm unicating
science through binge-w orth y animations.
v
Con ten ts
Dedication ii
A c kno wledgmen ts iii
List of Figures ix
List of Algorithms xv
Abstract xvi
Thesis o v erview xvii
1 Sync hronization in biology 1
1.1 Dance of the micr oswimmers . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
2 Flo w ph ysics at lo w Reynolds n um b er 7
2.1 The Na vier Stok es equa tions . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 The Reynolds n um b er . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.3 Stok es equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.4 The Oseen tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
vi
2.5 The Blak e tensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
3 Mathematical mo dels of ciliary sync hronization 16
3.1 Literature surv ey . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.2 Discrete cilia mo d el . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
4 A simple algorithm for large scale sim ulations of ciliary carp ets 22
4.1 Algorithm for ev aluating the Doubly-p erio dic Blak e T ensor . . . . . . . . . . 24
4.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.2.1 Emergence of metac hronal w a v es . . . . . . . . . . . . . . . . . . . . 32
4.2.2 Flo w pumping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
5 A Con tin uum mo del of ciliary carp ets 38
5.1 Calculating cilia-i nduced fluid flo ws . . . . . . . . . . . . . . . . . . . . . . . 39
5.2 Order parameter field s for sync hronization and their equations of motion . . 43
5.2.1 Calculating cilia induced forces a nd fluid flo ws . . . . . . . . . . . . . 48
5.3 Stabilit y of isotropic a nd sync hronized states . . . . . . . . . . . . . . . . . . 50
5.3.1 Stabilit y of the isotropic state . . . . . . . . . . . . . . . . . . . . . . 51
5.3.2 Stabilit y of the sync hronized state . . . . . . . . . . . . . . . . . . . 53
5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
6 Comparison of the n umerical and theoretical approac hes 60
6.1 Commen t on the Kuram oto order parameter . . . . . . . . . . . . . . . . . . 61
6.2 The Kuramoto ellipse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
vii
6.3 Comparing gro wth rates of the isotropic state . . . . . . . . . . . . . . . . . 64
6.4 Comparing gro wth rates of the fully sync hronized state . . . . . . . . . . . . 66
6.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
7 Conclusions 70
Bibliograph y 72
App endices 89
A Ev aluating the b oundary co nditions to obtain the fluid in teraction Kernel . . 90
B An alternativ e metho d to analyze linear stabilit y of the fully sync hronized state 95
C Comparison of cilia-i nduced flo wfields . . . . . . . . . . . . . . . . . . . . . . 100
viii
List of Figures
1.1 Illustration of the exp erimen t with p endulums, b y Huygens himself. Image
obtained from [1] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Illustration of cilia in the lung epithelium. Image credit: [2] . . . . . . . . . 4
1.3 Ciliary co ordination across length scales. Image credit: [3] Repro duced under
the fair use p olicy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4 Sync hronization of cilia and flagella: (a) Sp erm flagella b eat in-phase in close
pro ximit y [4]. (b) Chlam ydomonas cell swims in a breaststrok e fashion [5].
(c) E. c oli flagella bundle up [6]. (d) Cilia self-organize in to metac hronal
w a v es [7]. Image credit: [8] . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.1 An organism with a single degree of freedom lik e a scallop cannot ac hiev e net
lo comotion at lo w Reynolds n um b ers. Image credit: Wikip edia . . . . . . . . 11
2.2 Sc hematic represen tation of a force monop ole ab o v e a no-slip w all. Without
loss of generalit y , the p oin t force c an b e placed at (0,0,h) . . . . . . . . . . . 12
2.3 Flo w induced b y a unit strength stok eslet in free space (left), a stok eslet +
an tistok eslet (middle) and a stok eslet b ound b y a w all atz = 0 , i.e. a Blakelet
(righ t). The flo wfields are sho wn in the top ro w, and v elo cit y profiles in the
z -direction, at (x,y) = (0.5,0.5) are sho wn in the b ottom ro w. . . . . . . . 15
3.1 Simplified represen tation of b eating cilium: A oscillating filamen t is appro xi-
mated as a rigid b ead mo ving in a c ircular tra jectory . . . . . . . . . . . . . . 18
ix
3.2 (a) Exp erimen tally measured flo wfield around a c hlam ydomonas, (b) Flo w-
field induced b y a 3-stok eslet mo del. This sho ws that lo w er order mo dels are
go o d appro ximations to the flo wfields seen in the real w orld. This image has
b een repro duced from [9], repro duced h ere under the F air use p olicy . . . . . 19
3.3 (A) Sc h ematic represen tation of the ciliary carp et arranged in a doubly p eri-
o dic square lattice, ab o v e a no-slip b oundary (sho wn primary b o x only). (B)
First (blac k) and second (blue) force h armonics of the driving force. . . . . . 19
4.1 Sc hematic represen tation of the computational domain, consisting of the fun-
damen tal domain (red), images (grey) and in teraction stencil (blue). The
zo omed in view of the in teraction stencil sho ws the fo cal cilium at the cen ter
of the lattice (dark blue) and the neigh b oring cilia (red). W e compute the in-
duced v elo cit y at the fo cal cilium b y considering the influence of all neigh b ors
within the in teraction stencil. . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.2 Con v ergence of fluid v elo cit y . A : Cosine fields of a system of 151×151
cilia, in three sample states. B : Flo w sp eed (blac k) induced at the cen ter of the
lattice approac hes the result obtained using the F ast Multip ole Metho d [10]
(gra y) for increasing shells of in teracting neigh b ors. C : The relativ e p ercen t-
age difference in the induced flo w sp eed, as compared to the FMM result
(gra y) and due to increasing n um b er of s hells (blac k) as sho wn in (4.1) . . . 27
4.3 (A) Induced v elo cities at the lo cation of the fo cal cilium (blue), b y ev ery
other cilium inside the in teraction stencil (red) are pre-computed and stored.
(B) The tra jectory of eac h cilium is discretized in to n
θ
angular lo cations.
v
p,q
i←j
is the induced v elo cit y at angular lo cation p of the fo cal cilium i , due a
neigh b oring cilium, j , at angular lo cation q , and is computed using the Blak e
tensor. (C) Sc hematic represen tation of the v elo cit y lo okup table. Eac h en try
is a 2-dimensional v elo cit y v ector p opula ted as p er the description in (B). . . 28
4.4 The induced v elo cit y on ciliumj , due to a neigh b oring ciliumi is obtained b y
lo oking up the four en tries from the table, and linearly in terp olating b et w een
them. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
x
4.5 Log scale plots of (A) relativ e error, and (B) computational time as a function
of n um b er of shells. Solid and dashed lines indicate sim ulations with and
without lo okup, resp ectiv ely . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
4.6 W e p erform sim ulations using a carp et of 151×151 cilia, arranged in a square
lattice with constan t spacing d = 1 , at a heigh t h = 0.5 ab o v e a no-slip w all,
for100 time units. W a v e-lik e spatial patterns emerge when the individual cilia
are driv en b y the first (left) and second (righ t) force harmonics. Sim ulation
parameters used: Isotropic initial state, b/d = 0.2, a/d = 0.05 , 10 shells of
in teracting neigh b ors, and using n
θ
= 36 discrete p oin ts along the tra jectory
of eac h cilium. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
4.7 Cosine field obtained from sim ulations at t = 100 units, for increasing v al-
ues of cilium heigh t. P o c k ets of larger w a v elength emerge for longer cilia.
Sim ulations parameters used: Isotropic initial state, second force harmonic,
b/d = 0.2, a/d = 0.05 and 10 shells of in teracting neigh b ors, and using
n
θ
= 36 discrete p oin ts along the tra jectory of eac h cilium. . . . . . . . . . 35
4.8 Cum ulativ e directional flux induced b y the 151×151 carp et of cilia, o v er 100
units of time, for the first (blac k) and second (blue) force harmonics. . . . . 36
5.1 Ciliary La y er: A doubly-infinite arra y of cilia, wherein eac h cilium mo deled as
a b ead rotating on a fixed circular tra jectory under the influence of a p ositiv e
phase-dep enden t force F(θ) is appro ximated b y a force densit y la y er in the
limit of infinite n um b er of cilia. . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.2 In a system ofN phase oscillators (cilia in our case), eac h momen t of the order
parameter Y
n
, captures a differen t kind of phase order around the circle. The
oscillators ha v e b een sho wn on the sam e circle for ease of illustration. . . . . 44
xi
5.3 T op ro w: Sample force profiles, F(θ) , corresp onding to the first through the
fourth harmonics in (3.1). Middle ro w the transformation from angular phase
θ to a new phase ϕ in whic h the in trinsic angular sp eed, Ω 0
= 2π/T
0
, of a
single cilium is constan t (sho wn in inset). Here, T
0
is the in trinsic time p e-
rio d defined b y
R
2π
0
d
¯
θ/Ω(
¯
θ) and is in v arian t with the force harmonic, since
the in tegrals of the sine and cosine terms drop out. Bottom ro w: Cilia that
are uniformly distributed in ϕ (ϕ ∈ Unif[0,2π) ) exhibit a non-uniform dis-
tribution in θ that reflect main features of the force harmonic: single-p eak
distribution for the first force harmonic, double-p eak distribution for the sec-
ond force harmonic, and so on. . . . . . . . . . . . . . . . . . . . . . . . . . 46
5.4 Stabilit y of the isotropic state. The gro wth rate of the instabilit y in the
F ourier space, γ(k) , computed from (5.41), is sho wn for first four force har-
monics. In all four cases, almost all mo des are linearly unstable, with the
instabilit y b eing particularly strong for t he second force harmonic. . . . . . . 53
5.5 Stabilit y of the fully sync hronized state. The Ly apuno v exp onen t, µ (k) , com-
puted from (B.14), is sho wn for first four force harmonics. Unlik e the isotropic
state, in all four cases, some mo des are linearly stable while others are not.
But the instabilit y is again particularl y strong for the second force harmonic. 57
xii
6.1 (A) Numerical sim ulations of a 151×151 arra y of cilia sho w the emergence of
tra v eling w a v e patterns. Snapshots tak en at t = 100 for h/d = 0.5 (left) and
h/d = 1.5 (righ t). (B) Time ev olution of the Kuramoto p olar order parameter
P sho ws no clear distinction b et w een isotropic and tra v eling w a v e patterns.
(C) W e define P
α
for α∈ [0,2π] and plot (α,P
α
) in p olar co ordinates at t w o
snapshots tak en from the isotropic state and tra v eling w a v e pattern. Data
p oin ts are sho wn in blac k and elliptic b est-fit in blue. (D) Time ev olution
of the eccen tricit y and angle of the Kuramoto ellipse sho ws clear correlations
b et w een i ts eccen tricit y with the gro wth of the tra v eling w a v e instabilit y and
angle with dominan t direction of w a v e fron t (see Supplemen tal Mo vie). Sim-
ulations w ere p erformed for 100 time units using dt = 0.1 , a = 0.05 , b = 0.2 ,
d = 1 . First force harmonic: F = 1 + 0.5cosθ + 0.5sinθ , Second force
harmonic: F = 1+0.5cos2θ +0.5sin2θ . . . . . . . . . . . . . . . . . . . . . 62
6.2 Illustration of Kuramoto ellipse using snapshots of an isotropic state (top)
and a syn thetically created w a v e state (b ottom). Left column sho ws the
cosine field, obtained b y taking the cosine of the phase of the cilium. Cen ter
column sho ws the corresp onding Kuramoto order field– Real (Y
1
(x)) – whic h
is calculated b y coarse graining o v er one shell of neigh b ors. The righ tmost
column sho ws the p oin ts (α,P(α)) in green, and the corresp onding Kuramoto
ellipse in gra y . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
6.3 Isotropic state. (A) Gro wth rate γ(k) of the instabilit y predicted b y the con-
tin uum theory (colormap) and discrete sim ulations (blac k dots) in the F ourier
space. (B) Starting from a random Kuramoto field Y
1
(x) , map to F ourier
space
ˆ
Y
1
(k) and use γ(k) from panel A to step forw ard in time. Snapshots of
Re|Y
1
(x,t)| sho w the emergence of tra v eling w a v es similar to those obtained
from direct n umerical sim ulations of Eq. (3.5). (C) Kuramoto ellipse in tro-
duced in Fig. 1E obtained from the sample sim ulation in panel B (blac k dots)
and theory corresp onding to 200 realizations (grey), with a v erage sho wn in
dark grey . (D) Time ev olution of eccen tricit y and angle of the Kuramoto
ellipse; sample sim ulation (blac k line) a nd theory (grey). . . . . . . . . . . . 68
xiii
6.4 Sync hronized state: Ly apuno v exp onen t µ (k) obtained from theory and
sim ulations (A) o v er the en tire w a v e space (colormap), and (B) at cross-
sectionsk
y
d =−π/4 andk
y
=−π/2 with theory in solid lines and sim ulations
in dashed lines. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
A.1 Stabilit y of the fully sync hronized state. The Ly apuno v exp onen t, µ (k) , com-
puted from (B.14), is sho wn for first fo ur force harmonics. . . . . . . . . . . 99
A.1 Comparison of flo w fields for a unit stok eslet in half space, orien ted in the p os-
itiv e x -direction. (A) Flo w fields in the ciliary plane based on the con tin uum
theory (left) and particle mo del (righ t), resp ectiv ely . V elo cit y magnitude is
sho wn as a colormap and streamlines in white lines. In the particle mo del,
w e use the Blak e tensor in the half-space together with a truncated image
system to em ulate the effect of the doubly-p erio dic domain (B) Flo ws in a
plane p erp endicular to the ciliary plane. (C) Absolute v alue of the difference
in the x -comp onen t of the v elo cit y along y = 0 (solid) and x = 0 (dashed) at
z = h . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
A.2 (A) Sc hematic depiction of a lattice of cilia at random phase. (B) Flo w field
in the ciliary plane, flo w magnitude (colormap) and streamlines (white lines),
based on the con tin uum theory (left) and particle mo del (righ t), resp ectiv ely .
The red arro ws indicate the relativ e strength and orien tation of eac h Stok eslet.
(C) Flo w field in a plane p erp endicula r to the ciliary plane. . . . . . . . . . 102
A.3 (A) Sc hematic depiction of a lattice of fully sync hronized cilia. (B) Flo w field
in the ciliary plane, flo w magnitude (colormap) and streamlines (white lines),
based on the con tin uum theory (left) and particle mo del (righ t), resp ectiv ely .
The red arro ws indicate the relativ e strength and orien tation of eac h stok eslet.
(C) Flo w field in a plane p erp endicula r to the ciliary plane. . . . . . . . . . 103
xiv
List of Algorithms
1 Building the v elo cit y l o okup table . . . . . . . . . . . . . . . . . . . . . . . . 29
2 Computing induced v elo ci ties using lo okup . . . . . . . . . . . . . . . . . . . 32
xv
Abstract
Motile cilia are fundamen tal micro-actuators in cellular biology that b eat cyclically to driv e
flo ws. In dense ciliary carp ets, suc h as those in the mammalian lungs, brains, and repro duc-
tiv e tracts, fluid transp ort emerges from the co ordinated activit y of thousands of m ulticiliated
cells, eac h con taining h undreds of cilia. Despite progress in analyzing cilia co ordination in
certain mo dels, a general theory for co ordination in the limit of large arra ys o f cilia remains
lac king. Starting from discrete arra ys of cilia wherein eac h cilium is represen ted b y a w ell-
kno wn oscillator mo del, w e devise a fast n umerical algorithm for in v estigating the dynamics
of thousands of h ydro dynamically in teracting cilia. W e then dev elop a con tin uum theory
in the limit of infinitely-man y indep enden tly b eating cilia b y com bining to ols from activ e
matter and classical Stok es flo w. Isotropic and sync hronized co ordination states are unstable
in b oth discrete and con tin uum mo dels. T ra v eling w a v es emerge in b oth theories regard-
less of initial conditions, indicating that these emergen t w a v es are stable global attractors,
with w a v elength and direction determined b y the microstructure and activit y at the cilium
lev el. Our theory op ens up the prosp ect of establishing causal explanations of the effect of
individual cilia oscillations on large scale co ordination and fluid transp ort in ciliary carp ets.
xvi
Thesis o v erview
The con ten t in this thesis has b een organized in to the follo wing c hapters:
• W e motiv ate the problem in Chapter 1, with examples of a wide v ariet y of sync hroniza-
tion phenomena observ ed in nature. W e discuss the historical dev elopmen ts p ertaining
to sync hronization at lo w Reynolds n um b ers, and p oin t out the role of fluid-mediated
sync hronization.
• Chapter 2 is a primer on lo w Reynolds n um b er fluid dynamics. W e discuss the p e-
culiarit y of flo ws in viscosit y-dominated regimes, some useful in terpretations of the
Reynolds n um b e r, and presen t the go v erning equations of fluid flo w.
• In Chapter 3, w e pro vide an accoun t of the mathematical mo dels that ha v e previously
b een used to study ciliary sync hronization. W e also discuss eac h of their pros and cons,
and mak e a c hoice for the mo del to use for the purp oses of this thesis.
• In Chapter 4, w e devise an algorithm that enables large scale n umerical sim ulations of
carp ets of cilia, and presen t results of these sim ulations. W e rep ort the emergence of
w a v e-lik e patterns emerge in long time from and in tro duce no v el metrics to quan tify
them.
xvii
• In Chapter 5, dev elop a rigorous con tin uum theory to prob e the b eha vior of infinitely
large ciliary carp e ts. Sp ecifically , w e deriv e an expression for the v elo cit y induced b y
a carp e t of cilia, and use it to prob e the stabilit y of t w o p erio dic steady states.
• In Chapter 6 w e compare the n umerical results from c hapter 4, against the theoretical
predictions in c hapter 5 b oth qualitativ ely and quan titativ ely . W e p erform quan titativ e
comparison b y using the gro wth rate of the w a v e instabilit y as the metric, and sho w
that they a re in excellen t agreemen t with eac h other.
• W e conclude this thesis in Chapter 7. W e p oin t out the significance and implications of
our results and discuss with some ideas for near-term extensions of our w ork. Finally ,
w e offer some closing though ts on long-term op en questions left in the understanding
of ciliar y co ordination.
The author w ould lik e to declare here that significan t parts of c hapters 4 , 5 and 6
ha v e b een repro duced from three indep enden t man uscripts under preparation, for whic h the
author of this thesis is the lead author.
xviii
Chapter 1
Sync hronization in biology
S
ynchr oniza tion originates from the Greek w ords “c hronos” meaning time, and “syn”
meaning same. It underlies a strikingly large n um b er of phenomena around us– natural
and man-made. F or instance, gra vitational sync hron y causes the p erfectly cyclic seasonal
v ariations on our planet. A group of m usical instrumen ts pla y ed randomly generate noise,
but pla y ed in an orderly fashion they create m usic. F rom the momen t w e are b orn, the
neurons in our brain start firing sync hronously and the pacemak ers cells k eep our heart
running. Soldiers are ask ed to break lo c kstep while crossing bridges. Birds are kno wn to
flo c k together to fly efficien tly and fish sc ho ol collectiv ely to escap e predators. These are,
but a few examples of sync hron y w e are all familiar with. As scien tists, the natural question
to ask is, what causes this sync hron y? Its origin is not ob vious, esp ecially in view of the
second la w of thermo dynamics whic h dictates that the w orld m ust mo v e to w ards a state of
higher disorder!
The first record of scien tific inquiry in to the mec hanism b ehind sync hronization, is due
to the Dutc h ph ysicist Christian Huygens. He observ ed that a pair of p endulum clo c ks
1
Figure 1.1: Illustration of the exp erimen t with p endulums, b y Huygens himself. Image
obtained from [ 1 ]
moun ted on a w all tended to sp on taneously sync hronize after some time. He ev en tually
concluded it w as due to elastic coupling b y the w all, thereb y , sho wing for the first time,
that inanimate ob jects could sp on taneously sync hronize using simple ph ysical mec hanisms.
Mo dern realizations of this setup ha v e b een created using metronomes placed on a w o o den
plank, whic h in turn, is moun ted on empt y so da cans; a quic k searc h of y outub e turns up
man y videos demonstrations. Elastic coupling from the plank and the freedom to mo v e due
to the cok e cans causes the metronomes to sync hronize o v er time. Using biological oscilla-
tors suc h as fireflies, the cardiac pacemak er, etc. as the inspiration, rigorous mathematical
foundations for the theory of sync hronization w ere laid in the latter half of the 20th cen-
tury [ 11 , 12 , 13 , 14 ]. F or more information ab out the ric h b o dy of literature in this field,
see [ 1 ] for a textb o ok treatmen t and the review article [ 15 ]. F or a p opular science treatmen t
of this sub ject, the author of this thesis highly recommends the p opular science b o ok Sync
b y Prof. Stev en Strogatz [ 16 ].
The examples men tioned so far all refer to sync hronization among macr osc opic en tities,
2
but there are plen t y of examples of sync hron y in the micr osc opic w orld as w ell. This, in
fact, is the cen tral driving theme of this thesis. In the section that follo ws, w e w alk through
examples of sync hronization observ ed at v ery tin y , micro-metric, length scales induced b y
fluid coupling.
1.1 Dance of the microswimmers
A t the same time as Huygens w as observing sync hronizing p endula, his fello w Dutc h scien-
tist, An toni V an Leeu w enho ek w as w atc hing micro organisms under his microscop e, making
observ ations that con tin ue to riddle bioph ysicists cen turies later, to this da y in 2021!
Micro organisms use cilia and flagella, whic h are tin y hair lik e organelles that protrude
from cell b o dies, to lo comote. They basically use them as oars; w e will address the rele-
v an t flo w ph ysics in c hapter 2 . Cilia and fla gella
1
are found ubiquitously in biology– from
single-celled micro organisms lik e chlamydomonas , to higher mammals lik e h umans, and their
structure is remarkably conserv ed across ev olutionary time scales. They are ab out 5−10µm
in length and serv e the purp ose of lo comotion and pumping at micrometric length scales b y
a whip-lik e bac k-and-forth motion.
But wh y should they b e of in terest? Apart from pure, unadulterated scien tific curiosit y ,
it is b ecause they pla y a k ey role in executing ph ysiological pro cesses in v ertebrates. A few
examples of their role are: pushing dirt-laden m ucus out of the respiratory tract [ 17 ], driving
cerebrospinal fluid in the nerv ous system [ 18 ], assisting fertilization b y helping the sp erm cell
1
The difference b et w een cilia and flagella seems situational, rather than functional. When flagella o ccur
in higher n um b ers, lik e in carp ets, they seem to b e referred to as cilia. W e will adopt the ph ysicist’s p oin t
of view and use them in terc hangably throughout this thesis.
3
Figure 1.2: Illustration of cilia in the lung epithelium. Image credit: [ 2 ]
(whic h swims using a flagellum) swim to o v a in the F allopian tub e of the female repro ductiv e
system [ 19 , 20 ]. Dysfunctional cilia ha v e b een link ed to c hronic lung infections [ 21 ], infertilit y
and man y other diseases[ 22 ]. See figure 1.2 for an illustration of the role of cilia in the h uman
respiratory tract. Scien tists [ 23 ] studying em bry ogenesis disco v ered that the cilia inside the
Hensen’s no de collectiv ely induce a flo w, kno wn as the no dal flo w, whic h is resp onsible for
breaking the left-righ t symmetry in a mouse em bry o. In an astonishing result, they found
that mec hanically altering the flo w to go in the opp osite direction lead to the mouse’s en tire
b o dy plan b eing flipp ed!
A k ey feature of cilia is that their in ternal structure has b een conserv ed across ev olution-
ary time scales. The cilia in lo w er organisms lik e V olvox and p ar ame cium ha v e an in ternal
structure that is v ery similar to the cilia found in mammals. An upshot of this feature is
that it offers the ease of studying them exp erimen tally using lo w er, simpler organisms, and
4
Figure 1.3: Ciliary co ordination across length scales. Image credit: [ 3 ] Repro duced under
the fair use p olicy .
translating the findings to higher organisms. A more imp ortan t b enefit, ho w ev er, is that
a thorough understanding of cilia and their role in biology could shed ligh t on one of the
deep est outstanding questions in dev elopmen tal biology– the origin of multic el lularity .
Cilia displa y co ordinating b eating patterns across a wide range of length scales [ 3], as
illustrated in figure 1.3 . They ha v e an in tricate in ternal structure p o w ered b y motor proteins
that co ordinate to pro duce oscillatory motion at the lev el of an individual cilium. In addition
to this, a large n um b er of cilia on tissue can sp on taneously co ordinate their b eat to pro duce
co ordination across length scales m uc h larger than that of an individual cilium itself. A
few examples are sho wn in Fig. 1.4 . The flagella of sp erm cells are kno wn to sync hronize
in-phase in close pro ximit y [ 4 ], the flagella of the alga Chlamydomonas are kno wn to displa y
in-phase and an ti-phase gaits of swimming [ 5 ], and cilia on paramecia b eat in metac hronal
co ordination to prop el the organism in its living en vironmen t.
In a classic w ork, T a ylor [ 24 ] appro ximated sp erm flagella as infinite w a ving s heets and
pro v ed that sync hronizing in-phase minimized energy lost in the form of viscous dissipation.
5
Figure 1.4: Sync hronization of cilia and flagella: (a) Sp erm flagella b eat in-phase in close
pro ximit y [ 4 ]. (b) Chlam ydomonas cell swims i n a breaststrok e fashion [ 5 ]. (c) E. c oli flagella
bundle up [ 6 ]. (d) Cilia self-organize in to metac hronal w a v es [ 7 ]. Image credit: [ 8 ]
Ho w ev er, minimizing energy cannot b e a generic dynamic al principle b ehind sync hronization.
This naturally leads to the question– ho w, then, do cilia sync hronize at all? Since most
cilia found in biology are not connected to the cen tral nerv ous system, there is no cen tral
“leader” causing them to sync hronize. Some of the mec hanisms prop osed include bio c hemical
signalling, basal coupling and flo w-induced coupling. The emerging picture seems to suggest
that h ydro dynamic pla ys a k ey role in ciliary sync hronization. T o b etter understand this, the
next c hapter will b e dedicated to w ards understanding the flo w ph ysics on tin y , micrometric
scales in whic h cilia op erate.
6
Chapter 2
Flo w ph ysics at lo w Reynolds n um b er
F
luid dynamics is that branc h of ph ysics that deals with the mo v emen t of fluids. It has
a div erse range of applications, from b oth basic science and engineering p ersp ectiv es.
F o r example, it is used b y engineers in designing ro c k et propulsion systems, airplanes and
turbines, among man y other applications. On the other hand, it is also used to understand
flo w of glaciers, the swimming of fish and lo comotion of bacteria. The Na vier-Stok es equa-
tions are a set of non-linear partial differen tial equations that go v ern the flo w of fluids. This
c hapter serv es as a primer on the fluid dynamical principles at the length scales relev an t to
cilia and flagella, wherein, the b eha vior of fluids can b e coun ter to our in tuition from daily
life.
7
2.1 The Na vier Stok es equations
An incompressible, Newtonian fluid of densit y ρ and viscosit yη , mo ving with v elo cit yu(x,t)
satisfies the Na vier Stok es equations:
ρ
∂u
∂t
|{z}
unsteady term
+ (u·∇)u
| {z }
inertial term
!
= −∇p
| {z }
pressure term
+ η∇
2
u
| {z }
viscous term
+f, ∇·u = 0, (2.1)
where p(x,t) is the pressure field, and f includes all the external forces acting on the fluid
(gra vit y , for instance). The ph ysical significance has b een lab elled under of ea c h term in the
equation ( 2.1 ).
2.2 The Reynolds n um b er
The Reynolds n um b er, arguably the most imp ortan t dimensionless group in all of fluid
dynamics, is a n um b er that captures the relativ e imp ortance of inertial and viscous terms
in ( 2.1 ). Consider a b o dy of length L , immersed in a fluid of densit y ρ and viscosit y η ,
swimming at a sp eed U . The inertial term then has the scaling∼ ρU
2
/L , while the viscous
term scales as∼ ηU/L
2
. The Reynolds n um b er can then b e expressed mathematically as,
Re =
F
inertial
F
viscous
=
ρUL
η
. (2.2)
A h uman swimming in w ater t ypically swims at ∼ 1m/s and has a b o dy length of ∼ 2m .
Using ρ = 1000kg/m
3
and η = 10
−3
Pa·s . This giv es Re∼ 10
6
. The Re for submarines and
airplanes is ev en higher. In con trast, the bacterium E. c oli has a b o dy length of∼ 2µm and
8
swims at a sp eed of∼ 30µm/s . This giv es a v ery small Re of∼ 10
−6
.
But what do these n um b ers r e al ly me an ? In order to answ er this question, w e follo w the
in terpretation laid out b eautifully in [ 25 ] of the Reynolds n um b er as a dimensionless coasting
distance. T o see wh y , w e ask the follo wing question: if the swi mmer stops the swimming
motion at time t = 0 , then ho w long do es it coast for b efore coming to a complete stop?
W e no w try to estimate the order of magnitude of this time, sa y τ , in eac h of inertial (eg.
h uman swimming) and viscosit y-dominan t (eg. bacterium swimming) flo w regimes.
In inertial regimes, the swimmer m ust decelerate follo wing Newton’s la ws of motion.
Assuming the swimmer has a densit y ρ
s
, the deceleration is of the order ∼ ( mass) ×
( acceleration) = ρ
s
L
3
U/T . The drag force has an inertial scaling of ρU
2
L
2
. A balance
b e t w een these t w o forces yields τ ∼ Lρ
s
/ρU , or equiv alen tly , the distance coasted b y the
swimmer b efore coming to a stop is ∼ UT ∼ ρ
s
/ρ . Since h umans ha v e a densit y that is of
the same order of magnitude, w e ma y conclude that the coasting distance for us is of the
same order of magnitude as our b o dy length.
In viscosit y-dominated regimes, the fundamen tal difference is that the fluid drag scales
linearly with the v elo cit y∼ ηLU . Rep eating the analysis from the previous paragraph yields
τ ∼ ρ
s
L
2
/η , and a normalized coasting distance of∼ ρ
s
LU/η = Reρ
s
/ρ . Hence, in absolute
terms, a bacterium coasts for a distance less than the atomic length scale. This also leads
us to another useful ph ysical insigh t– that the bacterium m ust con tin uously exp end energy
to k eep mo ving.
W e end this section with a rumination to put things in p ersp ectiv e– in order to exp erience
what it feels lik e to b e a swimming bacterium, a h uman w ould ha v e to swim in tar !
9
2.3 Stok es equations
Normalizing eac h quan tit y with the appropriate c haracteristic scale, equation ( 2.1 ) can
rewritten in terms of the Reynolds n um b er as
Re
∂u
∂t
+(u·∇)u
!
=−∇p+η∇
2
u+f, ∇·u = 0. (2.3)
As discussed in the previous section, Re≪ 1 at small length scales. T aking this to the limit
of Re→ 0 , w e arriv e at the Stokes e quations :
−∇p+η∇
2
u+f = 0, ∇·u = 0, (2.4)
whic h are the appropriate equations to use to mo del bacterial scale flo ws. W e will return to
equation ( 2.4 ) frequen tly in the future c hapters of this thesis.
Notice that terms that dropp ed out of the full Na vier-Stok es equations, b y taking the lim-
its Re→ 0 , w ere the non-linear and time-dep enden t terms. This results in Stok es equations
ha ving t w o unique prop erties:
1. Line arity : whic h implies that solutions can b e constructed b y linear sup erp osition in
a v ariet y o f geometries; and
2. Time indep endenc e : time has no effect at lo w Reynolds n um b ers.
T ak en together, these prop erties lead to what is termed kinematic r eversibility . A direct
consequence of this on lo comotion at lo w Re is encapsulated in the Sc al lop the or em , coined
b y Ed Purcell in his seminal talk [ 26 ]. It states that a single degree of freedom swimmer,
10
Figure 2.1: An organism with a single degree of freedom lik e a scallop cannot ac hiev e net
lo comotion at lo w Reynolds n um b ers. Image credit: Wikip edia
lik e a scallop, cannot ac hiev e net displacemen t at lo w Reynolds n um b ers. T o ac hiev e net
displacemen t, a lo w Re swimmer m ust ha v e a non-recipro cal swimming gait.
In fact, this b eha vior is quite w ell kno wn in con text of cilia and flagella– they are kno wn to
break symmetry in a b eat cycle b y using a p o w er strok e, when the cilium pushes forw ard lik e
a rigid bar, and a reco v ery strok e, in whic h the cilium mo v es bac k to the starting p osition in
a flexible manner. Moreo v er, cilia also break symmetry globally at a group lev el, b y forming
metac hronal w a v es [ 27 ], to ac hiev e net pumping.
2.4 The Oseen tensor
W e are going to close out this c hapter with a discussion of the impulse resp onse of Stok es
equations. Consider a p oin t force f = (f
x
,f
y
,f
z
) concen trated at p osition x
S
= (x
S
,y
S
,h)
in an infinite fluid domain. The flo w induced b y it can b e found using the Stok es equation
−∇p+η∇
2
u+fδ(x−x
S
) = 0, ∇·u = 0. (2.5)
11
Figure 2.2: Sc hematic represen tation of a force monop ole ab o v e a no-slip w all. Without loss
of generalit y , the p oin t force can b e placed at (0,0,h) .
This can b e solv ed exactly [ 28 ] to obtain the Green’s function of the Stok es equations as
u(x) =S(x,x
S
)·f, (2.6)
wherex is the p oin t of in terest in space, and
S(x,x
S
) =
1
8πη|x−x
S
|
I+
(x−x
S
)⊗(x−x
S
)
|x−x
S
|
2
!
(2.7)
is a second order tensor kno wn as the Ose en tensor. Suc h a v elo cit y solution is called a
stokeslet . It is w orth noting here that higher order singularities can b e obtained b y taking
deriv ativ es of the Stok eslet solution.
2.5 The Blak e tensor
The solution in ( 2.7 ) is v alid in un b ounded domains. Ho w ev er, most phenomena of in terest
in biology o ccur close to b oundaries. In suc h scenarios, enforcing the no-slip condition at
12
the b oundary required adding an image system. Fig. 2.2 sho ws the sc hematic represen tation
of a singularities– a p oin t forcef lo cated at a p oin tx
S
in the domain at a heigh t h ab o v e a
plane w all. The corresp onding image system is placed at the lo cation x
I
, to ensure no-slip
condition at the w all.
The image system can b e constructed in t w o w a ys– the original one due to Blak e [ 29 ]
and a more recen t, elegan t reform ulation due to [ 30 ]. W e briefly review the final result of
b oth approac hes here.
Blak e form ulation: Blak e [ 29 ] constructed the image system b y placing a com bination
of higher order singularities at a reflected p oin t ab out the z-axis, x
I
= (x
S
,y
S
,−h) . The
resulting induced v elo cit y at a p oin t x = (x,y,z) , kno wn as a Blak elet, is giv en b y u(x) =
B(x,x
S
)·f , whereB is the Blak e tensor giv en b y
B(x,x
S
) =S(x,x
S
)−S(x,x
I
)+2h
2
S
D
(x,x
I
)−2hS
SD
(x,x
I
). (2.8)
Here,S
D
andS
SD
are the mo dified source and Stok es doublets, resp ectiv ely:
S
D
(x,x
I
) =
1
8πη
P·∇
x−x
I
|x−x
I
|
,
S
SD
(x,x
I
) =P·∇(S(x,x
I
)·e
z
),
(2.9)
whereP is a 3×3 matrix with diagonal elemen ts [1,1,−1] and all the non-diagonal elemen ts
equal to zero.
13
Gim butas form ulation: Gim butas [ 30 ] constructed the image system b y using a reflected
force in addition to lo cation of the p oin t force. The induced v elo cit y at a p oin t x in the
domain is giv en b y
u(x) =u
A
(x)−u
B
(x)−u
C
(x). (2.10)
where u
A
and u
B
are the v elo cit y fields induced b y the free space stok eslet, and an an ti-
stok eslet placed at the image lo cation, resp ectiv ely:
u
A
(x) =S(x,x
S
)·f, u
B
(x) =S(x,x
I
)·f
I
. (2.11)
where, imp ortan tly ,f
I
= (f
x
,f
y
,−f
z
) . u
C
is a correction term giv en b y
u
C
(x) = z∇ϕ (x)−[0,0,ϕ (x)] (2.12)
and ϕ (x) is a harmonic p oten tial giv en b y
ϕ (x) =−f
z
G
s
(x,x
I
)+hG
D
[f
I
](x,x
I
) (2.13)
with
G
S
(x,x
I
) =
1
4π|x−x
I
|
, (2.14)
G
D
[f
I
](x,x
I
) =∇G
S
(x,x
I
)·f
I
. (2.15)
14
Figure 2.3: Flo w induced b y a unit strength stok eslet in free space (left), a stok eslet +
an tistok eslet (middle) and a stok eslet b ound b y a w all at z = 0 , i.e. a Blakelet (righ t).
The flo wfields are sho wn in the top ro w, and v elo cit y profiles in the z -direction, at (x,y) =
(0.5,0.5) are sho wn in the b ottom ro w.
Flo wfields: The top ro w in Fig. 2.3 sho ws the flo ws induced b y a unit strength stok eslet
(left), a stok eslet + image system consisting only of an an ti-stok eslet (cen ter), and a stok eslet
+ full image system. The b ottom ro w sho ws the u and v v elo cities along z -direction, at
(x,y) = (0.5,0.5) . As indicated b y the middle column, just an an ti-stok eslet is not enough
to satisfy no-slip condition at the w all, and higher order singularities are needed.
Instead of ha ving the force concen trated at a single p oin t in the domain, it is p ossible
to w ork out solutions to the Stok es equations b y spreading out the dirac delta function in to
blob functions; this is referred to as regularized stok eslet; see [ 31 ] for details. Ho w ev er, for
the purp oses of this thesis, w e are going to confine ourselv es to the singular solution only .
15
Chapter 3
Mathematical mo dels of ciliary
sync hronization
I
n this c hapter, w e review the curren t state of the literature related to sync hronization
of cilia. In the pro cess, w e briefly discuss the v arious mathematical descriptions that
ha v e b een used to mo del cilia, and the salien t results that ha v e b een obtained using them.
W eighing on the pros and cons of these approac hes, w e mak e a c hoice of the mo del to use
for the rest of this thesis, and end with a detailed description of the same.
3.1 Literature surv ey
Cilia co ordination has b een analyzed using a wide v ariet y of mo dels. In the most bio-
ph ysically complete mo dels, cilia are represen ted as elastic filamen ts [ 32 , 33 , 34 , 35 , 36 ].
Metac hronal w a v es w ere found to emerge in 1D [ 37 ] and 2D [ 33 ] arra ys using suc h mo dels,
16
but small n um b er of cilia w ere considered due to limitations in computational p o w er. An-
other class of mo dels in v olv e coarse-graining the cilium as a 1D or 2D phase oscillator, kno wn
as r owers and r otors , resp ectiv ely . These mo dels are esp ecially p opular due to their analyt-
ical tractabilit y and exp erimen tal accessibilit y . Metac hronal w a v es ha v e b een observ ed in
c hains of ro w ers [ 38 , 39 ]. W e p oin t the reader to a thorough review [ 40 ] for further details.
F or the purp oses of this thesis, w e use the r otor description of the cilium as it offers
a go o d balance b et w een the simplicit y of a ro w er mo del and the complexit y of a filamen t
mo del. [ 41 ] in tro duced the mo del of cilia as b eads rotating in elliptical tra jectories, and it
w as sho wn that a pair of suc h cilia sync hronize if their tra jectories are tilted with resp ect to a
no-slip b oundary . Phase-lo c king b eha vior is also observ ed if an additional degree of freedom
is in tro duced, in the form of compliance along the normal direction to the tra jectory , ev en
for b eads rotating in circular orbits [ 42 ]. P airwise sync hronization of limit cycle rotors has
also b een observ ed exp erimen tally via activ ely driv en colloidal particles [ 43 , 44 , 45 ]. In sum,
breaking the time-rev ersal symmetry is crucial to sync hronization [ 46 ].
Using this observ ation, [ 47 ] generalized the rotor mo del to arbitrary tra jectories and
sho w ed that a non-constan t phase-dep enden t driving force can driv e h ydro dynamic sync hro-
nization. In particular, sync hronization in 1-d arra y using this minimal mo del giv es rise
to metac hronal w a v es from random initial conditions with b oth free and p erio dic b ound-
aries [ 42 , 48 ], and carp ets of cilia w ere found to exhibit either full phase lo c king, completely
disordered or spiral w a v e states [ 49 ]. Recen tly , Meng et al. sho w ed that the stabilit y of
metac hronal w a v es is determined b y the form of the tangen tial driving force for eac h ro-
tor [ 50 ], and [ 51 ] demonstrated the stabilit y of metac hronal w a v es using a more complicated
phase oscillator description of the cilium.
17
Figure 3.1: Simplified represen tation of b eating cilium: A oscillating filamen t is appro xi-
mated as a rigid b ead mo ving in a circular tra jectory .
3.2 Discrete cilia mo del
In the rotor mo del of cilia, eac h cilium is mo deled as a rigid b ead of radius a rotating along
a circular tra jectory of radius b , kno wn as the r otor mo del [ 46 , 41 , 42 , 47 , 50 ]. Fig. 3.1
illustrates the concept pictorially . A t this stage, it is w orth p oin ting out this is not just a
spheric al c ow appro ximation of ph ysics. It w as sho wn using exp erimen ts [ 9 ] that lo w er order
mo dels using singularities do offer a go o d represen tation of the flo ws induced. A result from
their w ork has b een repro duced in figure 3.2 , under the fair use p olicy .
W e consider an infinite t w o-dimensional arra y of suc h cilia where the cen ter x of the rotor
is lo cated on a square lattice of spacing d in a plane at a distance h from a flat substrate
(see Fig. 3.3 A). W e parameterize the phase of eac h cilium b y the angle θ along its circular
tra jectory , and in tro duce the (x,y,z) co o rdinates where (x,y) span the ciliary plane and z
describ es distance from the substrate. W e also in tro duce the normal and tangen tial unit
v ectors n = (−sinθ,cosθ,0) and t = (cosθ,sinθ,0) , suc h that the b ead lo cation can b e
describ ed b y its planar p osition v ector r =x +bn and the b ead v elo cit y b y ˙ r = b
˙
θt . Eac h
18
Figure 3.2: (a) Exp erimen tally measured flo wfield around a c hlam ydomonas, (b) Flo wfield
induced b y a 3-stok eslet mo del. This sho ws that lo w er order mo dels are go o d appro ximations
to the flo wfields seen in the real w orld. This image has b een repro duced from [ 9 ], repro duced
here under the F air use p olicy .
Figure 3.3: (A) Sc hematic represen tation of the ciliary carp et arranged in a doubly p erio dic
square latti ce, ab o v e a no-slip b oundary (sho wn primary b o x only). (B) First (blac k) and
second (blue) force harmonics of the driving force.
cilium is driv en indep enden tly b y a tangen tial force F(θ) = Ft , whose magnitude F(θ) is
phase-dep enden t and strictly p ositiv e; F(θ) can b e represen ted via its harmonics as
F(θ) = F
0
(1+
∞
X
i=1
α
i
cosiθ +β
i
siniθ), (3.1)
where α
i
and β
i
are co efficien ts of the F ourier series expansion of F in terms of θ . T w o
sample driving force profile s, corresp onding to the first and second harmonics, are depicted
in Fig. 3.3 B.
19
The ciliary plane is submerged in an semi-infinite three-dimensional (3D) fluid domain,
with a flat b ounding substrate at z = 0 . Considering t ypical cilia length and time scales, vis-
cous drag forces are dominan t and the fluid motion is b est mo deled using the incompressible
Stok es equations ( 2.4 ), rewritten here for con v enience.
η∆ u =−∇p, ∇·u = 0. (3.2)
Here,u(x,y,z)≡ (u,v,w) is the 3D fluid v elo cit y , p(x,y,z) the pressure, andη the viscosit y .
F o r ease of notation later (in c hapter 5 ), it is con v enien t at this stage to in tro duce the 2D
fluid v elo cit y v(x,y)≡ (u,v) ev aluated at z = h in the ciliary plane.
Balance of forces on a c ilium j in the doubly-infinite ciliary lattice dictates that
F
j
+N
j
−ζ
b
˙
θ
j
t
j
−v(r
j
)
= 0, (3.3)
where the subscript j is used to denote the quan tit y p ertaining to cilium j . P articularly , the
first term in ( 3.3 ) is the activ e tangen tial force F
j
= F
j
t
j
= F(θ
j
)t
j
, and the second term
N
j
= N
j
n
j
is a normal constrain t force that guaran tees the b ead remains on the desired
circular tra jectory . The last term denotes the drag force (with constan t drag co efficien t
ζ = 6πηa ) on the j
th
cilium and accoun ts for h ydro dynamic coupling among all cilia. F or
a single cilium, the constrain t force is iden tically zero and the in trinsic angular sp eed
˙
θ =
Ω( θ) = F(θ)/ζb is phase-dep enden t. F or the h ydro dynamically-coupled lattice of cilia, the
20
constrain t force is giv en b y
N
j
=−ζ(n
j
⊗n
j
)·v(r
j
) =−ζ (n
j
·v(r
j
))n
j
, (3.4)
where ⊗ denotes the outer tensor pro duct, and the time-ev olution of the cilia phase is
go v erned b y a system of coupled differen tial equations
˙
θ
j
= Ω j
+
1
b
t
j
·v(r
j
). (3.5)
Here, Ω j
= Ω( θ
j
) = F(θ
j
)/ζb refers to the in trinsic angular sp eed of cilium j .
21
Chapter 4
A simple algorithm for large scale
sim ulations of ciliary carp ets
M
any in teresting phenomena at lo w Reynolds n um b ers in v olv e micro-swimmers in ter-
acting with surfaces, b oth in natural and lab oratory settings [ 52 ]. A few examples
include mo v emen t of sp ermatozoa [ 53 ], formation of bacterial colonies [ 54 ] and motion of
Bro wnian particles. These organisms ha v e cilia or flagella to assist them to self-prop el in
fluid en vironmen ts. Theoretical mo dels used to study in teracting en tities often in v olv e mo d-
eling either the ciliated organism itself, or the individual cilia themselv es, as p oin t forces in
confined fluid domains [ 55 ]. While the sync hronization b et w een pairs of cilia is reasonably
w ell understo o d using theoretical mo dels, inquiries in to the sync hronization among large
n um b e r of cilia ha v e b een limited to a few h undred to a thousand cilia o wing to computa-
tional limitations. T o put this in p ersp ectiv e, this n um b er is ab out an order of magnitude
smaller than the n um b er of cilia on microb es lik e V olvox or p ar ame cium [ 56 , 57 ].
In this c hapter, w e are in terested in sim ulating the collectiv e dynamics in large arra ys
22
of cilia. Sp ecifically , w e use the rotor mo del of cilia describ ed in the previous c hapter. A
p opular and con v enien t approac h, to sim ulate largeN -b o dy systems, is to use doubly p erio dic
b ound ary conditions. Ho w ev er this presen ts computational c hallenges for an y reasonably
large n um b er of cilia, b ecause the presence of infinite n um b er of p erio dic images coupled
with the slo w deca y of h ydro dynamic in teraction leads to the conditionally con v ergen t nature
of the summation to compute fluid v elo cities. Previous studies ha v e addressed the flo w
generated b y force monop oles in doubly p erio dic domains [ 58 , 59 , 10 ]. Sp ecifically , [ 58 , 59 ]
used Ew ald decomp osition to split the summation in to t w o parts– one ev aluated in the
ph ysical space, and another in the F ourier space. This ho w ev er, requires a careful c hoice
of parameters for optimal p erformance. In [ 10 ], the authors reform ulated the summation in
suc h a w a y as to ensure zero net force in the primary domain, whic h then allo ws for using the
F a st Multip ole Metho d (FMM). While these metho ds offer high accuracy and computational
efficiency , they can b e tec hnically more in v olv ed and c hallenging to implemen t in co de.
In this c hapter, w e con tribute to this b o dy of literature b y devising a fast, accessible
algorithm that is easy to implemen t from scratc h. This c hapter is structured as follo ws: In
section 4.1 , w e discuss the c hallenges asso ciated with p erforming sim ulations with a large
n um b e r of cilia, and outline an algorithm to o v ercome them. In section 4.2 , w e presen t the
results of our n umerical sim ulations whic h sho w noisy tra v elling w a v e-lik e patterns emerging
from random initial conditions. Finally , w e conclude with a summary of the con tributions
in section 4.3 .
23
Figure 4.1: Sc hematic represen tation of the computational domain, consisting of the funda-
men tal domain (red), images (grey) and in teraction stencil (blue). The zo omed in view of
the in teraction stencil sho ws the fo cal cilium at the cen ter of the lattice (dark blue) and the
neigh b oring cilia (red). W e compute the induced v elo cit y at the fo cal cilium b y considering
the influence of all neigh b ors within the in teraction stencil.
4.1 Algorithm for ev aluating the Doubly-p erio dic Blak e
T ensor
T o close the mo del in ( 3.5 ), w e need to ev aluate the fluid v elo cit y v at cilium j induced b y
the infinite planar arra y of cilia. T o this end, w e first consider that eac h cilium is a force
monop ole of strength F
j
+N
j
and mak e use of the singularit y solution near a rigid w all
represen ted b y the Blak e tensor B ; see [ 29 , 60 , 30 ]. W e find that the con tribution of the
constrain t forceN
j
is negligible suc h that, for a finite n um b er of cilia, u
j
=
P
i
B(r
j
−r
i
)·F
i
.
Then, to accoun t for the infinite ciliary lattice, w e represen t the lattice b y a finite arra y of
n
c
×n
c
cilia in a doubly-p erio dic square domain of sizeL×L , suc h that the flo w generated b y
a ciliumi is due to its o wn motion in the fundamen tal domain as w ell as to its doubly-infinite
set of planar images; see left panel of figure 4.1 . W e get thatu
j
=
P
nc×nc
i
B
lattice
(r
j
−r
i
)·F
i
,
where B
lattice
·F
i
= lim
M→∞
P
M
−M
P
M
−M
B·F
i
, where M is a domain index. This double
sum is only conditionally-con v ergen t giv en that, due to the presence of the flat w all at
z = 0 , the tensor B deca ys as 1/r
2
in the far field, r b eing the distance from the force
24
singularit y [ 29 ]. Numerical ev aluation of this doubly-infinite sum requires emplo ying Ew ald
sum tec hniques [ 61 , 62 , 63 , 59 , 58 ]. This in v olv es, as alluded to previously , splitting the
sum in to t w o parts– one part to b e ev aluated in the real space and the other in the F ourier
space. F or maxim um computationa l efficiency , the splitting parameter has to b e c hosen
judiciously . Another approac h to ev aluate the doubly-infinite sum is to use the F ast Multip ole
Metho d [ 10 ]. These metho ds are in v olv e tec hnical c hallenges ev en is used lik e blac k b o xes,
as installing sp ecialized libraries can b e lab orious and not as p ortable as one w ould lik e.
Instead, w e prop ose to obtain a go o d appro ximation of the v elo cit y field without resorting
to a full-fledged doubly-p erio dic sum, b y using the ph ysics of the problem itself.
The fi rst k ey feature to note is the screening of h ydro dynamic in teraction among cilia
due to the no-slip b oundary , whic h ensures the con v ergence of the conditionally-con v ergen t
sum. That is, the flo wfield v generated b y all cilia at the lo cation of a “fo cal cilium”, can
b e represen ted as a finite summation. This finite summation consists of con tributions from
a finite, y et large, n um b er of nearest neigh b ors, whic h are lo cated within the “in teraction
stencil” around it. This approac h has previously b een used in, for example, [ 64 ]. The
in teraction stencil consists of a n um b er of concen tric shells, n
s
, around the fo cal cilium, as
illustrated in the righ t panel of figure 4.1 . n
s
shells corresp onds to a total of (2n
s
+1)
2
cilia
inside the stencil, including the fo cal cilium. Cilia lo cated outside the outermost shell are
assumed to ha v e negligible effect on the fo cal cilium of the stencil. Note that while using this
approac h, the set of in teracting neigh b ors need not necessarily con tain cilia from the image
b o x – while the cilia closer to the edge of the primary b o x in teract with p erio dic images,
the ones closer to the cen ter ma y not. In orde r to quan tify the error in tro duced b y finite
neigh b or truncation, w e compare the truncated v elo cit y to the result obtained from the F ast
25
Multip ole Metho d [ 10 ], as sho wn in figure 4.2 .
In figure 4.2 A, w e sho w the cosine field of a 151× 151 lattice of cilia for three sample
states– isotropic, syn thetic and noisy w a v e states. In figure 4.2 B, the blac k line sho ws the
truncated v elo cit y at the fo cal cilium, i.e., the cilium at the cen ter of the domain, as a
function of increasing n um b er of shells, n
s
. The gra y line is obtained b y ev aluating the same
v elo cit y using FMM [ 10 ], whic h tak es in to accoun t the infinite images, and is hence in v arian t
with n
s
. W e then calculate the relativ e p ercen tage difference b et w een the truncated and
FMM solutions as follo ws:
Relativ e difference =
∥v
ns+1
−v
ns
∥
∥v
ns+1
+v
ns
∥
×100. (4.1)
As sho wn in figure 4.2 C, the truncated solution appro ximates the true solution for finite
n um b e r of shells to a reasonable error.
The second feature of use in our n umerical algorithm is the regular arrangemen t of all
cilia on a square lattice, and their constrained motion along fixed circular tra jectories. W e
pre-compute the effects of all cilia in the stencil on the fo cal cilium in a pairwise fashion.
W e first discretize eac h circular tra jectory in to n
θ
angular lo cations as sho wn in figure 4.3 B.
The v elo cit y at the fo cal cilium j at angular lo cation p , induced due to a neigh b or cilium i
at angular lo cation q is computed using the Blak e tensor, and stored in a lo ok-up table for
future reference. The en tries in the lo ok-up table are 2-dimensional v elo cit y v ectors indexed
using the relativ e p osition of the neigh b oring cilium with resp ect to the fo cal cilium on the
2D lattice, and their discrete angular lo cations along the tra jectories (p,q )–v
p,q
j←i
. The en tries
corresp onding to self-in teraction, i.e. v
p,q
j←i
, are all set to zero, resulting in a lo okup table
26
Figure 4.2: Con v ergence of fluid v elo cit y . A : Cosine fields of a system of 151× 151
cilia, in three sample states. B : Flo w sp eed (blac k) induced at the cen ter of the lattice
approac hes the result obtained using the F ast Multip ole Metho d [ 10 ] (gra y) for increasing
shells of in teracting neigh b ors. C : The relativ e p ercen tage difference in the induced flo w
sp eed, as compared to the FMM result (gra y) and due to increasing n um b er of shells (blac k)
as sho wn in ( 4.1 ) .
27
Figure 4.3: (A) Induced v elo cities at the lo cation of the fo cal cilium (blue), b y ev ery other
cilium inside the in teraction stencil (red) are pre-computed and stored. (B) The tra jectory
of eac h cilium is discretized in ton
θ
angular lo cations. v
p,q
i←j
is the induced v elo cit y at angular
lo c ation p of the fo cal cilium i , due a neigh b oring cilium, j , at angular lo cation q , and is
computed using the Blak e tensor. (C) Sc hematic represen tation of the v elo cit y lo okup table.
Eac h e n try is a 2-dimensional v elo cit y v ector p opulated as p er the description in (B).
with (2n
s
+ 1)
2
×n
2
θ
× 2 en tries. This brings do wn the cost of the algorithm from O(n
4
c
)
to O(n
2
c
n
2
s
) , whic h is no w linear in terms of n um b er of cilia n
2
c
. A step-b y-step alg orithmic
description to construct the v elo cit y lo okup table is giv en in 1 . A dditional sp eed gains are
ac hiev ed b y v ectorizing sections of the co de, wherev er p ossible. It is p ossible to use symmetry
of the square lattice to pre-compute only a subset of the cilia within the stencil explicitly ,
namely , cilia whose cen ters lie within the mark ed o ctan t of the in teraction stencil. But giv en
that the computational cost of building the lo okup table is extremely lo w compared to total
sim ulation time, w e forgo this tric k.
28
Algorithm 1 Building the v elo cit y lo okup table
Input: n
θ
: n um b er of discrete angular p oin ts
Input: F(θ) : Driving force
Input: Stencil: (2n
s
+1)×(2n
s
+1) lattice of cilia
Input: j : index of fo cal cilium
Output: Induced v elo cit yv
p,q
j←i
1: for i = 1 : (2n
s
+1)
2
do
2: θ
p
i
, p = 1,2,...n
θ
: discretized lo cation of cilium i in to n
θ
angular lo cations—
3: end for
4: for i = 1 : (2n
s
+1)
2
do
5: if i̸= j then
6: for p = 1 : n
θ
do
7: Calculater
p
j
: p osition v ector of fo cal cilium j , at angular lo cation p
8: for q = 1 : n
θ
do
9: Calculater
q
i
: p osition v ector of neigh b oring cilium i , at angular lo cation q .
10: Calculate stok eslet strengthF
q
i
.
11: Compute induced v elo cit yv
p,q
j←i
=B(r
q
i
−r
p
j
)·F
q
i
12: end for
13: end for
14: else
15: Set self-induced v elo cit y to zero: v
j←i
=0∀p,q .
16: end if
17: end for
29
Figure 4.4: The induced v elo cit y on cilium j , due to a neigh b oring cilium i is obtained b y
lo oking up the four en tries from the table, and linearly in terp olating b et w een them.
During run time while in tegrating equation ( 3.5 ), the pre-computed v alues of induced
v elo cities are simply lo ok ed-up from the tabl e based on the relativ e p ositions of cilia. Re-
mem b er that the relativ e p osition of a giv en pair of cilia tak es in to accoun t their relativ e
p osition on the 2D grid, as w ell as their discrete angular lo cations along their tra jectories.
A second order accurate in terp olation sc heme is used to obtain the induced v elo cit y on an y
cilium; see algorithm 2 and figure 4.4 for a pictorial description. This sp eeds up our program
b y sidestepping the ev aluation of the Blak e tensor at eac h time step, whic h is computation-
ally more exp ensiv e than lo okup. This further bring the computational cost b y reducing the
constan t factor that do es not app ear in theO notation.
Finally , w e v alidate our algorithm b y p erforming c hec ks on con v ergence and s p eed gains.
W e run sim ulations using a 151×151 system of cilia with h/d = 0.5 , b/d = 0.2 , a/d = 0.05 ,
n
θ
= 36 , and v ary the n um b er of shells in p o w ers of 2, i.e. n
s
= 1,2,4,...,32 . W e n umerically
in tegrate equation ( 3.5 ) for 5 time units, using the R unge-Kutta metho d with a fixed time-
step of dt = 0.1 . W e record the v ector con taining the angular phases of all 151
2
= 22801
30
cilia, denoted b y θ , at the end of the sim ulation, and c ompute the relativ e error as
Relativ e error = (∥θ(2n
s
)−θ(n
s
)∥)/∥θ(2n
s
)∥. (4.2)
A dditionally , w e also record the w all clo c k time for these sim ulations, and the results are
summarized in figure 4.5 . W e find that our sim ulation has an order of con v ergence of ≈ 1.66 ,
and a time complexit y of ≈ 1.95 with resp ect to the n um b er of shells. Note that this is
close to the predicted quadratic scaling. W e rep eat these sim ulations without discretizing
the circular tra jectories, i.e. b y discarding the lo okup approac h and ev aluating the Blak e
tensor at eac h time step instead. As indicated b y figure 4.5 , w e find that while the order of
con v ergence of b oth approac hes remains comparable, the lo okup approac h results in sp eed
gains of close to an order of magnitude. Again, note that the scaling of time with resp ect
to the n um b er of shells is close to the predicted v alue of 2 . Based on the results s ho wn in
Figs. 4.2 and 4.5 , w e c ho ose the v alue of n
s
= 10 for all the results rep orted in the rest of
this thesis, unless noted otherwise.
31
Algorithm 2 Computing induced v elo cities using lo okup
Input: θ
i
, i = 1,2,...N : angular phase of eac h cilium
Input: N(i), i = 1,2,...N : list of in teracting neigh b ors for eac h cilium i
Output: v
i
, i = 1,2,...N : Induced v elo cit y on eac h cilium i , in the lattice
1: for eac h cilium i = 1 : N do
2: Discr etize : Using θ
i
, find the t w o nearest a ngular grid p oin ts θ
p
i
, θ
p+1
i
3: Inter action stencil : Find neigh b ors inside the in teraction stencil cen tered on cilium i –
N(i)
4: for all j∈N(i) do
5: Discr etize : Using θ
j
, find the t w o nearest angular grid p oin ts θ
q
j
, θ
q+1
j
6: L o okup : v
p,q
i←j
,v
p,q+1
i←j
,v
p+1,q
i←j
,v
p+1,q+1
i←j
7: Line ar interp olation : to ev aluatev
i←j
8: end for
9: end for
10: v
i
=
P
j
v
i←j
: total v elo cit y induced on cilium i
4.2 Results
4.2.1 Emergence of metac hronal w a v es
Using the algorithm describ ed in the previous section, w e in tegrate equation ( 3.5 ) forw ard in
time for 100 time units and study the emergen t dynamics in a system of 22801 cilia, arranged
on a 151×151 square lattice, starting from an isotropic state. Our main finding is that the
32
Figure 4.5: Log scale plots of (A) relativ e error, and (B) computational time as a function
of n um b er of shells. Solid and dashed lines indicate sim ulations with and without lo okup,
resp ectiv ely .
isotropic state is unstable and b ey ond an initial transien t phase, the cilia self-organize to
giv e rise to w a v e-lik e patterns, suc h as the ones sho wn in figure 4.6 .
W e find that these spatial w a v es are robust to initial conditions, and are dep enden t
largely on the activ e force profile of an individual cilium. F or instance, first and second
force harmonics lead to w a v e-lik e patterns in long-time, while others do not. In terestingly ,
this is consisten t with the exp erimen tally measured force profile of the flagellar b eat [ 65 ],
thereb y , reemphasizing the fact that the rotor mo del of cilia do es indeed capture the es-
sen tial ph ysical mec hanism of metac hronal sync hronization. W e find that the phase-defects
from the transien t phase in the sim ulation p ersist to long times, similar to those observ ed
in [ 33 ]. Curiously , they o ccur despite the absence of an explicit noise term in our go v erning
equation ( 3.5 ) and w e hop e to in v estigate this in future studies.
Apart from the force harmonic, w e find that the cilium heigh t h/d also affects the nature
of the emergen t w a v es. Sp ecifically , w e find that p o c k ets of larger w a v elength start to emerge
for longer cilia lengths, as sho wn in figure 4.7 .
33
Figure 4.6: W e p erform sim ulations using a carp et of 151×151 cilia, arranged in a square
lattice with constan t spacing d = 1 , at a heigh t h = 0.5 ab o v e a no-slip w all, for 100 time
units. W a v e-lik e spatial patterns emerge when the individual cilia are driv en b y the first
(left) and second (righ t) force harmonics. Sim ulation parameters used: Isotropic initial
state, b/d = 0.2, a/d = 0.05 , 10 shells of in teracting neigh b ors, and using n
θ
= 36 discrete
p oin ts along the tra jectory of eac h cilium.
4.2.2 Flo w pumping
Metac hronal w a v es in ciliary carp ets serv e to enhance fluid transp ort at lo w Reynolds n um-
b ers [ 66 , 32 , 64 ]. It is w orth making the distinction that in [ 66 , 32 , 64 ] metac hronal sync hron y
w as imp osed on the ciliary carp et, as opp osed to the presen t w ork where the phase co or-
dination in to w a v e patterns is an emergen t phenomenon. Our goal is to quan tify the fluid
transp ort in our system as the w a v e patterns emerge.
T o measure fluid pumping in our system, w e recall that a force monop ole of strength F ,
orien ted parallel to a no-slip w all at heigh t h , pushes fluids at a rate giv en b y [ 67 , 68 ]
q = hF/πη. (4.3)
34
Figure 4.7: Cosine field obtained from sim ulations at t = 100 units, for increasing v alues of
cilium heigh t. P o c k ets of larger w a v elength emerge for longer cilia. Sim ulations parameters
used: Isotropic initial state, second force harmonic, b/d = 0.2, a/d = 0.05 and 10 shells of
in teracting neigh b or s, and using n
θ
= 36 discrete p oin ts along the tra jectory of eac h cilium.
The cycle a v eraged flux due to a single isolated cilium is
q =
h
πη
Z
T
0
Fdt. (4.4)
Using the transformation dt = dθ/
˙
θ , and noting that for a single, isolated cilium
˙
θ = Ω( θ) =
F(θ)/bζ , w e get
q =
hbζ
πη
Z
2π
0
tdθ = 0. (4.5)
Hence, a single cilium cannot ha v e an y net pumping o v er a cycle.
Using the linearit y of Stok es flo w, it is straigh tforw ard to see that the total instan taneous
flux due to the ciliary carp et is the sum of all cilia con tributions, i.e.
P
j
q
j
(t) , where j is
the cilia index. Th us, the cum ulativ e v olume of fluid pump ed b y the carp et un til a time t ,
is
Q =
h
πη
Z
X
j
F
j
(t)dt. (4.6)
In Figure 4.8 , w e sho w the cum ulativ e pumping due to the ciliary carp et in the direction
parallel to the w a v e propagation direction. Results for the fir st force harmonic are sho wn
35
Figure 4.8: Cum ulativ e directional flux induced b y the 151× 151 carp et of cilia, o v er 100
units of time, for the first (blac k) and second (blue) force harmonics.
in blac k and for the second force harmonic in blue. This directional flux w as obtained b y
pro jecting the flux in equation ( 4.6 ) on to the direction of w a v e tra v ersal. W e discuss the
tec hnical details of ho w to obtain the w a v e tra v ersal direction in a future c hapter (Chapter 6).
Numerically , w e find that the cilia driv en b y the first force harmonic pump fluid m uc h more
effectiv ely as compared to driving b y the second force harmonic. W e b eliev e that the reason
for this lies in the nature of w a v es that emerge from eac h t yp e of forcing; see accompan ying
supplemen tary mo vies in [ 69 ]. The first force harmonic has a unique direction of the effectiv e
strok e since it has exactly one maxima and one minima p er cycle. On the other hand, suc h a
direction is hard to define unam biguously for the second force harmonic, due to the presence
of m ultiple maxima and minima. This results in some patc hes of cilia pumping fluid in the
opp osite direction to the other patc hes and, in turn, resulting in a reduced net flo w.
4.3 Conclusions
In this c hapter, w e devised a simple, y et p o w erful, algorithm that enables large scale sim u-
lations to shed ligh t on the emergen t dynamics among h ydro dynamically coupled cilia. Our
algorithm mak es use of t w o k ey ideas to mak e the ev aluation of the p erio dic Blak e tensor
36
computationally tractable: (1) the h ydro dynamic screening from the w all, to truncate the
conditionally con v ergen t summation, and (2) the circular nature of the cilia tra jectories.
Using these sim ulations, w e demonstrated the emergence of symplectic metac hronal w a v es
in the ciliary carp et. Our algorithm, while less accurate than more sophisticated metho ds
in [ 58 , 59 , 10 ], has the adv an tage of b eing v ery simple to program from scratc h, and is ideally
suited for quic k-and-dirt y implemen tation to gain insigh ts on p oten tially in teresting ph ysics
in problems.
There are some assumptions whic h can b e readily relaxed to study a ric her parameter
space. One suc h assumption is that our algorithm relies on all cilia in the domain b eing
iden tical, and arranged on a regular lattice. This can b e remedied b y using an analytical
approac h w e dev elop in the next c hapter, to ev aluate induced v elo cities. This w ould in v olv e
taking a in v erse F ourier transform to obtain the in-plane v elo cit y in ph ysical space, and
subsequen tly in terp olating on to the lo cation of e ac h cilia. This w ould allo w us to prob e a
ric her parameter space including the lattice geometry , cilia placed on an irregular lattice
(lik e in [ 70 ]) and non-circular tra jectories. Another, more v ersatile approac h, w ould b e to
incorp orate the FMM to ol dev elop ed in [ 10 ] in to our sim ulation framew ork. This w ould ex-
pand the p ossible parameter space to include non-planar b eating of cilia, i.e., with tra jectory
tilted with resp ect to the w all, and to also test the effect of compliance along the normal
direction, as has b een realized exp erimen tally with optical t w eezers [ 71 ].
37
Chapter 5
A Con tin uum mo del of ciliary carp ets
C
iliar y carp ets found in nature are m uc h larger than the sim ulations discussed in
the previous c hapter. F or example, in higher organisms lik e rats, ciliary carp ets are
found to ha v e o v er 10
6
individual cilia [ 70 ]. This n um b er increases b y a few more orders
of magnitude in the case of h umans. As a consequence, hop es of understanding ciliary
co ordination in suc h large systems will remain b ey ond the realm of computational tractabilit y
in the near future, ev en with the use of more sophisticated computational metho ds. W e
seek to remedy this issue in the curren t c hapter. W e use the rotor description of cilia and
systematically coarse grain their b eha vior in to con tin uous fields.
W e presen t a theoretical framew ork for understanding the ph ysical conditions leading to
the co ordination of infinitely-man y , indep enden tly b eating cilia, immersed in a 3D fluid and
arranged on a 2D arra y ab o v e a substrate. F or this, w e deriv e collectiv e equations of motion
for the dynamics of h ydro dynamically coupled cilia ab o v e a no-slip w all, b y adapting the
argumen ts dev elop ed for bulk systems in [ 72 ]. W e then analyze the stabilit y of isotropic and
fully-sync hronized steady states of these equations and sho w that they are unstable to small
38
Figure 5.1: Ciliary La y er: A doubly-infinite arra y of cilia, wherein eac h cilium mo deled as a
b ead rotating on a fixed circular tra jectory under the influence of a p ositiv e phase-dep enden t
force F(θ) is appro ximated b y a force d ensit y la y er in the limit of infinite n um b er of cilia.
p erturbations. W e further sho w ho w the gro wth rate of these instabilities dep end on the
phase dep endence of the forcing at the lev el of an individual cilium.
This c hapter is structured as follo ws: in section 5.1 , w e deriv e the h ydro dynamic k ernel
that describ es the coupling b et w een rotors and the fluid flo ws that they generate collec-
tiv ely . In section 5.2 w e in tro duce the con tin uous fields that describ e sync hronization on
h ydro dynamic scales and coarse grain our rotor based mo del to obtain a con tin uum theory
for carp ets of b eating cilia. Then in section 5.3 , w e use the con tin uum equations to iden tify
the spatially homogeneous steady states of the ciliary carp et and c haracterize their stabilit y .
W e discuss and con textualize our results in section 5.4 .
5.1 Calculating cilia-induced fluid flo ws
The force densit y that the cilia exert on the fluid at an instan t of time t , is giv en b y
f(x,t) =
X
j
(F
j
+N
j
)δ(x−r
j
). (5.1)
39
where δ(x−x
j
) has units of in v erse area. Note that, for the remainder of this thesis, w e
assume that the dep endence of f on space x , and time t , is implicit, and suppress it for
notational con v enience. Here, all forces are exerted within the thin 2d la y er at a distance h
ab o v e the surface on whic h cilia are em b edded, see Fig. 5.1 . T o solv e for the fluid v elo cit y
field u that the force densit y f induces, w e adapt the metho d in [ 61 , 73 , 74 ] concerning the
flo ws generated b y a 2d force la y er in an infinite fluid to the geometry of our ciliary system,
as discussed next.
Giv en that the force densit y f(x) in ( 5.1 ) is confined to a la y er at z = h , w e partition
space in to t w o domains ab o v e and b elo w the la y er, whic h w e denote b y the ‘+ ’ and ‘− ’
sup erscripts resp ectiv ely . Rewriting the incompressible Stok es equations in eac h domain
giv es
η∇
2
u
+
−∇p
+
= 0, ∇·u
+
= 0, z∈ (h,∞),
η∇
2
u
−
−∇p
−
= 0, ∇·u
−
= 0, z∈ [0,h).
(5.2)
Due to the ciliary force la y er atz = h , the stress field σ =−pI+η
∇u+(∇u)
T
exp eriences
a jum p
Jσ·e
z
K
z=h
= σ·e
z
|
z=h
+
− σ·e
z
|
z=h
−
=f(x) =
X
j
(F
j
+N
j
)δ(x−x
j
). (5.3)
In addition, w e imp ose a no-slip condition at the substrate
u
−
(x,y,z = 0) = 0, (5.4)
40
and v elo cit y con tin uit y at the ciliary la y er
u
+
(x,y,z
+
→h) =u
−
(x,y,z
−
→h). (5.5)
W e also require the pressure and v elo cit y to b e b ounded in the far field
lim
z→∞
p
+
= 0, lim
z→∞
u
+
= constan t. (5.6)
This c hoice captures the ph ysics of v anishing w all shear stress at z→∞ , suc h as one w ould
exp ect, for instance, at a w ater air in terface far ab o v e the ciliary la y er.
T o solv e ( 5.3 )–( 5.6 ), w e tak e the F ourier transform with resp ect to the in-plane co ordi-
nates x = (x,y) [ 73 , 74 ]. Using k = (k
x
,k
y
) for the in-plane w a v en um b er co ordinates, the
momen tum equations in ( 5.2 ) b ecome
ik· ˆ v
±
+∂
z
ˆ w
±
= 0,
η(∂
zz
−k
2
)ˆ v
±
−ikˆ p
±
= 0,
η(∂
zz
−k
2
)ˆ w
±
−∂
z
ˆ p
±
= 0,
(5.7)
where hatted quan tities denote the in-plane F ourier transforms of their ph ysical space equiv-
alen ts, and k = ∥k∥ is the magnitude of the 2D w a v e n um b er v ector. Since the fluid is
incompressible, the pressure ob eys the P oisson equation. In F ourier space, this is giv en b y
(∂
zz
−k
2
)ˆ p
±
= 0, (5.8)
41
whic h is solv ed b y
ˆ p
±
= A
±
e
−k(z−h)
+B
±
e
k(z−h)
, (5.9)
where A
±
= A
±
(k) and B
±
= B
±
(k) are constan ts of in tegration to b e determined from the
b oundary conditions. Substituting ( 5.9 ) in to ( 5.7 ) and solving forv
±
and w
±
giv es
ˆ v
±
=
O
±
−
ikA
±
2ηk
(z−h)
!
e
−k(z−h)
+
Q
±
+
ikB
±
2ηk
(z−h)
!
e
k(z−h)
,
ˆ w
±
=
C
±
+
A
±
2η
(z−h)
!
e
−k(z−h)
+
D
±
+
B
±
2η
(z−h)
!
e
k(z−h)
.
(5.10)
Again,C
±
= C
±
(k), D
±
= D(k), O
±
=O
±
(k) andQ
±
=Q
±
(k) are to b e determined from
b oundary conditions as detailed in App endix A . Moreo v er, w e plot the flo wfields obtained
using this approac h, and compare it with the flo wfields obtained b y direct summation using
the Blak e tensor, in the app endix C .
Finally , using the b oundary condit ion ( 5.3 )–( 5.6 ) w e finally arriv e at
ˆ v =
ˆ
K·
ˆ
f, (5.11)
with
ˆ
K(k
x
,k
y
) enco ding the Stok es k ernel asso ciated with the infinite force densit y la y er at
z = h ,
ˆ
K =
1
4ηk
3
2k
2
(1−e
−2kh
)I+((1+2kh−2k
2
h
2
)e
−2kh
−1)k⊗k
, k̸= 0,
h
η
I, k = 0.
(5.12)
Next, w e use the equatio ns of motion ( 3.5 ), the force that cilia exert on the fluid ( 5.1 ),
42
and the resulting flo w field ( 5.12 ) as the basis from whic h w e will deriv e a con tin uum theory
for ciliated surfaces.
5.2 Order parameter fields for sync hronization and their
equations of motion
W e seek to dev elop a con tin uum theory for ciliary co ordination. T o do this w e first in tro duce
the order parameter fields that describ e sync hronization in the con tin uous limit and deriv e
their equations of motion. W e then exp ose ho w fluid flo ws can b e deriv ed from these order
parameters.
F ollo wing [ 72 ], w e prop ose to monitor the sync hronization order parameter fields
Y
n
(x,t) =
1
ρ
c
X
j
e
inϕ
j
δ(x−x
j
), (5.13)
where ϕ j
is the phase of the j -th cilium, and ρ
c
corresp onds to the cilia densit y , whic h is a
constan t b y construction. By definition, Y
1
is the Kuramoto order field, and Y
2
the nematic
order field. Again, as done in ( 5.1 ), w e suppressed the dep endence of Y
n
on space and time
in our notation. The spatial a v erage of Y
1
o v er the doubly-p erio dic domain leads to the
Kuramoto order parameter P = |⟨Y
1
⟩| [ 75 ]; v alues of P near zero indicate phase disorder
while v alues near one corresp ond to phase sync hron y . As illustrated in Fig. 5.2 , for eac h
v alue of n , the quan tit y |⟨Y
n
⟩| = 1 captures a differen t kind of phase order. F or instance,
|⟨Y
1
⟩| = 1 captures in-phase sync hron y ,|⟨Y
2
⟩| = 1 captures an ti-phase sync hron y ,|⟨Y
3
⟩| = 1
captures a Mercedes star-lik e arrangemen t of rotors around the phase circle, and so on.
43
Figure 5.2: In a system of N phase oscillators (cilia in our case), eac h momen t of the order
parameter Y
n
, captures a differen t kind of phase order around the circle. The oscillators
ha v e b een sho wn on the same circle for ease of illustration.
A priori the c hoice of the phase ϕ j
is not unique. F or example, it can b e θ
j
or an y
smo oth non-decreasing function of θ
j
. W e fix this freedom b y requiring that the phase of an
unp e rturb ed cilium, that is a cilium in a quiescen t bac kground fluid with v = 0 , ob eys
˙
ϕ j
= Ω 0
, (5.14)
where Ω 0
= 2π/T
0
is a constan t. Here, T
0
=
R
2π
0
d
¯
θ/Ω(
¯
θ) is the b eat p erio d of the cilium in
a quiescen t bac kground fluid. This definition ensures that ϕ (θ + 2π) = ϕ (θ) and implies a
44
mapping
ϕ j
(θ
j
) = Ω 0
Z
θ
j
0
d
¯
θ
Ω(
¯
θ)
, (5.15)
from the angular p ositionθ
j
of the cilium along its b eat cycle on to its phaseϕ j
. This mapping
is unique as long as the cilium do es not rev erse its direction of motion. In the case of phase
indep enden t driving, that is, when F(θ
j
) = constant , ϕ j
and θ
j
are iden tical. In all other
cases ϕ j
and θ
j
are distinct, as sho wn in Fig. 5.3 for four represen tativ e examples of the first
four force harmonics. Imp ortan tly , this phase definition implies that in a system where the
cilia phases ϕ j
are distributed uniformly on the unit circle, the total force cilia exert on the
fluid, and th us the induced flo w field, v anishes. W e will call this state phase isotropic.
W e substitute ( 5.15 ) in to ( 3.5 ); to arriv e at the equations of m otion written in terms of
ϕ j
in the general case where the flo w fie ld is not necessarily zero, i.e., v̸= 0 ,
˙
ϕ j
= Ω 0
+
Ω 0
bΩ j
t
j
·v(r
j
). (5.16)
T aking the time deriv ativ e of ( 5.13 ) and using ( 5.16 ), w e arriv e at the equations of motion
for the momen ts Y
n
˙
Y
n
(x) =
1
ρ
c
X
j
ine
inϕ
j
δ(x−x
j
)
Ω 0
+
Ω 0
bΩ j
t
j
·v(x
j
+bn
j
)
!
. (5.17)
W e next expand the induced v elo cit yv(x
j
+bn
j
) =v(x
j
)+b∇v(x
j
)·n
j
+O(b
2
) ab out the
tra jec tory cen terx
j
and arriv e at
˙
Y
n
(x) = inΩ 0
Y
n
(x)+
inΩ 0
ρ
c
X
j
e
inϕ
j
δ(x−x
j
)
t
j
bΩ j
·
v(x
j
)+b∇v(x
j
)·n
j
+O(b
2
)
!
. (5.18)
45
Figure 5.3: T op ro w: Sample force profiles, F(θ) , corresp onding to the first through the
fourth harmonics in ( 3.1 ). Middle ro w the transformation from angular phase θ to a new
phaseϕ in whic h the in trinsic angular sp eed, Ω 0
= 2π/T
0
, of a single cilium is constan t (sho wn
in inset). Here, T
0
is the in trinsic time p erio d defined b y
R
2π
0
d
¯
θ/Ω(
¯
θ) and is in v arian t with
the force harmonic, since the in tegrals of the sine and cosine terms drop out. Bottom ro w:
Cilia that are uniformly distributed inϕ (ϕ ∈ Unif[0,2π) ) exhibit a non-uniform distribution
in θ that reflect main features of the force harmonic: single-p eak distribution for the first
force harmonic, double-p eak distribution for the second force harmonic, and so on.
46
Considering b/d≪ 1 , dropping higher order terms and rearranging, w e get
˙
Y
n
(x) = inΩ 0
Y
n
(x)+
inΩ 0
ρ
c
X
j
e
inϕ
j
δ(x−x
j
)
"
v(x
j
)·
t
j
bΩ j
+∇v(x
j
) :
t
j
⊗n
j
Ω j
#
, (5.19)
where ‘:’ denotes tensor con traction. Next, it is useful to define
a(ϕ ) =
t(ϕ )
bΩ( ϕ )
, B(ϕ ) =
t(ϕ )⊗n(ϕ )
Ω( ϕ )
, (5.20)
where a(ϕ ) is a v ector, B(ϕ ) is a second order tensor. Their complex F ourier series ma y b e
written in terms of ϕ asa(ϕ ) =
∞
P
m=−∞
e
a
m
e
imϕ
,B(ϕ ) =
∞
P
m=−∞
e
B
m
e
imϕ
, where
e
a
m
=
1
2π
Z
2π
0
a(ϕ )e
−imϕ
dϕ,
e
B
m
=
1
2π
Z
2π
0
B(ϕ )e
−imϕ
dϕ.
(5.21)
Substituting this in to equation ( 5.19 ), w e get
˙
Y
n
(x) = inΩ 0
Y
n
(x)+
inΩ 0
ρ
c
X
j
e
inϕ
j
δ(x−x
j
)
"
v(x
j
)·
X
m
e
a
m
e
imϕ
j
+∇v(x
j
) :
X
m
e
B
m
e
imϕ
j
#
.
(5.22)
Exc hanging the order of the t w o summations and using the substitution prop ert y of the
Dirac delta,v(x
j
)δ(x−x
j
) =v(x)δ(x−x
j
) , w e finally arriv e at
˙
Y
n
(x) = inΩ 0
"
Y
n
(x)+
X
m
e
a
m
·v(x)+
e
B
m
:∇v(x)
Y
n+m
(x)
#
, (5.23)
whic h capture the phase dynamics of the system.
47
These equations pro vide a coupled set of linear partial differen tial equations that, together
with ( 5.11 ) and ( 5.12 ), describ e the dynamical ev olution of the macroscopic v ariables Y
n
(x) .
5.2.1 Calculating cilia induced forces and fluid flo ws
W e no w deriv e an expression forf(x) in terms of the coarse-grained v ariablesv(x) andY
n
(x) .
W e re define f(x) from ( 5.1 ):
f(x) =
X
j
δ(x−r
j
)(F
j
+N
j
), (5.24)
and use ( 3.3 ) and ( 3.5 ), to get
f(x) =
X
j
δ(x−r
j
)(ζbΩ j
t
j
−ζ(n
j
⊗n
j
)·v(r
j
)). (5.25)
Expanding this expression ab outx
j
and k eeping terms up toO(b
2
) yields
f(x) =F(x)+∇·Σ(x),
(5.26)
where the surface force densit y exerted b y the con tin uous ciliary la y er is
F(x) =
X
j
δ(x−x
j
)ζ(bΩ j
t
j
−(n
j
⊗n
j
)·v(x
j
)), (5.27)
and t he surface stress densit y exerted b y the ciliary la y er is
Σ(x) =
X
j
δ(x−x
j
)ζbn
j
⊗(bΩ j
t
j
−(n
j
⊗n
j
)·v(x
j
)). (5.28)
48
F o r notational con v enience, w e in tro duce,
c = bΩ j
t
j
, D =n
j
⊗n
j
, E = bΩ j
n
j
⊗t
j
, G = 2n
j
⊗n
j
⊗n
j
,
(5.29)
wherec =c(ϕ ) is a v ector,D =D(ϕ ) andE =E(ϕ ) are second order tensors, andG =G(ϕ )
is a third order tensor. Next, w e write their F ourier series expansions,
c =
∞
X
m=−∞
e
c
m
e
imϕ
j
, D =
∞
X
m=−∞
f
D
m
e
imϕ
j
, E =
∞
X
m=−∞
e
E
m
e
imϕ
j
, G =
∞
X
m=−∞
f
G
m
e
imϕ
j
.
(5.30)
Here
e
c
m
,
f
D
m
,
e
E
m
, and
f
G
m
are the asso ciated F ourier co efficien ts and follo w from form ulae
similar to those in ( 5.21 ). With this
F(x) = ρ
c
X
m
ζ
e
c
m
−
f
D
m
·v(x)
Y
m
, Σ(x) = ρ
c
X
m
ζb
e
E
m
−
f
G
m
·v(x)
Y
m
, (5.31)
and finally
f(x) = ρ
c
ζ
X
m
e
c
m
−
f
D
m
·v(x)
Y
m
+b∇·
(
e
E
m
−
f
G
m
·v(x))Y
m
!
, (5.32)
whic h is the cilia generated forcing on the system. The flo w v elo cit y v can b e calculated b y
solving ( 5.11 ), together with ( 5.12 ), whic h is linear i nf .
49
5.3 Stabilit y of isotropic and sync hronized states
W e use the con tin uum theory dev elop ed ab o v e to in v estigate the s tabilit y of t w o spatially
homogeneous states that do not pump fluid: (i) the isotropic steady state, and (ii) the fully
sync hronized steady state. W e are ultimately in terested in understanding the conditions
under whic h spatially non-homogeneous, patterned states whic h can pump fluid emerge
sp on taneously in these ciliary carp ets.
The spatially-homogeneous isotropic steady state of the ciliary system is defined b y
Y
m
= 0 , for all non-zero m . In this case, the force f that the cilia c ollectiv ely exert on
the fluid v anishes. Th us no fluid flo w is induced and the individual cilia act as if they w ere
indep enden t.
The fully sync hronized steady state is the state in whic h the phase of all cilia is the same,
i.e., ϕ j
= ϕ for all j . In this case,∇v = 0 andv is parallel to the instan taneous direction of
motion of all cilia suc h that ( 5.19 ) simplifies to
˙
Y
n
= in
1+
v(ϕ )
bΩ( ϕ )
!
Ω 0
Y
n
, (5.33)
where v = ∥v∥ = v·t . Th us, the fully sync hronized state is also a p erio dic steady state.
The same conclusion can also b e arriv ed at directly from ( 5.16 ). Note that in the case of an
M -fold symmetric forcing function, the system admits M additional spatially homogeneous
p erio dic steady states, as can b e directly seen from ( 5.16 ). Hereafter, w e mak e no further
commen t on these sp ecial cases.
The emergence of spatially anisotropic states that can pump fluid necessarily requires
50
that all spatially homogeneous steady states of the system are unstable. Next, w e in v estigate
the conditions under whic h the t w o spatially homogeneous steady states that w e iden tified
here go unstable, starting with the isotropic steady state.
5.3.1 Stabilit y of the isotropic state
W e in v estigate the linear stabilit y of the isotropic state sub ject to small p erturbations. F or
notational con v enience, w e use the X
∗
to indicate the steady state v alue of a quan tit y X
and δX to indicate a small p erturbation a w a y from the steady state. F or instance Y
n
(x) =
Y
∗
n
(x)+δY
n
(x) ,f(x) =f
∗
(x)+δf(x) , andv(x) =v
∗
(x)+δv(x) .
F or the isotropic state, w e ha v e Y
∗
0
= 1 , whic h is the constan t cilia surface densit y , and
Y
∗
n
= 0 for all n ̸= 0 . Also, f
∗
(x) and v
∗
(x) are iden tically zero. Substituting in to ( 5.23)
and ( 5.32 ), resp ectiv ely , w e get to linear order,
δ
˙
Y
n
(x) = inΩ 0
δY
n
+
e
a
−n
·δv(x)+
e
B
−n
:∇δv(x)
. (5.34)
Expanding ( 5.32 ), to linear order again, w e find
δf(x) = ρ
c
ζ
−
f
D
0
·δv−b∇·
f
G
0
·δv+
X
m
e
c
m
δY
m
+b∇·
e
E
m
δY
m
!
. (5.35)
W e next write equations ( 5.34 ) and ( 5.35 ) in F ourier space:
δ
ˆ
˙
Y
n
(k) = inΩ 0
δ
ˆ
Y
n
(k)+
e
a
−n
+i
e
B
−n
·k
·δˆ v(k), (5.36)
51
and
δ
ˆ
f(k) =−ρ
c
ζ(
f
D
0
+ibk·
f
G
0
)·δˆ v(k)+ρ
c
ζ
X
m
e
c
m
+ibk·
e
E
m
δ
ˆ
Y
m
(k), (5.37)
and substitute ( 5.37 ) in to ( 5.11 ). After some rearranging, w e arriv e at
h
I+ρ
c
ζ
ˆ
K·
f
D
0
+ibk·
f
G
0
i
·δˆ v(k) = ρ
c
ζ
ˆ
K·
X
m
e
c
m
+ibk·
e
E
m
δ
ˆ
Y
m
(k).
(5.38)
Without further appro ximations, to linear order, w e can write δˆ v(k) =
P
m
ˆ ν
m
(k)δ
ˆ
Y
m
(k)
where, ˆ ν
m
(k) is a v ector-v alued function,
ˆ ν
m
(k) =
h
I+ρ
c
ζ
ˆ
K·(
f
D
0
+ibk·
f
G
0
)
i
−1
·ρ
c
ζ
ˆ
K·
X
m
e
c
m
+ibk·
e
E
m
.
(5.39)
Substituting this expression for δˆ v in ( 5.36 ), w e arriv e at
δ
˙
ˆ
Y
n
(k) = inΩ 0
δ
ˆ
Y
n
(k)+(
e
a
−n
+i
e
B
−n
·k)·
X
m
ˆ ν
m
(k)δ
e
Y
m
(k)
!
=
X
m
ˆ
L
nm
(k)δ
ˆ
Y
m
(k),
(5.40)
where
ˆ
L
nm
(k) = inΩ 0
δ
nm
+(
e
a
−n
+i
e
B
−n
·k)·ν
m
(k)
, (5.41)
is a second order tensor that couples the p erturbations in the order parameter fields. The
largest eigen v alue γ(k) of
ˆ
L(k) is the gro wth rate of the fastes t gro wing eigen v ector of the
system. The sign of the real part of γ(k) determines the o v erall stabilit y of the isotropic
52
Figure 5.4: Stabilit y of the isotropic state. The gro wth rate of the instabilit y in the F ourier
space, γ(k) , computed from ( 5.41 ), is sho wn for first four force harmonics. In all four cases,
almost all mo des are linearly unstable, with the instabilit y b eing particularly strong for the
second force harmonic.
steady state.
T o ev aluate ( 5.41 ) n umerically , w e truncate the summation in ( 5.39 ) retaining only the
termsm =−5 throughm = +5 . This results in the coupling matrix
ˆ
L(k) (with en tries
ˆ
L
nm
)
ha ving dimension 11×11 . The largest eigen v alue giv es the gro wth rate γ(k) of the system.
The top ro w in Fig. 5.4 sho ws the surface plots of the gro wth rates for the first four force
harmonics, as a function of w a v e-n um b er pairsk = (k
x
,k
y
) . W e find that the isotropic state
is linearly unstable to p erturbations of almost all w a v en um b ers, with the instabilities b eing
strongest for the second force harmonic.
5.3.2 Stabilit y of the sync hronized state
In the fully sync hronized state, all cilia b eat in-phase and ϕ j
= ϕ ∗
for all j . This also implies
that θ
j
= θ
∗
and v(x
j
) = v
∗
for all j and th us the system is in steady state, see Eq ( 3.5 ).
The angular sp eed at steady state is
˙
θ
∗
= Ω ∗
+t
∗
·v
∗
/b . Note that the in trinsic angular
sp eed Ω ∗
= Ω( θ
∗
) is a p erio dic function of the angle θ
∗
. In this section, w e ask for the
linear stabilit y of the fully sync hronized state. W e p erform the analysis using a set of closure
53
conditions on Y
n
similar to those used for the isotropic state, but w e also sho w another w a y
to p erform this analysis b y in tro ducing the p erturbation directly in the angular phase θ , and
using the appropriate closure conditions in the App endix B .
As in the previous section, the sup erscript ()
∗
denotes quan tities at the steady state,
and the prefix δ indicates the p erturbation. W e start from the observ ation that in the fully
sync hronized state, the v elo cit y gradien t is iden tically zero
∇v|
∗
= 0, (5.42)
and the order parameters tak e the form
Y
∗
m
= e
imϕ
= (Y
∗
1
)
m
. (5.43)
Substituting ( 5.42 ) in to ( 5.32 ) giv es the force densit y in the fully sync hronized state
f
∗
= ρ
c
ζ
X
m
e
c
m
−
f
D
m
·v
∗
Y
∗
m
, (5.44)
and the fluid v elo cit y
v
∗
=
ρ
c
ζh
η
I+
ρ
c
ζh
η
X
m
f
D
m
Y
∗
m
!
−1
·
X
m
e
c
m
Y
∗
m
. (5.45)
W e next note that giv en Eq. ( 5.43 ), p erturbations δY
m
around the fully sync hronized steady
54
state can b e expressed in terms of δY
1
using the relation
δY
m
= mY
∗
m−1
δY
1
. (5.46)
With this, to linear order, ( 5.23 ) b ecomes
δ
˙
Y
1
= iΩ 0
"
δY
1
+
X
m
e
a
m
·δv+
e
B
m
:∇δv
Y
∗
m+1
+(m+1)
e
a
m
·v
∗
Y
∗
m
δY
1
#
. (5.47)
Similarly , taking v ariations of ( 5.32 ) giv es
δf = ρ
c
ζ
X
m
e
c
m
−
f
D
m
·v
∗
δY
m
−
f
D
m
·δvY
∗
m
+b∇·
e
E
m
δY
m
−b∇·
f
G
m
·δvY
∗
m
,
(5.48)
In F ourier space, ( 5.47 ) and ( 5.48 ) b ec ome
δ
ˆ
˙
Y
1
= iΩ 0
"
δ
ˆ
Y
1
+
X
m
e
a
m
+i
e
B
m
·k
·δˆ vY
∗
m+1
+(m+1)
e
a
m
·v
∗
Y
∗
m
δ
ˆ
Y
1
#
, (5.49)
and
δ
ˆ
f = ρ
c
ζ
X
m
e
c
m
−
f
D
m
·v
∗
mY
∗
m−1
δ
ˆ
Y
1
−
f
D
m
·δˆ vY
∗
m
+ibk·
e
E
m
mY
∗
m−1
δ
ˆ
Y
1
−ibk·
f
G
m
·δˆ vY
∗
m
.
(5.50)
Substituting ( 5.50 ) in to ( 5.11 ) w e obtain an expression for the v elo cit y p erturbation
δˆ v =νδ
ˆ
Y
1
(5.51)
55
where the v ector ν is giv en b y
ν =
"
I+ρ
c
ζ
ˆ
K·
X
m
f
D
m
+ibk·
f
G
m
Y
∗
m
#
−1
·
ρ
c
ζ
ˆ
K·
X
n
n
e
c
n
−
f
D
n
·v
∗
+ibk·
e
E
n
Y
∗
n−1
!
.
(5.52)
Finally substituting this bac k in to equation ( 5.49 ) w e arriv e at the linearized equation of
motion
δ
ˆ
˙
Y
1
=
ˆ
L(t)δ
ˆ
Y
1
(5.53)
where
ˆ
L(t) = iΩ 0
1+
X
m
e
a
m
+i
e
B
m
·k
·ν·Y
∗
m+1
+(m+1)
e
a
m
·v
∗
·Y
∗
m
!
. (5.54)
It is w orth noting that
ˆ
L here is a complex scalar, while it is a complex matrix in the
linearized equations ab out the isotropic state ( 5.40 ). This is a direct consequence of ( 5.43),
whic h states that the dynamics of the system in the fully sync hronized state can b e captured
using Y
1
alone without the need for higher momen ts. Eq. ( 5.53 ) is a first order ODE that
describ es the time ev olution of the p erturbation δ
ˆ
Y
1
(k,t) ; it admits the solution
δ
ˆ
Y
1
(k,t = τ) = δ
ˆ
Y
1
(k,t = 0)e
R
τ
0
ˆ
L(t)dt
. (5.55)
Equation ( 5.55 ) giv es the gro wth of t he p erturbation up to time t = τ . Since the fully
sync hronized state is a p erio dic steady state, w e are sp ecifically in terested in the gro wth
56
Figure 5.5: Stabilit y of the fully sync hronized state. The Ly apuno v exp onen t, µ (k) , com-
puted from ( B.14 ), is sho wn for first four force harmonics. Unlik e the isotropic state, in all
four cases, some mo des are linearly stable while others are not. But the instabilit y is again
particularly strong for the second force harmonic.
of the p erturbation o v er one b eat cycle, T
∗
=
R
2π
0
dϕ/ϕ
∗
. T o this end, w e calculate the
Ly apuno v exp onen t, µ (k) , whic h is giv en b y
µ (k) =
1
T
∗
Real
Z
T
∗
0
ˆ
L(t)dt
!
.
(5.56)
Finally , w e use the transformation dt = dϕ/
˙
ϕ ∗
to get
µ (k) =
1
T
∗
Real
Z
2π
0
ˆ
L(t)
˙
ϕ ∗
dϕ
!
. (5.57)
The steady state is linearly stable if the Ly apuno v exp onen t µ (k) < 0 , critically stable for
µ (k) = 0 , and linearly unstable when µ (k) > 0 .
W e ev aluate ( B.14 ) n umerically to obtain the Ly apuno v exp onen t at eac h w a v en um b er
pair k . Results are summarized in the b ottom ro w of Fig. 5.5 for four force harmonics.
As in the isotropic state, w e find that the stronger instabilit y o ccurs for the second force
harmonic. Ho w ev er, in con trast to the isotropic state where p erturbations are unstable for
all w a v en um b ers, here p erturbations of certain w a v en um b er pairs gro w, while others deca y .
57
5.4 Conclusions
In this w ork, w e deriv ed a con tin uum theory for ciliary surfaces based on a simplified rotor
mo del for b eating cilia. This theory enables the quan titativ e study of large ciliated systems,
whic h, giv en their sheer size, are outside the scop e of discrete approac hes.
W e used the theory to iden tify the spatially isotropic p erio dic steady states of the system
and to in v estigate their stabilit y . W e sho w ed that large ciliary carp ets ha v e t w o spatially
isotropic steady states. The first is the isotropic steady state, in whic h ciliary phases are
dra wn from a uniformly random distribution. The second is the fully sync hronized steady
state in whic h all cilia in the system ha v e the same phase. W e further sho w ed that b oth
of these steady states are generically unstable, whic h implies that these carp ets will ev olv e
in to patterned states, suc h as c himera states or w a v e states. Related w ork that studies the
stabilit y of w a v e states in a similar setup suggests in that metac hronal w a v es whic h pump
fluid are stable at least in parts of the phase space [ 50 ]. As alluded to in the previous c hapter,
our o wn n umerical in v estigations suggest that w a v es patterns can emerge sp on taneously from
isotropic or sync hronized states.
The theory that w e presen ted here iden tifies and trac ks the macroscopic v ariables Y
n
,
whic h c haracterize the state of sync hron y of a ciliary carp et on large time and length scales.
The structure of the asso ciated equations of motion ( 5.23 ) do ho w ev er presen t new c hallenges.
In particular the equations of motion couple all momen ts Y
n
to eac h other. Th us to mak e
the system tractable, closure assumptions need to b e made. Here, w e close the hierarc h y
of equations of motion around the t w o spatially isotropic states of in terest b y prop osing
t w o differen t closures. W e use Y
∗
n
= 0 for |n| ̸= 0 around the isotropic steady state and
58
Y
∗
n
= (Y
∗
1
)
n
around the fully sync hronized steady state. In [ 50 ], the authors face the same
c hallenge. They close the system b y p ostulating Y
n
= A
n
exp(inϕ (x,t)) where A
n
are fixed
constan ts, and only allo w the phase ϕ (x,t) to dynamically v ary , effectiv ely restricting their
analysis to systems with smo othly v arying phase angle fields. A general closure that captures
the full dynamics of the system is y et to b e found and will b e the sub ject of future researc h.
W e ac kno wledge that the theory presen ted here mak es imp ortan t simplifications. In
particular, the microscopic mo del that underlies our descriptions is form ulated suc h that all
forces are exerted along fixed circular tra jectories on a t w o dimensional plane. This enables
us to deriv e an analytic k ernel for deriving the flo w fields generated b y b eating cilia. More
w ork is needed to extend our approac h to allo w for forces that are applied all throughout
the cilia length and for describing b eat patterns whose shap es c hange under the effect of
collectiv ely generated fluid flo ws.
Despite these simplifications, our theory is mak es imp ortan t con tribution to analyzing
the emergence of large scale cilia co ordination. It can easily b e adapted to treat spatially
non-uniform systems suc h as the ones studied in [ 70 ] and more complex forcing functions,
suc h as the ones measured exp erimen tally in [ 65 , 51 ], for example. The w ork presen ted
here is th us an imp ortan t step to w ards dev eloping a strong understanding of the ph ysics of
sync hornization in large ciliary carp ets.
59
Chapter 6
Comparison of the n umerical and
theoretical approac hes
I
n this c hapter, w e quan titativ ely compare the n umerical and theoretical approac hes de-
v elop ed in the previous c hapters to study ciliary carp ets. W e start b y demonstrating
the inadequacy of the Kuramoto order parameter, whic h is obtained b y a v eraging o v er all
cilia, in capturing the app earance of spatial w a v e patterns. Then, w e in tro duce the “Ku-
ramoto ellipse” as a new metric that satisfies our requiremen t. This allo ws us to compare the
theoretical predictions of gro wth rates of instabilities from c hapter 5 against the n umerical
sim ulations from c hapter 4 for the isotropic initial state. W e then compare the Ly apuno v ex-
p onen ts for the fully sync hronized states; this turns out to b e more straigh tforw ard without
the need for an y new metrics.
60
6.1 Commen t on the Kuramoto order parameter
W e use the algorithm devised in c hapter 4 to sim ulate a 151× 151 , 2-dimensional arra y of
cilia for 100 time units. Figure 6.1 A sho ws the cosine fields of the emergen t state, for the
first t w o force harmonics for short and long cilia, starting from an isotropic state, i.e. the
sim ulations are initialized b y dra wing from a uniformly random distribution of ϕ , and the
angles θ are obtained using the b ac kw ard transformation ( 5.15 ).
Recollect that w e defined t he momen t fields
Y
n
(x,t) =
1
ρ
c
N
X
j
e
inϕ
j
δ(x−x
j
), (6.1)
where ρ
c
= N/L
2
is the cilia densit y , and taking the a v erage of Y
1
o v er the all cilia in the
primary domain giv es to the Kuramoto order parameter
P =|⟨Y
1
⟩| =
Z
L/2
−L/2
Z
L/2
−L/2
Y
1
(x,y)dxdy
. (6.2)
Figure 6.1 B sho ws the time series of P for the sim ulations of long and short cilia. The blac k
line corresp onds to the first force harmonic, and the blue line to the second force harmonic.
The v alue ofP at initial time is close to zero due to lac k of an y spatial order. Ho w ev er, notice
that P con tin ues to remain close to zero ev en at long time, long after the spatial patterns
are clearly visible in the cosine fields in 6.1 A. This is not surprising. A w a v eform, ev en a
noisy one, w ould con tain phase angles almost uniformly distributed in the range [0,2π) and
naturally , corresp onds to a P v alue close to zero. This demonstrates wh y the Kuramoto
order parameter, a v eraged o v er all cilia is not a useful metric to measure spatial patterns.
61
Figure 6.1: (A) Numerical sim ulations of a 151× 151 arra y of cilia sho w the emergence
of tra v eling w a v e patterns. Snapshots tak en at t = 100 for h/d = 0.5 (left) and h/d =
1.5 (righ t). (B) Time ev olution of the Kuramoto p olar order parameter P sho ws no clear
distinction b et w een isotropic and tra v eling w a v e patterns. (C) W e define P
α
for α∈ [0,2π]
and plot (α,P
α
) in p olar co ordinates at t w o snapshots tak en from the isotropic state and
tra v eling w a v e pattern. Data p oin ts are sho wn in blac k and elliptic b est-fit in blue. (D)
Time ev olution of the eccen tricit y and angle of the Kuramoto ellipse sho ws clear correlations
b et w een its eccen tricit y with the gro wth of the tra v eling w a v e instabilit y and angle with
dominan t direction of w a v e fron t (see Supplemen tal Mo vie). Sim ulations w ere p erformed
for 100 time units using dt = 0.1 , a = 0.05 , b = 0.2 , d = 1 . First force harmonic: F =
1+0.5cosθ +0.5sinθ , Second force harmonic: F = 1+0.5cos2θ +0.5sin2θ .
6.2 The Kuramoto ellipse
W e in tro duce the “Kuramoto ellipse” as a metric to quan tify the app earance of spatial w a v e
patterns. Starting from the discrete particle mo del, w e compute the Kuramoto order field
[ 75 ] at the lo cationx
j
of a cilium j b y a v eraging o v er all cilia in a lo cal neigh b orho o d of j ,
Y
1
(x
j
) =
1
ρ
c
X
ℓ∈N(j)
e
iϕ
ℓ
. (S1)
62
Here, w e c ho ose N(j) con taining one shell of neigh b ors around c ilium j . W e rep eat this
at the lo cation of eac h cilium on the grid to estimate the Kuramoto order field Y
1
(x) from
discrete cilia sim ulations.
Next, w e dra w straigh t lines starting from the origin of the fundam en tal domain at angles
α∈ [0,2π) , and compute the a v erage magnitude of order parameter along that direction:
P
α
=
Z
L/2
0
Y
1
(scosα,ssinα)ds
, (S3)
s b eing a dumm y in tegration parameter o v er all p oin ts that lie along the said line. Note that
P
α
is the absolute v alue of (a constrained v ersion of ) the Radon transform of the 2D field
Y
1
(x) [ 76 ]. W e then use a p olar co ordinate represen tation (α,P
α
) to obtain the data p oin ts
indicated b y the green dots in Figure 6.2 C. W e apply a Principal Comp onen t Analysis to
these data p oin ts to fit, what w e call, the K ur amoto el lipse , whose ma jor and minor axes
corresp ond to the first and second principle comp onen ts, resp ectiv ely .
F or concreteness, w e demonstrate the utilit y of the Kuramoto ellipse in Figure 6.2 , using
t w o sample states– an isotropic state (top ro w) and a syn thetically generated w a v e state
(b ottom ro w). In the isotropic state, (α,P
α
) are distributed in an almost circular fashion
around the origin, indicating lac k of spatial order. In the w a v e state, since spatial order is
maxim um along the direction of w a v e orien tation, the ellipse stretc hes out in that direction.
W e refer to the eccen tricit y of the ellipse as the K ur amoto e c c entricity , and the inclination
of the first principal comp onen t to the horizon tal as the K ur amoto angle . The Kuramoto
eccen tricit y is a measure of the degree of coherence of w a v es, while the inclinations of the first
and second principal comp onen t giv e us the directions of w a v e orien tation and propagation,
63
Figure 6.2: Illustration of Kuramoto ellipse using snapshots of an isotropic state (top) and
a syn thetically created w a v e state (b ottom). Left column sho ws the cosine field, obtained
b y taking the cosine of the phase of the cilium. Cen ter column sho ws the corresp onding
Kuramoto order field– Real (Y
1
(x)) – whic h is calculated b y coarse graining o v er one shell of
neigh b ors. The righ tmost column sho ws the p oin ts (α,P(α)) in green, and the corresp onding
Kuramoto ellipse in gra y .
resp ectiv ely .
W e no w apply this metric to data from our sim ulations. Figure 6.1 C sho ws the Kuramoto
ellipse at the initial and final time for the second force harmonic, and figure 6.1 D sho ws the
time series of the Kuramoto eccen tricit y and angle, clearly demonstrating the emergence of
spatial coherence.
6.3 Comparing gro wth rates of the isotropic state
In this section, w e compare the gro wth rates of the instabilit y in the isotropic state b et w een
con tin uum theory and n umerical sim ulations.
Recall that using the con tin uum mo del in c hapter 5 , w e arriv ed at the equation for the
64
time ev olution of the p erturbation in F ourier space
δ
˙
ˆ
Y
n
(k,t) =
∞
X
m=−∞
ˆ
L
nm
(k,t)δ
ˆ
Y
m
(k,t), (6.3)
and the maximal real part of the eigen-v alues of
ˆ
L ga v e us the gro wth rate γ(k) . In fig-
ure. 6.3 A, w e rep ort theγ(k) as a colormap o v er allk for the first and second force harmonics.
Note that so far this is the same information as in figure 6.3 . The maxim um w a v e n um b er
and resolution in w a v e space are c hosen to b e consisten t with a discrete 151× 151 system
of cilia with lattice spacing d = 1 ; th us, the maxim um w a v e-n um b er is 2π/2d = π/d and
the resolution 2π/151d . Next, w e sup erimp ose data from a sample discrete cilia sim ulation
(blac k dots) of ( 3.5 ); the dots represen t the w a v en um b ers at whic h |
ˆ
Y
1
(k,t)| is maximal
(exceeding a giv en threshold). Note the agreemen t b et w een the sim ulation result and the
maximal gro wth rate in the con tin uum theory .
T o further c hallenge the theory , starting from the initial state used in the sample sim u-
lation, w e calculate
ˆ
Y
1
(k,0) , and in tegrate forw ard in time
ˆ
Y
1
(k,∆ t) =
ˆ
Y
1
(k,0)(1+γ(k)∆ t)
according to the linear gro wth rate from the con tin uum mo del. It is w orth noting here that
while it is p ossible to go from from a giv enϕ j
distribution to a uniqueY
n
(x) distribution, the
con v erse is not true, i.e. a fixed Y
n
(x) distribution cannot b e uniquely mapp ed to a ϕ j
dis-
tribution. W e transform bac k to ph ysical space and compare Real[Y
1
(x,t)] to that obtained
from the sample sim ulation (Fig. 6.3 B) and again find remarkable agreemen t, indicating that
the w a v e patterns are go v erned b y the fastest gro wing linear mo de. T o further c hallenge the
theory , w e calculate the Kuramoto ellipse based on Mon te-Carlo computations using the
con tin uum mo del gro wth rate with 200 initial isotropic states; the a v erage is sho wn in dark
65
grey in Fig. 6.3 C and the range in ligh ter grey while results from the sample sim ulation are
sup erimp osed in blac k dots. The time ev olution of the ellipses’ eccen tricit y and angle are
sho wn in Fig. 6.3 D using the same color con v en tion. The ellipse eccen tricit y correlates with
the instabilit y gro wth rate and its angle in ph ysical space is orthogonal to the direction of
maxim um gro wth rate in F ourier space. W e conclude that the fastest gro wing w a v en um b ers
in the linear con tin uum theory are predictiv e of the long-term tra v eling w a v e patterns in the
nonlinear discrete mo del.
6.4 Comparing gro wth rates of the fully sync hronized
state
Again, recall that in c hapter 5 , w e ended up with the equation for the time ev olution of the
p erturbation in F ourier space
δ
˙
ˆ
Y
1
(k,t) =
ˆ
L(k,t) δ
ˆ
Y
1
(k,t). (6.4)
The corresp onding Ly apuno v exp onen ts µ (k) w ere calculate b y cycle-a v eraging the instan-
taneous gro wth rates
µ (k) =
1
T
Real
"
Z
T
∗
0
ˆ
L(k,t) dt
#
, (6.5)
with T
∗
=
R
2π
0
dϕ/
˙
ϕ ∗
b eing the time p erio d.
W e ev aluate µ (k) at eac h w a v en um b er k and compare to direct n umerical sim ulations
66
of ( 3.5 ) o v er one b e at cycle. W e initialize the sim ulation with w a v e p erturbations
δY
1
= ε exp(i(k
x
x
j
+k
y
y
j
)), (6.6)
whereε≪ 1 is the magnitude of the p erturbation,k
(·)
= 2πq
(·)
/d
√
N , andq
(·)
=−⌊
√
N/2⌋,...,⌊
√
N/2⌋ .
The phase of eac h rotor is then obtained usingδϕ = cos
−1
Real[δY
1
] . Note that this is a lo ose
appro ximation to coarse graining, as w e use only one cilium rather than a collection of them.
The resulting Ly apuno v exp onen ts µ (k) are extracted using the ratio of the magnitudes of
the p erturbation at initial and final times, i.e.
µ (k) = Real
"
ln
δ
ˆ
Y
1
(k,t = T
∗
)
δ
ˆ
Y
1
(k,t = 0)
!#
(6.7)
Results are sho wn in Fig. 6.4 A o v er the en tire w a v e space and in Fig. 6.4 B for select cross-
sections. W e again note the remarkable agreemen t in the Ly apuno v exp onen ts of the insta-
bilit y ab out the sync hronous states, obtained from theory and sim ulations.
6.5 Conclusion
The quan titativ e agreemen t b et w een the t w o approac hes is remarkable for b oth the p erio dic
steady states. This indicates that this theory could pro vide an analytical windo w in to the
mec hanisms of symmetry breaking, sp on taneous phase co ordinatio n, and self-organized fluid
transp ort in ciliary carp ets.
67
Figure 6.3: Isotropic state. (A) Gro wth rate γ(k) of the instabilit y predicted b y the con tin uum theory (colormap) and discrete
sim ulations (blac k dots) in the F ourier space. (B) Starting from a random Kuramoto field Y
1
(x) , map to F ourier space
ˆ
Y
1
(k)
and use γ(k) from panel A to step forw ard in time. Snapshots of Re|Y
1
(x,t)| sho w the emergence of tra v eling w a v es similar
to those obtained from direct n umerical sim ulations of Eq. ( 3.5 ). (C) Kuramoto ellipse in tro duced in Fig. 1E obtained from
the sample sim ulation in panel B (blac k dots) and theory corresp onding to 200 realizations (grey), with a v erage sho wn in dark
grey . (D) Time ev olution of eccen tricit y and angle of the Kuramoto ellipse; sample sim ulation (blac k line) and theory (grey).
68
Figure 6.4: Sync hronized state: Ly apuno v exp onen t µ (k) obtained from theory and sim-
ulations (A) o v er the en tire w a v e space (colormap), and (B) at cross-sections k
y
d = −π/4
and k
y
=−π/2 with theory in solid lines and sim ulations in dashed lines.
69
Chapter 7
Conclusions
C
iliar y co ordination is a complex problem and in v olv es in v oking to ols and insigh ts
from man y differen t fields of researc h including statistical ph ysics, applied mathe-
matics, biology to name a few.
In this thesis, w e started with a coarse-grained appro ximation, the so called r otor mo del
of a cilium, and dev elop ed an algorithm that enables large scale sim ulations of ciliary car-
p ets. W e sho w ed via n umerical sim ulations that the isotropic state is unstable and leads
to the emergence of metac hronal co ordination among cilia sp on taneously . T o quan tify their
w a v elength and direction, w e in tro duced a Kuramoto ellipse based on a simplified Radon
transform of the Kuramoto order field. This metric can p oten tially b e mo dified and ex-
pressed as a function of space to rev eal spatial heterogeneit y and defects in the tra v eling
w a v e patterns. W e sho w ed that while an individual cilium do es not pump an y fluid b y it-
self, the formation of spatial w a v es breaks time rev ersal symmetry to ac hiev e net pumping,
akin to those observ ed in ciliated tissues. The emergen t w a v es are indep enden t of initial
conditions, but dep end on the prop erties of individual cilia.
70
W e then dev elop ed a con tin uum theory framew ork for studying ciliary carp ets and con-
ducted a linear stabilit y analysis that predicts quan titativ ely the long-term tra v eling w a v e
patterns seen in the discrete cilia sim ulations. The gro wth rates from n umerics and theory
sho w ed remarkably go o d agreemen t with eac h other.
Both our n umerical and theoretical framew orks are amenable to future extensions, b y
relaxing a few simplifications w e made along the w a y . F or example, our theory can b e
readily extended to incorp orate non-circular limit-cycle represen tations of individual cilia
obtained from exp erimen tal data of cilia b eating patterns [ 65 ] and non-square, p oten tially
patc h y , lattice structures reconstructed from high-resolution images of ciliated tissues [ 70 ].
In addition, b oth framew orks can b e mo dified to study sync hronization in ciliary carp ets in
micro c hannels. This can b e done in theory b y adding a no-slip b oundary condition at the top
c hannel w all, and using the resulting v elo cit y expression in the linear stabilit y analysis. In
the computational framew ork, the p erio dic Blak e tensor can b e mo dified using appropriate
set of singularities, as p er the deriv ation in [ 77 , 78 ]. This setup could b e easier to access
and v alidate exp erimen tally . These extensions will lead to causal explanations of the effect
of individual cilia oscillations and their spatial distribution along ciliated tissues on the
emergence of large scale co ordination and fluid transp ort in ciliary carp ets.
Although w e restricted oursel v es to h ydro dynamic coupling alone in this w ork, recen t
literature has suggested basal coupling as another strong factor in causing sync hroniza-
tion [ 79 , 80 , 81 , 82 ]. The future of this line of inquiry could lead to ric h, m ulti-ph ysics
mo dels, pa ving w a y for a more holistic understanding of ciliary sync hronization. The author
of this thesis hop es to ha v e made a small con tribution to w ards this grand goal.
71
Bibliograph y
[1] Arkady Pik o vsky , Jurgen Kurths, Mic hael Rosen blum, and Jürgen Kurths. Synchr o-
nization: a universal c onc ept in nonline ar scienc es . Cam bridge univ ersit y press, 2003.
[2] Blausen Medical. Medical gallery of blausen medical 2014. WikiJournal of Me dicine ,
1(2):1–79, 2014.
[3] William Gilpin, Matthew Storm Bull, and Man u Prakash. The m ultiscale ph ysics of
cilia and flagella. Natur e R eviews Physics , 2(2):74–88, 2020.
[4] Da vid M W o olley , Rac hel F Cro c k ett, William DI Gro om, and Stuart G Rev ell. A study
of sync hronisation b et w een the flagella of bull sp ermatozoa, with related observ ations.
Journal of Exp erimental Biolo gy , 212 (14):2215–2223, 2009.
[5] Marco P olin, Idan T uv al, Kn ut Dresc her, Jerry P Gollub, and Ra ymond E Goldstein.
Chlam ydomonas swims with t w o “gears” in a eukary otic v ersion of run-and-tum ble
lo comotion. Scienc e , 325(5939):487–490, 2009.
[6] Linda T urner, William S R yu, and Ho w ard C Berg. Real-time imaging of fluorescen t
flagellar filamen ts. Journal of b acteriolo gy , 182(10):2793–2801, 2000.
72
[7] P an teleimon Romp olas, Ramila S P atel-King, and Stephen M King. An outer arm
dynein conformational switc h is required for metac hronal sync hron y of motil e cilia in
planaria. Mole cular biolo gy of the c el l , 21(21):3669–3679, 2010.
[8] Eric Lauga and Ra ymond E Goldstein. Dance of the microswimmers. Phys. T o day ,
65(9):30, 2012.
[9] Kn ut Dresc her, Ra ymond E Goldstein, Nicolas Mic hel, Marco P olin, and Idan T u-
v al. Direct measuremen t of the flo w field around swimming micro organisms. Physic al
R eview L etters , 105(16):168101, 2010.
[10] W en Y an and Mic hael Shelley . Univ ersal image systems for non-p erio dic and p erio dic
stok es flo ws ab o v e a no-slip w all. Journal of Computational Physics , 375:263–270,
2018.
[11] Arth ur T Winfree. Biological rh ythms and the b eha vior of p opulations of coupled
oscillators. Journal of the or etic al biolo gy , 16(1):15–42, 1967.
[12] Y oshiki Kuramoto. Self-en trainmen t of a p opulation of coupled non-linear oscillators.
In International symp osium on mathematic al pr oblems in the or etic al physics , pages
420–422. Springer, 1975.
[13] Charles S P eskin. Mathematical asp ects of heart ph ysiology . Cour ant Inst. Math , 1975.
[14] Renato E Mirollo and Stev en H Strogatz. Sync hronization of pulse-coupled biological
oscillators. SIAM Journal on A pplie d Mathematics , 50(6):1645–1662, 1990.
73
[15] Stev en H Strogatz. F rom kuramoto to cra wford: exploring the onset of sync hronization
in p opulations of coupled oscillators. Physic a D: Nonline ar Phenomena , 143(1-4):1–20,
2000.
[16] Stev en Strogatz. Sync: The emer ging scienc e of sp ontane ous or der . P enguin UK, 2004.
[17] Nathalie Spassky and Alice Meunier. The dev elopmen t and functions of m ulticiliated
epithelia. Natur e r eviews Mole cular c el l biolo gy , 18(7):423–436, 2017.
[18] Gregor Eic hele, Eb erhard Bo densc hatz, Zuzana Ditte, Ann-K athrin Gün ther, Shoba
Kap o or, Y ong W ang, and Christian W estendorf. Cilia-driv en flo ws in the brain third
v en tricle. Philosophic al T r ansactions of the R oyal So ciety B , 375(1792):20190154, 2020.
[19] RA Ly ons, E Saridogan, and O Djahan bakhc h. The repro ductiv e significance of h uman
fallopian tub e cilia. Human r epr o duction up date , 12(4):363–372, 2006.
[20] Björn A Afzelius, P er Camner, and Björn Mossb erg. On the function of cilia in the
female repro ductiv e tract. F ertility and sterility , 29(1):72–74, 1978.
[21] Ximena M Bustaman te-Marin and La wrence E Ostro wski. Cilia and m uco ciliary clear-
ance. Cold Spring Harb or p ersp e ctives in biolo gy , 9(4):a028241, 2017.
[22] BA Afzelius. Cilia-related diseases. The Journal of Patholo gy: A Journal of the
Patholo gic al So ciety of Gr e at Britain and Ir eland , 204(4):470–477, 2004.
74
[23] Shigenori Nonaka, Y osuk e T anaka, Y asushi Okada, Sen T ak eda, Akihiro Harada,
Y oshimitsu Kanai, Mizuho Kido, and Nobutaka Hiroka w a. Randomization of left–
righ t asymmetry due to loss of no dal cilia generating left w ard flo w of extraem bry onic
fluid in mice lac king kif3b motor protein. Cel l , 95(6):829–837, 1998.
[24] Geoffrey Ingram T a ylor. Analysis of the swimming of microscopic organisms. Pr o c e e d-
ings of the R oyal So ciety of L ondon. Series A. Mathematic al and Physic al Scienc es ,
209(1099):447–461, 1951.
[25] Eric Lauga. An in tro duction to the h ydro dynamics of lo comotion on small scales. Soft
Interfac es: L e ctur e Notes of the L es Houches Summer Scho ol: V olume 98, July 2012 ,
98:81, 2012.
[26] Edw ard M Purcell. Life at lo w reynolds n um b er. A meric an journal of physics , 45(1):3–
11, 1977.
[27] James Gra y . Ciliary movement . C am bridge Univ ersit y Press, 1928.
[28] Maciej Lisic ki. F our approac hes to h ydro dynamic green’s functions–the oseen tensors.
arXiv pr eprint arXiv:1312.6231 , 2013.
[29] JR Blak e. Infinite mo dels for ciliary propulsion. Journal of Fluid Me chanics , 49(2):209–
222, 1971.
[30] Zydrunas Gim butas, Leslie Greengard, and Shra v an V eerapaneni. Simple and efficien t
represen tations for the fundamen tal solutions of stok es flo w in a half-space. Journal of
Fluid Me chanics , 776, 2 015.
75
[31] Josephine Ainley , Sandra Durkin, Rafael Em bid, Priy a Boindala, and Ricardo Cortez.
The metho d of images for regularized stok eslets. Journal of Computational Physics ,
227(9):4600–4616, 2008.
[32] Natan Osterman and Andrej Vilfan. Finding the ciliary b eating pattern with optimal
efficiency . Pr o c e e dings of the National A c ademy of Scienc es , 108(38):15727–15732,
2011.
[33] Jens Elgeti and Gerhard Gompp er. Emergence of metac hronal w a v es in cilia arra ys.
Pr o c e e dings of the National A c ademy of Scienc es , 110(12):4470–4475, 2013.
[34] Hanliang Guo, Lisa F auci, Mic hael Shelley , and Ev a Kanso. Bistabilit y in the syn-
c hronization of actuated microfilamen ts. Journal of Fluid Me chanics , 836:304–323,
2018.
[35] Brato Chakrabarti and Da vid Sain tillan. Sp on taneous oscillations, b eating patterns,
and h ydro dynamics of activ e microfilamen ts. Physic al R eview Fluids , 4(4):043102,
2019.
[36] Yi Man and Ev a Kanso. Multisync hron y in activ e microfilamen ts. Physic al R eview
L etters , 125(14):148101, 2020.
[37] Sha y Gueron, K onstan tin Levit-Gurevic h, Nada v Liron, and Jacob J Blum. Cilia
in ternal mec hanism and metac hronal co ordination as the result of h ydro dynamical
coupling. Pr o c e e dings of the National A c ademy of Scienc es , 94(12):6001–6006, 1997.
[38] Christopher W ollin and Holger Stark. Metac hronal w a v es in a c hain of ro w ers with
h ydro dynamic in teractions. The Eur op e an Physic al Journal E , 34(4):1–10, 2011.
76
[39] Ev elyn Hamilton and Pietro Cicuta. Changes in geometrical asp ects of a simple mo del
of cilia sync hronization con trol the dynamical state, a p ossible mec hanism for switc hing
of swimming gaits in microswimmers. PloS one , 16(4):e0249060, 2021.
[40] Nicolas Bruot and Pietro Cicuta. Realizing the ph ysics of motile cilia sync hronization
with driv en colloids. A nnual R eview of Condense d Matter Physics , 7(1):323–348, 2016.
[41] Andrej Vilfan and F rank Jülic her. Hydro dynamic flo w patterns and sync hronization
of b eating cilia. Physic al r eview letters , 96(5): 058102, 2006.
[42] Thomas Niederma y er, Bruno Ec khardt, and P eter Lenz. Sync hronization, phase lo c k-
ing, and metac hronal w a v e formation in ciliary c hains. Chaos: A n Inter disciplinary
Journal of Nonline ar Scienc e , 18(3):037128, 2008.
[43] Jurij K otar, Luk e Deb ono, Nicolas Bruot, Stuart Bo x, Da vid Phillips, Stephen Simp-
son, Simon Hanna, and Pietro Cicuta. Optimal h ydro dynamic sync hronization of
colloidal rotors. Physic al r eview lette rs , 111(22):228103, 2013.
[44] Douglas R Brumley , Marco P olin, Timoth y J P edley , and Ra ymond E Goldstein.
Metac hronal w a v es in the flagellar b eating of v olv o x and their h ydro dynamic origin.
Journal of the R oyal So ciety Interfac e , 12(108):20141358, 2015.
[45] Armando Maestro, Nicolas Bruot, Jurij K otar, Nariy a Uc hida, Ramin Golestanian,
and Pietro Cicuta. Con trol of sync hronization in mo dels of h ydro dynamically coupled
motile cilia. Communic ations Physics , 1(1):1–8, 2018.
[46] P eter Lenz and Andrey R yskin. Collectiv e effects in ciliar arra ys. Physic al biolo gy ,
3(4):285, 2006.
77
[47] Nariy a Uc hida and Ramin Golestanian. Generic conditions for h ydro dynamic sync hro-
nization. Physic al R eview L etters , 106(5):058104, 2011.
[48] Douglas R Brumley , Nicolas Bruot, Jurij K otar, Ra ymond E Goldstein, Pietro Cicuta,
and Marco P olin. Long-range in teractions, w obbles, and phase defects in c hains of
mo del cilia. Physic al r eview fluids , 1(8):081201, 2016.
[49] Nariy a Uc hida and Ramin Golestanian. Sync hronization and collectiv e dynamics in a
carp et of microfluidic rotors. Physic al r eview letters , 104(17):178103, 2010.
[50] F anlong Meng, Rac hel R Bennett, Nariy a Uc hida, and Ramin Golestanian. Conditions
for metac hronal co ordination in arra ys of mo del cilia. Pr o c e e dings of the National
A c ademy of Scienc es , 118(32), 2021.
[51] An ton Solo v ev and Benjamin M F riedric h. Global metac hronal sync hronization and
activ e noise in cilia carp ets. arXiv pr eprint arXiv:2012. 11741 , 2020.
[52] Jens Elgeti and Gerhard Gompp er. Microswimmers near surfaces. The Eur op e an
Physic al Journal Sp e cial T opics , 225(11):2333–2352, 2016.
[53] P etr Denissenk o, V asily Kan tsler, Da vid J Smith, and Jac kson Kirkman-Bro wn. Hu-
man sp ermatozoa migration in micro c hannels rev eals b oundary-follo wing na vigation.
Pr o c e e dings of the National A c ademy of Scienc es , 109(21):8007–8010, 2012.
[54] He-P eng Zhang, A vraham Be’er, E-L Florin, and Harry L Swinney . Collectiv e motion
and densit y fluctuations in bacterial colonies. Pr o c e e dings of the National A c ademy of
Scienc es , 107(31):13626–13630, 2010.
78
[55] Sa v erio E Spagnolie and Eric Lauga. Hydro dynamics of self-propulsion near a b ound-
ary: predictions and accuracy of far-field appro ximations. Journal of Fluid Me chanics ,
700:105–147, 2012.
[56] Da vid L Kirk. V olvox: a se ar ch for the mole cular and genetic origins of multic el lularity
and c el lular differ entiation . Cam bridge Univ ersit y Press, 2005.
[57] Boris Martinac, Y oshiro Saimi, and Ching Kung. Ion c hannels in microb es. Physiolo g-
ic al r eviews , 88(4):1449–1490, 2008.
[58] Hoang-Ngan Nguy en and Karin Leiderman. Computation of the singular and regular-
ized image systems for doubly-p erio dic stok es flo w in the presence of a w all. Journal
of Computational Physics , 297:442–461, 2015 .
[59] Shriram Sriniv asan and Anna-Karin T orn b erg. F ast ew ald summation for green’s func-
tions of stok es flo w in a half-space. R ese ar ch in the Mathematic al Scienc es , 5(3):1–19,
2018.
[60] JR Blak e and A T Ch w ang. F undamen tal singularities of viscous flo w. Journal of
Engine ering Mathematics , 8(1):23–29, 1974.
[61] Hidenori Hasimoto. On the p erio dic fundamen tal solutions of the stok es equations
and their application to viscous flo w past a cubic arra y of spheres. J. Fluid Me ch ,
5(02):317–328, 1959.
[62] Da vid Sain tillan, Eric Darv e, and Eric SG Shaqfeh. A smo oth particle-mesh ew ald
algorithm for stok es susp ension sim ulations: The sedimen tation of fib ers. Physics of
Fluids , 17(3):033301, 2005.
79
[63] Ricardo Cortez and F ranz Hoffmann. A fast n umerical metho d for computing doubly-
p erio dic regularized stok es flo w in 3d. Journal of Computational Physics , 258:1–14,
2014.
[64] Y ang Ding, Janna C Na wroth, Margaret J McF all-Ngai, and Ev a Kanso. Mixing and
transp ort b y ciliary carp ets: a n umerical study . Journal of Fluid Me chanics , 743:124–
140, 2014.
[65] Douglas R Brumley , Kirst y Y W an, Marco P olin, and Ra ymond E Goldstein. Flagellar
sync hronization through direct h ydro dynamic in teractions. elife , 3:e02750, 2014.
[66] Boris Guirao and Jean-F rançois Joann y . Sp on taneous creation of macroscopic flo w and
metac hronal w a v es in an arra y of cilia. Biophysic al journal , 92(6):1900–1917, 2007.
[67] N Liron. Fluid transp ort b y cilia b et w een parallel plates. Journal of Fluid Me chanics ,
86(4):705–726, 1978.
[68] D.J Smith, J.R Blak e, and E.A Gaffney . Fluid mec hanics of no dal flo w due to em bry-
onic primary cilia. Journal of The R oyal So ciety Interfac e , 5(22):567–573, 2008.
[69] An up V Kanale, F eng Ling, Hanliang Guo, Mic hael Shelley , Sebastian F ürthauer, and
Ev a Kanso. Sp on taneous phase co ordination in ciliary carp ets. arXiv pr eprint arXiv: ,
2021.
[70] Guillermina R Ramirez-San Juan, Arnold JTM Mathijssen, Mu He, Lily Jan, W al-
lace Marshall, and Man u Prakash. Multi-scale spatial heterogeneit y enhances particle
clearance in airw a y ciliary arra ys. Natur e Physics , 16(9):958–964, 2020.
80
[71] Nicolas Bruot, Jurij K otar, Filipp o De Lillo, Marco Cosen tino Lagomarsino, and Pietro
Cicuta. Driving p oten tial and noise lev el det ermine the sync hronization state of h y-
dro dynamically coupled oscillators. Physic al r eview letters , 109(16):164103, 2012.
[72] Sebastian F ürthauer and Sriram Ramasw am y . Phase-sync hronized state of orien ted
activ e fluids. Phys. R ev. L ett. , 111:238102, Dec 2013.
[73] Hassan Masoud and Mic hael J Shelley . Collectiv e surfing of c hemically activ e particles.
Physic al r eview letters, 112(12):128304, 2014.
[74] T ong Gao, Rob ert Blac kw ell, Matthew A Glaser, MD Betterton, and Mic hael J Shelley .
Multiscale mo deling and sim ulation of microtubule–motor-protein assem blies. Physic al
R eview E , 92(6):062709, 2015.
[75] Y oshiki Kuramoto. Chemic al oscil lations, waves, and turbulenc e . Courier Corp oration,
2003.
[76] S Helgason and S Helgason. The R adon T r ansform . Springer, 1980.
[77] Nada v Liron and S Mo c hon. Stok es flo w for a stok eslet b et w een t w o parallel flat plates.
Journal of Engine ering Mathematics , 10(4):287–303, 1976.
[78] Iv an T anasijević and Eric Lauga. Hydro dynamic sync hronization in strong confine-
men t. Physic al R eview E , 103(2):02 2403, 2021.
[79] Kirst y Y W an and Ra ymond E Goldstein. Co ordinated b eating of algal flagella
is mediated b y basal coupling. Pr o c e e dings of the National A c ademy of Scienc es ,
113(20):E2784–E2793, 2016.
81
[80] Gary S Klindt, Christi an R uloff, Christian W agner, and Benjamin M F riedric h. In-
phase and an ti-phase flagellar sync hronization b y w a v eform compliance and basal cou-
pling. New Journal of Physics , 19(11): 113052, 2017.
[81] Y ujie Liu, Rory Cla ydon, Marco P olin, and Douglas R Brumley . T ransitions in syn-
c hronization states of mo del cilia through basal-connection coupling. Journal of The
R oyal So ciety Interfac e , 15(147):20180450, 2018.
[82] Hanliang Guo, Yi Man, Kirst y Y W an, and Ev a Kanso. In tracellular coupling mo du-
lates biflagellar sync hron y . Journal of the R oyal So ciety Interfac e , 18(174):20200660,
2021.
[83] Nicola P ellicciotta, Ev elyn Hamilton, Jurij K otar, Marion F aucourt, Nathalie Delgeh yr,
Nathalie Spassky , and Pietro Cicuta. En trainmen t of mammalian motile cilia in the
brain with h ydro dynamic forces. Pr o c e e dings of the National A c ademy of Scienc es ,
117(15):8315–8325, 2020.
[84] Martin B. Short, Cristian A. Solari, Sujo y Ganguly , Thomas R. P o w ers, John O.
Kessler, and Ra ymond E. Goldstein. Flo ws driv en b y flagella of m ulticellular organ-
isms enhance long-range molecular transp ort. Pr o c e e dings of the National A c ademy of
Scienc es , 103(22):8315–8319, 2006.
[85] Y ang Ding and Ev a Kanso. Selectiv e particle capture b y async hronously b eating cilia.
Physics of Fluids , 2 7(12):121902, 2015.
[86] Sébastien Mic helin and Eric Lauga. Unsteady feeding and optimal strok es of mo del
ciliates. Journal of Fluid Me chanics , 715:1–31, 2013.
82
[87] Hanliang Guo, Hai Zh u, R uo w en Liu, Marc Bonnet, and Shra v an V eerapaneni. Optimal
ciliary lo comotion of axisymmetric microswimmers, 2021.
[88] Sébastien Mic helin and Eric Lauga. Efficiency optimization and symmetry-breaking
in a mo del of ciliary lo comotion. Physics of Fluids , 22(11):111901, 2010.
[89] Christophe Elo y and Eric Lauga. Kinematics of the most efficien t cilium. Phys. R ev.
L ett. , 109:038101, Jul 2012.
[90] Douglas R Brumley , Marco P olin, Timoth y J P edley , and Ra ymond E Goldstein. Hy-
dro dynamic sync hronization and metac hronal w a v es on the surface of the colonial alga
v olv o x carteri. Physic al r eview letters , 109(26):268102, 2012.
[91] M. C. Marc hetti, J. F. Joann y , S. Ramasw am y , T. B. Liv erp o ol, J. Prost, Madan
Rao, and R. A diti Simha. Hydro dynamics of soft activ e matter. R ev. Mo d. Phys. ,
85:1143–1189, Jul 2013.
[92] Johanna Raidt, Claudius W erner, T ab ea Menc hen, Gerard W Doughert y , Heik e Ol-
bric h, Niki T Loges, Ralf Sc hmitz, P etra P ennekamp, and Heym ut Omran. Ciliary
function and motor protein comp osition of h uman fallopian tub es. Human R epr o duc-
tion , 30(12):2871–2880, 2015.
[93] Regina F aub el, Christian W estendorf, Eb erhard Bo densc hatz, and Gregor Eic hele.
Cilia-based flo w net w ork in the brain v en tricles. Scienc e , 353(6295): 176–178, 2016.
[94] Jurij K otar, Marco Leoni, Bruno Bassetti, Marco Cosen tino Lagomarsino, and Pietro
Cicuta. Hydro dynamic sync hronization of colloidal oscillators. Pr o c e e dings of the
National A c ademy of Scienc es , 107(17):7669–7673, 2010.
83
[95] Nicolas Bruot. Hydr o dynamic c oupling and synchr onization of c ol loidal oscil lators .
PhD thesis, Univ ersit y of Cam bridge, 2013.
[96] William Gilpin, Viv ek N Prakash, and Man u Prakash. V ortex arra ys and ciliary tangles
underlie the feeding–swimming trade-off in starfish larv ae. Natur e Physics , 13(4):380–
386, 2017.
[97] Laurence Irving. Ciliary curren ts in starfish. Journal of Exp erimental Zo olo gy ,
41(1):115–124, 1924.
[98] Janna C Na wroth, Hanliang Guo, Eric K o c h, Elizab eth A C Heath-Hec kman, John C
Hermanson, Edw ard G R ub y , John O Dabiri, Ev a Kanso, and Margaret McF all-Ngai.
Motile cilia create fluid-mec hanical microhabitats for the activ e recruitmen t of the host
microbiome. Pr o c e e dings of the National A c ademy of Scienc es , 114(36):9510–9516,
2017.
[99] Nariy a Uc hida and Ramin Golestanian. Hydro dynamic sync hronization b et w een ob-
jects with cyclic rigid tra jectories. The Eur op e an Physic al Journal E , 35(12):1–14,
2012.
[100] John Blak e. A mo del for the micro-structure in ciliated organisms. Journal of Fluid
Me chanics , 55(1):1–23, 197 2.
[101] An up V Kanale, Sebastian F ürthauer, F eng Ling, Hanliang Guo, Mic hael Shelley , and
Ev a Kanso. Con tin uum theory for carp ets of rotary cilia. arXiv pr eprint arXiv: , 2021.
[102] An up V Kanale, Hanliang Guo, Sebastian F ürthauer, F eng Ling, Mic hael Shelley , and
Ev a Kanso. Large scale sim ulations of ciliary w a v es. arXiv pr eprint arXiv: , 2021.
84
[103] Barath Ezhilan, Mic hael J Shelley , and Da vid Sain tillan. Instabilities and nonlinear
dynamics of concen trated activ e susp ensions. Physics of Fluids , 25(7):070 607, 2013.
[104] Eric R Bro oks and John B W allingford. Multiciliated cells. Curr ent Biolo gy ,
24(19):R973–R982, 2014.
[105] Arian S F orouhar, Mic hael Liebling, Anna Hic k erson, Abbas Nasiraei-Moghaddam,
Huai-Jen T sai, Ja y R Ho v e, Scott E F raser, Mary E Dic kinson, and Morteza
Gharib. The em bry onic v ertebrate heart tub e is a dynamic suction pump. Scienc e ,
312(5774):751–753, 2006.
[106] Keith T urton. 20 - pumps and pumping. In Dennis A. Sno w, editor, Plant Engine er’s
R efer enc e Bo ok (Se c ond Edition) , pages 20–1–20–32. Butterw orth-Heinemann, Oxford,
second edition edition, 2002.
[107] P atric k R Sears, W ei-Ning Yin, and La wrence E Ostro wski. Con tin uous m uco ciliary
transp ort b y primary h uman airw a y epithelial cells in vitro. A meric an Journal of
Physiolo gy-Lung Cel lular and Mole cular Physiolo gy , 309(2):L99–L108, 2015.
[108] Susyn Joan Kelly , V o jta Bro dec ky , Elizab eth M Skuza, PJ Berger, and Stanisla v
T atk o v. V ariabilit y in trac heal m uco ciliary transp ort is not con trolled b y b eating cilia
in lam bs in viv o during v en tilation with h umidified and nonh umidified air. A meri-
c an Journal of Physiolo gy-Lung Cel lular and Mole cular Physiolo gy , 320(4):L473–L485,
2021.
85
[109] Kakani Katija, Giancarlo T roni, Jo ost Daniels, Kelly Lance, Rob E Sherlo c k, Alana D
Sherman, and Bruce H Robison. Rev ealing enigmatic m ucus structures in the deep sea
using deeppiv. Natur e , 583(7814):78–82, 2020.
[110] JR Blak e. A note on the image system for a stok eslet in a no-slip b oundary . In Pr o c.
Camb. Phil. So c , v olume 70, pages 303–310, 1971.
[111] Hanliang Guo, Hai Zh u, and Shra v an V eerapaneni. Sim ulating cilia-driv en mixing and
transp ort in complex geometries. Physic al R eview Fluids , 5(5):053103, 2020.
[112] V eikk o F Gey er, F rank Jülic her, Jonathon Ho w ard, and Benjamin M F riedric h. Cell-
b o dy ro c king is a dominan t mec hanism for flagellar sync hronization in a swimming
alga. Pr o c e e dings of the National A c ademy of Scienc es , 110(45):18058–18063, 2013.
[113] F orest O Mannan, Miika Jarv ela, and Karin Leiderman. Minimal mo del of the h ydro-
dynamical coupling of flagella on a spherical b o dy with application to v olv o x. Physic al
R eview E , 102(3):033114, 2020.
[114] Hanliang Guo, Janna Na wroth, Y ang Ding, and Ev a Kanso. Cilia b eating patterns are
not h ydro dynamically optimal. Physics of Fluids , 26(9):091901, 2014.
[115] Christopher Brennen and Ho w ard Winet. Fluid mec hanics of propulsio n b y cilia and
flagella. A nnual R eview of Fluid Me chanics , 9(1):339–398, 1977.
[116] Eric Lauga and Thomas R P o w ers. The h ydro dynamics of swimming micro organisms.
R ep orts on Pr o gr ess in Physics , 72(9):096601, 2009.
86
[117] Julia M Y eomans, Dmitri O Pushkin, and Henry Sh um. An in tro duction to the h y-
dro dynamics of swimming micro organisms. The Eur op e an Physic al Journal Sp e cial
T opics , 223(9):177 1–1785, 2014.
[118] Eric Lauga. Bacterial h ydro dynamics. A nnual R eview of Fluid Me chanics , 48:105–130,
2016.
[119] Pietro Cicuta. The use of bioph ysical approac hes to understand ciliary b eating. Bio-
chemic al So ciety T r ansactions , 48(1):221–229 , 2020.
[120] M A Sleigh. The biolo gy of Cilia and Flagel la . Macmillan Co., New Y ork, 1962.
[121] Nobutaka Hiroka w a, Y asushi Okada, and Y osuk e T anaka. Fluid dynamic mec hanism
resp onsible for breaking the left-righ t symmetry of the h uman b o dy: the no dal flo w.
A nnual r eview of fluid me chanics , 41:53–72, 2009.
[122] James Ligh thill. Flagellar h ydro dynamics. SIAM r eview , 18(2):161–230, 1976.
[123] Da vid Sain tillan and Mic hael J Shelley . Instabilities and pattern formation in ac-
tiv e particle susp ensions: kinetic theory and con tin uum sim ulations. Physic al R eview
L etters , 100(17):178103, 2008.
[124] Ra ymond E Goldstein. Batc helor prize lecture fluid dynamics at the scale of the cell.
Journal of Fluid Me chanics , 807:1– 39, 2016.
[125] Sylv ain Chateau, Julien F a vier, Sébastien P oncet, and Um b erto d’Ortona. Wh y an-
tiplectic metac hronal cilia w a v es are optimal to transp ort bronc hial m ucus. Physic al
R eview E , 100(4):042405, 2019.
87
[126] John R Blak e. A spherical en v elop e approac h to ciliary propulsion. Journal of Fluid
Me chanics , 46(1):199–208, 1 971.
[127] Jens Elgeti, Roland G Winkler, and Gerhard Gompp er. Ph ysics of microswimmers—
single particle motion and collectiv e b eha vior: a review. R ep orts on pr o gr ess in physics ,
78(5):056601, 2015.
[128] Nariy a Uc hida and Ramin Golestanian. Sync hronization in a carp et of h ydro dynam-
ically coupled rotors with random in trinsic frequency . EPL (Eur ophysics L etters) ,
89(5):50011, 2010.
[129] Brato Chakrabarti, Sebastian F ürthauer, and Mic hael J Shelley . A m ultiscale bio-
ph ysical mo del giv es quan tized metac hronal w a v es in a lattice of cilia. arXiv pr eprint
arXiv:2108.01693 , 2021.
88
App endices
89
App endix A
Ev aluating the b oundary conditions
to obtain the fluid in teraction Kernel
T o solv e for the unkno wn co efficien ts A
±
, B
±
, C
±
, D
±
,E
±
,G
±
in ( 5.9 ) and ( 5.10 ), w e first
apply the b oundedness condition of ( 5.6 ) as z→∞ to eliminate the exp onen tially gro wing
terms in ( 5.9 ) and ( 5.10 ). Substituting the resulting expressions in to ( 5.7 ) and matc hing
co efficien ts leads to
ik·E
+
−kC
+
+
A
+
2η
= 0. (A.1)
Similarly , substituting expressions for ˆ p
−
, ˆ v
−
and w
−
in to equation ( 5.7 ) and matc hing
co efficien ts leads to
ik·E
−
−kC
−
+
A
−
2η
= 0, ik·G
−
+kD
−
+
B
−
2η
= 0. (A.2)
90
Next, w e enforce the no-slip condition ( 5.5 ) at z = 0 to obtain
E
−
+
ihkA
−
2ηk
!
e
kh
+
G
−
−
ihkB
−
2ηk
!
e
−kh
= 0,
C
−
−
hA
−
2η
!
e
kh
+
D
−
−
hB
−
2η
!
e
−kh
= 0.
(A.3)
By con tin uit y of v elo cit y at z = h , w e get
C
+
= C
−
+D
−
, E
+
=E
−
+G
−
(A.4)
Finally , w e turn to the jump condition ( 5.3 ) at the ciliary la y er. In comp onen t form,σ·e
z
can
b e expressed as σ·e
z
= (η(∂
z
u+∂
x
w),η(∂
z
v +∂
y
w),−p+2η∂
z
w), whic h can b e rewritten
in F ourier space as ˆ σ·e
z
= (η(∂
z
ˆ v+ikˆ w),−ˆ p+2η∂
z
ˆ w). Ev aluating ˆ σ·e
z
at z = h
+
and
z = h
−
and matc hing the jump in ( 5.3 ) giv es
ηk(E
+
−E
−
+G
−
)+
ik
2k
(A
+
−A
−
+B
−
) =
ˆ
f, 2ηk(C
+
−C
−
+D
−
) = 0. (A.5)
where
ˆ
f = (
ˆ
f
x
,
ˆ
f
y
) and the v ertical comp onen t off is iden tically zero b y construction.
Equations ( A.1 )-( A.5 ) form a closed system of algebraic equations that w e solv e analyt-
ically to obtain the co efficien ts A
±
, B
±
, C
±
, D
±
,E
±
,G
±
using the Mathematica noteb o ok
giv en b elo w.
91
Velocities and Pressure definitions
k = {k1, k2};
mk = Sqrt[k.k];
pp = apExp[-mk (z - h)];
pm = amExp[-mk (z - h)] + bmExp[mk (z - h)];
wp = (cp + ap / (2 μ) (z - h))Exp[-mk (z - h)];
wm = (cm + am / (2 μ) (z - h))Exp[-mk (z - h)] + (dm +
bm / (2 μ) (z - h))Exp[mk (z - h)];
ep = {ep1, ep2};
vp = (ep - Ik / (2 μ mk)ap (z - h))Exp[-mk (z - h)];
em = {em1, em2};
fm = {fm1, fm2};
vm = (em - Ik / (2 μ mk)am (z - h))Exp[-mk (z - h)] + (fm +
Ik / (2 μ mk)bm (z - h))Exp[mk (z - h)];
Define Incompressibility and boundary conditions
incomp1 = Ik.ep - mkcp + ap /2 / μ ⩵ 0;
incomp2 = Ik.em - mkcm + am /2 / μ ⩵ 0;
incomp3 = Ik.fm + mkdm + bm /2 / μ ⩵ 0;
sigma13p = μ (Ik1wp + D[vp 〚1 〛, z]);
sigma13m = μ (Ik1wm + D[vm 〚1 〛, z]);
sigma23p = μ (Ik2wp + D[vp 〚2 〛, z]);
sigma23m = μ (Ik2wm + D[vm 〚2 〛, z]);
sigma33p = 2 μ D[wp, z] - pp;
sigma33m = 2 μ D[wm, z] - pm;
Eqns = {
incomp1, incomp2, incomp3,
(vm 〚1 〛 /. {z → 0}) ⩵ 0,
(vm 〚2 〛 /. {z → 0}) ⩵ 0,
(wm ⩵ 0 ) /. {z → 0},
(vm 〚1 〛 ⩵ vp 〚1 〛) /. {z → h},
(vm 〚2 〛 ⩵ vp 〚2 〛) /. {z → h},
(wm ⩵ wp) /. {z → h},
(sigma13p ⩵ sigma13m - T1) /. {z → h},
(sigma23p ⩵ sigma23m - T2) /. {z → h},
(sigma33p ⩵ sigma33m - T3) /. {z → h}
};
2
Solve for all unknowns and derive kernel for z=h
sol = Solve[
Eqns, {ap, am, bm, ep1, ep2, em1, em2, fm1, fm2, cm, cp, dm}];
StandardForm [FullSimplify [
vp 〚1 〛 /. sol /. {z → h, T3 → 0}, { κ^2 ⩵ k1^2 + k2^2, κ > 0}]]
StandardForm [FullSimplify [
vp 〚2 〛 /. sol /. {z → h, T3 → 0}, { κ^2 ⩵ k1^2 + k2^2, κ > 0}]]
Out[161]//StandardForm=
1
4 κ
3
μ
ⅇ
-2 h κ
-T1 κ
2
1 - ⅇ
2 h κ
+2h κ (-1 +h κ) +
k2
2
T1 -1 + ⅇ
2 h κ
+2h κ (-1 +h κ) -k1k2T2 -1 + ⅇ
2 h κ
+2h κ (-1 +h κ)
Out[162]//StandardForm=
ⅇ
-2 h κ
-k1k2T1 -1 + ⅇ
2 h κ
+2h κ (-1 +h κ) +T2 2 × -1 + ⅇ
2 h κ
κ
2
-k2
2
-1 + ⅇ
2 h κ
+2h κ (-1 +h κ)
4 κ
3
μ
3
App endix B
An alternativ e metho d to analyze
linear stabilit y of the fully
sync hronized state
In this section, w e ev alute the linear stablit y of the fully sync hronized state using a differen t
set of closure conditions than the ones used in 5.3 .
W e seek to find an explicit expression for the angular sp eed
˙
θ
∗
in the sync hronized state,
whic h requires kno wledge of the coarse-grained cilia-induced v elo cit yv
∗
. The latter in turn
dep ends on the coarse-grained force densit yf
∗
. By virtue of ( 5.25 ), the force densit y in the
fully syc hronized state is giv en b y
f
∗
= ρ
c
(ζbΩ ∗
t
∗
−ζ(n
∗
·v
∗
)n
∗
) (B.1)
Substituting ( B.1 ) in to ( 5.11 ), using ( 5.12 ) and taking the in v erse F ourier transform, w e
95
arriv e at
v
∗
=
ρ
c
ζbh
η
Ω ∗
t
∗
, (B.2)
where w e used the fact that in the fully sync hronized state, the quan tities Ω ∗
, t
∗
, n
∗
and
v
∗
are spatially in v arian t and therefore, their spatial deriv ativ es m ust v anish lea ving b ehind
only the mean flo w con tribution corresp onding to the k = 0 term. Lastly , substituting ( B.2 )
in to the expression for
˙
θ
∗
leads to t he constan t angular sp eed in the sync hronized state
˙
θ
∗
=
1+
ρ
c
ζh
η
!
Ω ∗
. (B.3)
That is, the angular sp eed of eac h cilium is alw a ys larger than that of a single isolated cilium.
W e no w in tro duce a small p erturbation δθ ab out the sync hronized state and write θ =
θ
∗
+δθ andv =v
∗
+δv . Note that Ω( θ
∗
+δθ) = Ω( θ
∗
)+ (dΩ /dθ)|
θ
∗
δθ . Plugging in to ( 3.5 )
and simplifying using ( B.2 ) and ( B.3 ), w e arriv e at
δ
˙
θ =
dΩ dθ
θ
∗
δθ +
1
b
t
∗
·δv. (B.4)
Expressed in F ourier space, w e get
δ
˙
ˆ
θ =
dΩ dθ
θ
∗
δ
ˆ
θ +
1
b
t
∗
·δˆ v. (B.5)
T o ev aluate δˆ v , w e first substitute ( B.1 ) in to ( 5.11 ), use ( 5.12 ) and tak e the in v erse F ourier
96
transform to get
v = ρ
c
ζF
−1
[
ˆ
K][bΩ ∗
t
∗
−(n
∗
·v
∗
)n
∗
], (B.6)
where F
−1
[
ˆ
K] is an op erator resulting from the in v erse F ourier transformation of
ˆ
K , and
con tains spatial deriv ativ es (remem b er that ik⇌∇ ). Then, w e tak e v ariations,
δv = ζρ
c
F
−1
[
ˆ
K]
"
(b
dΩ dθ
θ
∗
t
∗
−bΩ ∗
n
∗
−(t
∗
·v
∗
)n
∗
−(n
∗
·v
∗
)t
∗
)δθ
#
−ζρ
c
F
−1
[
ˆ
K][n
∗
⊗n
∗
δv],
(B.7)
and express bac k in F ourier space to arriv e at
δˆ v = ˆ νδ
ˆ
θ, (B.8)
where
ˆ ν = ζρ
c
I+ζρ
c
ˆ
K·(n
∗
⊗n
∗
)
−1
·
ˆ
K·
b
dΩ dθ
θ
∗
t
∗
−bΩ ∗
n
∗
−(t
∗
·v
∗
)n
∗
!
. (B.9)
Substituting this result in to equation ( B.5 ), w e get
δ
˙
ˆ
θ =
ˆ
Lδ
ˆ
θ, (B.10)
where
ˆ
L =
Z
t
0
dΩ dθ
θ
∗
+
1
b
t
∗
· ˆ ν
!
. (B.11)
97
( B.10 ) describ es the time-ev olution of the p erturbation δ
ˆ
θ(k,t) in F ourier space. Note that
in this linear regime, p erturbations at differen t w a v en um b ers do not in teract with eac h other.
It admits the solution
δ
ˆ
θ(k,t) = δ
ˆ
θ(k,t = 0)e
R
τ
0
ˆ
Ldt
(B.12)
This equation giv es the gro wth of the p erturbation up to time t = τ . Since the fully
sync hronized state is a p erio dic steady state, w e are sp ecifically in terested in the gro wth
of the p erturbation o v er one b eat cycle, T
∗
=
R
2π
0
dϕ/ϕ
∗
. T o this end, w e calculate the
Ly apuno v exp onen t, µ (k) , whic h is giv en b y
µ (k) =
1
T
∗
Real
Z
T
∗
0
ˆ
L(t)dt
!
.
(B.13)
Finally , w e use the transformation dt = dϕ/
˙
ϕ ∗
to get
µ (k) =
1
T
∗
Real
Z
2π
0
ˆ
L(t)
˙
ϕ ∗
dϕ
!
. (B.14)
The equilibirum state is said to b e linearly stable if µ (k) < 0 , and unstable otherwise.
W e ev aluate ( B.14 ) n umerically to obtain the gro wth rate of p erturbations at eac h w a v e
n um b e r k and the results are summarized in Fi g. A.1 for four force harmonics. Note that
this is, somewhat exp ectedly , similar to the result in 5.5 , indicating that b oth sets of closure
conditions w ork w ell for the fully sync hronized state.
98
Figure A.1: Stabilit y of the fully sync hronized state. The Ly apuno v exp onen t, µ (k) , com-
puted from ( B.14 ), is sho wn for first four force h armonics.
99
App endix C
Comparison of cilia-induced flo wfields
W e compare the flo w-fields induced b y ciliary carp ets in the t w o mo dels: (1) the particle
mo del that uses direct summation of the Blak e tensor with truncation [ 102 ], and (2) the
con tin uum mo del that uses the expression ˆ v(k) =
ˆ
K(k)·
ˆ
f(k) [ 101 ]. Here,
ˆ
f(k) is the force
densit y , expressed in F ourier space with w a v e v ector k , and
ˆ
K(k) is the in tegral Kernel of
the Stok es equation in the half-space asso ciated with the discon tin uit y la y er due to the force
densit yf ; see Eq. (4) in the main do cumen t. Detailed deriv ation is presen ted in [ 101 ].
In the con tin uum theory , p erio d icit y is implicit b y virtue of expressing the mo del in
F o urier space. Sp ecifically , w e use w a v en um b er resolution dk = 2π/L , where L is the length
of the main square domain. The maxim um w a v en um b er is k
max
= 2π/dx , where dx is the
grid resolution in the ph ysical space. W e ev aluate ˆ v(k) at eac h w a v en um b erk and w e in v ert
it n umerically to get the v elo cit y field u(x,y,z) in ph ysical space. Note that u(x,y,z) can
b e decomp osed in to a horizon tal comp onen tv(x) =u(x,y,z = h) and a v ertical comp onen t.
Figure A.1 sho ws the flo w field generated b y a single force monop ole (cilium) placed ab o v e
a no-slip w all; Figure A.1 A sho ws the in-plane flo wfield, Figure A.1 B sho ws the flo wfields in
100
Figure A.1: Comparison of flo w fields for a unit stok eslet in half space, orien ted in the
p ositiv e x -direction. (A) Flo w fields in the ciliary plane based on the con tin uum theory
(left) and particle mo del (righ t), resp ectiv ely . V elo cit y magnitude is sho wn as a colormap
and streamlines in white lines. In the particle mo del, w e use the Blak e tensor in the half-
space together with a truncated image system to em ulate the effect of the doubly-p erio dic
domain (B) Flo ws in a plane p erp endicular to the ciliary plane. (C) Absolute v alue of the
difference in the x -comp onen t of the v elo cit y along y = 0 (solid) and x = 0 (dashed) at
z = h .
a plane p erp endicular to the plane of the cilium, and Figure A.1C sho ws the absolute v alue
of the difference in the x-v elo cit y b et w een the calculations using the t w o approac hes, along
the lines x = 0 and y = 0 .
Figure A.2 compares the flo w fields generated b y a 5× 5 lattice of cilia in the t w o
configurations of in terest – the isotropic (Fig. A.2 ) and the fully sync hronized (Fig. A.3 )
states.
101
Figure A.2: (A) Sc hematic depiction of a lattice of cilia at random phase. (B) Flo w field
in the ciliary plane, flo w magnitude (colormap) and streamlines (white lines), based on the
con tin uum theory (left) and particle mo del (righ t), resp ectiv ely . The red arro ws indicate the
relativ e strength and orien tation of eac h Stok eslet. (C) Flo w field in a plane p erp endicular
to the ciliary plane.
102
Figure A.3: (A) Sc hematic depiction of a lattice of fully sync hronized cilia. (B) Flo w field
in the ciliary plane, flo w magnitude (colormap) and streamlines (white lines), based on the
con tin uum theory (left) and particle mo del (righ t), resp ectiv ely . The red arro ws indicate the
relativ e strength and orien tation of eac h stok eslet. (C) Flo w field in a plane p erp endicular
to the ciliary plane.
103
Abstract (if available)
Abstract
Motile cilia are fundamental micro-actuators in cellular biology that beat cyclically to drive flows. In dense ciliary carpets, such as those in the mammalian lungs, brains, and reproductive tracts, fluid transport emerges from the coordinated activity of thousands of multi-ciliated cells, each containing hundreds of cilia. Despite progress in analyzing cilia coordination in certain models, a general theory for coordination in the limit of large arrays of cilia remains lacking. Starting from discrete arrays of cilia wherein each cilium is represented by a well-known oscillator model, we devise a fast numerical algorithm for investigating the dynamics of thousands of hydrodynamically interacting cilia. We then develop a continuum theory in the limit of infinitely-many independently beating cilia by combining tools from active matter and classical Stokes flow. Isotropic and synchronized coordination states are unstable in both discrete and continuum models. Traveling waves emerge in both theories regardless of initial conditions, indicating that these emergent waves are stable global attractors, with wavelength and direction determined by the microstructure and activity at the cilium level. Our theory opens up the prospect of establishing causal explanations of the effect of individual cilia oscillations on large scale coordination and fluid transport in ciliary carpets.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Multiscale modeling of cilia mechanics and function
PDF
Evaluating sensing and control in underwater animal behaviors
PDF
Mechanical and flow interactions facilitate cooperative transport and collective locomotion in animal groups
PDF
Transitions in the collective behavior of microswimmers confined to microfluidic chips
PDF
Physics-based and data-driven models for bio-inspired flow sensing and motion planning
PDF
Dynamic modeling and simulation of flapping-wing micro air vehicles
PDF
Passive flight in density-stratified fluid environments
PDF
Traveling sea stars: hydrodynamic interactions and radially-symmetric motion strategies for biomimetic robot design
PDF
An experimental study of the elastic theory of granular flows
PDF
Hydrodynamic stability of two fluid columns of different densities and viscosities subject to gravity
PDF
Modeling and simulation of circulating tumor cells in flow. Part I: Low-dimensional deformation models for circulating tumor cells in flow; Part II: Procoagulant circulating tumor cells in flow
PDF
Modeling of multiscale continuum-atomistic systems using homogenization theory
PDF
Coordination in multi-muscle systems: from control theory to cortex
PDF
A stochastic Markov chain approach for tennis: Monte Carlo simulation and modeling
PDF
Investigating the debris motion during extreme coastal events: experimental and numerical study
PDF
Supersonic flow study of a capsule geometry using large-eddy simulation
PDF
Continuum modeling techniques and their application to the physics of soil liquefaction and dissipative vibrations
PDF
Continuum modeling of reservoir permeability enhancement and rock degradation during pressurized injection
PDF
Reconstruction and estimation of shear flows using physics-based and data-driven models
PDF
Development and applications of a body-force propulsor model for high-fidelity CFD
Asset Metadata
Creator
Kanale, Anup V.
(author)
Core Title
Synchronization in carpets of model cilia: numerical simulations and continuum theory
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Mechanical Engineering
Degree Conferral Date
2021-12
Publication Date
02/15/2022
Defense Date
07/28/2021
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
cilia,collective behavior,continuum theory,fluid mechanics,OAI-PMH Harvest,self organization,synchronization
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Kanso, Eva (
committee chair
), Luhar, Mitul (
committee member
), Nakano, Aiichiro (
committee member
), Newton, Paul (
committee member
)
Creator Email
anupvkanale@gmail.com,kanale@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-oUC16208270
Unique identifier
UC16208270
Legacy Identifier
etd-KanaleAnup-10170
Document Type
Dissertation
Rights
Kanale, Anup V.
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright. The original signature page accompanying the original submission of the work to the USC Libraries is retained by the USC Libraries and a copy of it may be obtained by authorized requesters contacting the repository e-mail address given.
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Repository Email
cisadmin@lib.usc.edu
Tags
cilia
collective behavior
continuum theory
fluid mechanics
self organization
synchronization