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
/
Multiscale modeling of cilia mechanics and function
(USC Thesis Other)
Multiscale modeling of cilia mechanics and function
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
MULTISCALEMODELINGOFCILIA MECHANICSANDFUNCTIONS Thesisby FengLing InPartialFulfillmentoftheRequirementsforthe Degreeof DoctorofPhilosophyinMechanicalEngineering UNIVERSITYOFSOUTHERNCALIFORNIA LosAngeles,California 2022 DefendedFebruary18,2022 © 2022 FengLing ORCID:0000-0002-1766-073X Allrightsreservedexceptwhereotherwisenoted Tothefutureofopenaccessknowledge. ii Acknowledgements First,IwouldliketothankmyadvisorDr.EvaKansoforintroducingmetothebeautifulmathematics and fluid dynamics of cilia motion and the wonderful world of biological physics, for her many lessons on communication and presentation, for her trust in letting me help out with writing grant proposals,andlastlyforthefinancialsupportthatenabledmycontinuationinremainingachildish adultforyetanotherhalfadecade. Next, I would like to thank Dr. Hanliang Guo for being my advisor-before-an-advisor during Eva’sinitialsabbatical,aswellasforbeingtheleaderbyexampleinourpack. Hismentorshipand encouragementshavebeeninstrumentalinshapingmyprofessionalisminscientificresearch. I would also like to thank Dr. Janna Nawroth for her patience, her biological insights and experimentalexpertisethatareproventobecrucialforourstillongoingcollaborations. I also need to thank all members and visitors of our lab1, for being continuous sources of intellectual exchange and support, for being the first ones to receive my never-ending complains, andfororganizingallthefunactivitiesespeciallyduringthelockdown. 1Innowayacompletelist(alphabeticalorder),MohamadAlsalman,BrendanColvert,HaotianHang,SinaHeydari, ChenchenHuang,YushengJiao,MorganJones,AnupKanale,JingyiLiu,Dr.YiMan,Dr.BasileRadisson,Tommaso Redaelli,WinfriedSchmidt,Dr.LionelVincent, etc. etc. iii Iamalsoforeverindebtedtomyadorableparentswhohadsupportedmeallthewayfromwhen I was just an mischievous baby, and all the great roommates (especially Hussein Hammoud and MichaelBarton)thathadtraveledwithmefordailychoresaswellasambitiousroadtrips. I wish to extend special thanks to all the professors that I have encountered during my stay at USC, especially my committee members (Dr. Paul Newton, Dr. Christoph Haselwandter, Dr. Ivan Bermejo-Moreno, Dr. Assad Oberai) and instructors that I was a teaching assistant for (Dr. Anita Penkova,Dr.JulianDomaradzki,Dr.TakahiroSakai). Theyhavehelpedmesolidifyoldknowledge and taught me many new concepts in diverse area of physics, mechanics, and numerical methods, while also trained me to become a more observant researcher and a more patient and effective teacher. Inclosing,Iwouldliketoexpressmygratitudetoallotherfriendships(gymbuddysandGKC) and higher-order relationships that I have inevitably failed to mention above, funding agencies that supported our projects (NIH, NSF, ONR, ARO)2, along with all the convenient modern tools (Google, YouTube, Stack Exchange, Github, ...) that have been indispensable for steering my professionalandcreativeinspirations. 2Research performed in this thesis are supported through the following grants: NIH R01 Grant 1R01HL153622- 01A1,NSFINSPIREGrant1608744,ONRGrant N00014-17-1-2062, ARO Grant W911NF-16-1-0074 iv TableofContents Acknowledgements iii ListofTables vii ListofFigures viii Abstract xi Chapter1: Introduction 1 1.1 Internalstructureofasinglecilium . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Hydrodynamicinteractionforciliacoordination . . . . . . . . . . . . . . . . . . . 6 1.3 Ciliainflowgeneratingorgans . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Chapter2: Mathematicalmodelsofcilia 11 2.1 Elastohydrodynamicsofasinglemicrofilament . . . . . . . . . . . . . . . . . . . 11 2.2 Couplingmolecularmotordynamicstociliaoscillations . . . . . . . . . . . . . . 15 2.2.1 Linearizationaboutthestraightequilibrium . . . . . . . . . . . . . . . . . 22 2.3 Phasedoscillatormodelofciliarycarpet . . . . . . . . . . . . . . . . . . . . . . . 24 2.3.1 Modelinganheterogeneousciliatedtissue . . . . . . . . . . . . . . . . . . 28 2.3.2 Characterizationoftheemergentwavepatterns . . . . . . . . . . . . . . . 30 2.4 Aunifiedmodeloffluidtransportinciliatedducts . . . . . . . . . . . . . . . . . . 33 2.4.1 Brinkmanpermeabilitycoefficient . . . . . . . . . . . . . . . . . . . . . . 38 2.4.2 Constraintonciliarymaterialandconservationofinputpower . . . . . . . 39 2.4.3 LinkingBrinkmanporositytoexperimentalobservations . . . . . . . . . . 40 2.4.4 Pumpingmetrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 2.4.5 AnalyticalsolutionsoftheBrinkmanmodelintwolimits . . . . . . . . . . 43 2.4.6 Non-dimensionalizationandparametervalues . . . . . . . . . . . . . . . . 45 Chapter3: 3Dspinningand2Dwhippinginstabilityofagenericmicrofilament 47 3.1 Linearstabilityanalysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 3.2 Spinningversusflappingoscillations . . . . . . . . . . . . . . . . . . . . . . . . . 51 3.3 Axialforceprofiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.4 Bead-springmodel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 Chapter4: Beatdiversityfrommolecularmotordynamics 75 4.1 Effectofchangingboundaryconditions . . . . . . . . . . . . . . . . . . . . . . . 75 4.2 Effectofforward-aftmotoractivityasymmetry . . . . . . . . . . . . . . . . . . . 80 v 4.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 Chapter5: Coordinationacrossheterogeneoustissue 89 5.1 Transportfirst,coordinatesecond . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 5.2 Patchycarpetcanenhancemetachronalorder . . . . . . . . . . . . . . . . . . . . 93 5.3 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 Chapter6: Universaldesignrulesforciliarypumps 98 6.1 PerformanceofCiliaryPumps . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 6.2 Highercilia-to-lumenratiogiveshigherpumpingpressure . . . . . . . . . . . . . 103 6.3 Optimalductdesigndependsonpressurerequirements . . . . . . . . . . . . . . . 105 6.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 Chapter7: Conclusionandfuture outlook 110 7.1 End-to-endmodelofciliarytissuefunctions . . . . . . . . . . . . . . . . . . . . . 111 7.2 Lifebeyondcilia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 Bibliography 113 AppendixA:3Dfilamentelasto-hydrodynamics 126 AppendixB:Two-state,two-filament2Dciliummodel 133 AppendixC:RotorhydrodynamicsusingregularizedStokeslets 137 AppendixD:Activeporousmediadrivenchannelflow 142 D.1 CFDsetup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142 D.2 Optimizationprocedures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 D.3 Optimalductdesigns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 vi ListofTables 2.1 Characteristicscalesofciliaandflagella . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Characteristicscalesofdyneinkinetics . . . . . . . . . . . . . . . . . . . . . . . . 19 2.3 Characteristicscalesofrotor-basedciliarycarpet . . . . . . . . . . . . . . . . . . 28 2.4 Characteristicscalesofinternalciliatedducts . . . . . . . . . . . . . . . . . . . . 45 vii ListofFigures 1 Scopeofthepresentthesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 1.1 Axonemestructure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.2 Themorphospaceofciliatedductsinnature . . . . . . . . . . . . . . . . . . . . . 9 2.1 Modelingciliumbeatsasphasedoscillators. . . . . . . . . . . . . . . . . . . . . . 25 2.2 Ciliarycarpetsandmetachronalorderparameters . . . . . . . . . . . . . . . . . . 27 2.3 Kuramotoellipse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 2.4 Porousciliamodelina2-Dchannel . . . . . . . . . . . . . . . . . . . . . . . . . 35 2.5 Ciliarymaterialconstraintandestimatesbasedoncollecteddata . . . . . . . . . . 43 3.1 Three-dimensionalbuckling duetotipfollowerforce . . . . . . . . . . . . . . . . 51 3.2 Transitionfrom3Dspinningto2Dflapping . . . . . . . . . . . . . . . . . . . . . 52 3.3 Frequencyofoscillationsversuslinearstabilityanalysis . . . . . . . . . . . . . . . 53 3.4 Bendingenergyversusaxial force . . . . . . . . . . . . . . . . . . . . . . . . . . 55 3.5 Axialforceprofiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.6 Galleryoffilamentoscillations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.7 Linearstabilityanalysisofaxiallyloadedfilament . . . . . . . . . . . . . . . . . . 61 3.8 Phasediagrambasedonpower-lawforcingprofiles . . . . . . . . . . . . . . . . . 63 3.9 Two-linkchain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.10 Linearstabilityoftwo-linkandthree-linkchains . . . . . . . . . . . . . . . . . . . 65 viii 3.11 Transitiontonon-spinningsolutionsofthebead-springmodels . . . . . . . . . . . 66 3.12 Asymmetricbeatingpatterns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 4.1 Galleryofciliumwaveformsduetouniformmotoractivity . . . . . . . . . . . . . 76 4.2 Motorswitch-inhibitionforclampedboundarycondition . . . . . . . . . . . . . . 77 4.3 Motorswitch-inhibitionforhingedboundarycondition . . . . . . . . . . . . . . . 78 4.4 Analysisofslidingfeedback controlofdyneinkinetics . . . . . . . . . . . . . . . 79 4.5 Linear stability analysis for sliding control mechanism under asymmetric dynein kineticsprofiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.6 Linearstabilityanalysisforcurvaturecontrolmechanismunderasymmetricdynein kineticsprofiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.7 Wavetravelingdirectionsunderasymmetricdyneinkineticsprofiles. . . . . . . . . 83 4.8 Waveforms observed in experiments and accessible under asymmetric dynein kineticsprofiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.9 Excitationofflagellawavereversalviabistability . . . . . . . . . . . . . . . . . . 85 5.1 Spontaneousemergenceofmetachronalwavestates . . . . . . . . . . . . . . . . . 90 5.2 Cumulativefluxfordifferentharmonicforcing . . . . . . . . . . . . . . . . . . . . 92 5.3 Patchinesscanenhancepumpingandmetachronalorderformation . . . . . . . . . 94 6.1 Porousmediapumpgeneralizestheenvelopemodel . . . . . . . . . . . . . . . . . 99 6.2 Ciliadensity,ciliatedfraction,andthematerialconstraint . . . . . . . . . . . . . . 100 6.3 Meanflowspeedasafunctionofadversepressureandlumenaspectratio . . . . . 101 6.4 Designinganoptimalciliarypump . . . . . . . . . . . . . . . . . . . . . . . . . . 102 6.5 ComputationalmodelthatmapsCDmorphologytofluidpumping . . . . . . . . . 103 6.6 OptimalCDdesigns . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 A.1 Geometry-preservingdiscretizationofanelasticfilament . . . . . . . . . . . . . . 127 ix A.2 CantileverdeflectioninStokes’regime . . . . . . . . . . . . . . . . . . . . . . . . 129 A.3 EulerbucklinginStokes’regime . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 A.4 HelicalbucklingandMitchellinstabilityinStokes’regime . . . . . . . . . . . . . 130 B.1 Robustnesstoinitialconditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 B.2 Convergenceofnumericalsharpwavefrontsolution . . . . . . . . . . . . . . . . . 136 C.1 Schematicofthequasi-periodiccomputationaldomain . . . . . . . . . . . . . . . 138 C.2 Convergenceofquasi-periodicfluidvelocity . . . . . . . . . . . . . . . . . . . . . 139 D.1 ValidationoftheBrinkman solver . . . . . . . . . . . . . . . . . . . . . . . . . . 143 D.2 Performancescoreatdifferentadversepressuregradient . . . . . . . . . . . . . . . 146 D.3 Optimaldesignsthatmaximizepumpingefficiency . . . . . . . . . . . . . . . . . 147 D.4 Pumpingefficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 D.5 Optimalperformancescoreandpumpingefficiencyovertheentiremorphospace . . 149 x Abstract Figure 1: Scope of the present thesis. This thesis attempts to answer several questions on the multiscale physics of cilia motion and their flow functions using a cascade of mathematical and computationalmodels. Cilia are micron-size, hair-like protrusions present in almost all Eukaryotic cells. Their motile variety-alsocalledEukaryoticflagellumwhenappearinginsmallnumbers-enablescellstodirect fluidflowaroundthemselvesbymovingandcoordinatinginwave-likefashions. Theoriginofcilia activityarisesfromacomplexbutrelativelyconservednanoscalestructurecomposedofmanysemi- rigidmicrofilaments(microtubuledoublets)thatslideduetothousandsofmolecularmotorproteins (axonemal dyneins). Notably, cilia can work individually as propellers for animal spermatozoa or xi single-celled flagellates, move together in single or a few pairs like in algae Chlamydomonas in someways analogous to appendages of mammals, or show up in large groups on microorganisms like Paramecium and Stentors or on epithelial tissue of mammalian brain ventricle, airway, and reproductivetracts(e.g.,airwayciliacanhaveadensityofuptoonemillionpermm 2 !). The physics behind ciliary motion have long fascinated biologists and fluid dynamicists since their discovery more than 300 years ago. Nevertheless, there remain debates on the precise mechanisms responsible for the spontaneous oscillation of individual cilium and the long-range coordinationofmanycilia. Perhapsmoreimportantly,ourunderstandingofciliaphysicsareoften segregated by the cascading length scales inherent for any ciliated tissue. There is thus a pressing needtodevelopmultiscaleanalysesthatcanincorporatekeyinsightsatthenano-andmicro-scales intofunctionalpredictionsatthetissue-andorgan-level. Inthisthesis,threesuchcontributionswillbeexpoundedindetail: (i)Usingtwo-statemolecular motormodelsandthecontinuousfilamentrepresentationfortheelasto-hydrodynamicsofcilia,we demonstrate how boundary conditions and distribution of activities along a cilium can change the geometry and wave traveling direction of spontaneously emerging cilia oscillations [1–3]. (ii) Reducingindividualciliamotionintophasedoscillators,wetheoreticallyandnumericallyquantify how key attributes of cilium beat (effective and recovery strokes, planarity of motion) combined withstructuralheterogeneityparametersatthetissue-levelalterthequalityandformoftheemergent metachronal waves and associated flow functions [4]. (iii) Coarse-graining cilia activity as active porous-media, we illustrate how the morphology of a ciliated duct predicts its relative ability to performbulktransportversusfiltration, andestablishuniversaldesignrulesforciliarypumps[5]. Frustratingasitmaybe,thisworkmarksonlytheexcitingbeginningofmodelciliatedsystems fromthenanometerscaleallthewaytomacroscopicfunctions. Combinedwithfurtheradvancesin xii imagingtechniquesatdifferentlengthandtimescales,itistheauthor’shopethatthenumericaland theoreticaltechniquesdevelopedherewillserveasgroundworkfortheresolutionoftheremaining mysteriesaboutcilia,andperhapseventhedevelopmentofdiagnostictoolsforcilia-relateddiseases inthenear future. MainContributions 1. Man, Y., Ling, F. & Kanso, E. Cilia oscillations. Philosophical Transactions of the Royal SocietyB375,20190157(2019). 2. Ling, F., Man, Y. & Kanso, E. Flagellar wave reversal via forward-aft molecular motor asymmetry.inpreparation(2022). 3. Ling, F., Guo, H. & Kanso, E. Instability-driven oscillations of elastic microfilaments. JournaloftheRoyalSocietyInterface15, 20180594(2018). 4. Kanale,A.V.,Ling,F.,Guo,H.,Fürthauer,S.&Kanso,E.Spontaneousphasecoordination inmodelciliarycarpets.inpreparation(2022). 5. Nawroth, J. C. et al. Ciliated Duct Morphology Determines Fluid Pumping Function. sub- mitted (2022). xiii Chapter1 Introduction This chapter introduces some background information on the physics and biology of cilia me- chanics from these perspectives: (i) what is known about the mechanisms driving the individual oscillation of an cilium, (ii) evidence that cilia coordinate via hydrodynamic interactions, and (iii) theimportanceofciliaryorgansintheanimalkingdom. 1.1 Internalstructureof a single cilium Cilia are are driven into oscillatory motion by an intricate internal structure, referred to as the central axoneme, composed of microtubule doublets and dynein molecular motors. Despite the highly-conserved structure of the central axoneme across eukaryotic cells, cilia and flagella of different cells exhibit wildly different beating patterns, from non-planar, cone-like, motions to planar, wave-like, deformations [6, 7]. The mechanisms that regulate the activity of the dynein motors,causingthemtoproduceoscillatorymotions,remainelusive,andevenlessisknownabout the mechanisms leading to this diversity in beating patterns [8–11]. Several theoretical studies support the hypothesis that dynein motors are regulated by a geometric feedback from mechanical deformationstomolecularactivity[12–15]. Brokawwasthefirsttoshowthatdelayedfeedbackfromcurvaturetodyneinactivitycouldlead to oscillations [12, 16–18]. Hines and Blum developed detailed models of elastic filaments that 1 generate sustained oscillations with curvature feedback control [19]. Murase, Jülicher and others demonstrated the existence of oscillatory modes in sliding-control models [20–23]. Lindemann proposed a ‘geometric clutch’ hypothesis, where the dynein activity changes as a function of the spacing between doublet pairs leading to axoneme oscillations [24–26]. A comparison of these various feedback mechanisms to experimental observations seems to favor the hypothesis of regulation bycurvaturefeedbackcontrol[15]. Although these feedback models are appealing, Bayly and Dutcher argued convincingly that evidence supporting the hypothesis that dynein regulation is required for bending oscillations remains circumstantial [27]. They then proposed an alternative mechanism that does not require regulation of dynein activity to generate oscillatory motions [27]. The mechanism relies on a dynamicbucklinginstabilityinducedbytheinternaldyneinforces,whichapplyaxialstressestothe axoneme as it deforms. This instability is reminiscent to the classic Euler buckling instability, but instead of the familiar static instability under fixed forces, these follower forces lead to a dynamic buckling instability, known as a ‘flutter’ instability in aeroelasticity [28] or, more informally, as the ‘garden hose’ instability [29, 30]. To demonstrate that this dynamic instability could lead to sustainedoscillationsinciliaandflagella,BaylyandDutcherusedelaboratemodelsoftheaxoneme, including finite element analysis and a model of two filaments, representing pairs of microtubule doublets, connected via passive elements of elastic spring and viscous damping. Han and Peskin also proposed an elaborate modelof the axoneme structure that lead to sustainedoscillations via a dynamicinstability[31]. DeCanio,Lauga,andGoldsteinexploredtheconceptofinstability-driven oscillations in the context of a simpler model, consisting of a single microfilament, confined to planarmotion,clampedatoneendandactedonbyafollowerforceatthefreeend[32]. 2 Eukaryoticciliaandflagellaaredrivenintooscillatorymotionbyanintricateinternalstructure, referred to as the axoneme (Fig. 1.1). The axoneme structure, composed of microtubule doublets and dynein molecular motors, is highly conserved across evolutionary time and cell type [33, 34]. However, the mechanisms that regulate the activity of the molecular motors, causing them to produce oscillatory motions, remain elusive. At present, there is no universal, experimentally- tested, theory for describing the active forces and moments generated in the axoneme that lead to sustainedciliaoscillations. Details of the structure of the axoneme of motile cilia were first delineated by transmission electron microscopy in the middle of the twentieth century [34]. The axoneme consists of nine microtubuledoublets,connectedtoacentralpairofmicrotubulesviaradialspokes,whichisknown as the “9+2” structure as shown in Fig. 1.1. Nexin links connect the outer microtubule doublets. Duringtheirpowerstroke,outerandinnerdyneinarmsbindtoneighboringmicrotubules,generating equalandoppositeforcesonadjacentmicrotubuledoublets. Unbindingofthedyneinstalkrequires energy – the energy released by the hydrolysis of ATP. The key unresolved problem is the spatial and temporal regulation of the binding and unbinding of molecular motors and its bearing on the oscillationsofthecilium. Several experimental approaches have been proposed to address this problem [35–44]. Rapid freezingofliveciliasamplesbysuddenremovalofATPresultinrigorwaves[36],andreactivation of these waves with reintroduction of ATP show wave propagation towards the flagellum tip as if oscillationshadnotbeeninterrupted[37,45]. Theseresults,inconjunctionwithotherexperiments on the functions of dynein motors [46], suggest that cilia oscillations are in tight control of the molecular machinery. In contrast, during experimentation with low concentrations of ATP in the algae model system Chlamydomonas reinhardtii, cilia were observed to bend at almost 3 Bridge Principal Reverse Dynein 1 2 3 4 5 6 7 8 9 Δ F + F − a + − e z e y e y e x Figure1.1: AxonemestructureLeft: thecrosssectionofatypicalcilium. Dyneinmotorsbetween microtubule doublets induce sliding, which gets resisted by elastic links. A symmetry-breaking bridge between two doublets separates the dynein motors into principal (P, red) and reverse (R, blue) sides. Right: model cilium consists of two elastic filaments ( + and−) coupled by elastic linkeranddamper,withdyneinmotorsexertingslidingforces( + and − )ondifferentsidesofthe cilium. ThebendingoftheciliumresultsinarelativeslidingΔbetween+and−filaments. constant curvature, suggesting a static mode that is distinct from the dynamic beating mode [43, 47]. Reactivation of the dynamic mode was possible with gradual increase of ATP. The relaxed configurationsobservedatlowATPandthetransitiontooscillatorymotionswithincreasinglevels of ATPsuggest thatoscillations could beinduced by amechanical instabilitythat does notrequire finespatiotemporalregulationofthemolecularmotors. Rapid freezing methods were recently coupled to powerful cryogenic electron microscopy to provide detailed structures of the conformations of the molecular motors [39–41, 44]. Snapshots of the dynein conformations at various instances of the cilium beating cycle were then correlated withtheoverallwaveformoftheciliumtodrawexperimentally-basedhypothesesontheregulation of the molecular motors in relation to cilia oscillations [44]. Asymmetrical activation of the dynein motors correlates with local curvature, suggesting that reciprocal inhibition is likely the main mechanism of dynein control at full waveform beating. These results support the notion that 4 oscillationsarefinelycontrolledbythemolecularmotors[44],buttheydonotruleouttheexistence ofan‘open-loop’instabilitymechanismthattriggerstheonsetofoscillations. Existingtheoriesdesignedtoshedlightonthemechanismscontrollingthemolecularmachinery either assume geometric feedback control from the cilium configuration to the molecular motor activity[21,23–26,31,48–55],or,morerecently,relyonsteadydistributedaxialforceswherethe dyneinactivityneednotbeaffectedbytheciliumconfiguration[3,27,56]. Geometric feedback theories come in three flavors: sliding control, curvature control, or geometricclutch. Thekeyideaintheslidingcontroltheoryisthattheactivityofmolecularmotors, and thus the active forces they generate, is regulated by the tangential sliding distance between two adjacent microtubule doublets caused by the bending motion of the axoneme [12, 22, 23]. The curvature control theory considers that the magnitude of the active forces is proportional to the curvature of the centerline of the axoneme [49, 50], whereas the geometric clutch approach considerstheactiveforcestobegovernedbytheseparationdistancesbetweenadjacentdoublets[25, 26, 57]. Comparison of these three feedback mechanisms suggests that curvature control gives best-fitstoexperimentally-observedwaveforms[58,59]. However, regulation of dynein activity may not be required to generate oscillatory motions [3, 27, 32, 56]. Oscillations arise as a result of a dynamic buckling instability – a Hopf bifurcation – induced by the axial stresses applied by the dynein motors on the axoneme. These mechanisms areparticularlyappealingbecausetheyprovideasimplerexplanationforsustainedoscillationsthat doesnotrequirefinetuningofthemolecularmotoractivity,buttheydonotexplainthedifferential dyneinbindingobservedin[44]. 5 Mathematicalmodels,whetherinsupportoffeedbackcontrolorinstability-drivenoscillations, areanabstractionoftheaxonemestructure,ignoringseveraldetailsinfavorofanalyticalrepresen- tations of the cilium centerline. These models are motivated by the fact that the “9+2” axoneme is characterized by a bridge that connects two adjacent doublets, labeled 5 and 6 in Fig. 1.1, thus dividing the dynein motors into principal (P) and reverse (R) sides relative to the bridge [15, 53]. Motors on the two sides operate antagonistically resulting in cilia oscillations. This “tug-of-war” aspect of the motor activity is minimally captured in the context of two elastic filaments ( +) and (−) of length , representing opposite sides of the axonemal bridge and separated by a constant distance at their base. The (±) filaments are subject to an active force doublet that results in an activemomentdensity onthecenterline,allowingittodeformintheplaneofmotion,asshown inFig.1.1. 1.2 Hydrodynamicinteraction for cilia coordination In multi-ciliated cells and tissues, neighboring cilia coordinate their phase to beat in the form of a metachronalwave [60–64]. Ciliacoordinationandmicroscaleflowscanbeobservedinmicrodis- sected ex-vivo epithelia [65–67] and engineered in-vitro tissues [68]. Despite the fundamental importance of cilia coordination for the function of major organs, the emergence of this striking metachronalcoordinationremainsopaque,asexperimentsthatprobetheemergenceandrobustness ofciliacoordinationatthetissuelevelareinherentlychallenging[68]. Itisoftenhypothesizedthathydrodynamicinteractionsarethedominantdrivingforcesforthis coordination [69–71]. Namely, as cilia each beat individually, they exert mechanical forces on the surrounding fluid that get transmitted to other cilia via the fluid medium. The importance of 6 this hydrodynamic interactions in the coordinated beating of cilia has been observed experimen- tally in pairs of cilia isolated from the somatic cells of Volvox carteri [69] and Chlamydomonas reinhardtii [70] and captured quantitatively in-silico [72] using a simplified representation of the cilium internal machinery. We thus posit that the exact details driving cilia oscillations matter little to multi-cilia coordination, and we seek a model that captures the dominant physics at the level of the single-cilium beat kinematics and mechanical interactions between cilia. The latter, in addition to the fluid medium, could be mediated by elastic coupling through the cytoskeleton [70, 73–76]. Thesesubstratemediatedinteractionsseemtopredominateciliacoordinationinunicellu- lar organisms [70, 73, 75], however, their importance for cilia coordination in multi-ciliated cells and tissues, besides their role in organizing and maintaining the geometric integrity of the ciliary lattice,remainsunderdebate[68]. However,todate,evenwithanabundanceofexperimental[68, 70] and theoretical [63, 67, 77–79] work, we still lack a coherent framework for establishing structure-to-functionmapsfromsingle-ciliumcharacteristicstotissue-levelcoordinationandfluid pumping. 1.3 Ciliainflowgenerating organs Ducts that drive luminal flows by the coordinated beat of internal cilia are ubiquitous in animal biology. In humans, ciliated ducts (CDs) pump fluids in the airways [80], brain ventricles and spinal canal [65], and oviduct and efferent ductules of the reproductive system [81, 82]. Their failure is directly linked to major pathologies, including bronchiectasis [83], hydrocephalus [84], and ectopic pregnancy [85]. Cilia in these ducts are short relative to the lumen diameter and oriented perpendicular to the epithelial surface, in a ciliary carpet design [86]. Many animals also 7 feature ducts with a strikingly different cilia arrangement, the so-called ciliary flame [87], where tightly packed, comparatively long cilia beat longitudinally along a narrow lumen. Ciliary flames that are thought to pump fluid for the purpose of excretion provide a model system for human kidneydisease[88,89]. Despite the fundamental importance of CDs in animal physiology, functional studies of intact CDsarerare[82,90],andtherelationshipbetweenCDmorphologyandtheirability,orfailure,to pump fluid remains largely unexplored. Ciliary flames and carpets have not been compared func- tionally, nor have other CD morphologies been characterized in relation to these two fundamental designs. Existing studies have focused on exposed ciliated surfaces and externally ciliated organisms because it is difficult to measure ciliary beat and fluid flow in intact internal ducts. Cilia os- cillations, metachronal coordination, and microscale flows have been observed in microdissected ex-vivo epithelia [67], engineered in-vitro tissues [91, 92], and, more recently, in organ-on-chip models [93, 94]. Additionally, leveraging the remarkable conservation of the ultrastructure of motileciliaamongeukaryotes,diverseprotistandanimalmodelsystemshaveemergedforprobing the functional spectrum of cilia, from locomotion and feeding [70, 95] to symbiotic host-microbe partnerships [66, 96]. However, we have yet to relate the architecture of intact CDs inside of or- ganismstotheirpumpingperformance. Establishingthisconnectionwouldelucidatetheapparent dichotomybetweencarpetsandflames,informourunderstandingofhumanciliopathies[97],aidin comparing different internal fluid transport mechanisms in animals [98, 99], and help in assessing thehypothesizedrolesofCDsinexcretion[88,89]andhost-microbeinteractions[100]. Using morphometric data collected from CDs in all 34 animal phyla [101], a total of 58 duct systems representing 26 animal phyla are found and plotted on the morphological design space 8 Figure 1.2: The morphospace of ciliated ducts in nature. A. CDs plotted based on cilia-to- lumenratioℎ/ andlumendiameter andclassifiedbytheirfunctionsandphyla. Here,triangles representsystemswhoseprimaryfunctionsarethoughttobe(ultra-)filtration,excretion,orpressure generation. Circles denote systems that transport and mix food particles, liquids, or cells. The squares represent systems whose fluid transport functions are currently unknown or disputed. The color code indicates the phylogenetic clades listed in the color bar. The dotted line is a linear fit indicating robust correlation between and ℎ/ in CDs (R-square goodness of fit = 0.68). B. PhylogenetictreeshowingthephylaknowntofeatureCDsandincludedintheanalysis(boldtype) and the phyla where CDs are absent, not documented, or shaped in a way that cannot be captured by the morphospace in (a).The color code is the same as in (a). Tree design adapted from [101, 102]; of lumen diameter and cilia-to-lumen ratio ℎ/, see Fig. 1.2. The remaining 8 phyla either lackedevidenceofmotileciliaorCDsor,inthecaseofsponges(phylumPorifera),wereexcluded because of their complex duct morphologies [103, 104]. Plotting ℎ/ versus for all 58 duct systems showed a well-defined 2D morphospace with a strong correlation between and ℎ/ (Fig.1.2A,bottom). Highcilia-to-lumenratios(flames)tendedtoassociatewithnarrowducts,and lowcilia-to-lumenratios(carpets)withwiderducts,whereassystemswithlarge andlargeℎ/, for example, were markedly absent. Mapping known biological functions onto the morphometric distribution showed a clear association of flame designs with filtration, and of carpet designs with 9 bulkfluidtransport,clearance,andmixing. Asthedistributionofsystemsinthemorphospacewas not reflected by their phylogenetic distance, this suggested a convergent evolution of CD designs basedonfunction. Intuitively,thedesignofCDsforfluidpumpingisconstrainedbythemetaboliccostofbuilding andactuatingthecilia[105,106],limitingthenumberandlengthofciliaonagivenciliatedsurface. Ahighercilia-to-lumenratioℎ/ necessarilyalsoimplieslessspaceintheductlumenthatcanbe taken up by fluid. Minimizing cilia cost while maximizing fluid volume seems to favor a ciliary carpet design at large lumen diameter. However, the ciliary flame designs observed in purported filtration systems suggest that other functional considerations are at play. Specifically, structural data suggest that flame cells achieve filtration by pumping fluid through a flow-resisting sieve [87, 107, 108], indicating that flame designs may be specialized to pump fluid in the presence of high resistance to flow, or, equivalently, adverse pressure gradients that cause counter-currents relative tothedesiredflowdirection. In order to test the validity of this hypothesis using fluid dynamics, we will formulate a mathematical model that is capable of representing both extremes of CD designs, from ciliary carpetstociliaryflames. 10 Chapter2 Mathematicalmodelsofcilia In this chapter, a set of mathematical models for cilia will be developed in detail. Each model areeffectiveabstractionssuitableforanalyzingciliaactivityatacascadeofdifferentlength-scales. Westartwithaquickdescriptionoftheelastohydrodynamicsofasinglemicrofilamentsubmerged in stokes fluid. This model can effectively describe the motion of an isolated microtubule or microtubule doublet under the sliding forces of molecular motors. Then we move on to coupling it to a two-state model of the dynein motors responsible for the spontaneous oscillations of cilia. Nextwediscusshowtoreducefilamentbasedmodelintosimplifiedphasedoscillatorssuitablefor simulationandanalysisatthescaleofthousandstomillionsofcilia. Finallyacoarsegrainedactive porous media model of cilia is proposed that enables unified analysis of all internal ciliated ducts fromciliarycarpetstociliaryflameswhereciliafillsthelumenofacavity. 2.1 Elastohydrodynamics of a single microfilament Consideraninextensibleelasticfilamentoflength andcircularcross-sectionofradius,clamped at its base normally to a fixed plane (,); see Fig. A.1(a). Let{e 1 ,e 2 ,e 3 } be an inertial frame, attachedatthefilamentbase,suchthat {e 1 ,e 2 }spanthe(,)-plane. Thecenterlineofthefilament is described by the vector r(,) where is arc-length and is time. The local orientation of the 11 Table2.1: Characteristicscalesofciliaandflagella. Parameter Symbol Dimensionalvalue Bendingrigidity 800 pN·m 2 [110] Filamentlength 20m Fluidviscosity ≈ water,20 ◦ 10 −3 Pa·s Timescale 4 ⊥ / 0.546 s Frequencyscale /( 4 ⊥ ) 1.83 Hz Forcescale / 2 2 pN Forcedensity / 3 0.1 pN·m −1 Slendernessratio / 1/100 Normaldrag ⊥ = 4 log(/) 2.73× 10 −3 Pa·s filament is described by an orthonormal material frame {d 1 ,d 2 ,d 3 }, with d 3 = r/ = t the tangentialunitvector. Thechangeof{d 1 ,d 2 ,d 3 } alongthefilament,atsometime ,isgivenby d =κ ×d , = 1, 2, 3, (2.1) whereκ is a generalized curvature vector recording the infinitesimal change in orientation along the filament; see, e.g., [109]. Specifically, given the rotation matrix Q(,) that maps vectors expressed in the inertial frame to their counterparts in the filament material frame, κ is the vector representation,inthematerialframe,oftheskew-symmetricmatrixκ × =Q(Q ′ ) T ,wheretheprime denotesdifferentiationwithrespectto and() T isthetransposeoperator. ThebalanceofforcesandmomentsonacrosssectionofthefilamentaregivenbytheKirchhoff equations[111],subjecttotheclamped-freeboundaryconditions, N +f ℎ +f = 0, M +t×N= 0. (2.2) 12 Here,NandMaretheinternalelasticforceandbendingmoment,respectively,f ℎ thehydrodynamic forceperunitlength,andf theaxialforceperunitlengthduetothemolecularmotors. Takingthe cross-productofthetangentvectortwiththebalanceofmomentsin(2.2),theinternalelasticforce canbereadilyrewrittenas N=t×(M/)+Λt, (2.3) whereΛ is a Lagrange multiplier that enforces the inextensiblity constraintt·t= 1. Physically,Λ representstheaxialtensionalongthe filament. We take advantage of the small aspect ratio of the filament ( /≪ 1) and use the resistive- force theory that linearly relates f ℎ to the instantaneous distribution of velocities v =¤ r along the filament centerline, where the dot denotes differentiation with respect to time . Specifically, f ℎ = − ∥ tt+ ⊥ (d 1 d 1 +d 2 d 2 ) visrelatedtovviatheanisotropicdragcoefficients ⊥ = 4/log(/) and ∥ = ⊥ /,with= 2,(see,e.g.,[112]),equivalently, f ℎ =− ⊥ I− − 1 tt v. (2.4) Toarriveat (2.4),weintroducedtheidentitymatrixI=d 1 d 1 +d 2 d 2 +tt. The dynein motors generate tangential forces along the filament. To model these effects, we consider a compressive force density (force per unit length) f =− t acting along the centerline of the filament, reminiscent to the models presented in [27, 32]. Here, the distribution of the force density may vary along the filament, = (), such that the total compressive force exerted by thedyneinmotorsisgivenby = ∫ 0 d. Bendingandtwistmomentsgenerateddirectlybythe dyneinmotorsareignored. 13 We consider a linear constitutive relation between the internal elastic moment M andκ . For anaxisymmetricfilamentwithcircularcross-sectionofarea andYoung’smodulus,thislinear relation simplifies to M= Bκ , whereB= diag(1, 1, 2) is the bending rigidity tensor, expressed in the filament material frame, and = 2 is the stiffness coefficient. It is worth noting that the energy needed to produce twist is of order / larger than that to produce bending, therefore we not only ignore the twist moments due to dynein activity, but also the tangential component of the internalmoment,thatis,wesetM·t= 0,whichyieldsnotwistofthefilament;see[113]. In the absence of twist, the generalized curvature κ can be expressed as κ = 1 d 1 + 2 d 2 . Here, 1 (,) and 2 (,) arescalarcurvaturefunctions,relatedtothecurvature(,) andtorsion (,) intheFrenet-Serretformulation,suchthat= √︃ 2 1 + 2 2 and=[tan −1 ( 2 / 1 )]/. Itis importanttodistinguishherebetweentwistandtorsion. ItisalsoworthnotingthattheFrenet-Serret frame(n,b,t),wherenandbarethenormalandbinormalunitvectors,isnotamaterialframe. In particular,{d 1 ,d 2 ,d 3 ≡t} isrelatedto(n,b,t) througharigidrotationaboutthet-directionbyan angleequalto tan −1 ( 2 / 1 ). Thedimensionalparametersandcharacteristicscalesrelevantforciliaandflagellaaresumma- rized in Table 2.1. Here, we rewrite the system of equations in non-dimensional form using the filament length as the characteristic length scale, / 2 as the force scale, and 4 ⊥ / as the characteristic time scale. The non-dimensional system of equations is characterized by the drag anisotropy parameter= ⊥ / ∥ = 2, as well as by the parameters that describe the distribution of axialforcesalongthefilamentcenterline. Inthefollowing,weconsiderallvariablesandparameters tobenon-dimensional. 14 2.2 Couplingmolecular motor dynamics to cilia oscillations Wecouplesthecenterlinedynamicsofciliatoamicroscopicmodelofthemotoractivityfollowing thehierarchicalmodeldevelopedin[114]andanalyzedin[115]. Considerasimplifiedrepresentationoftheaxoneme[15]. The“9+2”axonemeischaracterized by a bridge that connects two adjacent doublets, labeled 5 and 6 in Fig. 1.1, thus dividing the dynein motors into principal (P) and reverse (R) sides relative to the bridge. Motors on the two sides operate antagonistically resulting in cilia oscillations. This “tug-of-war” aspect of the motor activity can be minimally captured in the context of two elastic filaments ( +) and (−) of length, representing opposite sides of the axonemal bridge and separated by a constant distance at their base. The (±) filaments are subject to an active force doublet that results in an active moment density onthecenterline,allowingittodeformintheplaneofmotion. The cilium centerline is described byr(,), where is the arc length measured from the base and is time. The positions of the (±) filaments are given by r ± (,)= r(,)± 2 n, wheren(,) is the unit normal along the centerline. We also introduce the unit tangentt(,) to the centerline. In a Cartesian coordinate, say(,) whose origin is located at the base of the centerline, we write t = [cos, sin], and n = [− sin, cos], where tan(,) = / is the local slope of the centerline,andr=[(,),(,)]. 15 The equations of motion of the centerline that arise from balancing the tangential and normal forces and bending moments, including the moment due to motor activity, are given by [23, 57, 116] F +f h =0, (2.5) ++ = 0, (2.6) wheref h isthehydrodynamicforceperunitlength,F=t+nand aretheinternalforcesand momentactingalongthecenterline,with and denotingthetangentialandnormalcomponents of the internal force. The notation F = F/ and = / is used to represent the spatial derivatives withrespectto. Due to the slenderness of the filament ( ≪ ), we model the hydrodynamic force using the resistive force theory at low Reynolds number. That is, we consider f ℎ to be proportional to the localvelocitywithanisotropicdragcoefficients, f h =−( ∥ tt+ ⊥ nn)·r , (2.7) where r = r/ represents the time derivative. The drag coefficients satisfy ⊥ = 2 ∥ with ∥ ∼ 4/ln(/). A more accurate representation of the hydrodynamic forces in terms the slender-body theory, which includes the logarithmic corrections with non-local hydrodynamic interactions,wasconsideredin[115]. 16 Equation (2.5) can be written in scalar form as follows. First, multiply both sides of (2.5) with ( −1 ∥ tt+ −1 ⊥ nn) toobtainthefollowingrepresentationofthefilamentvelocity r = −1 ∥ ( − )t+ −1 ⊥ ( + )n. (2.8) Here, we have substitutedt = n andn =− t. Now, assume the filament is not extensible and applytheconstraintr ·r = 0 toget, uponintroducing= ∥ / ⊥ , −(1+) − − 2 = 0, (2.9) − −1 2 +(1+ −1 ) + = ⊥ . (2.10) Equation(2.6)issimplifiedfurtherbyconsideringalinearconstitutiverelationforthebending moment = ,where isthebendingrigidity. SubstitutingintoEq.(2.6),wehave ++ = 0. (2.11) Equations(2.9-2.11)provideasetofthreecoupledpartialdifferentialequationsthatweusetosolve for(,),(,),and(,),subjecttoproperlychosenboundaryconditions. Here,weconsider one end of the filament to be clamped to the cell body at = 0, and the other end to be free. The force-andmoment-freeboundaryconditionsat= aregivenby | = = 0, | = = 0, | = = 0. (2.12) 17 Atthebase,fromr (= 0,)= 0,weget ( − )| =0 = 0, ( + )| =0 = 0. (2.13) Theclampedendadmitsalsofilamentisclamped,wehave | =0 = 0. (2.14) Ifthefilamentishinged,thetotalmomentat = 0vanishes,leadingto | =0 − ∫ 0 ( ′ ,) ′ = 0. (2.15) ThefilamentdynamicscanthusbeobtainedbysolvingEqs.(2.9-2.11)subjecttothesixboundary conditions in Eqs. (2.12), (2.13) and (2.14). either (2.14) for the clamped case or (2.15) for the hingedcase. In the limit of small deformations(,), we see from (2.9) and (2.10) that ∼ () while ∼( 2 ). Therefore,thetensionisnegligible. Wecancancelthenormalforcebysubstituting Eq.(2.11)into(2.10),leadingtoasinglegoverningequationoftheform ⊥ =− −( ) . (2.16) 18 Table2.2: Characteristicscalesofdyneinkinetics Dimensional Dimensionless 50m lengthscale 0.9–1.7nN·m 2 forcescale 50ms timescale ⊥ 10 −3 – 1Pa·s Sp∼( ⊥ / ) 1/4 0.1–10 1–5pN ∼ 2 / 100–10000 0.5–2.5pN ∗ ∼ / 1.5–2 2×10 3 pN·m −2 ∼ 2 2 / 50–100 5–7m/s ∼/ ∼0.5 200nm = 0.1–0.3 10 3 m −1 = 1–100 withtheboundaryconditions | = = 0, | = = 0, (2.17) | =0 = 0, and | =0 = 0, clamped, | =0 = ∫ 0 ( ′ , ′ ) d ′ , hinged. (2.18) The active moment per unit length is generated by the longitudinal force doublets (,) perunitlengththatrepresenttheeffectsofthemolecularmotoractivity, (,)=(,). (2.19) 19 Asthecenterlinedeformsandbendsundertheinfluenceof ,itinducesrelativeslidingΔbetween thetwo(±) filamentswheretheinternalforces (,) areapplied;seeFig.1.1. Therelativesliding Δ= ∫ 0 [|(r − ) |−|(r + ) |] ˜ isgivenby Δ(,)=[(,)−(0,)]. (2.20) This sliding is resisted by nexin protein cross-linkers that act as a linear spring of stiffness . In otherwords, the internalforce density consistsof activeand passiveparts, arising fromboth the motoractivityandthepassiveresponseofthenexincross-linkers[114,115], (,)=( + + + − − )−Δ. (2.21) Here, is the average density of motors along both filaments, ± are the fractions of motors on the(±) filaments that are in the bound state, ± is the load exerted by a single motor. To close the model in (2.21), we must model the binding kinetics of the molecular motors as well as the motor loads. Forthebindingkinetics,weuseacommontwo-statemechanochemicalmodelconsistingof boundandunboundmolecularmotors,withaconstanttotal(boundorunbound)numberofmotors onbothfilaments. Asinglemotorcanbindtotheoppositefilamentatarate andunbindatarate . Weassumethatthemotorsswitchbetweenthesetwostatesstochastically;then,thefractions ± ofattachedanddetachedmotorsfollowtheFokker-Planckevolutionequations ( ± ) =(1− ± )− ± . (2.22) 20 In [23, 48, 117], the exchange rates and depend on periodic potential landscapes that govern theinteractionofthemolecularmotorswiththefilaments. Asimplermodelwasproposedrecently in [114] based on experimental measurements; in this empirical model, the binding rate= is constant and the unbinding rate exponentially increases with the load ± exerted by the motor and local curvature , i.e. = exp(± ± / ± / ), where is a constant and , are characteristic load or curvature for detachment. To model the load ± , we use a linear force- velocity relationship ± =± 0 (1∓Δ / 0 ), whereΔ is the sliding velocity and and are the stall force and associated zero load velocity at which the motors are at complete rest; see [114, 115]andreferencesthereinformoredetails. Here, we introduce an additional relation= , where() is a normalized dimensionless duty ratio that could depend on the arclength variable. This way we can control the likelihood of dynein attachment differently as along the flagellum. This additional complication is motivated based on the knowledge that dynein motor complexes can be either assembled starting from the basal direction or from the distal tip with the components transported via inter-flagella transport (IFT) [118]. Put together, Eqs. (2.22) that govern the fractions ± of motors bound to the(±) filamentsbecome ( ± ) = (1− ± )− ± exp 1∓ Δ ± , (2.23) andtheresultingforce (,) perunitlengthcanberewrittenas (,)= ( + − − )−( + + − ) Δ −Δ. (2.24) 21 Equations (2.19), (2.23) and (2.24) need to be coupled to the filament equations, Eqs. (2.9-2.11), via(2.20)toobtainaclosedsystemofequations. Non-dimensionalequationsareobtainedbycon- sideringthelengthscaleofthefilament andthetimescale = 1/( + ) ofthemotorkinetics. Specifically, we define the sperm number Sp 4 = ℎ / , where ℎ = ⊥ / is the hydrodynamic force and = / 3 is the elastic force due to bending. We also define the dimensionless active moment = / , where = / 2 is the elastic bending moment. Additional parameters include the stall to critical molecular force ratio ★ = / , the dimensionless sliding moment = 2 / , the ratio = / of the cilium diameter to the characteristic displacement due to motor activity, and the duty ratio = . Table 1 summarizes all dimensional and non-dimensionalparameters. 2.2.1 Linearizationaboutthestraightequilibrium Thedimensionlessnonlinearequationsareoftheform(using= ∥ / ⊥ = 1/2) 0= − 3 2 − − 1 2 2 , (2.25) S 4 = − 2 2 + 3 + , (2.26) 0= ++ , (2.27) =( + − − )− ( + + − )− , (2.28) ( ± ) =(1− ± ) (2.29) −(1−) ± exp[ ∗ (1∓ )]. 22 Whenthefilamentisclampedatthebase,thenonlinearboundaryconditionsare | =0 = 0, ( + )| =0 = 0, ( − )| =0 = 0, (2.30) | =1 = 0, | =1 = 0, | =1 = 0. (2.31) Assuming small displacement around the straight filament configuration ( = ), we can linearizeEq.(2.25)to(2.27)andobtainthefollowinglinearequationsforthefilament S 4 = , (2.32) =− (2.33) − [( + − − )− ( + + − )]. Sincetheequilibriumsolutionof ± is =/(+(1−) ∗ ),wecanobtaintheleadingorder equationformotorstatesbysubstituting ± = ± intoEq.(2.29)andget =−[+(1−) ∗ ]+(1−) ∗ ∗ . (2.34) Eq. (2.16) of main text can be thus obtained by substituting Eq. (2.33) into (2.32). But for the purpose of computing dominant eigenvalues shown in Fig. 4.4, we can linearize the boundary conditionsandget | =0 = 0, | =0 = 0, (2.35) | =1 = 0, | =1 = 0, (2.36) 23 fortheclampedcase,andsolveEq.(2.32)and(2.33)asamatrixeigenvalueproblemdirectlyafter appropriatediscretizationofthespatialderivatives. Once we obtain the set of eigenvalue and eigenvector pairs(,) , it is possible to extract the average bending wave propagation direction for each mode by testing the sign of the mean phase gradient ∠= ∫ 1 0 tan −1 I() R() d. (2.37) Namely,wavewouldtendtotravelsinthe+ base-to-tipdirectionifandonlyif ∠< 0[53]. 2.3 Phasedoscillatormodel of ciliary carpet To pump fluid at the micron scale where viscosity is dominant [119], individual cilia beat in a non-reversible manner. Specifically, the two key characteristics of an asymmetric cilium beat is that they usually (i) have an effective and recovery stroke each capable of generating different amount of hydrodynamic forces, and (ii) beat in a nearly planar fashion that are not symmetric whenlookingatatop-downview;seeFig.2.1A,D.Inordertofullyunderstandhowlargearraysof hydrodynamically interacting cilia self-organize into states that collectively break symmetry and pump flows, it is then convenient and even necessary to distill such key more microscopic, local descriptions into as few of variables as possible so that we can effective isolate the influence of all thepotentialcompoundingvariables. Tothisend,westartfromawell-knownmodelfortheindividualciliumasanonlinearoscillator (Fig. 2.1B,E, [121–123]) where the effective strokes and beat eccentricities are captured literally as Stokeslet rotors with tilted or eccentric trajectories, then try to fit their resultant flow field with those generated by Stokeslets that move in planar, circular trajectories (Fig. 2.1C,F). Specifically, 24 h effective h recovery <F effective > <F recovery > 2 1 0 0 π 2π F 2 1 0 0 π 2π F 2 1 0 0 π 2π F 2 1 0 0 π 2π F h effective h recovery b minor b major b minor b major A B C D E F b h Figure 2.1: Modeling cilium beats as phased oscillators. A. Individual cilia beats have effective andrecoverystrokesthatgeneratesasymmetricflowviaforcedifferenceandasymmetricbending, ignoring viscosity and complex fluid effects [120]. B. This symmetry-breaking can be effectively capturedbyatiltedrotor[61,69]. C.Toleadingorderusingaplanarcircularrotor,thefar-fieldflow effects of this asymmetry is captured by our first harmonic forcing. Note that despite this fitting, the planar circular rotor do not generate net flow when left alone. D. Individual cilia beats often have a dominant beating direction/plane. E. This asymmetry can be represented by allowing an rotor following an elliptic trajectory. F. To leading order using a planar circular rotor, the far-field floweffectsofthisasymmetryiscapturedbyoursecondharmonicforcing. wenumericallycalculatetheflowfielddueatilted(andeccentric,respectively)rotornearano-slip wall at a large set of Cartesian coordinate points, then perform a nonlinear least squares fit of that withtheflowfieldgeneratedbyaplanarcircularrotorwithanundeterminedforcemagnitudealong the trajectory. In other word, we fit the induced velocity (u,v,w) generated by a planar circular rotor with variable speed to velocities generated by the tilted or eccentric rotor by changing the force magnitude∥F∥(). As a result, we see that, to leading order, the effective/recovery strokes can be captured by a equivalent planar, circular rotor with a variable forcing that oscillates once along its rotation (Fig.2.1C), while the beat eccentricity can be captured by a variable forcing that oscillates twice along its rotation on a planar, circular rotor (Fig.2.1F). Thus, it is plausible and 25 effective to investigate the emergent coordination in ciliary carpets by studying the hydrodynamic interactionsduetothousandsofinteractingplanar,circularrotors. To this end, consider a periodic two-dimensional array of cilia of box length , wherein each cilium is a spherical bead of radius moving on a circular trajectory of radius , centered at x ( = 1, 2,...) on a planar square lattice of base distance located at a height ℎ above a wall where the Cartesian coordinates, span the ciliary plane and the perpendicular direction. The instantaneous position of the ℎ cilium is determined by its phase , or, in vector form, by r = x +n , where n ≡(cos , sin , 0) is the normal unit vector in the ciliary plane. Each ciliumisdrivenindependentlybyanactiveforceF = t actingalongthetangentunitvectort and whosemagnitude =( )isphase-dependentandstrictlypositive. Generically,thisforcingcan berepresentedbyaninfinitesumofforceharmonics ()= 0 1+ Í ∞ =1 cos+ sin . Thissimplerepresentationoftheciliumasalimitcycleexertingaphase-dependentforce() on thefluidmediumcanbeadaptedtoaccommodatearbitrary,includingexperimentally-derived,cilia beatforms [69]. Theciliaryplaneissubmergedina3Dfluiddomainofviscosity ,boundedat= 0,withfluid velocityu(,,) driven by the collective beating of all cilia and governed by the incompressible Stokes equations. For notational convenience, we introduce the 2D velocity v = u(,, = ℎ), whichistheprojectionofu ontotheciliaryplane. Balanceofforcesoncilium intheciliaryplanedictatesthat− ¤ t −v(r ) +F +N = 0, wherethefirsttermdenotesthedragforce(withconstantdragcoefficient = 6)andaccounts forhydrodynamiccouplingamongcilia;thelasttermisaconstraintforcethatguaranteesthebead remainsonthedesiredcirculartrajectory. Projectingtheforcebalanceonthenormalandtangential 26 Figure 2.2: Ciliary carpets and metachronal order parameters. (A) individual cilia are repre- sented by a force monopole driven to move along a circular trajectory by a phase-dependent force (); first and second force harmonics are shown below. An doubly-periodic array of cilia on a square lattice in the(,) plane above a wall is represented in the continuum limit by a force discontinuity layer (shown in blue). (B) Cosine of the phase field in synchronized and isotropic states.(C)Numericalsimulationsof 151× 151ciliashowtheemergenceoftravelingwavepatterns. Snapshots taken at = 100 for ℎ/ = 0.5 (left) and ℎ/ = 1.5 (right). (D) Time evolution of the Kuramoto polar order parameter shows no clear distinction between isotropic and traveling wave patterns. First force harmonic is shown in black and second harmonic in blue. (E) We define () for∈[0, 2] and plot(,()) in polar coordinates at two snapshots taken from the isotropic state (panel B) and traveling wave pattern (bottom left of panel C). Data points are shown in black and elliptic best-fit in blue. (F) Time evolution of the eccentricity and angle of the Kuramotoellipseshowsclearcorrelationsbetweenitseccentricitywiththegrowthofthetraveling waveinstabilityandanglewithdominantdirectionofwavefront. Firstandsecondforceharmonics are shown in black and blue respectively, and solid and dashed lines represent eccentricities and angles respectively. Here and in all figures, simulations are performed for 100 time units using = 0.1, = 0.05, = 0.2, = 1. First force harmonic: = 1+ 0.5 cos+ 0.5 sin, Second forceharmonic: = 1+ 0.5 cos 2+ 0.5 sin 2. Adaptedfrom[4]. directions,respectively,leadstoanexpressionfortheconstraintforce,N =−(n ·v(r ))n ,anda setofcoupleddifferentialequationsthatgovernthetime-evolutionoftheciliaphase, ¤ =Ω + 1 t ·v(r ). (2.38) Here, Ω = / is a phase-dependent intrinsic angular speed. In this model, isolated cilia do not pump fluid since they do not break time-reversal symmetry ((2.38) is invariant under the 27 Table2.3: Characteristicscalesofrotor-basedciliarycarpet Parameter Symbol Value Rotorspacing ∼ 10m Rotorfrequency Ω ∼ 30 Hz Fluidviscosity ≈ water,20 ◦ 10 −3 Pa·s Forcingcoefficient , ∥α +β ∥≤ 1 Rotor(minor)radius 0.4∼ 1 Rotorheight ℎ 0.1∼ 2 Periodicboxsize ≲ 151 Coverageratio 0∼ 1 Patchwavenumber ∈N Patchdisorder 0∼ 1 transformation→− for v= 0). This symmetry can be broken by tilting the ciliary trajectory withrespecttotheboundingwall;see, e.g.,[77,120]. Notethattheequationsofmotioncanbeintegratedaslongasanexpressionofhydrodynamically induced velocities due to other Stokeslets can be obtained. For computational efficiency, this will beobtainviatheregularizedStokesletmethodaboveawall. Formoredetails,pleaseseeC. 2.3.1 Modelinganheterogeneousciliatedtissue In order to study the effect of tissue-level heterogeneity on the coordination of ciliary carpet, it is also important to introduce abstract representation of tissue organization parameters in our rotor-based model. To this end, we first introduce the coverage fraction , patch wavenumber , 28 and a variance parameter that controls the randomness built into our simulated patches [67]. Specifically,wefirstconstructascalarpotentialfunction (,)=F 1 (+ ) , 1 (+ ) , 1 ()= 1− 2 2 2 − + 1 2 , (2.39) , ∼N(0,) whereF(,) isaunitnormseparablefunctionthatcontrolstheboundaryshapeofeachpatch. In general,F couldtaketheformof min(,),,or(+)/2forsquare,circular,ordiamondshaped patches. For simplicity, we setF(,) = min(,).the Combine that with 1 () representing a trianglewavethatvariesperiodicallyfrom−1to 1withaunitperiod,wecaneasilydeterminethat exactly percentofthepotentialfunction(,) hasitsvalue≥(1− 2). Inotherwords,wecan construct a periodic box of ciliary patches by simply discarding rotorsr whose value of(r ) is belowthecoveragethreshold 1− 2. Another way to introduce heterogeneity is through the introduction of incoherent beating into our system. For example, we can make the intrinsic frequency Ω of each rotor to follow a stochastic distribution. In addition, it is also reasonable to introduce misalignment of effective stroke directions by varying the weighting between and such that the ‘bumps’ of each rotor forcing also follow a distribution. However, we left the exploration of such variations in future workwithaeyeonexperimentally-deriveddistributionparameters. 29 2.3.2 Characterizationoftheemergentwavepatterns We next seek to characterize the properties of the desired metachronal wave state. Introducing the momentfieldsdescribingphasecoordination, (x,)= 1 ∑︁ e i (x−x ), (2.40) where = / 2 is the cilia density and is constant by construction. 1 is the Kuramoto order field. The average of 1 over the doubly-periodic domain leads to the Kuramoto order parameter =|⟨ 1 ⟩| = ∬ 1 (,) ; values of near zero indicate phase disorder while values near one correspond to phase synchrony [124]. In Fig 2.2D, we plot versus time for the two cases shown in Fig. 2.2C (left column): converges to a steady-state value that is independent of initial conditions. Surprisingly, does not capture the change in spatial patterns that arise as the system evolvesfromisotropictolong-termtravelingwaves. To remedy this, we introduce what we call the Kuramoto ellipse. We first compute the direc- tional Kuramoto order parameter = ∫ 1 ( cos, sin) d , where is a dummy integration parameter along a line segment originating from the center of the doubly-periodic domain at an angle to the-direction, and of length equal to the size of the fundamental domain. Intuitively speaking, measurestheaverageorderparameteralongthedirection. Notethat canalsobe seen as the absolute value of the Radon transform of 1 (x) [125]. We arrive at a polar coordinate representation(, ); black dots in Fig. 2.2E and enlarged in 2.3C. Their distribution is roughly circularintheisotropicstateandellipticinthetravelingwavestate. 30 We use principal component analysis to fit (, ) to a Kuramoto ellipse (blue curves in Fig.2.2E).Fig.2.2FshowsthetimeevolutionoftheeccentricityandangleoftheKuramotoellipse starting from the isotropic state and subject to first and second force harmonics. The Kuramoto eccentricity and angle are good indicators of the traveling wave patterns that emerge at steady state: theKuramotoangleisorthogonaltothewavedirectionanditseccentricityislargerformore coherentwaves. Forconcreteness,wedemonstratetheutilityoftheKuramotoellipseinFig.2.3,onanisotropic state(toprow)andasyntheticallygeneratedwavestate(bottomrow),respectively. Intheisotropic state,(, ) are distributed in an almost circular fashion around the origin, indicating lack of spatial order. In the wave state, since spatial order is maximum along the direction of wave orientation, the ellipse stretches out in that direction. We refer to the eccentricity of the ellipse as the Kuramoto eccentricity, and the inclination of the first principal component to the horizontal as theKuramotoangle. TheKuramotoeccentricityisameasureofthedegreeofcoherenceofwaves, while the inclinations of the first and second principal component give us the directions of wave orientationandpropagation,respectively. Another tool that is helpful in determining the degree of spatial coordination of our rotor dynamicsisthetwo-point(orrathertwo-rotor)correlationfunction. Consider Corr ()= ∑︁ (r −r )< 1 cos( − ), (2.41) where are rotor phases, and r the rotor center locations. This function thus measures the closesness ofthephaseanglesofrotorswithinadistance fromeachother. 31 Figure 2.3: Kuramoto ellipse. Using snapshots of an isotropic state (top) and a synthetically created wave state (bottom) to illustrate the usage of Kuramoto ellipse. A. The cosine field, obtained by taking the cosine of the phase of the cilium. B. The corresponding Kuramoto order field– Real ( 1 (x))– which is calculated by coarse graining over one shell of neighbors. C. Polar coordinates(,()) ingreen,andthecorrespondingKuramotoellipseingray. In addition to measuring the degree of cilia coordination, it is also of utmost importance to studytheperformanceofourrotorsystemwithregardtoitsflowfunction. Sincewemodelciliary carpet as rotors inside a semi-infinite domain, recall that a force monopole of strength F, oriented paralleltoano-slipwallatheightℎ,pushesfluidsatarategivenby q=ℎF/ [120,126]. Using the linearity of Stokes flow, it is straightforward to see that the total instantaneous flux due to the ciliary carpet is the sum of all cilia contributions, i.e. Í q (), where is the cilia index. Thus, thecumulativevolumeoffluidpumpedbythecarpet,orthe cumulativeflux ,untilatime, is Q()= ℎ ∫ 0 ∑︁ F () d. (2.42) 32 Thus, if the direction of pumping is fixed, we can compare the average pumping capability of different rotor systems by evaluating the magnitude ∥Q()∥ of the full system divided by the numberofrotorspresentinourperiodicdomain(e.g.,seeFig.5.2). As a sidenote, [67] used a different way to measure fluid transport by simulating Brownian particle displacement as advected by the generated flow. Then a particle clearance time T can be calculated by measuring how long it takes for the set of particles to translate a fixed distance, e.g., for the periodic box size. We do not resort to this measure here as our model has only a single no-slip wall and net flux can be readily calculated as shown above, and our main focus for now is onunderstandingthesynchronizationcharacteristicsoftheciliarycarpet. 2.4 Aunifiedmodeloffluid transport in ciliated ducts Cilia-powered flows have a long history in fluid dynamics; They can be linked back to the famous Taylor’s swimming sheet [127] and Blake’s impermeable envelope model [128], that have been used to study cilia-powered locomotion [76, 129], and cilia-powered fluid transport [8, 67, 120, 130]. In these models, the tips of densely packed cilia are represented as an impermeable sheet undergoing transverse and longitudinal traveling waves, without accounting for the details of the oscillatory motions of individual cilia. Models that account for the mechanics of individual cilia also have a long history of development [63, 66, 126, 131–135]. Modifying these models to account for the dramatic diversity in cilia beating patterns and complex ciliated duct geometries is cumbersome because they suffer from having a large number of parameters to be fit to (often unavailable)experimentaldataandhighcomputationalcostforlargenumbersofcilia. 33 Recently, [136] proposed a coarse-grained model of active filamentous structures in viscous fluids that account for both the elastic response and driving activity of the filament assemblies. In thesemodels,oneiscontendwithaneffectiverepresentationofthefilaments(withoutrepresenting individualfilaments)asan activeporousmediumthatdrivesandrespondstothesurroundingfluid. The two-way coupling between the filaments and the surrounding fluid is described by modifying theBrinkmanequationsforflowspastfixedporousmedia[137,138]. The Reynolds number Re = ℓ/, whereℓ is the cilium length, the fluid velocity, the fluid density and its viscosity, reveals the relative importance of viscous versus inertial forces in governing the fluid behavior. Both ciliary carpets and flames have Reynolds number on the order of∼ 10 −3 –10 −2 , even using fluid viscosities as low as water. Therefore, in both systems inertial forces are negligible and viscous forces are dominant, and the fluid motion is best described in the context of the incompressible Stokes flow model. Manipulating the surrounding fluid in this drag-dominatedmicroscaleisnotintuitive. Forciliarycarpets,theasymmetricpowerandrecovery stroke and metachronal waves break the time-reversible symmetry of Stokes flow, resulting in net fluid transport and pumping [112]. For ciliary flames, symmetry is broken by traveling waves propagating along the cilia length. Given the similarities in the kinematics of the traveling waves andeffectiveflowregimesinthetwociliaryphenotyes,wepositthatbothcarpetsandflamescanbe analyzedusingaunifiedfluid-structureinteractionmodelbasedonamodificationoftheBrinkman equationsthatdescribeflowsinporousmedia. Here, since our goal is to relate the morphological parameters of the cilia and ciliated duct to thecharacteristicsofcilia-generatedflowsinsidetheseciliatedducts,wederivedasimplercoarse- grained model based on one-way coupling between the cilia; we ignore the elastic response of the porous medium to flows and represent its activity by a prescribed traveling wave motion. This 34 h/2 h/2 H Adverse Pressure Gradient ΔP Viscous Fluid Flow Speed U Porous Ciliary Layer L xc (n, t) ρc (x, t) a b c 0.15 0 10 15 20 25 30 0.1 0 -0.1 -0.15 x 0 L 0.2L 0.4L 0.6L 0.8L x 0 L 0.2L 0.4L 0.6L 0.8L v c (x, t) ζ c (x, t) t t 5 0.05 -0.05 Figure2.4: PorousCiliaModelina2-DChannel. a. DiscreteLagrangianparticlerepresentation of the quasi-1D longitudinal density wave inspired by cilia motion. Each column of particles are labeled as (,), and relates the particle formulation to the continuum model. As time progresses, the traveling wave moves to the left, but generates net fluid motion to the right. (b-c). Cilia velocity (,) and Brinkman coefficient (,) as interpolated on the Eulerian grid, respectively. Every row of particles inside the porous layer follow the shown velocity and density profile. Here increase in darkness indicates the passage of time. = 0.1, and ===== 1. yields a tractable minimalist computational framework capable of capturing the essential aspects ofbothshort(carpet)andlong(flame)ciliabeatinginsideaduct. Consider a periodic cross-section along the longitudinal-direction of the ciliated duct, with length in the-direction and width in the orthogonal-direction. The fluid inside the duct is driven by porous layers of cilia growing from the two side walls of the rectangular channel. The fluidmotionisgovernedbytheincompressibleBrinkmanequation −∇ 2 u+∇+(u−v)=−Δ/·e , ∇·u= 0, (2.43) where is the fluid viscosity, u(,,) the fluid velocity field, (,,) the pressure field, and (,,) the Brinkman coefficient that relates permeability to fluid drag forces. The velocity field v(,,) represents the motion of cilia in the ciliary layers, and Δ/ on the right hand side represents a uniform pressure gradient in the adverse direction to the desired flow direction. 35 We consider periodic boundary conditions in the longitudinal -direction and no-slip boundary conditionsinthe-direction, u(0,,)=u(,,), (0,,)= (,,), u(, 0,)=u(,,)=0. (2.44) The minimalist form of activity that pumps fluid in such a channel is to abstract cilia motion intoalongitudinalwavethattravelsinthe-directionbutremainsinsyncalongthe-direction;see Fig.2.4. Inthisquasione-dimensionalsetup,ciliaareonlypresentintheregions∈(0,ℎ/2] and ∈[−ℎ/2,),soitishelpfultodefinethefunction ()= 1−H(−ℎ/2)+H(−+ℎ/2), whereH is the Heaviside function (H() = 0 for < 0 and 1 otherwise), to decompose the cilia drivenvelocityfieldandBrinkmancoefficientas, v(,,)= (,) ()e , (,,)= (,) (), (2.45) where (,) and (,) are the one-dimensional velocity and Brinkman coefficient due to the longitudinalciliarywaves;seeFig.2.4a. To compute (,), let (,) denote the Lagrangian label of the ℎ -column of particles, where is an integer from 0 to /, and is the average spacing between each particle in the -direction. Consider that individual particles oscillate with amplitude and frequency, and, 36 with no loss of generality, they coordinate with wavelength that is identical to the channel length . ThustheLagrangiandescriptionoftheparticlesmotionisgivenby (,)=+ cos(2/+). (2.46) Inthecontinuumlimitof→,weget (,)=+ cos(2/+). (2.47) Information regarding inter-cilium spacing is preserved by defining a uniform packing density of ¯ (,)= 1/ insomereferenceconfiguration. As cilia oscillate, their total number inside the periodic channel is conserved implying that ∫ ¯ d = ∫ d . Therefore, the cilia packing density ( ,) in the current configuration shouldsatisfy ( ,)= ¯ −1 = ¯ 1− 2 sin(2/+)/ . (2.48) Snapshotsof / ¯ isshowninFig.2.4b,c. Similarly,theLagrangianvelocity ( ,) inducedby ciliaoscillationis ( ,)= =− sin(2/+). (2.49) with the understanding that = ( ,) on the right-hand side of both (2.48) and (2.49). This inversemapfromLagrangian toEuleriancoordinatescanbeobtainedbyinverting(2.47). Since (2.47) is transcendental in, closed-form expression of this inverse map is not readily available. 37 However, such map always exists whenever (2.48) is non-singular, which is possible as long as < /2. Practicallyspeaking,thisinversemapcanbeapproximatedeitherviaaseriesexpansion for small , or through numerical interpolation for finite . Considering the latter, we calculate (,) and ( ,) separately in a dense enough grid and interpolate for (,) at given (in MATLAB,weusev=interp1(x_c,v_c,x));seeFig.2.4b. 2.4.1 Brinkmanpermeabilitycoefficient WederiveanexpressionfortheBrinkmancoefficient (,) asfollows. Wefirstestimateitsvalue in the reference configuration (Fig. 2.4). Here, each cilium is represented as a stationary array of equally-spaced particles of radius that span a distanceℎ/2 in the -direction with horizontal packingdensity ¯ = 1/;seeleftmostpanelofFig.2.4. ConsistentwiththeBrinkmanmodel,the total drag caused by the motion of all particles is equal to the superposition of drag forces due to eachindividualparticle. Thusthetotaldragforceinsidetheductduetociliamotionis =( ¯ )( ℎ 2 )(6 )= 3 ¯ ℎ . (2.50) Here, ¯ is the total number of cilia in the duct, andℎ/2 the number of "pores" or particles that represent each cilium. From (2.43), it is clear that the Brinkman coefficient converts velocity to pressuregradientandhasunits[(Pam −1 )·(ms −1 ) −1 ]. Inthereferenceconfiguration,wehave ¯ = 1 1 ℎ = ¯ (3 ¯ ℎ) ℎ = 3 ¯ 2 , (2.51) 38 where we first divide the drag force coefficient by the total ciliated area ℎ to get a pressure, then bytheinter-ciliumspacing= 1/ ¯ togetapressuregradientacrossthemeandistancewherethe dragforceisapplied. Given (2.51), the time varying Brinkman coefficient in the current configuration follows from theconservationlaw ∫ d = ∫ ¯ d, (,)= ¯ −1 = 3 ¯ (,). (2.52) Substitutinginto(2.45),wearriveatanexpressionfortheBrinkmancoefficient (,,). Snapshots of the Brinkman coefficient (,,) as well as the corresponding coarse-grained velocity vector field v(,,) areshowninthetworightmostpanelsofFig.2.4. 2.4.2 Constraintonciliarymaterialandconservationofinputpower The operation of CDs is accompanied by a metabolic cost associated with building and actuating thecilia. Here,wetranslatethisintuitionintoourmathematicalmodel. We emphasize that, in our quasi-1D model, the ciliary layers are treated as particles of radius ,fixedthroughouttobethesameasthetypicalradiusofacilium, = 0.125m,packedtightlyin the-direction and with center-to-center distance in the-direction in a reference configuration (Fig.2.4). Thus,thetotalactivematerialinaciliatedductisgivenby( ¯ )(ℎ/2). Bydefinition, the total material in the ciliated layers is conserved. Considering the radius of individual cilia is conserved,thisconstraintontheciliarymaterialisequivalentto ¯ ℎ= constant. (2.53) 39 A simple calculation shows that this constraint is analogous to a constraint on the ratio of cilia volume to duct surface area. The volume of an individual cilium is 2 ℎ, while ¯ represents the number of cilia per unit duct surface area (in our model, ¯ is the cilia count per unit duct length). Thus,alimitontheratiooftotalciliavolumetoductsurfaceareais ¯ ( 2 ℎ),which,considering constant,isequivalenttoalimitonthematerialconstraint ¯ ℎ. Tocompareciliatedductsthatdemandthesamegrowthandactuationcostsoftheciliarylayers, we fix the total ciliary material in the duct across simulations by setting ¯ ℎ constant, and we considerciliarywaveswithfixedamplitude ,wavelength= 1/,andfrequency. Thisimplies that the velocity of the ciliary waves in (2.49) is fixed and the total drag force in (2.50) generatedintheciliarylayersisalsofixed. Thisisequivalenttoconservationoftotalpowerinput, whichisdiscussedfurtherin§D.2. Itisinstructiveforthediscussioninthemaintexttovisualizetheeffectofthematerialconstraint ¯ ℎ= constant over the entire morphospace(,ℎ/); see Fig D.5A. Given a constant amount of ciliarymaterial( ¯ ℎ=constant),Ductswithhighercilia-to-lumenratioandlargerlumendiameter havelower longitudinalreferenceciliadensity ¯ comparedtonarrower,morecarpet-likeducts. To estimate the amount of ciliary material in the duct, that is, to estimate the value of the constantin(2.53),welinktheBrinkmanporositytoexperimentalobservations,asdiscussednext. 2.4.3 LinkingBrinkmanporositytoexperimentalobservations In the classic 2D Brinkman model, porosity is defined as the area ratio = fluid / total of the fluid fluid to the total area total . In our quasi-1D model, the area porosity is given by = 1−(ℎ/2)( 2 )/(ℎ)= 1−(/2)(/). Recallingthatbydefinitionthereferencedensityof 40 cilia ¯ = 1/, we get the following relationship between the Brinkman porosity and reference ciliadensity ¯ = 1− ¯ 2 , or, equivalently, ¯ = 2(1−) . (2.54) This is useful for expressing the constraint on ciliary material (2.53) in terms of . Noting that and are held constant in all simulations, we arrive at equivalent expressions of the material constraint, ¯ ℎ= constant, or, equivalently, (1−)ℎ= constant. (2.55) When comparing model and data, we match with the porosity obtained from measurements as follows. For ciliary carpets, the density of cilia within the ciliary layer is usually independent of the cilia-to-lumen ratio ℎ/ and can be spatially heterogeneous. We used a baseline estimate of= 50% to account for both inter-cilium spacing and patchiness of such systems, and we used the cilia radius = 0.125 m. For CDs with cilia traversing the entire lumen, which is common for longitudinally oriented cilia, the measured area ratio of cilia to total lumen is interpreted both as a direct estimate of the complement of porosity 1− as well as the dimensionally adjusted cilia-to-lumenratio(ℎ/) 2 . Thenusingciliumradius,weestimatethevaluesof ¯ andmaterial constraints ¯ ℎ (Fig.2.5B). 41 2.4.4 Pumpingmetrics The effectiveness of a pump is determined by both the quantity of flow it can generate and the efficiency at which it can operate. Since net flow is only possible in -direction of the channel, we firstdefinetheflowprofile ()= 1 ∫ 0 u(0,,)·e d= 1 ∫ 0 u(,, 0)·e d, (2.56) where = 2/ and the equality follows from the periodicity of the applied forcing and channel boundary conditions in. The net flow of the channel is thus =⟨ ⟩ , where⟨·⟩ symbolizes averagingalong-directionoverchannelheight. Fromthiswecandefinetheareaflowrate as =. (2.57) Next,wedefinetheefficiencyofthechannelpumpasthepowerratio = rateofchangeoffluidmechanicalenergy powerinputbyciliarylayers = 1 2 ⟨¤ ⟩ 2 ⟨ ⟩ ⟨ ⟩ = 2 3 ¯ ⟨ ⟩ 2 (2.58) where ¤ isthemassflowrate( isareadensityofthefluid)and ⟨·⟩ symbolizesaveragingalong -directionofthechannelasineq.(2.56). Theaveragespeed⟨ ⟩ oftheprescribedtravelingwave isgivenby ⟨ (, 0)⟩ = 1 ∫ 0 ( , 0) d = 2 /. (2.59) 42 Figure 2.5: Ciliary material constraint and estimates based on collected data. A. Effect of material constraint ¯ ℎ = constant on the longitudinal cilia density ¯ . The density ¯ of the ciliary material decreases as the lumen diameter increases. B. Histogram of estimated ¯ ℎ using data collected in Suppl. Data. Here we use(ℎ/) 2 the measured cilia to lumen area ratio to estimatereferenceciliadensitywheneverappropriate,andassumesabaselineciliadensityof50% forcarpet-likedesignsifnoimagesareavailable. Puttogether,thepowerefficiency becomes = 3 ¯ ℎ . (2.60) where = 6 2 2 4 −1 isaconstantcoefficientinalloursimulations. Notethatbyconstruction, the power input is conserved in all simulations that satisfy the material constraint ¯ ℎ= constant, thus,theefficiencyisproportionaltothepoweroutput;specifically, ∝ 3 . 2.4.5 AnalyticalsolutionsoftheBrinkmanmodelintwolimits Wederiveanalyticalsolutionsintwolimits: thelimitofinfinitesimallyshortciliarycarpet( ℎ→ 0 and ¯ →∞) and the limit of a duct completely filled with flame cilia ( ℎ= ) in the absence of adversepressuregradients. 43 As ℎ→ 0 and ¯ →∞, the ciliary layers collapse into infinitesimally thin ciliary carpets that move according to the prescribed longitudinal density wave, such that the governing equation reduces to Stokes flow with slip velocity boundary conditions at the top and bottom walls. This scenarioisconceptuallyequivalenttothestudyofciliaryflowbasedontheenvelopemodel[128]. Since we have an expression (2.59) for the average flow generated by , linearity of the Stokes equation implies that, with no adverse pressure gradients, the mean flow speed inside the entire channelwouldalsobeexactlythatoftheprescribedtravelingwave− 2 −1 . Ontheotherhand, whenℎ= 0 and ¯ is finite, no cilia is present and the no-slip boundary condition implies no flow canbegenerated. If we haveℎ=, the entire inner duct is subjected to Brinkman flow with no-slip boundaries at thetop and bottom. Thus we havev(,,)= (,)ˆ e for all∈(0,). It is convenient here to introduce the stream function such thatu=∇×. The incompressibility constraints is thus satisfied automatically by , and the governing equation can be rewritten in terms of by taking curlof (2.43)toarriveat −∇ 4 +∇ 2 +∇×v= 0, (2.61) subjecttotheboundaryconditions∇×= 0at= 0and=. Sinceourchoiceoflongitudinal wave gives∇×v = 0 when ℎ = , we get that the fluid velocity is identically zero, u = 0 everywhere,inabsenceofadversepressuregradients. Whenadversepressureispresent,itreduces tothepressuredrivenPoiseuilleflowproblemin2Dperiodicchannelwithanadditionalresistance contributedbytheporouscilialayer,whichadmitsaclosedformsolution. 44 Table2.4: Characteristicscalesofinternalciliatedducts Parameter Symbol Value Metachronalwavelength 100m Wavefrequency Ω 15∼ 30 Hz Waveamplitude /2 Fluidviscosity ≈ water,20 ◦ 10 −3 Pa·s Cilia-to-lumenratio ℎ/ 0∼ 1 Ductlumendiameter 5∼ 200m Adversepressure Δ 0∼ 100 Pa Materialconstraint ¯ ℎ 1∼ 1000 meanflowspeed = ∫ 0 u·e ≤/(4) 2.4.6 Non-dimensionalizationandparametervalues OurcontinuumBrinkmanmodelin(2.43)involveseightindependentdimensionalparameters: : viscosityofthefluidmedium, : widthofthechannel(lumen diameter), ℎ: combineddepthoftheBrinkmanlayers,ℎ∈(0,), ¯ : referenceciliadensity; ¯ = 1/,where isthelongitudinalinter-ciliumdistance, : lengthofthechannel;= 1/,where isthewavenumberoftheciliarywave, : waveamplitude,∈(0,/2), : wavefrequency, Δ: adversepressure. We write (2.43) in non-dimensional form by choosing as our characteristic scales the viscosity of the fluid medium , the duct length, which we set to be equal to the ciliary wavelength, and the ciliabeatfrequency(CBF). Wewritethepressureandvelocityindimensionlessform, e = , e u= u , e v= v . (2.62) 45 Wealsowritetheremainingparametersindimensionlessform, e = , e ℎ= ℎ , e = , e ¯ = ¯ , e = 2 , f Δ= Δ −1 . (2.63) Thenon-dimensionalcounterpartoftheBrinkmanequation(2.43)isgivenby −∇ 2 e u+∇e + e (e u− e v)=− f Δ·ˆ e , (2.64) We used (2.64) in our numerical simulations (§D). To convert between dimensionless and dimen- sional quantities, we used by default the viscosity of water = 10 −3 Pa·s, CBF = 15 Hz, and a ciliary wavelength of = 100 m unless noted otherwise. In all simulations, we set the wave amplitude to = /(2). We varied the cilia-to-lumen ratio ℎ/ in increments of 5%, and we examinedlumendiameters rangingfrom 1–1000m,andmaterialconstraints ¯ ℎrangingfrom ¯ ℎ= 1to 1000(Fig.2.5). 46 Chapter3 3D spinning and 2D whipping instability of a generic microfilament Inthisstudy,weconsideranelasticmicrofilamentinaviscousfluid,clampedatoneendandacted uponbyadistributionofaxialforces. Unlike[32],thefilamentisfreetoundergothree-dimensional (3D) motions. To investigate the filament deformations in 3D, we adapt the numerical framework advanced by [139, 140], as well as a combination of linear stability theory and a low-order bead- spring model. We show that the filament does indeed exhibit a Hopf bifurcation that leads to sustained oscillations, consistent with [32]. However, unlike [32], the filament undergoes non- planarspinningmotions,reminiscenttothecone-likebeatingmotionofcilia. Importantly,atlarger axial forces, we show the existence of another bifurcation, not reported previously, that causes the filament to transition from 3D spinning to planar wave-like deformations. We investigate the transitionfrom3Dspinningto2Dflappingundervariousaxialforceprofilesandinthebead-spring model. We conclude by commenting on the significance of these results to biological cilia and flagella. 47 3.1 Linearstabilityanalysis Considerastraightfilamentsubjecttoadistributionofaxialforces f =− ()t,alongthefilament centerline, where () is a general function of . In order to assess the stability of the straight configuration,welinearizetheequationsofmotionaboutthestraightequilibriumandsolveforthe dominanteigenvalue. To begin, we substitute (2.3) and (2.4) into the force balance in (2.2) to arrive at the vector equation ¤ r=(Λ ′ + ′ − )t+(Λ 2 − ′′ 2 )d 1 −(Λ 1 − ′′ 1 )d 2 , (3.1) where we used the fact that = √︃ 2 1 + 2 2 . Equation (3.1) leads, upon further simplifications (see details in the supplemental document), to three nonlinear scalar equations in terms of the filament tensionΛ(,) andfilamentshape: curvature (,) andtorsion(,). Expandingtheseequations about the straight equilibrium state = = 0, and linearizing with respect to and, it can be shownthatabovevectorequationleadstotwoscalarequationsonly, Λ ′′ = ′ (3.2) ′′′′ −Λ ′′ +¤ = 0. (3.3) The third equation is trivially satisfied at the linear level. Since these equations do not depend on torsion(,),thelinearanalysiscannotcapture3Ddeformationsofthefilament. Thatistosay,in thelinearregimenearthestraightequilibrium,thefilamentundergoes planardeformations. Thus, withoutlossofgenerality,wecanequatethearc-length withthe coordinateintheinertialframe 48 andset= ′′ ,wheretheprimenotationheredenotesdifferentiationwithrespectto . Substituting into (3.3) and integrating with respect to twice, we arrive at the fourth-order partial differential equation, ′′′′ −Λ ′′ +¤ = 0, (3.4) for the infinitesimal displacement (,) in the(,)-plane. To close this equation, we integrate (3.2)fromthefreetipofthefilamentanddefinetheinternaltractionforce Λ()= ∫ ( ˜ ) d˜ . (3.5) Assume the solution(,) is separable in the form of(,) =() , we obtain the boundary valueproblem ′′′′ −Λ ′′ + = 0, (0)= ′ (0)= 0, ′′ ()= ′′′ ()= 0. (3.6) The first two boundary conditions correspond to the clamped end and the latter two to the free end. Ingeneral,thedifferentialequationin(3.6)belongstotheclassoflineardifferentialequations with non-constant coefficients, given that Λ() is a function of for a general distribution of axial forces. Therefore, depending on the form of (), this boundary value problem may not admit an analytical solution. The boundary conditions welcome the use of the free vibration modes of 49 a clamped-free beam as test functions for a numerical solution of this boundary value problem. Namely,weexpand() ontotheeigenfunctionsoftheclamped-freelinearbeamequations ()= ∞ ∑︁ =1 Φ (), (3.7) where areconstantcoefficientsand Φ () aregivenby Φ ()= ∞ ∑︁ =0 cosh( )− cos( ) − [sinh( )− sin( )], (3.8) with =(cos + cosh )/(sin + sinh ) and the-throotofthetranscendentalcharacter- isticequation cos cosh + 1= 0. This“assumed-modes”method,anextensionoftheclassical Ritzmethod[141,142],wasusedtoanalyzethestabilityofflagelladynamicsin[27,53]. Finally,weproject (3.6)ontothesebasisfunctionsandrewriteitinweakformas ∫ 0 Φ ′′ Φ ′′ d− ∫ 0 Φ ΛΦ ′′ d + ∫ 0 Φ Φ d = 0. (3.9) This equation can be rewritten in matrix form as(K+C+M)Y= 0. Thus, the linear stability of the filament away from the straight configuration can be obtained via the eigenvalue problem K+C =−M. The real part of the resulting indicates the growth (or decay) of pertubations andtheimaginarypartcorrespondtothefrequencyofagivenmodalshape Φ (). 50 0 0.02 0.04 0.06 -0.5 0.0 0.5 1.0 Time Time Bending Energy (b) (c) (d) 0 0.02 0.04 0.06 0 1 2 3 4 F a s (a) Spinning (locked curvature) Position Figure3.1: Three-dimensionalbucklingduetotipfollowerforce. (a)Filamentsubjecttoanaxial followerforceof ,concentratedatitstip. (b)Startingfromastraightconfiguration,thefilament buckles into a locked-curvature configuration and undergoes 3D spinning motion for = 75. (c) The,,and positionsofthefilamentdistaltippointand(d)itsbendingenergyversustime. 3.2 Spinningversusflapping oscillations We first consider the set-up analyzed in [32] of an elastic filament subject to a follower axial force of magnitude concentrated at the distal tip of the filament; see figure 3.1(a). In [32], the filament is confined to undergo planar motions. Starting from a small perturbation away from the straightconfiguration,theauthorsidentifythreedynamicalbehaviorsdependingonthemagnitude of the axial force : (i) for small force, the filament returns monotonically to its original straight configuration, (ii) as the force value increases, the filament displays decaying oscillations back to the straight configuration, and (iii) for forces above a given threshold, the filament settles into finite-amplitudeperiodicdeformations,whichtheyreferredtoas planarflapping . Here,westartwithasmallnon-planarperturbationawayfromthestraightconfigurationandwe numericallysolveforthefilamentdynamics. Infigure3.1,weset = 75,whichliesintheplanar flapping regime of [32]. Unlike the flapping behavior reported in the latter, the filament buckles into a configuration with locked curvature and undergoes three-dimensional spinning about the - directionasillustratedinfigure3.1(b). Thetipdisplacementofthefilamentanditsbendingenergy 51 as a function of time are reported in figure 3.1(c) and (d), respectively. The (,) coordinates of the filament tip settle into periodic motions while the -coordinate goes to a constant value. The filament bending energy also goes to a constant value, emphasizing that, in its buckled state, the curvatureofthefilamentisconstant(seesupplementalmovieS1). -0.6 -0.4 -0.2 0 0.2 0.4 0.6 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0 0.05 0.1 0.15 0 1 2 3 4 5 6 7 8 9 Time Bending Energy 0 0.05 0.1 0.15 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 z Time (b) (c) (d) (a) F a = 190 F a = 180 x y Figure 3.2: Transition from 3D spinning to 2D flapping. (a) 3D spinning motion at = 180 versus 2D flapping motion at = 190. (b) when projected onto the(,)-plane, the tip of the spinningfilamenttracesacirculartrajectory,whilethetrajectoryofthetipoftheflappingfilament collapses to a fixed line. (c) and (d) shows that both the position of the distal tip and the total bending energy of the filament become constant in time for spinning, but remain oscillatory for flapping. Bysymmetryofthegeometryandmaterialpropertiesofthefilament,planarinitialperturbations lead to filament motions that remain confined to the perturbation plane. In particular, for planar initial perturbations, the filament in figure 3.1 undergoes 2D flapping motions as reported in [32] (seesupplementalmovieS2). However,non-planarperturbationsleadto3Dspinningmotions(see figure3.1andsupplementalmovieS1). Takentogether,theseresultsgobeyondthefindingsof[32] toshowthat,subjecttonon-planarperturbations,planarflappingcanalsobeunstable. We next examine the filament behavior for increasing values of . We observe a transition from3Dspinningto2Dwave-likeoscillationsasshowninfigure3.2(a). For = 180,thefilament settles into a buckled configuration with locked curvature and undergoes 3D spinning about the -direction, whereas for = 190, the filament converges in finite time to planar oscillations. The 52 0 50 100 150 100 200 300 Frequency F a Im[λ] 2D 100 0 20 40 60 80 F a 0 2̟ 1 (a) (b) 0 100 200 300 400 500 Spinning (locked curvature) Flapping F a 3D 2D Im[λ] 2̟ 1 Re[λ] 2̟ 1 Figure 3.3: Frequency of oscillations versus linear stability analysis. (a) real (dashed line) and imaginary (solid black line) parts of the dominant as a function of for the range of reported in [32]. A Hopf bifurcation at ≈ 37.7 marks the transition from a stable straight configurationtosustainedoscillations. Thenonlinearoscillationfrequenciesarecomputedfor non- planar perturbations (red line) and confined planar perturbations (blue line). Confined filaments oscillateatfrequenciesthatareinreasonableagreementwiththosepredictedbythelinearstability analysis,while3Dspinningmotionsoccurathigherfrequencies. Therelativefrequencydifference (1−flapping/spinning frequency) ranges from ∼ 15% to∼ 20%. (b) Same analysis extended to a larger range of – dashed box in lower left corner corresponds to the results in (a). We observe asecondtransitionat ≈ 188fromstable3Dspinningtostable2Dflappingmotions,evenwhen thefilamentissubjectedto non-planarperturbations. Atlargevaluesof ,theflappingfrequency deviatessignificantlyfromthelinearstabilityanalysis. (,) motionofthefilamenttipisshowninfigure3.2(b): theredlinecorrespondsto3Dspinning ( = 180) and the black line to 2D flapping ( = 190). The flapping motion occurs in an arbitrary plane, that depends on the non-planar initial perturbations. However, the frequency and 53 amplitude of these wave-like deformations are independent of initial conditions. For 3D spinning, the-coordinate of the filament tip and the filament bending energy go to constant values in finite time; see red lines in figure 3.2(c) and (d), respectively. For the 2D flapping motion, they oscillate with constant amplitudes (black lines). It is important to distinguish between the 2D flapping motion obtained here and the planar motions of [32]. The planar motions in [32] are unstable to non-planar perturbations and are obtained as a result of confining the filament in one plane. Here, the filament dynamically converges to planar deformations, even when subject to large non-planar initial perturbations (see supplemental movie S3). In other words, for this value of , the planar motionisrobusttonon-planarinitial perturbations. Weplotinfigure3.3(a)theoscillationfrequencyofthefilamentasafunctionof fortherange ofvaluesconsideredin[32]. Wecomparethesefrequenciestotherealandimaginaryvaluesofthe dominant eigenvalues obtained from the linear stability analysis about the straight configuration discussedin§3.1. Linearstabilityanalysisshowsthattheimaginarypartofthedominanteigenvalue (solidblackline)becomescomplexat ≈ 20whereastherealpart(dashedblackline),associated withtherateofgrowthoftheinitialperturbation,remainsnegativeuntil ≈ 37.7. Thatistosay, the straight filament configuration is monotonically stable for < 20 and stable with decaying oscillations for 20 < < 37.7. The filament undergoes a Hopf bifurcation at = 37.7, as the complex pair of eigenvalues crosses the imaginary axis, leading to unstable motions with growing oscillations. The linear oscillation frequencies are given by the imaginary part of the dominant eigenvalueIm[]/2. Thelinearstabilityresultsarequantitativelyconsistentwiththosein[32]. Thefrequencyofthenonlinearoscillationsiscomputedbylookingatcorrelationsofpeaksinthe timeevolutionoftheaveragepositionofthefilament. Nonlinearoscillationsareanalyzedfornon- planar and planar (confined) perturbations, depicted in red and blue, respectively. For > 37.7, 54 50 100 150 200 250 300 350 0 1 2 3 4 5 6 7 8 9 10 0 F a Bending Energy 3D 2D Flapping (stable) Spinning (stable) 2D Flapping (unstable) Figure 3.4: Bending energy versus axial force. The bending energy of spinning motion (red) increases logarithmically with following the curve b = ln( − 38.7) (solid black line). The envelope of the flapping bending energy, from minimum to maximum energy, is depicted for the unstable (confined) 2D flapping for < 188 (blue) and the stable 2D flapping for > 188 (redandblue). Theminimumoftheflappingbendingenergyofflappingstartstodecreasenearthe transition. non-planar perturbations give rise to 3D spinning motions at frequencies that are consistently largerthanthefrequenciesofthe2Ddeformationsassociatedwithplanarperturbations. Thelatter are comparable to the frequencies predicated by the linear stability analysis. Although the linear analysis correctly captures the transition from stable to unstable straight configuration, it only predictsplanardeformationsandcannotdistinguishbetween2Dand3Doscillations. Thefactthat thelinearanalysisonlycapturesplanardeformationsisduetothestructureofthelinearequations, which decouples torsion from bending curvature, as discussed in §3.1. Our nonlinear analysis indicates that, subject to non-planar perturbations, the Hopf bifurcation at = 37.7 leads to a buckled configuration with locked curvature that spins at a higher frequency (around 15% larger) thanthefrequencyofconfineddeformations(seefigure3.3a). 55 In figure 3.3(b), we extend this analysis to larger values of , with the lower left corner, highlighted by a narrow dashed line corresponding to figure 3.3(a). Our analysis indicates the presenceofanotherbifurcationat ≈ 188thatcausesthefilamenttotransitionfrom3Dspinning to2Dflappingmotions,evenwhensubjectedtonon-planarperturbations. Attheselargevaluesof ,theflappingfrequencydeviatessignificantlyfromthelinearstabilityanalysis;thelatterfailsto capture the second transition from 3D spinning to 2D flapping. We emphasize that these flapping motions are fundamentally distinct from the flapping motions reported in [32] for < 100. The latter are unstable to non-planar perturbations. For > 188, the filament dynamically converges toplanarflapping. Before we proceed, a few comments on the scaling of the dimensional force and instability threshold with the filament length are in order. The dimensional force is equal to / 2 . Thus, in dimensional form, for a filament of length = 20 m (see Table 2.1), the first transition from straight configuration to 3D spinning occurs at / 2 = 37.7/ 2 = 75.4 pN and the second transition from 3D spinning to 2D flapping occurs at / 2 = 188/ 2 = 376 pN. Shorter filamentsrequirelargervaluesofthedimensionalforcetotriggertheseinstabilities. We examine more closely the robustness of the 3D spinning and 2D flapping motions with respect to the amount of initial perturbation away from the straight equilibrium. Our numerical experiments show that, for 37.7 < < 188, any minute perturbation to the planar symmetry puts the filament into a 3D spinning state, whereas > 188 always leads to planar flapping even for dramatically non-planar initial conditions. These results indicate that the second transition is not sensitive to initial perturbations. They further suggest that for 37.7 < < 188, the basin of attractionofthe3Dspinningmodeisthesetofallnon-planarinitialconditions,whilefor > 188, 56 thebasinofattractionofthe2Dflappingisthefullspaceofinitialconditions,planarandnon-planar (SeesupplementalmoviesS3andS4). Toelucidatethephysicalmechanismsatplayintheseregimes,wenotethattheworkdonebythe axial force is balanced by both the elastic bending energy of the filament and the work dissipated via viscous drag due to the filament motion. In 3D spinning regime, the filament maintains a locked curvature, characterized by a constant bending energy. In figure 3.4, we plot the bending energy of the filament versus the force value for both 3D spinning (red) and planar flapping (blue) which are unstable for < 188. In its locked-curvature configuration, the bending energy grows logarithmically with . When undergoing 2D flapping deformations, the bending energy accesseslargerrangesofenergyvalues. Thisexplainswhythespinningfrequenciesarelargerthan the frequencies of confined 2D flapping: since a smaller amount of the work done by the axial force is absorbed by the bending energy of the locked filament than by the confined 2D flapping filament,moreworkisavailablefor3Dspinning. Inthe2Dflappingregime,theworkdonebythe axial force is too large; the filament cannot settle on a locked-curvature configuration that stores a sufficient amount of bending energy to allow the excess work to be dissipated via viscous drag by arigidrotationofthelockedfilament. Thefilamenthastodeform. Tosummarize,theresultsinthissectionshowthatthenonlineardynamicsofanelasticfilament inaviscousfluid,clampedatoneendandsubjecttoafollowercompressiveforceatitsfreeend,is farricherthanpreviouslyrecognized. Asthemagnitudeoftheaxialforce increases,fourdistinct regimesofdynamicalbehaviorsareobserved: (i)amonotonicallystableregime,(ii)aregimewith decaying oscillations, (iii) 3D spinning motions at a buckled configuration with locked curvature, and(iv)2Dflappingmotionswherethe filamentdynamicallyconvergetoplanaroscillations. 57 f a s a s αf a f a s 1 (a) (b) Figure 3.5: Axial force profiles. Schematic showing the distribution of the axial forces along the filament centerline: (a) an axial force density that varies linearly from the base of the filament to itstip,and(b)theaxialforceisconcentratedatadistanceof fromthefilamentbase. 3.3 Axialforceprofiles We examine the robustness of the four regimes identified in §3.2 to the distribution of axial forces alongthefilamentcenterline. Specifically,weconsidertwoforceprofiles. Thefirstprofileconsists ofalinear forcedensity,withtwoparameters and , f =− [1+(2− 1)]t. (3.10) Here, istheintensityoftheforcedensity,rangesfrom−1to 1and determinestheslopeof thelinearforcedistribution;negative correspondstoaforcedensitythatdecreasesas increases towardsthetipofthefilamentandpositive toaforcedensitythatincreasestowardsthefilamenttip. Thetotalcompressiveforceexertedonthefilamentisgivenby = ∫ 1 0 [1+(2− 1)] d= . The second profile consists of a concentrated force at a distance from the base point of the filament f =− (− )t, (3.11) 58 s a = 0.3 s a = 0.5 s a = 0.9 (a) (b) (c) (d) α = – 0.5 α = 0 α = 1 F a F a α s a Frequency Frequency Figure3.6: Galleryoffilamentoscillations. In(a)and(b),filamentssubjecttolineardistribution of axial forces exhibit buckling and 3D spinning motion only despite changes in parameter . But the threshold of instability decreases and the radius of the circle formed by the filament tips increasesforhigher. In(c)and(d),filamentsaresubjectedtoaconcentratedaxialforceatvarying location seesqualitativelydistinctregimesofmotion. For < 0.5,thedistaltipofthefilament remain taller than the largest circle formed by the spinning motion. When forces are concentrated closer to the tip, we observe the 3D spinning to 2D flapping transition of figure 3.2. Note that in both (b) and (d), changes in the threshold of initial instability of static straight equilibrium are well-matched by the linear stability analysis (solid black lines). = 125 and 300 respectively for snapshotsshownin(a)and(c). where(− ) is the Dirac delta function and ranges between 0 to 1; = 1 corresponds to the force at the tip of the filament. The total compressive force exerted on the filament is given by = ∫ 1 0 (− )d= . Aschematicdepictingthesetwoprofilesisshowninfigure3.5. In figure 3.6(a), we set = 125 and consider three values of. For=−0.5 the axial force decreasestowardsthefilamenttip, = 0correspondstoaconstantforcedensityalongthefilament, and for = 1 the axial force increases towards the filament tip. In all three cases, the filament bucklesintoalockedconfigurationandundergoes3Dspinningmotion. Wesystematicallyanalyze 59 the behavior of the filament as a function of the two parameters and . The results are shown in the parameter space(, ) in figure 3.6(b): black dots correspond to cases where the straight configuration is stable, the solid black line denotes the Hopf bifurcation marking the transition to growing oscillations based on the linear stability analysis of § 3.1. The colormap denotes the frequency of the nonlinear oscillations of the filament. The filament motion is represented next to each point(, ) by superimposing the filament configuration at various time steps within the oscillationcycle. Forcedensitiesthatdecreasetowardsthetipofthefilament( < 0)requirelarger forcevalues toreachthethresholdfortheHopfbifurcation. InallcasespasttheHopfbifurcation, thefilamentbucklesandundergoes3Dspinningatfrequencieslargerthanthefrequenciespredicted bythelinearstabilityanalysis. In figure 3.6(c), we consider a concentrated force = 300 and vary its location along the filament. At = 0.3, thefilament bucklesinto alocked curvatureand spinsabout the -direction. The filament exhibits large curvatures near the point of application of the axial force, with flatter curvaturenearthetip. For = 0.5,thefilamentalsobuckles,exhibitingsimilarbehaviorbutmore gradual variations in curvature. For = 0.9 the filament undergoes planar flapping motion, even when subject to non-planar perturbations. A gallery of the filament behavior as a function of the two parameters and is depicted in figure 3.6(d). The filament first transitions from a straight configurationto3Dspinningmotion. Asecondtransitionoccursfrom3Dspinningto2Dflapping motion as and increase further, as highlighted by the dashed black line in the upper right corneroftheparameterspace. Indimensionalform,theforcethreshold / 2 requiredtotrigger theseinstabilitiesincreasesasthefilamentlengthdecreases. For both profiles of axial forces, the first transition, marking the destabilization of the straight configuration under the distribution of axial forces , is captured by the linear stability analysis 60 -1 -0.5 0 0.5 1 50 100 150 200 0 1 2 3 4 0 0.5 1 1.5 2 0.2 0.4 0.6 0.8 0 0.5 1 1.5 2 0.2 0.4 0.6 0.8 150 250 350 450 0 1 2 3 4 α ×10 3 -1 -0.5 0 0.5 1 α Re[λ] Im[λ]/2π s a s a ×10 4 Re[λ] Im[λ]/2π ×10 3 F a F a (b) (a) ×10 2 Figure3.7: Linearstabilityanalysisofaxiallyloadedfilament. Weplotthefullresultsoflinear stability analysis for (a) the linear forcing profiles and (b) the concentrated forcing profiles. The real parts of indicate the stability of straight equilibrium, while the imaginary parts predict the frequencyofoscillations. Similartofigure3.3,theimaginarypartsbecomepositivebeforethereal parts, indicating the existence of decaying oscillations. The frequencies predicted here by linear stabilityareinroughagreementwiththoseshowninfigure3.6. of §3.1. As we vary the two parameters corresponding to each forcing profile, (, ) and ( , ), respectively,linear stability analysisindicates a transition from stableto unstable straight configuration through a Hopf bifurcation as shown in figure 3.7. However, linear stability analysis predicts growing planar oscillations at frequencies less than the frequencies of the 3D spinning motionsobtainedinthenonlinearsystem. Further,thislinearanalysisdoesnotcapturethesecond transition,fromnon-planarspinningtoplanarflapping,observedintheconcentratedforceprofile. Based on the gallery of filament oscillations reported in figure 3.6, the following observations are in order: (1) The first buckling instability, captured by the linear stability stability analysis, always leads to 3D spinning motions. (2) The second transition from 3D spinning to 2D flapping seems to occur when the force distribution is biased towards the tip of the filament; we explored thisobservationinmoredetailnext. (3)Awayfromthetransitionthresholds,thefilamentbehavior is insensitive to changes in the parameters of the axial force, even to changes in the form of the axial force profile. These observations could have important implications on the relevance of this instability-drivenmechanismtobiologicalciliaandflagellaasdiscussedin3.5. 61 To better understand the effect of the distribution of axial forces along the filament on this transition from 3D spinning to 2D flapping, we introduce a forcing profile that allows us to continuouslychangefromalinearly-distributedforceprofiletoaforceconcentratedatthefilament tip. Tothisend,weconsider f =− (+ 1) t, (3.12) where controls how much force is concentrated near the tip in a gradual way. At = 1, the force expression in (3.12) is identical to that in (3.10) for= 1. As →∞, this force converges to(3.11)for = 1. Again,wehave = ∫ 1 0 (+ 1) d= Figure 3.8 is a phase diagram that summarizes the behavior of the filament as a function of and . We findthat for abig range ofaxial forces, potentiallyexceeding the biologically-relevant range, when = 1, that is, for a linear force distribution, there is no transition from 3D spinning to 2D flapping motions. As exceeds 1, that is, for force distributions that are larger towards the filament tip, we observe a transition from spinning to flapping motions at biologically-relevant valuesoftheaxialforce (asdiscussedin§3.5). Forevenlargerforces,anothertransitionoccurs fromplanarflappingtotoroidalflappingmotions,andeventochaotic-likebehavior. Representative trajectoriesofthemeanpositionofthefilamentaresuperimposedontothephasespacetoillustrate the complex nature of these motions. These results indicate that higher concentrations of axial forcesnearthedistalportionofthemicrofilamentacceleratethedevelopmentofcomplexbehaviors. It is worth noting here that we calculated the Lyapunov exponents for large values of concentrated at the tip of the filament leading to chaotic-like behavior. Namely, we calculated the Lyapunovexponentsforarbitraryinitialconditionsbasedonfiniteseparation inthemeanposition of the filament (results not shown). Upon the removal of a rotational symmetry due to the initial 62 0 0.5 1 1.5 2 2.5 3 3.5 4 250 500 1000 2000 4000 Flapping Spinning Multi-Periodic + Chaotic F a p Figure3.8: Phasediagrambasedonpower-lawforcingprofiles. Weexpandfindingsoffigure3.6 byshowingthetransitionbehaviorsatmuchhigher usingthepower-lawforcingprofilesof(3.12). We see that only when the axial forces increase at a power larger than≈ 1.5 would 3D spinning transitions to 2D flapping. Moreover, at very large forces, we can also obtain much more exotic trajectories. The patterns shown above are formed by the average positions of the filament over its arclength(greyscalerepresentspassageoftime). Euler buckling, trajectories from two nearby initial conditions diverge initially, but saturate at a large distance, giving rise to Lyapunov exponents that are consistent with chaos. Similar chaotic trajectories are reported in recent models of active biological filaments [143]. It is not clear if thischaoticbehaviorisbiologically-relevantforabeatingflagellumbecausetheselarge maylie beyondthecapabilitiesofthedyneinmolecularmotors. However,itisnotcompletelysurprisingto obtainchaoticbehaviorin thisdissipativesystemduetothe resemblanceofthefilamentequations of motion to the Kuramoto-Sivashinsky equation that models diffusive instabilities in a laminar flamefront. To conclude this section, we note that the spinning to flapping transition is insensitive to the anisotropy in the fluid drag force. Our numerical experiments show that when changing , the filament continues to transition from 3D spinning to planar flapping then to 3D motion, albeit 63 requiringadifferentamountofaxialforce. Thisobservationwillallowustoignoredraganisotropy inourreduced-order,bead-springmodeldiscussednext. Time Bending Energy (b) (c) (d) (a) Spinning (Stable) F a Flapping (Unstable) ψ 2 ψ 1 100 150 200 0 0.1 0.2 0.3 0.4 0 50 Figure 3.9: Two-link chain. (a) A two-link chain subject to an axial follower force = 4, concentrated at its tip. (b) Starting from a straight configuration, non-planar perturbations lead to 3D spinning motions. (c) By symmetry, Planar perturbations lead to flapping motions, but these solutionsareunstabletothethree-dimensionalperturbations. (d)themaximumbendingenergyof planarflappingislargerthantheconstantbendingenergyofthespinningmotion. 3.4 Bead-springmodel The elastic filament, subject to an axial force at the free tip, buckles and undergoes 3D rigid spinning at a locked curvature before it transitions to 2D flapping motions at large values of . These nonlinear effects are not captured by the linear stability analysis. Here, we consider a bead- springmodelofthefilamentwiththeaimofreproducingthesenonlineareffectswiththeminimum numberofdegreesoffreedom. We model the elastic filament as a chain of beads located at r , where = 1,...,, and connected successively via inextensible rods of length / from the origin r =(0, 0, 0) T ; see figure 3.9 for a two-link example. To emulate the filament elasticity we attach to each bead ≠ atorsionalspring,thatproducesaspringmoment , =− where isthespringconstantand := cos −1 (t −1 ·t ) is the the angle between successive links, with t :=(r −r −1 )/∥r −r −1 ∥. 64 We enforce the clamped boundary condition at r weakly by including a spring moment , = − cos −1 (e 3 ·t 1 ), proportional to the angle between the equilibrium axis and the first link. The chainissubjecttoanaxialforceactingatitstip,onbead,oftheformF =− t . Eachbeadis subjecttoanisotropicdragforceF ℎ, =−v . 1 2 3 4 5 6 -0.2 0 0.2 0.4 0.6 0.8 1 2 3 4 5 6 -0.2 0 0.2 0.4 0.6 0.8 3D 2D F a 3D 2D F a F a Im[λ] 2̟ 1 Im[λ] 2̟ 1 (a) (b) Re[λ] 2̟ 1 Re[λ] 2̟ 1 Figure3.10: Linearstabilityoftwo-linkandthree-linkchains. Nonlinearfrequenciesofthe(a) two- and (b) three-link chains compared to the dominant eigenvalue from linear stability analysis. Similar to figure 3.3, the frequency of 3D spinning is consistently larger than the frequency of planarflapping. Thelatterisunstabletonon-planarperturbations. We rewrite these quantities in non-dimensional form using the total chain length as the characteristic length scale and 3 / as the characteristic time scale. To this end, the non- dimensional force is given by / as opposed to the non-dimensional force in the filament 65 2 4 6 8 10 12 14 16 18 20 F a N=2 Flapping N=5 Spinning Bending Energy 0 1 2 3 4 5 6 7 8 9 Windmilling N=4 N=6 N=3 Number of Links F a 2 3 4 5 6 20 15 7 5 3 (a) (b) Figure 3.11: Transition to non-spinning solutions of the bead-spring models. We perform a continuation type analysis for different number of links, where we let the chain settle onto a stable solution at a given and successively increase to obtain new stable solutions. (a) The maximum and minimum envelope of the resultant bending energy as a function of . Spinning motionsbecomeunstableatlarge for≥ 3. (b)Chainbehaviorasafunctionofthenumberof linksand . flappingsolutionsbecomestableforafixedrangeof onlyat = 6. model 2 /. Thus, the range of relevant values of the dimensionless force in the chain model may differ from that of the filament model. In the following, we consider all quantities are non-dimensional. We write the position vector to bead as r = r −1 +t /, where the vector t is expressed in spherical coordinates and relative to bead, namely, t ≡(cos cos , cos sin , sin ) T . Wecanthenexplicitlyexpressthevariationsr andusetheprincipleofvirtualwork, F ·r + ∑︁ =1 F ℎ, ·r + −1 ∑︁ =0 , = 0, (3.13) to obtain 2 governing equations for an-link chain in terms of and . Note that after proper nondimensionalization, the planar two-link equations of [32] can be recovered by setting = 0, 66 = 0, 1. For = 1,inextensibilitydictatesthatthetipfollowerforcehasnoeffect,thusthestraight equilibrium configuration is globally asymptotically stable. The system admits both spinning and flapping trajectories for all ≥ 2. However, since torsion is related to 3 r/ 3 , to account for the effectsoftorsion,weneed ≥ 3. Infigure3.9,weshowthebehaviorofatwo-linkmodelsubjecttoanaxialfollowerforce = 4, concentrated at its tip. Starting from a straight configuration, non-planar perturbations bring the chain into 3D spinning motion. Planar flapping motions can only be obtained for planar initial perturbations,andareunstabletonon-planarperturbations,consistentwiththefullfilamentmodel. The bending energy of spinning is constant over time while the bending energy of the flapping motion oscillates as a function of time, accessing larger energy values than the spinning motion, similartothefullfilamentmodel. Welinearizetheequationsofmotionforatwo-linkchainandathree-linkchainandexaminethe dominanteigenvaluesofthelinearsystem. Thestatevariables from decoupleatthelinearlevel, similartothedecouplingofcurvatureandtorsioninthefullfilamentmodel. Figure3.10showsthe nonlinear frequencies of the two- and three-link chains, subject to 3D and 2D perturbations, and compares them to the results of the linear stability analysis. As in the full filament model, a Hopf bifurcation leads to sustained oscillations of the chains, and as in figure 3.3, the frequencies of 3D spinning are consistently larger than the frequencies of planar flapping of the confined chain. The lattersolutionsareunstabletonon-planarperturbations. Inordertoexaminethetransitionfrom3Dspinningto2Dflappingmotionsinthechainmodel, weperformacontinuationtypeanalysisforfivedifferentlinknumbers = 2to 6. Westopat = 6 because it corresponds to the chain with the smallest number of links that unambiguously exhibit the transition from 3D spinning to 2D flapping. For each , we start with = 0 and increase 67 its value systematically as follows: we let the chain settle into a stable solution at a given value of , we then increase the value of to obtain a new stable solution, and repeat successively. In figure3.11(a)weshowthemaximumandminimumenvelopeofthebendingenergyasafunctionof . Sincethe3Dspinningmotionscorrespondtorigidrotationsatalockedcurvatureofthebuckled filament, the bending energy is constant, and its maximum and minimum are equal. Observe that spinning become unstable at large values of only for ≥ 3. In figure 3.11(b), we illustrate the chain behavior by superimposing several snapshots of the same chain during one period of its oscillations. We show the filament behavior as a function of the number of links and the tip force . For 3≤ < 6 we observe a transition from spinning to ‘windmilling’ motions, where the chain or part of the chain rotates (windmills) roughly about an axis in the(,) plane. For = 6, planar flapping solutions become stable for a range of . This is also true for > 6 (results not shown for brevity). In the limit of large , the discrete chain model exhibits similar results to the discrete filament model of Section A, albeit without drag anisotropy. That is to say, instability-driven oscillations occur under the assumption of local isotropic drag, indicating that they are insensitive to the hydrodynamic force model. This is in contrast to systems that require hydrodynamicinteractionsalongthechaintotriggersustainedoscillations[144]. The transition away from the 3D spinning motion is consistent with our physical intuition in §3.2. When spinning, the work done by the tip force is balanced by the constant elastic energy stored in the locked configuration of the chain and the work dissipated through viscous drag due to the “rigid” rotation about the -axis. As the tip force increases, to maintain this balance, the three-link chain opts to minimize time-dependent deformations by rigidly rotating the whole chain about an axis in the(,) plane. This rotation is permitted because the clamped boundary condition in the chain model is imposed weakly, via an elastic spring at. Since these rotations 68 activate the elastic spring at the base, the corresponding elastic energy is non-zero. Also, these rotations lead to larger linear velocities and thus larger drag forces (that is, larger dissipation) than rotationsaboutthe-axis. Forfour-andfive-linkchains,theupperlinkofthechainsuccumbsfirst undertheaxialforceandbeginstorotateinawindmill-likemotion,asillustratedinfigure3.11(b). However, as the number of links increase, these hybrid solutions – consisting of rigid spinning aboutthe-axiswithwindmill-likerotationsthatarespatially-localizedatthedistallinks–become unstable. At = 6,thechainclearlyexhibitswave-likedeformations. 3.5 Discussion We considered an elastic microfilament, clamped at one end in a viscous fluid, subject to a distribution of axial forces. We identified two major transitions in the filament behavior. First, the filamenttransitionsfromastraightconfigurationtoabuckledconfiguration,withlockedcurvature, and undergoes 3D spinning. The second transition occurs at larger values of the axial force and marks a destabilization of the 3D spinning motions, giving rise to 2D flapping oscillations. These periodic motions, both 3D spinning and 2D flapping, are robust to large perturbations in the filament configuration. They are also robust, away from the bifurcation values, to perturbations in the parameters of the axial forces, and even to changes in the force profile altogether (as long as theforcedensityincreasessufficientlynonlinearlytowardsthetipofthefilament). Thesefindings, althoughinthecontextofasimplifiedmodel,supporttheideathatanopen-loop,instability-driven mechanism could explain not only the origin of sustained oscillations but also the wide variety of periodicbeatingpatternsobservedinciliaandflagella. 69 (c) (d) (a) (b) F a s Time 0.13 0.138 0.146 0 1 -4 -2 0 2 4 0.356 0.36 0.364 Time 0 0.5 1 -15 -10 -5 0 5 0.32 0.33 0.34 0.35 0.36 Time 0 0.5 1 -15 -10 -5 0 5 θ Figure 3.12: Asymmetric beating patterns. (a) A naturally-curved filament with a reference curvature ofκ = 2d 1 is subjected to an axial force at its tip. We define the rotation angle in the(,)-plane associated with the vertical projection of the tip point. (b) At = 100, an asymmetric 3D spinning pattern develops: monotonically decreases and the vertical tip position variesperiodically. (c)At = 200,thefilamentengagesinasymmetricflappingwhileitslowly precesses about the vertical axis: drifts on average. (d) At = 250, the asymmetric flapping is nearly planar but switches between flapping on the left and right sides of that plane: is periodic. Thefilamentsnapshotsaretakenfromtheregionshighlightedingrey. Cilia and flagella, despite their highly conserved internal axoneme structure across cell types, exhibitdistinctbeatingpatternsdependingonthecelltype,environmentandfunction. Forexample, sperm cells propel themselves by planar wave-like deformations of their flagellum [145]. The biflagellate algal cell Chlamydomonas swim using synchronized flapping oscillations of their flagella, referred to as ‘breast stroke,’ and can transition to unsynchronized wave-like oscillations, referred to as ‘free style,’ that allows them to turn and reorient themselves [146, 147]. In densely- packedciliatedsurface,motileciliaoftenbeatin3Dpatterns,withanasymmetryintheireffective andrecoverystroke. Agalleryofvariousplanarandnon-planarciliabeatingpatternscanbefound in[6,7]. The model presented in this paper does not accurately represent the intricate internal structure of cilia. It rather assumes that its net effect on the cilium centerline can be approximated by a distribution of axial forces of constant magnitude. An accurate representation of the effect of the 70 axoneme structure on the centerline motion would give rise to a distribution of axial and normal forces that are coupled to the centerline geometry. These limitations notwithstanding, we use the current model to roughly gauge whether the instability-driven oscillations and the transitions reported in this study are relevant to the mechanics of cilia and flagella. To this end, we provide a rough estimate of the total axial force that is biologically attainable. We estimate the total amount offorcesthatthedyneinmotorsprovideinsideatypicalcilium. Theaxonemeconsiststypicallyof ‘9+2’ microtubule doublets, with nine outer doublets and a central pair (see figure 1.1). Dynein molecules form an array of cross-bridges between pairs of outer microtubule doublets. The force exerted on microtubules by a single molecular motors is known to be of the order of a few pico- newtons. So assuming 8 dynein arms per 100 nm and 3-8 pN of force per dynein arm [148], we would get a per-doublet force density estimate of 200-600 pN·m −1 .This gives us about 2-5 nN·m −1 per cilium assuming all 9 doublets are activated. In reality, due to the internal structural dampingandotherreactionforces(mostdoubletsexperiencedyneinforcesinbothcompressiveand extensiledirections),thetotalcompressiveforcedensityislikelytobemuchsmaller. Inourmodel, a force density on the order of 20 pN·m −1 can lead to sustained oscillations and is sufficient to eveninducethetransitionfrom3Dspinningtoplanarflappingregime. Thiscalculationisbasedon the cilia length and bending rigidity estimate of Table 2.1, where we have chosen a typical length of 20m as cilia and flagella can vary in length between 10-50 m, and a high stiffness value of 800pN·m 2 comparingtovaluesof25-800pN·m 2 reportedinliterature[55,110,113,134,149]. These numbers suggest that the axial forces due to dynein activity in cilia and flagella should be large enough to support these instability-driven oscillations, as noted in more elaborate models of [27,31,56]. 71 Another important consideration for exploring the relevance of these instability-driven oscil- lations to biological cilia is that cilia often exhibit asymmetric beating patterns. To account for asymmetric oscillations in our filament model, we incorporate a non-zero reference curvature of the filament. Experimental evidence from thealgal cells Chlamydomonas supports the notion that their flagella are naturally curved [58]. In our model, we set M = B(κ −κ ), where κ is a non-zero reference curvature. We consider κ = 2d 1 , as depicted schematically in the top of Figure3.12(a). Figure3.12showsthefilamentdynamicsunderthreevaluesoftheaxialforce at the free tip of the filament. We identify three distinct behaviors of the naturally curved filament: (i) “asymmetric” spinning (figure 3.12b), reminiscent to the spinning behavior of the naturally straight filament, but here the spinning axis is not normal to the base wall to which the filament is clamped; (ii) “asymmetric” precession while flapping (figure 3.12c), where the flapping motion is combined with a slow precession about the vertical axis, and flapping frequency about five times faster than the precession frequency; and (iii) “asymmetric” flapping (figure 3.12d) that is nearly, butnotcompletelyplanar;thefilamentswitchesbetweentheleftandrightsidesofaflappingplane in an alternating fashion, as evident from the shadows in the inset of Figure 3.12d. These regimes can be cleanly delineated from the time evolution of a cumulative rotation angle, defined in the (,)-plane as the angle to the vertical projection of the filament tip, see the bottom schematic in figure3.12(a). Inthetoprowoffigures3.12(b-d),weplotthechangein afterthefilamentreaches itslimitcyclebehavior. At = 100(figure3.12b), changemonotonicintime,indicatingthatthe filament is simply spinning. At = 200 (figure 3.12c), the fast periodic fluctuation of reflects the flapping motion, and the slow drift of its average value captures the filament precession. At = 250 (figure 3.12c), varies in(−,] indicating that the motion is purely flapping. Due to the asymmetry incorporated in the system, the vertical position for the filament tip changes with 72 time in all three behaviors; see bottom row of figures 3.12(b-d). In particular, over one period of oscillation(somehighlightedingreyinfigure3.12),thefilamenttipperiodicallygetsclosertothen furtherawayfromthebasewall,inawayreminiscenttotheeffectiveandrecoverystrokesinciliary beatingpatterns. Theworkinthisstudyisafirststeptowardsunderstandingtheinterplaybetweenthemolecular motor activity, the elastic properties of biological filaments, and the fluid mechanical forces in generatingthewidevarietyofbeatingpatternsobservedacrossflagellatedcellsandinthetransitions between multiple beating modes in the same cell as the level of molecular motor activity changes. Thisworkfocusedonthecaseofasinglefilamentinaviscousfluid,usingtheresistiveforcetheory. One of the main advantages of this simple model is that it avoids the requirement of numerous materialparameters,unknownexperimentally,suchasthoseneededinafull-fledgedfiniteelement simulation involving more details of the internal structure of flagella [42, 56]. This simple model also allows us to focus on the fundamental mechanisms underlying flagellar oscillations and the varioustransitionsthatcouldpotentiallyberesponsibleforthevarietyofbeatingpatternsobserved innature. However,afewremarksonthelimitationsofthemodelareinorder. Ourmodeldoesnot fullyreflecttheinternalmechanicsofcilia. Inlaterchapters,wearederivingtheequationsofmotion and effective forces acting on the centerline from more elaborate, coupled, multi-filament models of the axoneme structure. The present model does not account for the long range hydrodynamic effects of viscous fluids. Future extensions of this work will account for the full hydrodynamics, usingacombinationofslenderbodytheoryandregularizedStokesletformulation,inthepresence of the cell body or the cell wall [55, 150]. Future work will also focus on the effects of anisotropy inthebendingrigiditymatrixB,whichencodesinformationabouttheinternalaxonemestructure. Experimental evidence indicates that the presence or absence of certain structures within the 73 flagellum could lead to 3D versus planar beating patterns in the Chlamydomonas flagella [11]. Thissuggeststhatanisotropyofbendingrigiditycouldbeanothercrucialeffectresponsibleforthe varied motion of cilia and flagella. It is also important to account for the effects of the intrinsic oscillations of the dynein motor forces due to the time scales of their binding and unbinding to microtubules and due to thermal fluctuations [148]. Another interesting avenue for exploration is to treat this problem as a reformulation of the classic anti-phase vs in-phase synchronization of two planar oscillating filaments [151]. Then our discovery can be recast as two filaments each representing oscillation in xz- or yz- plane with nonlinear geometric coupling would lead to a transition from anti-phase synchronization, a spinning motion, to in-phase synchronization, a motion with a fixed plane! Lastly, this instability-driven mechanism could have implications on thecoordinationsofmultipleactivefilamentswithnon-localhydrodynamicinteractions[55]. The presence of neighboring filaments and other solid boundaries may affect the transition thresholds andbiasplanarversusnon-planarbeatingpatterns. 74 Chapter4 Beatdiversityfrommolecularmotordynamics In this chapter, we examine what kind of ciliary waveform diversity can be achieved using the simplified two-state, two-filament model of cilium as explained in Section 2.2. In particular, we willdemonstratehow: (i)theboundaryconditionofthebasalendofthecilium,(ii)theforward-aft molecular motor activity asymmetry, and (ii) the type of geometric feedback mechanism change thewavetravelingdirectiononthecilium. 4.1 Effectofchangingboundary conditions Extendingtheanalysisofclampedfilamentsin[115],weshowthathingedfilamentsalwaysexhibit wave propagation from base-to-tip, in contrast with the fact that clamped filaments exhibit both base-to-tip and tip-to-base oscillations. In hinged filaments, we observe that molecular activities form sharp propagation fronts while the filament oscillates smoothly. For clamped filaments, in additiontotheanalyzingdirectionofwavepropagationovertheparameterspaceofmotoractivity and Sperm number as done in [115], we analyze in detail the change in amplitude and frequency ofoscillationsoverthesameparameterspace. Weobservesharptransitionsinbothamplitudeand frequencyacrossthelinesmarkingthechangeinwavedirection. Weconcludebycomparingthese trendstoexperimentalobservationsofciliabehavior. 75 We first consider the model cilium with clamped boundary conditions. We vary the sperm number S and activity and examine the resulting beating patterns. Fig. 4.1 shows that, for S = 10, as increases, the bending waves change from tip-to-base propagation to base-to-tip propagation,asreportedin[115]. Allresultsareinsensitivetoinitialconditions(seeAppendixC). For = 6800, as S increases, the bending wave switches from base-to-tip propagation back to tip-to-base, and eventually remain stationary for large values of S. Increasing S is equivalent to increasing the cilium length, which, when holding all other parameters constant, leads to cases whereactivityisinsufficienttotriggeroscillations. Wenextexaminetherelationshipbetweenthefilamentgeometryandthedyneinmotordynamics. According to (2.23), the dynein detachment rates ± are governed by∓Δ . Fig. 4.2 shows that μ a = 6320 Sp = 10 μ a = 6800 Sp = 5 μ a = 6800 Sp = 10 μ a = 6800 Sp = 15 μ a = 10160 Sp = 10 μ a = 14000 Sp = 10 Figure 4.1: Gallery of cilium waveforms due to uniform motor activity. Top: as motor activity increases, the wave form changes its propagating direction from tip-to-base, to base-to-tip. Bottom: as sperm number Sp increases, the wave form changes back and eventually become stationary for large Sp. In all simulations, we used a mesh sizeΔ= 0.01 and an integration time stepofΔ= 3.125× 10 −4 . 76 μ a = 14000 Sp = 10 0 0.02 0.04 0.06 0.08 0.1 -1.5 -1 -0.5 0 0.5 1 1.5 -3 -2 -1 0 1 2 3 s 0 0.02 0.04 0.06 0.08 0.1 0 0.2 0.6 0.4 0.8 1 n + n − Δ + − Δ n ± s 0 0.2 0.6 0.4 0.8 1 0 0.02 0.04 0.06 0.08 0.1 -1.5 -1 -0.5 0 0.5 1 1.5 -3 -2 -1 0 1 2 3 0 0.02 0.04 0.06 0.08 0.1 0 0.2 0.6 0.4 0.8 1 n + n − Δ + − Δ n ± 0 0.2 0.6 0.4 0.8 1 μ a = 6320 Sp = 10 Figure 4.2: Motor switch-inhibition for clamped boundary condition. Snapshots of sliding magnitudeΔ and dynein motor activities on+ and− filaments, as well as the values of Δ, + , − along arclength at different times (opacity indicates passage of time). Top row shows example with tip-to-base waves at a small value of , while bottom shows that with base-to-tip waves after the second transition. Both examples have a clamped boundary condition at the base. In all simulations,weusedameshsizeΔ= 0.01andanintegrationtimestepofΔ= 3.125× 10 −4 . changes in the slidingΔ correspond to the activation and inhibition of molecular motors along the +and−filaments. Specifically,thetip-to-basebendingwavecorrespondstoatip-to-basetraveling wave in the motor activation dynamics ± , and when the bending wave travels from base to tip, so does the motor activations. This is in accord with the experimental findings of [44], where a ‘switch-inhibition’mechanismwasproposed. For hinged boundary conditions (equation (2.15)), the model cilium produces a bending wave travelingfrombasetotip,consistentwiththelinearizedanalysisof[23],forallSpermnumbersSp andactivity tested. Notethatforsmallvaluesof ,themotorbindingandunbindingdynamics switch between the+ and− filaments with no prominent propagation of activity from base to tip; see top row of Fig. 4.3. This implies that the base-to-tip bending waves are mostly due to elasto- hydrodynamics. For larger values of (bottom row of Fig. 4.3), the motor binding dynamics 77 μ a = 3500 Sp = 5 0 0.02 0.04 0.06 0.08 0.1 -1.5 -1 -0.5 0 0.5 1 1.5 -3 -2 -1 0 1 2 3 s 0 0.02 0.04 0.06 0.08 0.1 0 0.2 0.6 0.4 0.8 1 n + n − Δ + − Δ n ± s 0 0.2 0.6 0.4 0.8 1 0 0.02 0.04 0.06 0.08 0.1 -1.5 -1 -0.5 0 0.5 1 1.5 -3 -2 -1 0 1 2 3 0 0.02 0.04 0.06 0.08 0.1 0 0.2 0.6 0.4 0.8 1 n + n − Δ + − Δ n ± 0 0.2 0.6 0.4 0.8 1 μ a = 1000 Sp = 5 Figure 4.3: Motor switch-inhibition for hinged boundary condition. Snapshots of sliding magnitude Δ and motor activities on+ and− filaments, as well as the values of Δ, + , − as functionsofarclengthatdifferenttimes(opacityindicatespassageoftime). Toprowshowsresults forasmallervalueofmotoractivity . Bothexampleshavebendingwavestravelingfrombaseto tip,butatsmaller ,thereislittlepropagationin ± ,thusthefilamentbendingwavesareattributed toelastohydrodynamics. Inallsimulations,weusedameshsizeΔ= 0.01andanintegrationtime stepofΔ= 3.125× 10 −4 . form wave patterns that travel from the base to its tip. Fig. 4.3 also shows the formation of sharp propagation fronts (or shocks) in ± . We numerically verified the existence of these solutions for distincttimesteps(seeAppendixC).Thesesharpfrontsimplythatthemolecularmotordynamics followanon-linearwaveequation. We analyze the behavior of the clamped filament as a function of the parameter space Sp and . Analysisofhingedfilamentisomittedherebecausethedirectionofpropagationofthebending wavesdoesnotchangeaswevarytheseparameters. Specifically,inFig.4.4,wereportthebeating amplitude and frequency of oscillations for Sp = 5− 15, = 2000− 14000. We also report the evolution of the dominant eigenvalue based on linearized analysis (2.16) as a function of for Sp= 10. The linear stability results show an initial Hopf bifurcation leading to spontaneous 78 5 10 15 5 10 15 Sp Sp 0.2 0.4 0.6 0.8 1 1.2 1.4 0 0.05 0.1 0.15 0.2 0.25 0.3 ×10 4 0.2 0.4 0.6 0.8 1 1.2 1.4 ×10 4 Amplitude Frequency 0 0.1 0.2 0.3 0.4 0.5 0.6 14 0 2 4 6 8 10 12 μ a μ a ×10 -3 0 1 2 3 Sp = 10 Re Im Dominant Eigenvalue Figure4.4: Analysisofslidingfeedbackcontrolofdyneinkinetics. Amplitudeandfrequencyof oscillations. Amplitudeisdefinedinthetransverse( )direction. Fundamentalfrequencyobtained through autocorrelation of the time histories of the nonlinear waveform. Color correspond to amplitude and frequency values as indicated by the colorbars to the right. Black squares indicate that multiple traveling or standing waves compete and a steady state behavior was not observed. Dominanteigenvaluesfromlinearstabilityanalysisareshownwithrespecttoactivity atSp= 10. The light gray region highlight the regime post the first bifurcation, where the filament undergoes spontaneous oscillation with tip-to-base bending waves. Dark gray indicates the regime post the second bifurcation where the filament the bending waves travel from base-to-tip. In both regions, frequencydecreasesandamplitudeincreaseswithrespecttoincreaseinactivity. NotethatSperm numbershownisnormalizedwithrespecttodragcoefficient ⊥ (insteadoffluidviscosityin[115]). Inallsimulations,weusedameshsizeΔ= 0.01andanintegrationtimestepofΔ= 3.125×10 −4 . oscillations from tip-to-base, and a second transition leading to base-to-tip traveling waves, as reportedin[115]. Acloseexaminationoftheamplitudeandfrequencyofoscillationsacrossthese transitions (delineated by solid lines) in the left two panels of Fig. 4.4, show the following: (i) there is a jump or discontinuity in amplitude and frequency at the second transition when bending waves change the direction of propagation, and (ii) in both regimes, the amplitude of oscillations increase with increasing activity and the frequency decreases. This is in contrast to experimental observations where an increase in ATP concentration (for which is a proxy) correspond to saturatedgrowthinbothbeatingamplitudesandfrequencies[42]. 79 4.2 Effectofforward-aftmotor activity asymmetry Inprevioussection,dyneinmotordutyratio(dimensionlessattachmentrate)iskeptasaconstant, independentofthearclengthvariable. Takingasimilarapproachweusedinthecaseoftangential follower force profiles [3], it is again instructive from a reduced-order modeling perspective to see how the resultant cilium beats would change if we varying linearly as a function of arclength while keeping the total integrated activity constant (/∈ [−1, 1], ∫ 1 0 = const.); also see Section3.3. Recall that under sliding control, as the total integrated activity increases, the clamped-free ciliumfirstundergoaHopfbifurcationintooscillationswithbendingwavestravelfromtip-to-base, then go through a second transition where the flagella waves reverses its direction to go from its proximalbasetoitsdistaltip(Fig.4.4). Now,introducingthefollowingdutyratioprofiles ()= 1− (2− 1) , (4.1) we show that the there continue to be two potential bifurcating transitions based on the linear stabilityanalysisawayfromthestraightequilibrium,seeFig.4.5. Specifically, when activity is concentrated near the base region ( / =−1, left column of Fig. 4.5), the first Hopf bifurcation leads directly to base-to-tip bending waves, as indicated by the negative phase gradient ∠(,) < 0 (see Eq. (2.37)). Moreover, these linearly unstable oscillatory modes also loses their dominance at a larger value of , indicating a potential second transitionwherewavedirectioncouldpotentiallyreverse. 80 Figure 4.5: Linear stability analysis for sliding control mechanism under asymmetric dynein kinetics profiles. Left. When motor activity profile is base-heavy, dominant solutions to the linearized problem indicates two transitions from no sustained oscillations, to base-to-tip waves, then to a non-linear instability regime. Center. Comparing with the uniform case shown in Fig. 4.4 above, dominant solution first go through a Hopf bifurcation to tip-to-base waves, then a second transition zone. Right. As activity shifts to be more tip-heavy, Hopf bifurcation wave directions remain tip-to-base while both transition thresholds decreases. Colorbar represents the meanphasegradientofthedominantmode ∠,whichindicatesflagellawavetravelingdirection; seeEq.(2.37). As activity shift to be more concentrated near the distal end (/ = 1, right column of Fig. 4.5), the first Hopf bifurcation regains its tip-to-base bending waves ( ∠(,) > 0), with both first and second transition thresholds (black lines) shifts downwards to lower values of over all activity . This makes intuitive sense as more activity near the free end should make it easier toexcitetheunstablemodesofthecilium. Interesting to note, however, if we change the geometric feedback mechanism from sliding- to curvature-basedruleswheremotorsdetachifandonlythelocalciliumcurvaturereachesathreshold, a similar parametric study reveals that only base-to-tip waves ( ∠(,) < 0) are achievable and no the linear oscillatory modes remain dominant for all tested ranges of without a possible secondtransitionzone;seeFig.4.6. 81 Figure 4.6: Linear stability analysis for curvature control mechanism under asymmetric dynein kinetics profiles. Under curvature-based geometric feedback mechanism, tip-heavy ac- tivity cannot force the dominant Hopf bifurcation modes away from base-to-tip waves. Colorbar represents the mean phase gradient of the dominant mode ∠, which indicates flagella wave travelingdirection;seeEq.(2.37). Inordertoverifyourlinearanalysis,weperformacontinuationstylenonlinearsimulationakin to what was done in Fig. 3.3 but using our two-state, two-filament cilium method as described in SectionB,withtheadditionofincorporatingavariablemotordutyratio()accordingtoEq.(4.1). In Fig. 4.7, we see that indeed the second transition zones revealed in Fig. 4.5 materialized in flagella wave reversals: when activity profile is base-heavy, simulated wave forms go from no sustainedoscillations,tobase-to-tipwaves,thentotip-to-basewaves(leftcolumn);otherwise,cilia oscillations conform to the finding of the previous section where we first have tip-to-base bending waves, then move on to base-to-tip waves as overall activity parameter increases (center and right column). Moreover, in all cases, there exists a bi-stable zone where both wave traveling directionsareachievabledependingontheinitialconditions. However, when comparing across columns of Fig. 4.7, the most striking difference is that the bi-stable zone increases in size as activity shift towards base-heavy profiles and come to nearly cover the entire tested parameter space in the limiting case of/ =−1. This is an intriguing 82 Figure4.7: Wavetravelingdirectionsunderasymmetricdyneinkineticsprofiles. Weshowthe flagella wave traveling direction with the thrust generated for base-heavy, uniform and tip-heavy activity profiles. It is evident that bistable region is maximized when activity is concentrated near thebasalregion,andtheregionisminimizedatthetip-heavyprofile. Inaddition,weseethatwhen activity is stronger near the base, the net thrust force is more prograde then retrograde. Colorbar represents the net tangential force⟨f ·t⟩, which indicates flagella wave traveling direction; see Eq.(B.9). finding as it suggests that base-heavy acitivity profiles can in theory maintain the most amount of freedom in ‘choosing’ which wave direction is desired. In addition, when activity profile are base-heavy (/< 0), the resulting simulated waveform can also be more diverse: (i) for small Sperm numbers, we recover wide beatings that are qualitatively similar to stereotyped cilia beats, (ii) for larger Sperm numbers, we could get both base-to-tip and tip-to-base flagella waves; see Fig.4.8A. 83 Figure 4.8: Waveforms observed in experiments and accessible under asymmetric dynein kinetics profiles. A. Simulated waveforms starting from different initial conditions and various Sperm number and three values /. = 8000 in all cases. B. Experimentally observed waveforms of Leishmania Mexicana extracted from the data of [118]. The resulting waveforms havebeenrecenteredandalignedatthebasetoremovetherelativemotionofthecellbodybetween each frames. Colorbar represents the net tangential force⟨f ·t⟩, which indicates flagella wave travelingdirection;seeEq.(B.9). Inaddition,whencomparingsimulatedwaveformswithexperimentalobservations(Fig.4.8B) extracted from Leishmania Mexicana [118], a parasite that can control its wave directions, it is evident that regardless of their wave traveling directions, the tip of their flagellum do not move nearly as much as their basal counterparts; a characteristic that is shared with our base-heavy profile results at intermediate Sperm numbers. This hints that even for tip-to-base waveforms observedin-vivo,theflagellumcouldstillbeoperatinginaregimewheremolecularmotoractivity isbase-heavy. Infact,apossiblemechanismforflagellawavereversalcouldoccurinsequence: (i)startfrom, (ii) tip gets excited (iii) continue beating with. To confirm or refute this hypothesis might require liveimagingofthemolecularmotoractivityduringsuchareversalevent. 84 Figure 4.9: Excitation of flagella wave reversal via bistability. Simulated wave reversal event using time-varying motor activities. Increasing shades of gray indicates passage of time, and the reversaleventcyclesthroughpartA→B→D→C.A.Initially,motoractivityconcentratednearthe proximal base of the cilium (/ < 0) drive the bending waves to go from its base towards the distal tip. During this stage, bending wave amplitudes are larger near the basal end. B. As activity moves to become more uniform, base-to-tip waves continue to dominate and propagates throughout the entire flagellum. D. When molecular motor near the distal ends are selectively activated (/ > 0), even if only for a brief moment, tip-to-base waves quickly dominates the entire flagellum. C. Bistability regims shown in the left column of Fig. 4.7 ensures that bending waves will continue to propagate from tip-to-base even if motor activities switch back to a more uniform or even slightly base-heavy profiles. Now the bending wave amplitudes again are larger nearthebasalend,comformingthetheexperimentalobservationsextractedinFig.4.8B.However, if either the total activity level drops below the bi-stable threshold or activity in the distal half oftheflagellumstopscompletely,wewillrestartthecycleandgobacktopanelAwithbase-to-tip waves. Here = 8000, Sp= 8 forallpanels. 4.3 Discussion We examined a model of cilia oscillations with sliding and curvature feedback control [114, 115], and examined differences in filament dynamics with respect to different boundary conditions and activity asymmetries along the cilium length. As the activity level increases, the model exhibits a 85 Hopf bifurcation that leads to sustained oscillations. The direction of deformation waves depends ontheboundaryconditions. Forhingedboundaryconditions,thewavespropagatefrombase-to-tip only. For clamped boundary conditions, the waves propagate from tip-to-base following this Hopf bifurcation, and a second bifurcation exists at higher activity levels leading to reversal in wave propagationfromtip-to-basetobase-to-tip. As cilia emerge from a basal body whose details differ among cilia types [152], the exact mechanical conditions at the base of cilia and flagella are very complex. As such, the physical boundaryconditionsofactualciliaareprobablyinbetweenclampedandhingedandlikelytovary betweenorganismsandduringthedifferentstagesofdevelopment. Still,itisveryintriguingtonote that for certain species in the order Trypanosomatida known to control the direction of flagellar beat propagation [153–156], flagellum are sometimes clamped to the cell body at multiple points. This might suggest that micro-organisms capable of flagella wave reversals could admit boundary conditionsthatarebeclosertoaclampedfilamentthanahingedone. Moreimportantly,weshowedthatbase-heavymotoractivityprofilesnaturallyproducebase-to- tipbendingwaveswhiletip-heavyactivitygeneratestip-to-basewavesmoreeasilyforourclamped model cilia. This finding is in strong agreement with the experimental observations of [118], where a lack of outer dynein arm expressions was seen to correlate with dramatic decrease in tip-to-basewaveforms. Andtheexistenceofabi-stablezoneforbothwavedirectionalsoprovides amechanisticexplanationforhowflagellawavereversaleventscanoccur(seeFig.4.9). While our model generates bending waves that are correlated with switch-inhibition of motor activation similar to experimental observations [44], other evidence points to the incompleteness of this sliding control model. In sliding-control model, the oscillation amplitude increases with increasing motor activity and the frequency decreases. These trends are in striking contrast to 86 experimentalobservationsofcilia,wheretheoscillationamplitudeandfrequencyundergosaturated growth with increasing ATP levels [42]. Moreover, in this sliding feedback model, the waveforms differfromthebeatingpatternsobservedincilia. Toproducecilia-likeoscillations,[115]introduced curvature control and bias in the molecular motor kinetics. The mechano-chemical processes that would form the basis of such bias are not clear. Based on these discrepancies, we are inclined to conclude that the sliding feedback control, although may play a role in the molecular motor kinetics,isunlikelytobetheonlymechanismdrivingciliaoscillations. One possible solution is to combine curvature-control and sliding-control simultaneously to get the best of both worlds. With minimal parameter tuning, e.g., a linear weight between two control mechanisms, it is possible to produce base-to-tip waves with a better frequency response as a function of total activity for smaller Sperm number, while maintain the direction switching ability for larger Sperm numbers. This would also qualitatively conform to the observation of “static” modes for short cilium [43]. Moreover, such a combination can also be motivated from a non-ad-hoc, mechanistic perspective, as the wildly accepted geometric clutch theory is shown to beacombinationofslidingandcurvaturecontrolinalinearizedanalysis[53]. To date, mathematical models of cilia have generated a gallery of potential mechanisms that canleadtosustainedoscillations. Recentstudieshavebeguntoaddresstherelativemeritsofthese candidate mechanisms in comparison to experiments [13, 14, 58, 157–160]. However, at present, thereisnoexperimentalbenchmarkofciliaoscillationsunderagreedupon‘nominal’and‘altered’ conditions that can be used as a testbed for theoretical models. To evaluate mathematical models, we need to compare their behavior to experimental results, but, more importantly, we need to test their predictive ability under altered conditions. The analysis presented in [157] of oscillation frequency versus cilium length in the Chlamydomonas is an example. Four theoretical models 87 were tested in [157], none of which was able to reproduce the experimentally-observed trends in frequency versus length. Until we, as a community, establish clear experimental benchmarks that canbeusedtoscrutinizetheoreticalmodelsfortheirabilitytoreproduceciliabeatingpatternsand to predict changes in these patterns under altered conditions, the regulation of molecular motor activityinciliaoscillationsremainsanopenproblem. 88 Chapter5 Coordinationacrossheterogeneoustissue This chapter explores how the phased oscillator model of cilia informs us about the spontaneous emergence of metachronal waves, especially when they are organized in a heterogenous, patchy fashion just like how real ciliary carpets appear in-vivo. We systematically demonstrates how individualciliabeatsandtissue-leveldisorderaffectsthelong-termsynchronizationdynamicsand resultingflowpatternsinciliarycarpets. 5.1 Transportfirst,coordinate second In [4], we used both large scale numerical simulation and a coarse-grained continuum theory to predict the instability of the noisy isotropic states, and showed that metachronal wave states will spontaneously emerge in a very robust manner; see Fig. 5.1. Specifically, since the isotropic state expressed in the moment field coordinates (Eq. (2.40)) satisfies 0 (x) = 1 and = 0 for all other , it is easy to show that it is a periodic steady state or an equilibrium state; see [4, 161]. Thus, linearstabilityanalysisaboutsmallperturbationsabouttheisotropicstateresultsinlinearevolution equationsthat,whenexpressedinFourierspace,giveaneigenvalueproblemoftheform ¤ ˆ = ∞ ∑︁ =−∞ ˆ ˆ , (5.1) 89 Figure 5.1: Spontaneous emergence of metachronal wave states. A. Growth rate of different Fourier modes of the Kuramoto order field (Eq. (2.40)) as predicted by the linearized continuum theory (colormap) and discrete simulations (black dots). B. Starting from a random Kuramoto order field, successive (linear) application of the Fourier mode growth rates produces results that are qualitatively indistinguishable with those produced using the direct numerical simulations of Eq.(2.38). C.TheKuramotoellipse(Fig.2.3C)obtainedfrom200realizations(lightgrey)andthe mean (dark gray) of direct numerical simulations, with the corresponding results of part B shown inblack. D. TimeevolutionofeccentricityandangleoftheKuramotoellipse. Seealso[4]. where(ˆ·)indicatesFouriertransform,(·)markstheperturbationvariable,and ˆ dependsonthe FouriertransformoftheStokesletkernel;see[4,161]. Moreover,since decaysinsizequicklyas grows for states near the isotropic condition, Eq. (5.1) can be closed and computed numerically viatruncation. In Fig. 5.1A, we report the growth rate(k) (maximal real part of the eigenvalues of ˆ (k)) as a colormap over all wavenumbersk for the first and second force harmonics shown in Fig. 2.1. Evidently,theisotropicstateisunstableforalmostallwavenumbers( > 0),withstrongerinstability exhibitedbythesecondforceharmonic. Wesuperimposedatafromasinglerealizationofdiscrete cilia simulation (black dots) using (2.38); the dots represent the wavenumbers at which| ˆ 1 (k,)| 90 is maximal (top 10%). Note the agreement between the simulation result and the maximal growth rateinthe continuumtheory. To further challenge the theory, starting from the initial state used in the sample simulation, wecalculate ˆ 1 (k, 0),andintegrateforwardintime ˆ 1 (k,)= ˆ 1 (k, 0)exp((k)) accordingtothe linear growth rate from the continuum model. We transform back to physical space and compare Real[ 1 (x,)] to that obtained from the sample simulation (Fig. 5.1B) and again find remarkable agreement, indicating that the wave patterns are governed by the fastest growing linear mode. Lastly,wecalculatetheKuramotoellipsebasedonMonte-Carlocomputationsusingthecontinuum model growth rate with 200 initial isotropic states; the average is shown in dark grey in Fig. 5.1C and the range in lighter grey while results from the sample simulation are superimposed in black dots. The time evolution of the ellipses’ eccentricity and angle are shown in Fig. 5.1D using the same color convention. The ellipse eccentricity correlates with the instability growth rate and its angle in physical space is orthogonal to the direction of maximum growth rate in Fourier space. We conclude that the fastest growing wavenumbers in the linear continuum theory are predictive ofthelong-termtravelingwavepatternsinthenonlineardiscretemodel. A carful examination of the above results shows that second forcing harmonics, an abstract representation of eccentricity of cilium beats, produces a stronger metachronal wave instability; also see Fig. 5.2A. Indeed, this difference could be why [77]’s incomplete report that only second harmonicforcesleadtometachronalwaves. However,astrongercoordinationdoesnotnecessarilytranslatetobetterfluidtransportfunction. AsshownbyFig.5.2B,firstharmonicforcing,whichcapturestheessentialfeatureofeffectiveand recoverystrokes,dramaticallyoutperformssecondharmonicforcingintermsofbulkfluidtransport. This is because second harmonic forces possesses an antipodal symmetry, which should not be 91 Figure 5.2: Cumulative flux for different harmonic forcing. A. Snapshot of 151× 151 steady state rotor configurations for (I) first harmonic forcing ( 1 = 1, = 0), (II) mixed forcing with larger first harmonics ( 1 = 2/3, 2 = 1/3, = 0), (III) mixed with larger second harmonics ( 1 = 1/3, 2 = 2/3, = 0), and (IV) pure second harmonic forcing ( 2 = 1, = 0). Tinted portionin(IV)indicatesmetachronalwavegoinginopposingdirections. Colorbarindicatescosine of phase angles; see Fig. 2.2. B. Per cilium cumulative flux induced by the rotor carpet over 2000 units of time; see Eq. 2.42. Clearly, first harmonic components are responsible for pumping. C. Two-pointspatialcorrelationofthephaseanglesasafunctionofdistanceinunitsofrotorsspacing. able to distinguish metachronal waves that travel in exactly the opposite directions by symmetry argumentsandlinearityofStokes’flow! Infact,wavesseparatedbythegrainboundarieshighlighted inpanel (IV)ofFig. 5.2Atravelpreciselyagainst wavesinother partsofthe periodicdomain,and explainsthelowcumulativefluxreportedinFig.5.2B.Nevertheless,asconfirmedbothqualitatively and quantitatively in Fig. 5.2A and C, any additional second harmonic forcing components will increasethecorrelationbetweenneighboringrotors. 92 5.2 Patchycarpetcanenhance metachronal order Depending on the exact initial conditions, defects and grain boundaries emerge spontaneously alongside metachronal order for a particular realization of our rotor carpets. As defect formation dependsstronglyontheinitialvariationoftherotorphases,itisintuitivetobelievethatasmallerset ofrotorscanleadtomuchmorecoherentmetachronalwaves. Infact,comparingourresultsbased on151x151rotorsandwavesreportedin[77]formedby11×11rotors,wecanseeaclearincrease in defect complexity. Therefore, it is reasonable to hypothesize that the typical patchy layout of a actualciliarycarpetcouldalsohelpreducetheformationofdefectsandgrainboundaries. Totestthishypothesis,weapplythetissue-levelorganizationparametersintroducedinSection 2.3.1 to our simulated rotor carpet. Specifically, we reduce coverage ratio and vary patch wavenumber ,butkeepthepatchvariance to0forsimplicity. InFig.5.3AandB,weillustratehowthemetachronalorderchangeinthepresenceofpatchiness as the rotors move according to first or second harmonic force profiles, respectively. While it is difficult to note clear differences for the first harmonic forcing results, it is very clear that the introduced ‘air-gaps’can stronglyshape thegrain shapesdue toopposing travelingwaves: defects are much more likely to emerge diagonal to the patches so that waves tend to originate from the edges and produce more predictable patterns. The tinted region signal waves that travel in the opposingdirection(totheleft). Quantitatively, the two-point correlation function reveals that meta-order also emerged for the caseoffirstharmonicforcing;seeFig.5.3C.Ascoverageratiodecreases,weseeperiodicpatterns thatindicatedistantpatchesarecoordinateorevensynchronizedwitheachother. Asaresult,wesee 93 thatthecumulativefluxgeneratedperciliumcanactuallyincreasewhilepatchinessaredecreasing the total number of active rotors present in the periodic domain (Fig. 5.3D)! However, this result should be interpreted with care as our plain circular rotors, being minimalist abstractions of real cilia,cannotpumpanyfluidindividually;whatwearemeasuringis,albeitimportant,essentiallya higher order synchronization effect. Moreover, changes in the asymptotic ∥Q∥ value might not be realized by the trajectory of real Brownian particles; e.g., while this net flux measured at infinite height is always two-dimensional, a group of Brownian particles located near the cilia/rotors can gainaslow driftinthe+ directionduetotheno-slipwallatthebase. Figure 5.3: Patchiness can enhance pumping and metachronal order formation. A. Visual comparison between fully covered ( = 1) and patchy ( = 0.5, = 4) for rotors under first harmonic forcing. B. Visual comparison between fully covered (= 1) and patchy (= 0.5, = 4)forrotorsundersecondharmonicforcing. Highlightedrotorsformmetachronalwavesthatmove in opposite directions. C. The distant patches can coordinate with each other across the ‘air-gaps’ as shown by the two-point correlation of the rotor phases. D. As coverage ratio decreases, the cumulativefluxgeneratedperciliumfirstdropsbelowthefullycoveredcarpetthenincreasesabove it. E. For illustrative purposes, the same phenomenon can also be observed for second harmonic forcing. It is evident that the introduced patchiness speeds up the formation of metachronal order (notethechangeinsimulationtimescale). 94 5.3 Discussion We developed numerical algorithms and continuum theory for analyzing phase coordination in large scale ciliary carpets. Existing simulations have been limited to a few hundred cilia [61, 62], which is an order of magnitude smaller than the number of cilia on microbes like the Volvox or Paramecium[162,163],andmanyordersofmagnitudesmallercomparedtothenumberciliafound, for instance, in mammalian airways [67]. Our work provides a framework to begin to explore the structure-to-function maps from individual cilia oscillations to emergent coordination and fluid pumpinginlargescaleciliarysystems. Weshowedinbothnumericsandtheorythat,inlargescaleciliarysystems,spatially-homoge- neous isotropic states that do not pump fluids are unstable when cilium are subjected to periodic forcing adjustments. In particular, we showed forcing that represent the effective and recovery stroke of cilium beat and vary only once along the beating cycle (first harmonic) lead to the emer- gence of wave patterns that break time reversal symmetry and produce net pumping, a result that was not reported in [77]. We also confirmed the results of [77] and see forcing that encodes the eccentricity of cilium beat produces much stronger coherent waves. Notably, in all cases, the bulk property of the emergent waves are independent of initial conditions. Our theory in the limitofinfinitenumberofcilia predictsquantitativelythelong-termtravelingwavepatternsinthe nonlinear numerical simulations, indicating that linear instabilities dominate the emergent waves. Our findings emphasize the importance of the single-cilium beat kinematics to both tissue-level coordinationandfluidpumping. 95 We further demonstrated how patchy tissue organizations can surprisingly accelerate the for- mation of global metachronal order and outperform the uniformly covered tissues. Our findings are consistent with recent suggestion that spatial heterogeneity in cilia organization and tissue architecture enhance particle clearance, and go beyond these results to simultaneously probe the effectofspatialheterogeneityonciliacoordinationandcilia-inducedflows. Ciliatedtissuesinmajororgans,suchasinthemammalianlungs,brains,andreproductivetracts, havespecializedfluidpumpingfunctions. Specializationcanbecompromisedifmultiplestatesof coordination co-exist. In contrast, ciliated organisms that rely on cilia for locomotion require the existence of multiple coordination states to perform distinct swimming and turning motions [70, 164–166]. Multiple coordination states have been reported in cilia models consisting of filaments and rowers that beat in a plane perpendicular to the substrate [55, 72, 78, 79, 167], reminiscent of cilia with planar beat kinematics. In addition to the mere size of the ciliary system, the cilium beatkinematicsisamarkeddifferencebetweenthelatterstudiesandourstudy. Inourmodel,each rotor beats cyclically in a plane parallel to the substrate, akin to cilia with out-of-plane effective and recovery strokes. We found that metachronal coordination is a global attractor that is robust to large perturbations, as long as the single-cilium beat and cilia organization within the tissue remainintact. Ourresultssuggestthatfunctionalspecificityoftheciliatedtissueisinterlacedwith the architecture of the ciliary system. These findings corroborate existing assays that examine the cilium beat kinematics to distinguish between diseased and healthy cilia [168–170], and suggest complementarytissue-levelparameters(e.g.,cilialatticestructure)tobegintoevaluatemetachronal coordinationandlinkittofluidpumpinginhealthyanddysfunctionalciliatedtissues [171,172]. Wemadeafewsimplificationsthatcanbereadilyrelaxedinfutureextensionsofthisstudy. For example,ourtheorycanbereadilyextendedtoincorporatenon-circularlimit-cyclerepresentations 96 of individual cilia obtained from experimental data of cilia beating patterns [69] and non-uniform lattice structures reconstructed from high-resolution images of ciliated tissues [67]. These exten- sions will allow us to quantitatively calibrate the model for specific ciliated tissues, and to predict howperturbationstothesingle-ciliumkinematicsandciliaorganizationwithinthetissueaffectthe emergenceoflargescalecoordinationandfluidpumpingfunction. 97 Chapter6 Universaldesignrulesforciliarypumps 6.1 PerformanceofCiliary Pumps Our porous media model generalizes the classical impermeable sheet model of ciliated sur- faces [127]. In the limit of ℎ→ 0 while keeping ¯ ℎ constant, the ciliary layers collapse into infinitesimallythinciliarycarpetsthatmoveaccordingtotheprescribedlongitudinaldensitywave, such that the governing equation reduces to Stokes flow with slip velocity boundary conditions at thetopandbottomwalls. Thisscenarioisconceptuallyequivalenttothestudyofciliaryflowbased on the envelope model [128]. In this limit, mean flow speed is maximized at the steady state in the absence of adverse pressure since the flow profile can remain constant for all. However, such pumps are powerless in face of steady adverse pressure; the flow profile reduces to that of Poiseuille’sflowwithaconstantoffset;seeFig.6.1a-b. Similarly, if cilia, now represented as active Brinkman layers, remain restricted in height with respectto thechannellumen height(smallℎ/),mean flowspeed canremain competitivewith respect to the theoretical limit proposed by the envelope model in absence of adverse pressure (Fig.6.1c). Whilethesepumpscanperformbetterinfaceofadversepressuregradient(Fig.6.1d), thenetflowdirectionstillreversesforlargeenough |Δ|. Incontrast,whenciliaheightℎapproaches lumenheight,pumpingperformanceundernoadversepressureissacrificedforacontinuedability 98 ΔP U x U x a b e f 0.3 30 c d Figure 6.1: Porous media pump generalizes the envelope model. (a-b). Instantaneous velocity field u(,,) and flow profile () produced by the envelope model of ciliated duct, where fluid is pumped by slip velocities at the impermeable boundaries. Pumping is maximized without adverse pressure, but even with a relatively small adverse pressure gradient ofΔ/ −1 =−1, netflowreversesdirection. (c-d). Whenpumpingwithaporouslayerofcilia(ℎ/= 0.2),adversepressuregradientcanonly reverse net flow direction in regions where cilia is absent. Density distribution of cilia divided bythereferencedensity ¯ isshowningrayscale. (e-f). theporousmediamodelcanmodelciliated ducts with long cilia (ℎ/ = 0.9) that can also pump fluid against large adverse pressure gradient Δ. ℎ= 1000inc-f,Δ/ −1 =−50ind,f,and = 0.159≈/2 forallpanels. topumpagainstverylargeadversepressuregradient(Fig.6.1e-f). Thisisbecausewhiletheadded adversepressureforcescaneasilypushtheflowintheStokesflowlayerbackward,theyarerelatively ineffectiveagainsttheactiveresistanceinsidetheBrinkmanlayers. This ability to withstand incoming adverse pressure is mediated by the value of average cilia density ¯ . InFig.6.2a,weshowhowmeanflowspeed changesasafunctionofciliatedfraction ℎ/ and ¯ . On the left panel where adverse pressure is set to 0, we see that the flow speed 99 Flow Speed U [μm/s] ΔP/μωL -1 = −50 -1000 -100 0 100 0 20 40 60 80 100 120 Ciliated Fraction h/H Ciliated Fraction h/H 0 1 0.2 0.4 0.6 0.8 0 1 0.2 0.4 0.6 0.8 ρ c h = 10 ρ c h = 10 2e 2d 2f ρ c = 1000 ρ c = 100 ρ c = 10 Flow Speed U [μm/s] 0 20 40 60 80 100 120 ρ c h = 1000 1 5 10 100 1 5 10 100 ΔP = 0 2c 10 -10 -1000 -100 0 100 10 -10 ρ c h = 1000 a b Figure 6.2: Cilia density, ciliated fraction, and the material constraint. a. Net flow speed as a function of ciliated fraction ℎ/ for increasing values of ¯ at zero (left) and finite adverse pressure ofΔ/ −1 =−50 (right) in gray scale. Dashed line shows the theoretical flow speed limitof⟨ ⟩ . b.(ℎ/) withdifferentvaluesofconstant ¯ ℎ indarkbluegradient. Under large adversepressure,positiveflowisonlypossibleathighciliatedfractionforlarge or ℎ. Square marks shows the corresponding flow speed for Fig. 6.1c-f. = 0.159 and/ = 1 for all panels. Notethescalechangesin-axis(bi-directionallogmappingisusedwhereconvenient[173]). increases with higher and higher cilia density, but remain bounded above by the analytical speed limit of⟨ ⟩ . In particular, note that even though we increase ¯ by one order of magnitude from 10 to 100, the resultant speed increments in a much slower fashion. Similarly on the right panel, we show that flow speed again sway towards positive flow direction as ¯ increases from 10 to 1000 even in face of an adverse pressure gradient of Δ/ −1 =−50. Note that only higher ciliated fraction combined with large enough cilia density can maintain flow in the desired positive direction. However, fully-covered quasi-1D channel do not pump fluids due to continuity constraints, and only act as additional passive resistance when pressure gradients are present; see also App. 2.4.5. More importantly, if we wish to compare pumps with different ciliated fractions 100 0 −20 −5 Adverse Pressure ΔP/μωL -1 −10 −15 −200 0 −1e3 −600 200 h/H = 0.2 h/H = 0.4 h/H = 0.6 h/H = 0.8 Flow Speed U [μm/s] -2000 0 -1000 0 0.2 0.4 0.6 0.8 1 ∆P/μωL -1 = −20 ∆P = 0 H/L = 1 Ciliated Fraction h/H 500 ∆P/μωL -1 = −20 ∆P = 0 h/H = 0.5 Aspect Ratio H/L 100 0 -100 -10000 0.1 1 10 2 1 -1000 a b c Figure6.3: Meanflowspeedasafunctionofadversepressureandlumenaspectratio. a . Flow speed versusΔ for various ℎ/; is linear inΔ but nonlinear in ℎ/. b. When adverse pressure is present, the flow direction at lower ciliated fractions ℎ/ reverses, and the magnitude of the reverse flow is larger as the lumen diameter increases. c. For sufficiently tall channel (≫ ), negativeΔ produces a backward flow speed that is quadratic with respect to lumen diameter,inlinewithPoiseuilleflowinarectangularchannel. ¯ ℎ= 1000and = 0.159forall panels. fairly by keeping input power constant, we need to keep the material constraint ¯ ℎ as a constant. ShowninbluelinesinFig.6.2,weseethatforsufficientlylargeconstraintconstant,theflowspeed performance approaches that of large cilia density. Together these results indicate that there exists apointofdiminishingreturnforhighciliadensityormaterialconstraintconstant. Next we systematically explore the effect of the ciliated fraction ℎ/, lumen aspect ratio/ and adverse pressure has on mean flow speed under the material constraint ¯ ℎ = 1000. In Fig. 6.3a, we confirm that Darcy’s law applies and mean flow speed is linear with respect to Δ. However, since is non-linear in terms of ℎ/, at larger ℎ/, decays more slowly as Δ increases, hence the ability to produce positive flow at higher adverse pressure. Thus even thoughlargerciliatedfractionissub-optimalinabsenceofadversepressuregradients,theyexhibit better performance for sufficiently large Δ; see Fig. 6.3b. In addition, as lumen aspect ratio / increases, it becomes harder to pump againstΔ as pressure-driven flow is quadratic in (Fig.6.3c). Asaresult,toconstructaciliarypumpthatremainfunctionalagainstverylargeadverse 101 h/H = 0.99 H/L = 0.05 h/H = 0.01 H/L = 10 -10 8 -10 4 -10 0 -10 -4 10 -6 10 -4 10 -2 10 0 -10 6 -10 2 -10 -2 Adverse Presssure ΔP/μωL -1 Area Flow Rate Q [mm 2 /s] 10 -3 10 -1 10 -5 h/H = 0.5 H/L = 1 Figure6.4: Designinganoptimalciliarypump. Foranyfixedgeometricaldesign(ciliatedfraction ℎ/ andlumenaspectratio/),flowrateandadversepressurehasatrade-offakintoallpumps in general. Thus the most efficient operating point for each geometric design is at the corner of eachcurve. ¯ ℎ= 1000and = 0.159. pressure,bothnarrowductopenings(relativetowavelength)andhighciliatedfractionarerequired simultaneously. And indeed such trend is evident in the morphospace of (,ℎ/) among ciliated pumpsobservedinbiology[5]. Furthermore, we can demonstrate this design criterion by examining our cilia-inspired pumps in the functional space of flow rate vs. operating pressure, as commonly done in the analysis of engineered pumps (Fig. 6.4) [90]. Specifically, for a set of fixed geometric parameters such as ciliated fraction, lumen aspect ratio, and material constraint, the pump model produces smaller (area) flow rate as adverse pressureΔ increases in a sharply nonlinear fashion. And to pump at maximum efficiency , pumps of all design need to operate near the corner points as seen in Fig.6.4. 102 Figure 6.5: Computational model that maps CD morphology to fluid pumping. A. Cilia are modeled collectively as an active porous medium with prescribed unidirectional traveling density waves (,) that produce instantaneous flows u(,,) with cross-sectional flow profile = ∫ 0 u·ˆ e and net flow speed = ∫ 0 . Ciliary carpet (ℎ/ = 0.2) and flame (ℎ/= 0.9)designsatΔ= 0(IandII)andΔ=−50 −1 (IIIandIV)atwhichonlytheflame design pumps fluid in the desired direction. B. Flow speed versus ℎ/: under large adverse pressure,positiveflowisonlypossibleathighcilia-to-lumenratio. C. versusΔ: atlargerℎ/, decays more slowly asΔ increases, hence the ability to produce positive flow at larger adverse pressure. D. Flow speed over the entire morphospace (,ℎ/) atΔ= 0. E. Adverse pressure gradientsΔ for which a net flow speed of = 100m s −1 is achievable. Contour lines indicate that flame designs continue to pump as Δ increases. In all panels, ¯ ℎ = 1000, = 20 Hz and = 100m,andinpanelsa–c,= 100m. 6.2 Highercilia-to-lumen ratio gives higher pumping pressure We first considered two systems with ciliary carpet ( ℎ/ = 0.2) and flame ( ℎ/ = 0.9) designs and same lumen diameter = 100m, and we set ¯ ℎ= 1000 as an upper bound of total ciliary 103 Figure6.6: OptimalCDdesigns. A.PumpingperformanceoftwoextremeCDdesign(solidblack lines, ¯ ℎ= 1000),andthelarvaceanciliatedfunnel(blue, ¯ ℎ= 200)andesophaguscarpet(pink, ¯ ℎ = 40). Flow rates and adverse pressures that maximize pumping efficiency ∝ 3 / ¯ ℎ at different ciliary material constants ¯ ℎ∈{1, 5, 10, 100, 1000} are superimposed (dotted curves). B.OptimalCDmorphologiesinthe(,ℎ/)morphospace(whitelines)exhibitashiftfromcarpet toflamedesignsatlargervaluesofmaximumadversepressure(orangecolormap),accompaniedby adecreaseinoptimalflowrate(greencolormap). Highestflowratesareachievedbycarpetdesigns with high lumen diameter and low cilia-to-lumen ratio, and highest pressure is generated by flame designswithlowerlumendiameterandhighercilia-to-lumenratio. BiologicaldatafromFig.1.2A are superimposed with biological function, as known, indicated by the shape of the corresponding marker. Gray-outregionindicateslowperformancescore(< 50%),seeFig.D.5. material in all surveyed CDs (Fig. 2.5). The carpet design generated higher flow speeds at zero adversepressure,butfailedandallowedreversalofluminalflowinthepresenceofadversepressure Δ(Fig.6.5aIandII),whereaspositiveluminalflowremainedrobustintheflamedesign(Fig.6.5a III and IV). In fact, at zero adverse pressure, positive flows were generated by all cilia-to-lumen ratiosℎ/ (Fig.6.5b),exceptinthetwolimitingcasesℎ/→ 0andℎ/→ 1(Methods§2.4.5). Maximal flow speeds on the order of 100 m s −1 were achieved at small ℎ/, consistent with experimentallyobservedflowspeedsinciliarycarpets,suchasinairwayepithelia [174,175]. The flame design, by contrast, was more effective in maintaining positive flow speeds under adverse pressure for a wide range of Δ values (Fig. 6.5c); this agrees with the hypothesized filtration functionofciliaryflameswhichinvolvesfluidtransportinthepresenceofadversepressures[88]. 104 6.3 Optimalductdesign depends on pressure requirements This trade-off between net flows and sustained adverse pressure is best appreciated in the pump functionspace(,Δ),where=isthenetcross-sectionalareaflowrate(Fig.6.6a): aciliary pump(,ℎ/) operates along a characteristic curve connecting its maximal generated flow rate atΔ = 0 to the maximal adverse pressureΔ sustained without causing flow reversal = 0. Each ciliary pump functions along its own characteristic curve, depicted in solid lines in Fig. 6.6a for representative carpet and flame designs. Intuitively, a ciliary pump (,ℎ/) performs best when operating near the transition from maximal flow rate to maximal pressure (depicted as stars forthelarvacean). Forthelarvaceanfunnel,ourresultsatthesetransitionalvaluessuggestthat,whenpumpingin water with a CBF of = 60 Hz and ciliary wavelength of = 50 m, it can provide a maximal flowrateofabout0.2 lhr −1 atΔ≈ 4Pam −1 ,whichamountstogenerating640Paofpressure for the 160m long duct. In contrast, a similar calculation estimates that the larvacean esophagus (L = 100 m, = 20 Hz) can provide a flow of nearly 40 l hr −1 against a pressure of less than 0.1Paover its2500mlength. To probe optimal ciliary pump designs, we introduced a pumping efficiency defined as the ratio of the rate of change of the fluid mechanical energy to the power input by the ciliary layers (Methods §D.2). We maximized efficiency over the 3D space(,ℎ/,Δ) for fixed power input. We repeated the same process for different levels of power input, corresponding to ciliary constants ¯ ℎ ranging from 1 to 1000 to reflect the entire span found in our morphometric survey (Fig. 2.5). Results from our optimization routine (dotted lines in Fig. 6.6A) are indeed located at 105 thetransitionsfrommaximalflowrate tomaximaladversepressureΔ. Displayingtheseresults on the morphospace(,ℎ/) (Fig. 6.6B,C), we found that the most efficient designs (solid white lines labeled by the corresponding value of ¯ ℎ) exhibited a shift from carpet to flame designs at larger values of maximum adverse pressure (orange colormap). This shift is accompanied by a decreaseintheoptimalflowrate(greencolormap). Specifically,optimaldesignsoverthefullrange of empirically-derived values of ¯ ℎ span a band from carpets with large lumen diameter and small cilia-to-lumen ratioℎ/ that can generate high volume flows against weak adverse pressure gradients,toflamedesignswithsmalllumendiameter andhighcilia-to-lumenratioℎ/ thatcan pumpfluidsatlowflowratesagainststrongeradversepressuregradients. Designsthathavelow andlowℎ/ aresub-optimalbecause,atlowℎ/ andΔ (Fig.6.6B),foraconstantpowerinput, ducts with larger are always more efficient in transporting flows (Fig. 6.6C). Designs that have high andhighℎ/ arealsosub-optimalbecause,athighℎ/ andΔ,foraconstantpowerinput, thelumendiametercannotincreasewithoutloweringtheciliadensityand,inturn,theinputpower density, thus limiting the value ofΔ that can be sustained, potentially leading to flow reversal. Putdifferently,althoughtheoreticallyincreasing athighℎ/ isviable,therequiredciliadensity andinputpowergrowdisproportionatelyandlikelybeyondbiologicallimits. 6.4 Discussion We presented and analyzed the performance of a cilia-inspired active porous-media pump model. Our results demonstrate mathematically the importance of morphological parameters for ciliary pumps, and has strong implications for the understanding ciliated ducts and tissue in biological settings [5]. We show that permeability is not only natural but necessary for cilia to function as 106 pumps, especially when faced with adverse flow conditions. Specifically, only ciliary pumps with largeenoughciliadensity,highenoughciliatedfraction,andsmallenoughlumensizecanproduce netpositiveflowagainstlargeadversepressuregradients. Inreality,thissituationfrequentlyoccurs as cilia are responsible for transporting material against forces such as gravity or filter material throughhighlyresistivepathways[87,89]. Thisworkestablishesaquantitativelinkbetweenthemorphologyandthefluidpumpingfunction ofciliatedducts. Asaconsequence,wecannowusestructuralfeaturesofagivenofCDtopredictits relativeabilitytoeithertransportfluidsinbulk,ortoperformfiltration. Caseinpoint,ourfindings support the decades-old hypothesis that the larvacean’s ciliated funnel (Fig. 1.2a) helps maintain the internal fluid volume through unidirectional injection of filtered seawater into the blood [176], instead of being limited to sensory and endocrine functions [177–179]. In the Hawaian Bobtail Squid(Euprymnasocolopes),wecanpredictthattheciliatedconduittothelightorganthathoststhe bacterialsymbiont(VibrioFischeri)transportsfluidsinductandantechamber,whereasfiltrationor pressurizationmayoccurinthebottleneckregion[5]. Furthermore,thebottleneckdiametervaries dynamicallyduringdifferentstagesofcolonization[100],suggestingaswitchbetweenfiltrationand clearance functions that could help maintain the host-symbiont association. Similarly, the ciliated pore canals of the madreporite in sea urchins appear to change dynamically between a carpet and flame-like design. Our results support the proposed functional switch between fluid pump and valve that maintains the inner hydrostatic pressure [180]. We also note that many systems fall into intermediaterangesofthemorphospace,whichneithermaximizeflowratenorpressuregeneration, and tend to balance these trade-offs in a way that optimizes pumping efficiency. It is possible that these hybrid morphologies facilitate additional functionalities, such as chemosensing [181], mixing[82],orsignalling[65]. 107 Our findings may provide general insights into the evolution of waste excretion strategies in animalsandofferafunctionalexplanationwhycilia-basedexcretoryorgansareabundantinanimals notexceedingafewmilligramsinbodymass,whereaslargeranimalsrelyonbloodpressuredriven systems[182]. Ourstudysuggeststwomechanisms. First,oursurveyshowsthemarkedabsenceof CDs with large andℎ/. Since cilia have limited length and emanate at finite density from the tissue surface, the lumen diameter cannot increase without eventually lowering the cilia-to-lumen ratio. Thismaterialconstraintcorrespondsisanalogoustoaconstraintoninputpowerandprovides anupperboundarytothepressureandfiltrationflowrateanindividualflamepumpcangenerate. AlongwithassessingthepumpingfunctionofhealthyCDs,ourstudymayaidinthefunctional interpretation of disease phenotypes [183] and inform the development of drugs that disrupt ex- cretory functions of parasites [184]. Moreover, our approach provides a rational design tool for tissue-engineeredandsyntheticciliateddevicesthatperformfiltrationandflowcontrol,inaddition toexistingcarpet-likeconfigurationsforbulktransport[42]. On the other hand, insights based on our model are also useful for the design of microfluidic pumps inspired by cilia [185–187]. Our results form maps the geometric design parameters that controlthedistributionofactiveelementstofunctionalparametersthatareusefulforperformance and cost-benefit analysis. In addition, despite our choice of symmetric beating profile for each Lagrangian particle,thissimplificationcouldbedesirableforengineeredsystems. A few words on generalizations and limitations of our model. In addition to our bulk flow analysis, it would be useful to study particle clearance or residence time in such active porous media framework, as for many ciliated systems, the ultimate goal is to understand how particles embedded in the fluid medium moves so that we can understand the detailed transport and even sensing capability of such pumps. In order to limit the degree of freedom for analysis, we opted 108 for a prescribed, symmetric, and quasi-1D model of cilia activity. It would be straightforward to incorporate asymmetry in beating pattern and two-dimensional traveling waves, or even introduce recordedciliamotionasinputtoourmodel. Thiscouldalsoresolvethecounterintuitivebehaviorof fullchannelstallat(ℎ=)duetocontinuityequation. Additionally,generalizingourformulation to an axisymmetric form or one that allows different boundary shapes could also bring additional geometric insight on how ciliated pump perform in more realistic conditions. Moreover, the unsteadytermcanalsobeintroducedtostudythetransientbehaviorofsuchciliarypumpsandtake careofthephaselagthatcouldbeintroducedduetohighfrequencyofciliaoscillations[188]. And lastly,ashydrodynamicinteractionandsynchronizationbetweenciliaisalsoofgeneralinterestfor thestudyofciliatedtissue,ciliaactivitycanalsobeintroducedsuchthatonecanstudythetwo-way couplingbetweenciliaandfluidmotion[79,136,161,189,190]. Our survey might be missing alternative CD designs in the “hidden biology” of less studied animals [191]. We also omitted ciliated ducts whose geometry is not well described by our 2D morphospace, such as the ciliated chamber of sponges [103, 104, 192] and some of the unique ciliated spaces found in ctenophores [193–196]. No distinction was made between different types of fluids, such as mucus and water, or between unidirectional and mixing flows. We considered a Newtonian fluid model, heuristic cilia waveforms and material constraints, and fundamental measuresofpumpingperformanceinordertolimitthenumberofvariablesandrevealthedominant structure-functionrelationshipsofciliatedducts. Takentogether,our2Dmorphospaceandfunctionalmaparerepresentativeandrelevantforthe majority of ciliated ducts known to biology and provide a tool for assessing the physiological role andintegrityofciliatedorgans. 109 Chapter7 Conclusionandfutureoutlook Inthisthesis,acascadeofmathematicalandcomputationalmodelsofciliaandciliatedtissuewere constructedinhopesofsolvingsomelongstandingmysteriesofciliaatdramaticallyvaryinglength scales. Fourmainresultsstandasfollows: (i) We showed for the first time that a clamped microfilament subject to Stokesian drag and active follower forces undergoes a fully nonlinear, 3D spinning Hopf bifurcation before a second instability transitioning it into 2D planar beating motions. The results have already inspired or servedasfoundationforseveralworkinthecommunity[1,151,197]. (ii) Building on top of an existing two-state, two-filament model of 2D beating cilia, we demonstrate that sliding feedback control and front-back asymmetry of molecular motor activity can reproduce flagella wave reversal behaviors as seen in experiments. Compounded on our analysisonflagellawaveondifferentboundaryconditionsandgeometricfeedbackmechanism,this frameworkofanalysishaspromisingpotentialtoexplainhowthegreatdiversityofciliaoscillation naturallyoccurdespitetheirhighlyconservedinternalstructures. (iii)Capturingkeycharacteristicsofindividualciliabeatasharmonicforcingonphasedoscilla- tors,wenumericallyandtheoreticallyshowedthatmetachronalwavesspontaneouslyemergefrom isotropic initial conditions, and quantified how effective/recovery strokes governs fluid pumping, andtheeccentricityofciliumbeatdriveslongrangemetachronalcoordinationofciliarycarpet. In 110 addition, we showed that structural heterogeneity at the tissue-level can surprisingly improve the qualityoftheglobalmetachronalorderandassociatedflowfunctions. (iv) Coarse-graining cilia activity in internal ciliated ducts as active porous-media, we showed that the spatial organization of cilia inside the duct and the size of the duct can predict its relative ability to perform bulk transport versus filtration, and thus established a universal design rule for allciliarypumps. 7.1 End-to-endmodelof ciliary tissue functions As we have separted created models that span the nanoscale dynamics of molecular motors to milimeterandcentimeterciliatedtissuefunctions,thenaturalnextstepistocombinedifferentlevel ofabstractionintoacoherentframeworkofanalysis. Inordertoformulatesuchapipeline,thereareseveralindividualchallengesforeachpresented model: (i) for the two-state, two-filament model of single cilium, an extension to three-dimension with arbitrary boundary conditions could greatly increase its utility. The challenge here is a balance between details and performance: while it is straightforward to expand the model back to include the 9-filament structure, it might not be able to provide as much insight as it might cost to simulate. A potential comproise could be to develop a general theory where 3d filament motions aredecomposedassynchronizationoftwoplanarfilaments(see3.5). (ii) Our current approach when faced with large sets of tunable parameters is to reduce quam maxime for the cleanest possible analysis. A different, perhaps more practically minded approach, would be to incorporate experimental data at the appropriate level of this pipeline. For example, 111 thetissue-levelparametersinourphasedoscillatormodelcouldserveasgreatcandidatesforimage basedmeasurementinputs. (iii) Throughout this study, the role of complex fluid and viscoelasticity have been largely ignored. In addition,ciliaare knowntobeimportantsensing organelles,thusincorporatingmulti- physics components such as advection-diffusion of signaling compounds could also be greatly desirable. 7.2 Lifebeyondcilia Our line of research can also provide insights for the broader field of active matter beyond the confines of cilia and flagella. For example, our phased oscillator model of ciliary carpet can be seen as “swarmalators,” as the coordination of the oscillators is strongly coupled to their spatial position. In fact,italsomakesaveryuniquemusicalinstrument(MATLABApplet). Lastly, during the course of this PhD, the following contributions were also made: (i) we studiedtheinversekinematicsofoctopusarmmovementsforthedevelopmentofdecentralizedsoft underwater robotic appendages [198], (ii) we applied modern reinforcement learning techniques (PPO)tobetterunderstandthenavigationandcontrolofthree-linkfish[199],and(iii)weanalyzed the stochastic dynamics of fish schools driven by behavior rules and hydrodynamic interactions (MATLAB Applet). Lifeindeedcontinuestoflourishbeyondcilia! 112 Bibliography 6. Sleigh, M. A. Patterns of ciliary beating. Symposia of the Society for Experimental Biology 22,131–150(1968). 7. Brennen,C.&Winet,H.Fluidmechanicsofpropulsionbyciliaandflagella. AnnualReview ofFluidMechanics9,339–398(1977). 8. Blake,J.R.Amodelforthemicro-structureinciliatedorganisms.JournalofFluidMechanics 55,1–23(1972). 9. Chwang, A. & Wu, T. Y. A note on the helical movement of micro-organisms. Proceedings oftheRoyalSocietyofLondonB:BiologicalSciences178,327–346(1971). 10. Hirokawa,N.,Tanaka,Y.,Okada,Y.&Takeda,S.Nodalflowandthegenerationofleft-right asymmetry.Cell 125,33–45(2006). 11. Meng, D., Cao, M., Oda, T. & Pan, J. The conserved ciliary protein Bug22 controls planar beatingofChlamydomonasflagella. JCellSci127,281–287(2014). 12. Brokaw, C. J. Bend propagation by a sliding filament model for flagella. Journal of Experi- mentalBiology55,289–304(1971). 13. Brokaw, C. J. Thinking about flagellar oscillation. Cell Motility and the Cytoskeleton 66, 425–436(2009). 14. Riedel-Kruse,I.H.,Müller,C.&Oates,A.C.Synchronydynamicsduringinitiation,failure, andrescueofthesegmentationclock.Science317, 1911–1915(2007). 15. Sartori,P.,Geyer,V.F.,Scholich,A.,Jülicher,F.&Howard,J.Dynamiccurvatureregulation accounts for the symmetric and asymmetric beats of Chlamydomonas flagella. eLife 5, e13258(2016). 16. Brokaw, C. J. Computer simulation of flagellar movement: I. Demonstration of stable bend propagation and bend initiation by the Sliding Filament Model. Biophysical Journal 12, 564–586(1972). 17. Brokaw, C. J. & Rintala, D. R. Computer simulation of flagellar movement. III. Models incorporating cross-bridge kinetics. Journal of mechanochemistry & cell motility 3, 77–86 (1975). 18. Brokaw, C. J. Computer simulation of flagellar movement. VI. Simple curvature-controlled modelsareincompletelyspecified. BiophysicalJournal 48,633–642(1985). 19. Hines, M. & Blum, J. J. Three-dimensional mechanics of eukaryotic flagella. Biophysical Journal 41, 67–79(1983). 113 20. Murase, M., Hines, M. & Blum, J. J. Properties of an excitable dynein model for bend propagationinciliaandflagella. JournalofTheoreticalBiology139, 413–430(1989). 21. Hilfinger,A.,Riedel-Kruse,I.,Howard,J.&Jülicher,F.HowMolecularMotorsShapeThe FlagellarBeat.BiophysicalJournal 96, 196a(2009). 22. Hilfinger, A., Chattopadhyay, A. K. & Jülicher, F. Nonlinear dynamics of cilia and flagella. PhysicalReviewE 79, 051918(2009). 23. Camalet, S. & Jülicher, F. Generic aspects of axonemal beating. New Journal of Physics 2, 24(2000). 24. Lindemann, C. B. A model of flagellar and ciliary functioning which uses the forces trans- versetotheaxonemeastheregulatorofdyneinactivation.CellMotilityandtheCytoskeleton 29,141–154(1994). 25. Lindemann,C.B.A"geometricclutch"hypothesistoexplainoscillationsoftheaxonemeof ciliaandflagella. JournalofTheoreticalBiology168,175–189(1994). 26. Lindemann, C. B. Geometric clutch model version 3: The role of the inner and outer arm dyneinsintheciliarybeat.CellMotilityandtheCytoskeleton52, 242–254(2002). 27. Bayly, P. V. & Dutcher, S. Steady dynein forces induce flutter instability and propagat- ing waves in mathematical models of flagella. Journal of the Royal Society Interface 13, 20160523. issn:1742-5689(2016). 28. Dowell, E. H., Curtiss, H. C., Scanlan, R. H. & Sisto, F. A modern course in aeroelasticity (Springer,1989). 29. Paidoussis, M. P. Fluid-structure interactions: slender structures and axial flow (Academic press,1998). 30. Paidoussis, M. P. The canonical problem of the fluid-conveying pipe and radiation of the knowledge gained to other dynamics problems across applied mechanics. Journal of Sound andVibration310,462–492(2008). 31. Han, J. & Peskin, C. S. Spontaneous oscillation and fluid–structure interaction of cilia. ProceedingsoftheNationalAcademyofSciences115, 4417–4422(2018). 32. De Canio, G., Lauga, E. & Goldstein, R. E. Spontaneous oscillations of elastic filaments induced by molecular motors. Journal of the Royal Society Interface 14, 20170491. issn: 1742-5689(2017). 33. Satir,P.&Sleigh,M.A.Thephysiologyofciliaandmucociliaryinteractions.AnnualReview ofPhysiology52, 137–155(1990). 34. Satir,P.&Christensen,S.T.Overviewofstructureandfunctionofmammaliancilia.Annual ReviewofPhysiology69, 377–400(2007). 114 35. Woolley,D.Theratspermtailafterfixationbyfreeze-substitution inTheFunctionalAnatomy of the Spermatozoon: Proceedings of the Second International Symposium, Wenner-Gren Center,Stockholm,August197323(1975),177. 36. Gibbons,B.H.&Gibbons,I.Propertiesofflagellar"rigorwaves"formedbyabruptremoval of adenosine triphosphate from actively swimming sea urchin sperm. The Journal of Cell Biology63, 970–985(1974). 37. Goldstein,S.F.Startingtransientsinseaurchinspermflagella. TheJournalofCellBiology 80,61–68(1979). 38. Gibbons,I.Ciliaandflagellaof eukaryotes. JCellBiol 91, 107s–124s(1981). 39. Movassagh, T., Bui, K. H., Sakakibara, H., Oiwa, K. & Ishikawa, T. Nucleotide-induced global conformational changes of flagellar dynein arms revealed by in situ analysis. Nature Structural&MolecularBiology17, 761(2010). 40. Lin, J. et al. Building blocks of the nexin-dynein regulatory complex in Chlamydomonas flagella. JournalofBiologicalChemistry286, 29175–29191(2011). 41. Lin, J., Okada, K., Raytchev, M., Smith, M. C. & Nicastro, D. Structural mechanism of the dyneinpowerstroke.Naturecellbiology16, 479(2014). 42. Chen,D.&Zhong,Y.Acomputationalmodelofdyneinactivationpatternsthatcanexplain nodalciliarotation.BiophysicalJournal 109,35–48(2015). 43. Geyer, V. F., Sartori, P., Friedrich, B. M., Jülicher, F. & Howard, J. Independent control of the static and dynamic components of the Chlamydomonas flagellar beat. Current Biology 26,1098–1103(2016). 44. Lin, J. & Nicastro, D. Asymmetric distribution and spatial switching of dynein activity generatesciliarymotility.Science 360,eaar1968(2018). 45. Tani,T.&Kamimura,S.Reactivationofsea-urchinspermflagellainducedbyrapidphotol- ysisofcagedATP.JournalofExperimentalBiology 201,1493–1503(1998). 46. Brokaw, C. J. Control of flagellar bending: a new agenda based on dynein diversity. Cell MotilityandtheCytoskeleton28, 199–204(1994). 47. Mukundan, V., Sartori, P., Geyer, V., Jülicher, F. & Howard, J. Motor regulation results in distal forces that bend partially disintegrated Chlamydomonas axonemes into circular arcs. BiophysicalJournal 106, 2434–2442(2014). 48. Jülicher, F. & Prost, J. Cooperative molecular motors. Physical Review Letters 75, 2618 (1995). 49. Brokaw,C.J.Bendpropagationalongflagella. Nature 209,161(1966). 115 50. Hines,M.&Blum,J.Bendpropagationinflagella.I.Derivationofequationsofmotionand theirsimulation.BiophysicalJournal 23, 41–57(1978). 51. Dillon,R.H.&Fauci,L.J.Anintegrativemodelofinternalaxonememechanicsandexternal fluid dynamicsinciliarybeating. JournalofTheoreticalBiology207,415–430(2000). 52. Dillon, R. H., Fauci, L. J. & Omoto, C. Mathematical modeling of axoneme mechanics and fluid dynamics in ciliary and sperm motility. Dynamics of Continuous Discrete and ImpulsiveSystemsSeriesA10, 745–758(2003). 53. Bayly, P. & Wilson, K. Analysis of unstable modes distinguishes mathematical models of flagellarmotion. JournaloftheRoyalSocietyInterface12, 20150124(2015). 54. Goldstein,R.E.,Lauga,E.,Pesci,A.I.&Proctor,M.R.Elastohydrodynamicsynchroniza- tionofadjacentbeatingflagella. PhysicalReviewFluids 1,073201(2016). 55. Guo,H.,Fauci,L.,Shelley,M.J.&Kanso,E.Bistabilityinthesynchronizationofactuated microfilaments. JournalofFluidMechanics836, 304–323(2018). 56. Hu, T. & Bayly, P. V. Finite element models of flagella with sliding radial spokes and interdoublet links exhibit propagating waves under steady dynein loading. Cytoskeleton 75, 185–200(2018). 57. Bayly, P. V. & Wilson, K. S. Equations of interdoublet separation during flagella motion revealmechanismsofwavepropagationandinstability.BiophysicalJournal107,1756–1772 (2014). 58. Sartori,P.,Geyer,V.F.,Scholich,A.,Jülicher,F.&Howard,J.Dynamiccurvatureregulation accounts for the symmetric and asymmetric beats of Chlamydomonas flagella. eLife 5, e13258(2016). 59. Sartori, P., Geyer, V. F., Howard, J. & Jülicher, F. Curvature regulation of the ciliary beat throughaxonemaltwist.PhysicalReviewE 94, 042426(2016). 60. Michelin, S. & Lauga, E. Efficiency optimization and symmetry-breaking in a model of ciliarylocomotion. PhysicsofFluids22, 111901(2010). 61. Osterman, N. & Vilfan, A. Finding the ciliary beating pattern with optimal efficiency. ProceedingsoftheNationalAcademyofSciences108, 15727–15732(2011). 62. Elgeti, J. & Gompper, G. Emergence of metachronal waves in cilia arrays. Proceedings of theNationalAcademyofSciences110,4470–4475(2013). 63. Ding, Y., Nawroth, J. C., McFall-Ngai, M. J. & Kanso, E. Mixing and transport by ciliary carpets:anumericalstudy.JournalofFluidMechanics743,124–140(2014). 64. Guo, H., Zhu, H., Liu, R., Bonnet, M. & Veerapaneni, S. Optimal Ciliary Locomotion of AxisymmetricMicroswimmers.arXivpreprintarXiv:2103.15642927,A22(2021). 116 65. Faubel, R., Westendorf, C., Bodenschatz, E. & Eichele, G. Cilia-based flow network in the brainventricles.Science 353,176–178(2016). 66. Nawroth, J. C. et al. Motile cilia create fluid-mechanical microhabitats for the active re- cruitment of the host microbiome. Proceedings of the National Academy of Sciences 114, 9510–9516(2017). 67. Ramirez-SanJuan,G.R.etal.Multi-scalespatialheterogeneityenhancesparticleclearance inairwayciliaryarrays.NaturePhysics16, 1–7(2020). 68. Pellicciotta,N.etal.Entrainmentofmammalianmotileciliainthebrainwithhydrodynamic forces.ProceedingsoftheNationalAcademyofSciences117,8315–8325.issn:0027-8424 (2020). 69. Brumley,D.R.,Wan,K.Y.,Polin,M.&Goldstein,R.E.Flagellarsynchronizationthrough directhydrodynamicinteractions.eLife3,e02750(2014). 70. Wan, K. Y. & Goldstein, R. E. Coordinated beating of algal flagella is mediated by basal coupling.ProceedingsoftheNationalAcademyofSciences113,E2784–e2793(2016). 71. Hamilton, E. & Cicuta, P. Changes in geometrical aspects of a simple model of cilia syn- chronization control the dynamical state, a possible mechanism for switching of swimming gaitsinmicroswimmers.PLOSOne 16, e0249060(2021). 72. Man,Y.&Kanso,E.Multisynchronyinactivemicrofilaments. PhysicalReviewLetters125, 148101(2020). 73. Geyer, V. F., Jülicher, F., Howard, J. & Friedrich, B. M. Cell-body rocking is a dominant mechanism for flagellar synchronization in a swimming alga. Proceedings of the National AcademyofSciences110, 18058–18063(2013). 74. Klindt, G. S., Ruloff, C., Wagner, C. & Friedrich, B. M. In-phase and anti-phase flagellar synchronization by waveform compliance and basal coupling. New Journal of Physics 19, 113052(2017). 75. Liu, Y., Claydon, R., Polin, M. & Brumley, D. R. Transitions in synchronization states of model cilia through basal-connection coupling. Journal of the Royal Society Interface 15, 20180450(2018). 76. Guo, H., Man, Y., Wan, K. Y. & Kanso, E. Intracellular coupling modulates biflagellar synchrony. JournaloftheRoyalSocietyInterface18, 20200660(2021). 77. Meng, F., Bennett, R. R., Uchida, N. & Golestanian, R. Conditions for metachronal coor- dination in arrays of model cilia. Proceedings of the National Academy of Sciences 118 (2021). 78. Solovev,A.&Friedrich,B.M.Synchronizationinciliacarpets:multiplemetachronalwaves are stable,butonewavedominates.NewJournalofPhysics(2021). 117 79. Chakrabarti, B., Fürthauer, S. & Shelley, M. J. A multiscale biophysical model gives quan- tizedmetachronalwavesinalatticeofcilia.arXivpreprintarXiv:2108.01693(2021). 80. Bustamante-Marin,X.M.&Ostrowski,L.E.Ciliaandmucociliaryclearance.ColdSpring HarborPerspectivesinBiology9,a028241(2017). 81. Raidt, J. et al. Ciliary function and motor protein composition of human fallopian tubes. HumanReproduction30, 2871–2880(2015). 82. Yuan, S. et al. Motile cilia of the male reproductive system require miR-34/miR-449 for development and function to generate luminal turbulence. Proceedings of the National AcademyofSciences116, 3584–3593(2019). 83. Tilley,A.E.,Walters,M.S.,Shaykhiev,R.&Crystal,R.G.Ciliadysfunctioninlungdisease. AnnualReviewofPhysiology77, 379–406(2015). 84. Carter,C.S.etal.AbnormaldevelopmentofNG2+PDGFR-+neuralprogenitorcellsleads to neonatal hydrocephalus in a ciliopathy mouse model. Nature Medicine 18, 1797–1804 (2012). 85. Blyth, M. & Wellesley, D. Ectopic pregnancy in primary ciliary dyskinesia. Journal of ObstetricsandGynaecology28, 358–358(2008). 86. Gilpin, W., Bull, M. S. & Prakash, M. The multiscale physics of cilia and flagella. Nature ReviewsPhysics2,74–88(2020). 87. McKanna, J. A. Fine structure of the protonephridial system in Planaria. Zeitschrift für ZellforschungundMikroskopischeAnatomie92, 509–523(1968). 88. Rink,J.C.,Vu,H.T.-K.&Alvarado,A.S.Themaintenanceandregenerationoftheplanarian excretorysystemareregulatedbyEGFRsignaling.Development 138,3769–3780(2011). 89. Vu,H.T.-K.etal.Stemcellsandfluidflowdrivecystformationinaninvertebrateexcretory organ.eLife4,e07405(2015). 90. Vogel, S. Living in a physical world X. Pumping fluids through conduits. Journal of bio- sciences32, 207–222(2007). 91. Pellicciotta, N. et al. Cilia density and flow velocity affect alignment of motile cilia from braincells.JournalofExperimentalBiology223 (2020). 92. Pellicciotta,N.etal.Entrainmentofmammalianmotileciliainthebrainwithhydrodynamic forces.ProceedingsoftheNationalAcademyofSciences117,8315–8325.issn:0027-8424 (2020). 93. Nawroth,J.C.etal.Amicroengineeredairwaylungchipmodelskeyfeaturesofviral-induced exacerbation of asthma. American Journal of Respiratory Cell and Molecular Biology 63, 591–600(2020). 118 94. Sone,N.etal.MulticellularmodelingofciliopathybycombiningiPScellsandmicrofluidic airway-on-a-chiptechnology.ScienceTranslationalMedicine13(2021). 95. Solari, C. A., Ganguly, S., Kessler, J. O., Michod, R. E. & Goldstein, R. E. Multicellularity and the functional interdependence of motility and molecular transport. Proceedings of the NationalAcademyofSciences103,1353–1358(2006). 96. Kanso, E. A.,Lopes, R. M.,Strickler, J. R., Dabiri,J. O. & Costello,J. H. Teamwork in the viscous oceanic microscale. Proceedings of the National Academy of Sciences 118. issn: 0027-8424(2021). 97. Hildebrandt,F.,Benzing,T.&Katsanis,N.Ciliopathies.NewEnglandJournalofMedicine 364,1533–1543(2011). 98. Andrikou, C., Thiel, D., Ruiz-Santiesteban, J. A. & Hejnol, A. Active mode of excretion acrossdigestivetissuespredatestheoriginofexcretoryorgans.PLOSBiology17,e3000408 (2019). 99. Ichimura, K. & Sakai, T. Evolutionary morphology of podocytes and primary urine- producingapparatus.AnatomicalScienceInternational 92, 161–172(2017). 100. Essock-Burns, T., Bongrand, C., Goldman, W. E., Ruby, E. G. & McFall-Ngai, M. J. Inter- actions of Symbiotic Partners Drive the Development of a Complex Biogeography in the Squid-VibrioSymbiosis.mBio11(2020). 101. Dunn, C. W., Giribet, G., Edgecombe, G. D. & Hejnol, A. Animal phylogeny and its evolutionary implications. Annual Review of Ecology, Evolution, and Systematics 45, 371– 395(2014). 102. Feuda, R. et al. Improved modeling of compositional heterogeneity supports sponges as sistertoallotheranimals.CurrentBiology 27,3864–3870(2017). 103. Asadzadeh, S. S., Larsen, P. S., Riisgård, H. U. & Walther, J. H. Hydrodynamics of the leuconspongepump.JournaloftheRoyalSocietyInterface16, 20180630(2019). 104. Leys, S. P. et al. The sponge pump: the role of current induced flow in the design of the spongebodyplan.PLOSOne 6,e27787(2011). 105. Marshall,W.F.,Qin,H.,Brenni,M.R.&Rosenbaum,J.L.Flagellarlengthcontrolsystem: testing a simple model based on intraflagellar transport and turnover. Molecular Biology of theCell 16, 270–278(2005). 106. Ishikawa, H. & Marshall, W. F. Ciliogenesis: building the cell’s antenna. Nature Reviews MolecularCellBiology12, 222–234(2011). 107. Scimone, M. L., Srivastava, M., Bell, G. W. & Reddien, P. W. A regulatory program for excretorysystemregeneration inplanarians.Development 138, 4387–4398(2011). 119 108. Ruppert,E.E.&Smith,P.R.Thefunctionalorganizationoffiltrationnephridia. Biological Reviews63, 231–258(1988). 109. Audoly,B.&Pomeau,Y.Elasticityandgeometry:fromhaircurlstothenon-linearresponse ofshells(OxfordUniversityPress,2010). 110. Xu, G. et al. Flexural rigidity and shear stiffness of flagella estimated from induced bends andcounterbends.Biophysical Journal 110,2759–2768(2016). 111. Kirchhoff, G. Ueber das Gleichgewicht und die Bewegung eines unendlich dünnen elastis- chenStabes.JournalfürdiereineundangewandteMathematik 56,285–313(1859). 112. Lauga, E. & Powers, T. R. The hydrodynamics of swimming microorganisms. Reports on ProgressinPhysics72, 096601 (2009). 113. Eloy, C. & Lauga, E. Kinematics of the most efficient cilium. Physical Review Letters 109, 038101(3July2012). 114. Oriola,D.,Gadêlha,H.&Casademunt,J.Nonlinearamplitudedynamicsinflagellarbeating. RoyalSocietyOpenScience4,160698(2017). 115. Chakrabarti, B. & Saintillan, D. Spontaneous oscillations, beating patterns, and hydrody- namicsofactivemicrofilaments. PhysicalReviewFluids4,043102(2019). 116. Machin,K.Wavepropagationalongflagella. J.exp.Biol 35, 796–806(1958). 117. Jülicher, F. & Prost, J. Spontaneous oscillations of collective molecular motors. Physical ReviewLetters 78,4510(1997). 118. Edwards, B. F. L. et al. Direction of flagellum beat propagation is controlled by proxi- mal/distal outer dynein arm asymmetry. Proceedings of the National Academy of Sciences 115, E7341–e7350(2018). 119. Purcell,E.M.LifeatlowReynoldsnumber.AmericanJournalofPhysics45,3–11(1977). 120. Smith, D. J., Gaffney, E. A. & Blake, J. R. Modelling mucociliary clearance. Respiratory Physiology&Neurobiology163,178–188(2008). 121. Uchida, N. & Golestanian, R. Synchronization and collective dynamics in a carpet of mi- crofluidicrotors. PhysicalReviewLetters 104,178103(2010). 122. Uchida, N. & Golestanian, R. Generic conditions for hydrodynamic synchronization. Phys- icalReviewLetters106, 058104(2011). 123. Uchida, N. & Golestanian, R. Hydrodynamic synchronization between objects with cyclic rigidtrajectories.TheEuropeanPhysicalJournal.E,Softmatter 35, 9813–9813(2012). 124. Kuramoto,Y. Chemicaloscillations,waves,andturbulence(CourierCorporation,2003). 125. Helgason,S.&Helgason,S. TheRadonTransform(Springer,1980). 120 126. Liron, N. Fluid transport by cilia between parallel plates. Journal of Fluid Mechanics 86, 705–726(1978). 127. Taylor, G. Analysis of the swimming of microscopic organisms. Proceedings of the Royal SocietyofLondon.SeriesA209,447–461(1951). 128. Blake, J. Infinite models for ciliary propulsion. Journal of Fluid Mechanics 49, 209–222 (1971). 129. Michelin, S. & Lauga, E. Optimal feeding is optimal swimming for all Péclet numbers. PhysicsofFluids23, 101901(2011). 130. Chrispell,J.C.,Fauci,L.J.&Shelley,M.Anactuatedelasticsheetinteractingwithpassive andactivestructuresinaviscoelasticfluid. PhysicsofFluids25, 013103(2013). 131. Liron, N. & Mochon, S. The discrete-cilia approach to propulsion of ciliated micro- organisms.JournalofFluidMechanics75, 593–607(1976). 132. Gueron, S. & Liron, N. Ciliary motion modeling, and dynamic multicilia interactions. BiophysicalJournal 63, 1045(1992). 133. Blake,J.,Liron,N.&Aldis,G.Flowpatternsaroundciliatedmicroorganismsandinciliated ducts.JournalofTheoreticalBiology98, 127–141(1982). 134. Guo,H.,Nawroth,J.C.,Ding,Y.&Kanso,E.Ciliabeatingpatternsarenothydrodynamically optimal.PhysicsofFluids26, 091901(2014). 135. Guo, H. & Kanso, E. Evaluating efficiency and robustness in cilia design. Physical Review E 93, 033119(2016). 136. Stein, D. B. & Shelley, M. J. Coarse graining the dynamics of immersed and driven fiber assemblies.PhysicalReviewFluids4, 073302(2019). 137. Wuttanachamsri,K.&Schreyer,L.EffectsofCiliaMovementonFluidVelocity:IINumer- icalSolutionsOveraFixedDomain.TransportinPorousMedia 134,471–489(2020). 138. Wuttanachamsri, K. & Schreyer, L. Effects of Cilia Movement on Fluid Velocity: I Model of Fluid Flow due to a Moving Solid in a Porous Media Framework. Transport in Porous Media136,699–714(2021). 139. Bergou, M., Wardetzky, M., Robinson, S., Audoly, B. & Grinspun, E. Discrete elastic rods inACMtransactionsongraphics(TOG)27(2008),63. 140. Gazzola,M.,Dudte,L.,McCormick,A.&Mahadevan,L.Forwardandinverseproblemsin themechanicsofsoftfilaments. RoyalSocietyOpenScience5,171628(2018). 141. Ritz,W.UbereineneueMethodezurLosunggewisserVariationsproblemedermathematis- chenPhysik.JournalfurMathematik 135,s–1(1909). 121 142. MacDonald, J. Successive approximations by the Rayleigh-Ritz variation method. Physical Review43,830(1933). 143. Pearce, S. P., Heil, M., Jensen, O. E., Jones, G. W. & Prokop, A. Curvature-sensitive kinesin binding can explain microtubule ring formation and reveals chaotic dynamics in a mathematical model. Bulletin of Mathematical Biology arXiv preprint arXiv:1803.07312 (2018). 144. Laskar, A. et al. Hydrodynamic instabilities provide a generic route to spontaneous bio- mimetic oscillations in chemomechanically active filaments. Scientific Reports 3, 1964 (2013). 145. Woolley, D. M., Crockett, R. F., Groom, W. D. & Revell, S. G. A study of synchronisation betweentheflagellaofbullspermatozoa,withrelatedobservations. JournalofExperimental Biology212,2215–2223(2009). 146. Leptos, K. C. et al. Antiphase synchronization in a flagellar-dominance mutant of Chlamy- domonas.PhysicalReviewLetters 111,158101(2013). 147. Wan, K. Y., Leptos, K. C. & Goldstein, R. E. Lag, lock, sync, slip: the many “phases” of coupledflagella. JournaloftheRoyalSocietyInterface11, 20131160(2014). 148. Shingyoji, C., Higuchi, H., Yoshimura, M., Katayama, E. & Yanagida, T. Dynein arms are oscillatingforcegenerators.Nature 393,711–714(1998). 149. Guo, H. & Kanso, E. Evaluating efficiency and robustness in cilia design. Physical Review E 93, 033119(2016). 150. Wróbel, J. K., Cortez, R., Varela, D. & Fauci, L. Regularized image system for Stokes flow outsideasolidsphere.JournalofComputationalPhysics317,165–184(2016). 151. Man,Y.&Kanso,E.Multisynchronyinactivemicrofilaments. PhysicalReviewLetters125, 148101(2020). 152. Fawcett, D. W. & Porter, K. R. A study of the fine structure of ciliated epithelia. Journal of Morphology94, 221–281(1954). 153. Gittleson, S. M. Flagellar activity and Reynolds number. Transactions of the American MicroscopicalSociety,272–276(1974). 154. Baron, D. M., Kabututu, Z. P. & Hill, K. L. Stuck in reverse: loss of LC1 in Trypanosoma bruceidisruptsouterdyneinarmsandleadstoreverseflagellarbeatandbackwardmovement. JournalofCellScience120,1513–1520(2007). 155. Gadelha,C.,Wickstead,B.&Gull,K.FlagellarandciliarybeatinginTrypanosomemotility. Cell MotilityandtheCytoskeleton64, 629–643(2007). 156. Wheeler, R. J. Use of chiral cell shape to ensure highly directional swimming in Try- panosomes.PLOScomputationalbiology13,e1005353(2017). 122 157. Bottier, M., Thomas, K. A., Dutcher, S. K. & Bayly, P. V. How does cilium length affect beating?BiophysicalJournal 116,1292–1304(72019). 158. Lindemann, C. B. & Lesich, K. A. The geometric clutch at 20: stripping gears or gaining traction?Reproduction 150,R45–r53(2015). 159. Woolley, D. M. Flagellar oscillation: a commentary on proposed mechanisms. Biological Reviews85, 453–470(2010). 160. Lindemann,C.B.&Lesich,K.A.Flagellarandciliarybeating:theprovenandthepossible. JCellSci123, 519–528. issn:0021-9533(Feb.2010). 161. Kanale,A.V.,Ling,F.,Shelley,M.,Fürthauer,S.&Kanso,E.Continuumtheoryforcarpets ofmodelcilia.arXivpreprintarXiv:(2021). 162. Kirk, D. L. Volvox: a search for the molecular and genetic origins of multicellularity and cellulardifferentiation 33(CambridgeUniversityPress,2005). 163. Martinac, B., Saimi, Y. & Kung, C. Ion channels in microbes. Physiological reviews 88, 1449–1490(2008). 164. Polin, M., Tuval, I., Drescher, K., Gollub, J. P. & Goldstein, R. E. Chlamydomonas swims with two ‘gears’ in a eukaryotic version of run-and-tumble locomotion. Science 325, 487– 490(2009). 165. Cortese,D.&Wan,K.Y.Controlofhelicalnavigationbythree-dimensionalflagellarbeating. PhysicalReviewLetters126,088003(2021). 166. Wan,K.Y.Flagella:Anewkindofbeat. eLife10, e67701(2021). 167. Guo,H.,Zhu,H.&Veerapaneni,S.Simulatingcilia-drivenmixingandtransportincomplex geometries.PhysicalReviewFluids5,053103(2020). 168. Chioccioli,M.etal.Quantitativehigh-speedvideoprofilingdiscriminatesbetweenDNAH11 and HYDIN variants of primary ciliary dyskinesia. American journal of respiratory and criticalcaremedicine199, 1436–1438(2019). 169. Chioccioli,M.etal.AHow-Toguideto:Quantitativehigh-speedvideoprofilingtodiscrim- inatebetweenvariantsofprimaryciliarydyskinesia.bioRxiv, 614966(2019). 170. Blanchon, S. et al. Deep phenotyping, including quantitative ciliary beating parameters, and extensive genotyping in primary ciliary dyskinesia. Journal of Medical Genetics 57, 237–244(2020). 171. Feriani, L. et al. Assessing the collective dynamics of motile cilia in cultures of human airwaycellsbymultiscaleDDM.BiophysicalJournal 113,109–119(2017). 172. Cicuta, P. The use of biophysical approaches to understand ciliary beating. Biochemical SocietyTransactions48, 221–229(2020). 123 173. Webber, J. B. W. A bi-symmetric log transformation for wide-range data. Measurement ScienceandTechnology24, 027001(2012). 174. Sears, P. R., Yin, W.-N. & Ostrowski, L. E. Continuous mucociliary transport by primary human airway epithelial cells in vitro. American Journal of Physiology-Lung Cellular and MolecularPhysiology309,L99–l108(2015). 175. Kelly, S. J., Brodecky, V., Skuza, E. M., Berger, P. & Tatkov, S. Variability in tracheal mucociliary transport is not controlled by beating cilia in lambs in vivo during ventilation withhumidifiedandnonhumidifiedair. AmericanJournalofPhysiology-LungCellularand MolecularPhysiology320,L473–l485(2021). 176. Ruppert,E.E.Structure,ultrastructureandfunctionoftheneuralglandcomplexofAscidia interrupta(Chordata,Ascidiacea):clarificationofhypothesesregardingtheevolutionofthe vertebrateanteriorpituitary.ActaZoologica71, 135–149(1990). 177. Bassham,S.&Postlethwait,J.H.Theevolutionaryhistoryofplacodes:amoleculargenetic investigationofthelarvaceanurochordateOikopleuradioica.Development132,4259–4272 (2005). 178. Deyts,C.,Casane,D.,Vernier,P.,Bourrat,F.&Joly,J.-S.Morphologicalandgeneexpression similarities suggest that the ascidian neural gland may be osmoregulatory and homologous to vertebrate peri-ventricular organs. European Journal of Neuroscience 24, 2299–2308 (2006). 179. Holmberg,K.TheciliatedbrainductofOikopleuradioica(Tunicata,Appendicularia).Acta Zoologica63,101–109(1982). 180. Tamori, M., Matsuno, A. & Takahashi, K. Structure and function of the pore canals of the sea urchin madreporite. Philosophical Transactions of the Royal Society of London. Series B:BiologicalSciences351, 659–676(1996). 181. Shah, A. S., Ben-Shahar, Y., Moninger, T. O., Kline, J. N. & Welsh, M. J. Motile cilia of humanairwayepitheliaarechemosensory.Science 325,1131–1134(2009). 182. Bartolomaeus,T.&Ax,P.ProtonephridiaandMetanephridia-theirrelationwithintheBila- teria.JournalofZoologicalSystematicsandEvolutionaryResearch30, 21–45(1992). 183. Ott, E. et al. Pronephric tubule morphogenesis in zebrafish depends on Mnx mediated repression of irx1b within the intermediate mesoderm. Developmental Biology 411, 101– 114(2016). 184. Valverde-Islas, L. E. et al. Visualization and 3D reconstruction of flame cells of Taenia solium(cestoda).PLOSOne 6,e14754(2011). 185. Khaderi, S. N. et al. Magnetically-actuated artificial cilia for microfluidic propulsion. Lab onaChip11,2002–2010(2011). 124 186. Hanasoge, S., Hesketh, P. J. & Alexeev, A. Microfluidic pumping using artificial magnetic cilia.Microsystems&nanoengineering4,1–9(2018). 187. UlIslam,T.,Bellouard,Y.&denToonder,J.M.Highlymotilenanoscalemagneticartificial cilia.ProceedingsoftheNationalAcademyofSciences118(2021). 188. Wei,D.,Dehnavi,P.G.,Aubin-Tam,M.-E.&Tam,D.IsthezeroReynoldsnumberapprox- imationvalidforciliaryflows? PhysicalReviewLetters122,124502(2019). 189. Kanale,A.V.etal.Spontaneousphasecoordinationinciliarycarpets.arXivpreprintarXiv: (2021). 190. Kanale, A. V. et al. A simpleO() algorithm for simulations of ciliary carpets. arXiv preprintarXiv:(2021). 191. Dunn,C.W.,Leys,S.P.&Haddock,S.H.Thehiddenbiologyofspongesandctenophores. TrendsinEcology&Evolution30, 282–291(2015). 192. Leys, S. P. & Eerkes-Medrano, D. I. Feeding in a calcareous sponge: particle uptake by pseudopodia.TheBiologicalBulletin211, 157–171(2006). 193. Norekian, T. P. & Moroz, L. L. Neural System and Receptor Diversity in the ctenophore Beroeabyssicola.JournalofComparativeNeurology527, 1986–2008(2019). 194. Tamm,S.L.Ciliaandthelifeofctenophores. InvertebrateBiology133,1–46(Jan.2014). 195. Gemmill, J. F. Ciliary action in the internal cavities of the ctenophore Pleurobrachia pileus Fabr.ProceedingsoftheZoologicalSocietyofLondon88, 263–265(1918). 196. Presnell, J. S. et al. The presence of a functionally tripartite through-gut in Ctenophora has implicationsformetazoancharactertraitevolution.CurrentBiology26,2814–2820(2016). 197. Westwood, T. A. & Keaveny, E. E. Coordinated motion of active filaments on spherical surfaces.PhysicalReviewFluids6,L121101(2021). 198. Ling, F. & Kanso, E. in Bioinspired Sensing, Actuation, and Control in Underwater Soft RoboticSystems213–228(Springer,2021). 199. Jiao,Y.etal.Learningtoswiminpotentialflow. PhysicalReviewFluids6,050505(2021). 125 AppendixA 3Dfilamentelasto-hydrodynamics To solve for the filament dynamics in three-dimensions, we numerically integrate the governing equations (2.2), together with the corresponding boundary conditions. To this end, we discretize the filament’s centerline into + 1 vertices r 1 ,...,r +1 and straight edgesℓ = r +1 −r , where = 1,...,, in the spirit of [139] and [140]. The unit tangent to edge is defined as t = ℓ /ℓ whereℓ =∥ℓ ∥. Thetranslationalvelocitiesv 1 ,...v +1 areassignedtovertices. Thebodyframes {d 1, ,d 2, ,d 3, =t } andtherotationmatricesQ arenaturallyassignedtoedges. To obtain a discrete representation of the generalize curvature vector κ , we start from the definition of the associated skew-symmetric matrix κ × = Q(Q ′ ) T . The solution to this first-order differential equation is of the form Q(+Δ) = exp(−κ × Δ)Q(), where Δ is a segment of constantcurvature. Thus,writingQ = exp(−κ × ℓ )Q −1 ,wedefinethediscretecurvaturematrixas κ × =−log Q Q T −1 /ℓ ,whereℓ = 1 2 (ℓ −1 +ℓ )isa‘Voronoi’integrationdomainfromthemidpoint ofthepreviousedgetothatofthenextedge. Givenκ × ,wecanreadilyevaluateκ andthediscrete elasticmomentsM inthefilamentmaterialframe. WesolveforthediscreteelasticforcesN inthematerialframeusing(2.3). Theinextensibility constraint is enforced weakly by considering a large tensile stiffness along the filament’s centerline. Next, we transform N to the inertial frame, and substitute into the balance of forces in (2.2) using the expression for the hydrodynamic force density f ℎ from (2.4) to get the inertial 126 frame velocities v . We use standard time integrators for stiff equations (MATLAB’s ode15s) to propagate the filament position r forward in time. To close the system, we enforce the clamped boundaryconditionbyfixing r 1 =0 andt 1 =e 3 ,whileleavingthefreeendunconstrained. Reference Configuration Current Configuration FigureA.1: Geometry-preservingdiscretizationofanelasticfilament. We discretize the filament’s centerline into + 1 vertices r 1 ,...,r +1 and straight edges ℓ = r +1 −r , where = 1,...,, in the spirit of [139] and [140]; see Fig. A.1. The unit tangent to edge is defined as t =ℓ /ℓ whereℓ =∥ℓ ∥. Inextensibility is enforced weakly by considering large tensile stiffness = along the filament’s centerline. That is to say, t is allowed to experience slight extension or compression, inducing a local axial strains =(ℓ / ˆ ℓ )t −d 3 , where ˆ ℓ is a strain-free reference configuration length. The linear velocities v 1 ,...v +1 are assigned to vertices. The body frames are naturally assigned to edges, so is the rotation matrix Q; we let {d 1, ,d 2, ,d 3, =t } andQ denotethediscreterepresentationofthesequantities. Here,wedefined the rotation matrix Q using the convention x = Qx, where x is the position vector expressed in body-frame coordinates (see appendix A). Throughout this work, the subscript will be used to emphasizewhenvectorsarerepresentedinthebodyframe. Theskew-symmetricmatrixκ × =(Q/) T Qcanbeexpressedinthebodyframebyachange ofbasisasκ × =Qκ × Q T =Q(Q/) T . Thesolutiontothedifferentialequation κ × =Q(Q/) T 127 is of the formQ(+Δ)= exp(−κ × Δ)Q(), given that the curvature stays constant in the region Δ. Wecanthusdefinethediscretecurvature κ × , suchthatQ = exp(−κ × , ℓ )Q −1 ,leadingto κ × , =− log Q Q T −1 ℓ . (A.1) Here,ℓ = 1 2 (ℓ −1 +ℓ )isa‘Voronoi’integrationdomainspanningfromthemidpointoftheprevious edgetothatofthenextedge. Rememberthatthisdefinitionofdiscretecurvatureisonlydefinedat aninterior vertex= 2,...,. To obtain the bending dynamics of the filament, we first use equation (3) of the main text to solve for the internal elastic force N in the body frame of the current geometric configuration of thefilament, N , =t , × D B , κ , +A (κ , ℓ )×(B , κ , ) | {z } frametransport +()s , . (A.2) Here, the finite difference operator D and discrete average operatorA take quantities on-th and (+ 1)-th vertices as inputs and output quantities on the-th edge (= 1,···,). Sinceℓ andκ , are only defined at interior vertices, we need to pad trivial (zero) quantities for boundary vertices (= 1and+ 1). Inotherwords,the -thoutputofD[(·) ] is(·) ,+1 −(·) . Similarly,the -th outputofA[(·) ] is[(·) ,+1 +(·) ]/2. Notethatinextensibilityisenforcedweaklybysetting theaxialcomponentofN , tobe()s , withalargestiffness . 128 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -0.035 -0.03 -0.025 -0.02 -0.015 -0.01 -0.005 0 0.005 11 vertices 21 vertices 41 vertices 81 vertices Euler-Bernoulli 20 40 60 80 Number of Vertices 10 -4 10 -3 10 -2 1st order FigureA.2: CantileverdeflectioninStokes’regime. Acantileverbeamsubjectedtoatiploadof 0.1/ 3 submergedinviscousfluidsissimulatedanditsstateat = 20 4 /isshowntoconverge toEuler-Bernoullitheory. 11 21 41 81 Number of Vertices 10 -2 10 -1 Relative Error (%) 1st order Figure A.3: Euler buckling in Stokes’ regime. Euler buckling of a clamped-clamped beam is simulated in viscous fluids under vertical compressive tip loads. The critical force thresholds convergetothetheoreticalvalueinthelimitofspatialrefinement. Snapshotsoftheactualbuckled shapesare shownforatiploadof 80/ 2 . Next, we invert the force balance in equation (1) using the expression for the fluid force from equation(2)ofthemaintexttogettheinertialframevelocityv ofeachvertex v = 1 I+(− 1)A t ⊗A t D N −f , ℓ . (A.3) Here the operatorsD andA actually convert edge quantities to+ 1 vertex quantities unlike in (A.2). 129 Twist Glue FigureA.4: HelicalbucklingandMitchellinstabilityinStokes’regime. Helicalbucklingoccurs when an elastic rod is twisted at its ends to a sufficient degree. Mitchell instability describes the buckling behavior when a sufficiently twisted rod is glued at its two ends. In both case, the time snapshotsfromoursimulationprovethatwerecoverthedesiredbehavior. Quantitativeverifications areomittedforbrevity. Inthenumericalimplementation,wesolve(A.2)and(A.3)inthestrain-freereferenceconfigu- ration. Therefore we need to adjust all quantities expressed in the current geometric configuration to the strain-free reference state (we use the hat notation ˆ (·)to denote the latter). Provided that cross-sectionsretaintheircircularshapesatalltimes,weonlyneedtoinsertalocaldilationscalar = ℓ / ˆ ℓ at appropriate places in the discretized equations to ensure consistency. Specifically, due to the conservation of volume of an infinitesimal volume element, a dilation in edge length translates to a shrinkage of cross sectional area; namely, we haveℓ = ˆ ℓ ˆ which implies that = ˆ / . Consequently, the discrete bending/twisting rigidity tensor ˆ B and tensile stiffness ˆ in the stain-free reference configuration are related to B and in the current configuration via B = ˆ B / 2 and = ˆ / . All length quantities need to include a factor of when converted to the strain-free reference configuration. In addition, because we used the integration domain ℓ in the momentum equation (A.2), we need to average all quantities that vary with arc-length over the 130 voronoi regionℓ . For example, the bending rigidity B and dilation factor are averaged overℓ using ˆ B = ˆ B ˆ ℓ + ˆ B −1 ˆ ℓ −1 2 ˆ ℓ , E = ˆ ℓ + −1 ˆ ℓ −1 2 ˆ ℓ . Whenadjustingthegeneralizedcurvaturewithrespecttothedilationfactor,wewriteκ = ˆ κ /E . We use standard time integrator for stiff equations (MATLAB’s ode15s) to solve (A.3) and propagate the filament position r forward in time. To closed the system, we enforce the clamped boundaryconditionbyfixing r 1 =0 andt 1 =e 3 ,whileleavingthefreeendunconstrained. We test our numerical method on the deflection and buckling behavior of a cantilever beam submerged in viscous fluid. Since the steady state shape of a cantilever beam under transverse tip load would remain the same in the presence of fluid drag, it is sound to compared our simulation results with the linear Euler-Bernoulli theory at a large (in the absence of local filament shear). Indeedweobservefirstorderconvergenceofourdeflectedfilamentat = 20 4 / asweincrease spatial resolution; See Fig. A.2. Similarly, since the onset of Euler buckling of a cantilever beam underaxialcompressionatoneendshouldalsoremainequalregardlessofthesurroundingmedium, we compare the threshold of buckling event with the linear theory result of = 4 2 (/ 2 ). Again we observe first order convergence in the relative error in the limit of spatial refinement. Moreover,ournumericalsimulationalsogivethenonlinear buckledshapeofthefilament. To conclude, we note that the actual algorithm developed is slightly more general in that we can allow twist in the filament following the methods in [139]. In the presence of properly-chosen applied twists, the algorithm correctly exhibit the localized helical buckling instability as well as theMitchellinstability[200–204];SeeFig.A.4. 131 References 200. Coyne,J.Analysisoftheformationandeliminationofloopsintwistedcable.IEEEJournal ofOceanicEngineering15, 72–83(1990). 201. VanderHeijden,G.&Thompson,J.Helicalandlocalisedbucklingintwistedrods:aunified analysisofthesymmetriccase.Nonlineardynamics21, 71–99(2000). 202. Neukirch, S., Van Der Heijden, G. & Thompson, J. Writhing instabilities of twisted rods: frominfinitetofinitelength. JournaloftheMechanicsandPhysicsofSolids50,1175–1191 (2002). 203. Van der Heijden, G., Neukirch, S., Goss, V. & Thompson, J. Instability and self-contact phenomena in the writhing of clamped rods. International Journal of Mechanical Sciences 45,161–196(2003). 204. Goriely, A. Twisted elastic rings and the rediscoveries of Michell’s instability. Journal of Elasticity84, 281–299(2006). 132 AppendixB Two-state,two-filament2Dciliummodel Following [115], we discretize the arc length uniformly with + 1 points, with mesh size Δ= 1/. Weusethesecondordercentraldifferenceformulaforthespatialderivatives, (·) , ≃ 1 2Δ [−(·) −1 +(·) +1 ] (·) , ≃ 1 Δ 2 [(·) −1 − 2(·) +(·) +1 ]. (B.1) Weusesecond-orderforwardandbackwarddifferenceformulatoenforcethebaseandtipboundary conditions,respectively. Toget nonlinearsolutionsfor,wefirstsubstituteEq.(2.26)and(2.28)into(2.27)andobtain S 4 [ ( + − − )−+ +]= ( + + − )[ − 2 2 + 3 + ]. (B.2) Thenateachtimestep,wefirstcalculatethenormalforce () andtension () usingEq.(2.25) and (B.2). We then enforce boundary conditions for () and () using (−1) . Next, we employ asecondorderbackwardschemetoEq.(2.26)andget S 4 2Δ h 3 (+1) − 4 () + (−1) i = () − 2 () ( () ) 2 + 3 () () + () () . (B.3) 133 We enforce the boundary conditions for () at this step. Finally, we use implicit Euler scheme on Eq.(2.29)toget 1 Δ h (+1) ± − () ± i = h 1− (+1) ± i −(1−) (+1) ± exp h ∗ 1∓ () i , (B.4) where () =Δ −1 [ (+1) − () ]. t = 0 steady orbit s 0 0.02 0.04 0.06 0.08 0.1 0 0.2 0.6 0.4 0.8 1 n + t = 0 steady orbit s 0 0.02 0.04 0.06 0.08 0.1 0 0.2 0.6 0.4 0.8 1 n + Figure B.1: Robustness to initial conditions. Initial valudes (red) for r and final steady orbit (grayscale) of r and + before = 25 are shown for two different initial conditions. Initial conditions for motor states are always chosen to be constant with respect to (not shown). In the toprow,westartfromanalmoststraightconfigurationexceptthatthemidpointisperturbedslightly, ( = 0.5, = 0) = 0.001. In the second row, the initial condition for is set to be 0.2 sin(3). It is clear that resultant steady orbit for both and + are identical, except that due to the hinged boundarycondition,themeanangularpositionofthefinalperiodicsolutionisdeterminedbythatof theinitialcondition. Inbothcases,themotorstate + exhibitidenticalsharpwavefrontpropagation from base to tip. Here we used Sp= 5 and = 3500, identical to second row of Fig. 4.3 in the maintext. WealsousedameshsizeΔ= 0.01andanintegrationtimestepofΔ= 3.125× 10 −4 . 134 Whensolvingforthehingedboundarycondition,wereplacethefirstformulainEq.(2.30)with | =0 − ∫ 1 0 ( ′ ,) ′ = 0. (B.5) SubstitutingEq.(2.27)into(B.5),weget | =0 − ∫ 1 0 ( ′ ,) ′ = | =0 + ∫ 1 0 ( +) ′ = | =0 − | =0 + ∫ 1 0 , = ∫ 1 0 , (B.6) whereweusedtheboundarycondition | =1 = 0. Wearriveat ∫ 1 0 = 0. (B.7) This constraint needs to be enforced when solving for the normal and tensile components of the forceviaEq.(2.25)and(B.2). Thediscretizedversionofthisintegralconstraintis 1 2 1 + ∑︁ =2 + 1 2 +1 = 0. (B.8) Allotherequationsfollowthesameformasintheclampedcase. Note that in order to determine the average flagella wave traveling direction for nonlinear simulations,allweneedistotakethetangentialcomponentofhydrodynamicdragforce f ℎ · ˆ =(− − ), (B.9) 135 compute its average over one period and the entire arclength (⟨f ℎ · ˆ ⟩), then check its sign. This works since by total force balance, it must equal to the total reaction force experienced by the clampedjointinthetangentialdirection. t = 0 steady orbit s 0 0.02 0.04 0.06 0.08 0.1 0 0.2 0.6 0.4 0.8 1 n + ∆t Relative Error 10 −3 10 −2 10 −1 10 0 10 −5 10 −4 10 −3 Figure B.2: Convergence of numerical sharp wave front solution. The almost straight initial condition (red) forr and final steady orbit (grayscale) of r and + are shown next to each other for Δ = 10 −4 . The relative error (L2 norm) of the final solution at = 25 at variousΔ and fixed mesh sizeΔ= 0.01 is computed againstΔ= 2× 10 −5 and shown on the right. Black crosses are relative error of position r, whereas red circles and blue squares represent that of motor state + and − ,respectively. Itisevidentthatthenumericalmethodisconvergingtoaconsistentsolution, whereasharpwavefrontpropagationfrombasetotipfor ± . HereweusedSp= 5and = 3500, identicaltosecondrowofFig.4.3inthemaintext. We have verified that our results for the clamped boundary conditions are similar to those reported in [114, 115]. Here, we present our validation results for the hinged boundary condition only. We choose the case shown in Fig. 4.3 of the main text as our test case. In Fig. B.1, we use two distinct initial conditions and show that the filament behavior is robust to variations in initial conditions. In Fig. B.2, we vary the integration time step and show that the filament behavior and molecular motor dynamics are robust to the integration time step. In all simulations in this work, wechooseameshsizeΔ= 0.01andanintegrationtimestepofΔ= 3.125× 10 −4 . 136 AppendixC RotorhydrodynamicsusingregularizedStokeslets We propose an accessible algorithm withO() complexity to simulate ciliary dynamics. The algorithm relies solely on familiarity with the singularity solutions to the Stokes equations, and exploits key geometric and dynamical features of the ciliary system. Specifically, we consider a quasi-periodicciliarysystem,whichwedefineheretoconsistofafundamentalrectangulardomain containing a 2D array of cilia, and eight image copies forming the immediate neighbors to the fundamental domain; see Fig. C.1. This simplification is consistent with ciliary systems that are dense but never infinite; mathematically, it eliminates the need for a specialized treatment of the doubly-periodic image system without reducing theO( 2 ) complexity of the problem. To further reduce the asymptotic computational complexity toO(), we exploit the fast decay of hydrodynamicinteractionnearano-slipwallandchoosetoignorepairwiseinteractionslargerthan a pre-determined radius. In effect, our approach can be seen as performing only the ‘easy’ part of a fast multipole method and ignoring the far-field interaction effects completely. Nevertheless, for our problem of interest, we demonstrate that this approximation would not introduce catastrophic errorwhendonecarefully;seeFig.C.2. 137 Figure C.1: Schematic of the quasi-periodic computational domain. A. The computational domain, consisting of the fundamental domain (red), images (gray) and interaction stencil (blue). Thezoomedinviewoftheinteractionstencilshowsthefocalciliumatthecenterofthelattice(dark blue)andtheneighboringcilia(red). B.Theciliainsidethefundamentaldomainarearrangedona squarelattice,withlatticespacing,atadistanceℎaboveano-slipboundary. Adaptedfrom[190]. Tocomputethepair-wiseinteractionofeachcilia/Stokeslet,weusethefollowingsemi-infinite domainregularizedStokesletformulationbasedonimagesystems[205]. Thisinvolvescalculating thefollowingregularizedsingularities: u Stokeslet =f 1 ()+(f·(x−x ))(x−x ) 2 (), u potentialdipole =q 1 ()+(q·(x−x ))(x−x ) 2 (), u Stokesdoublet =(g·b)(x−x ) 2 ()+(g·(x−x ))b 2 ()+ (b·(x−x ))g ′ 1 ()/+(g·(x−x )) (b·(x−x ))(x−x ) ′ 2 ()/, u rotlet =L×∇ (∥x−x ∥). (C.1) where 1 () = ′ /− , 2 () =( ′′ − ′ )/ 3 , 1 () = ′ /− , 2 () =( ′′ − ′ )/ 3 with being the chosen blob function and , , the corresponding smoothed Green’s and 138 Figure C.2: Convergence of quasi-periodic fluid velocity . A. Cosine fields of a system of 151× 151 cilia at heightℎ= 0.5 above a plane wall, in the isotropic (top row) and synthetic wave (bottom row) states. B. Flowfields ( v 9 , where subscript indicates the number of shell is used for pairwise interactions) in the (/ℎ = 1) and ( = 1) planes, computed inside a small 5× 5 patch at the center, marked by the orange box in A. C. Flow velocity (solid black) induced at the center of the lattice approaches the result obtained using the Fast Multipole Method [206] (dashed gray) for increasing shells of interacting neighbors. D. The relative error in the induced flow computed ( =∥v ∗ −v ∥/∥v ∗ ∥), using v ∗ = v 2 (solid black), and v ∗ = v 9 (solid gray), v ∗ =v 2P (dashedgray). Thebluedottedlinesshowthenumberofshellsneededforarelativeerror of 1%. Adaptedfrom[190]. biharmonic solution. Andf,q,g,L being the respective strengths of the singularities andb the directionofthedoublet. 139 Assume, without loss of generality, the no-slip wall is aligned to e 1 direction at distance away from origin. Then we can find induced velocity at location y due to a Stokeslet at x =(+ℎ , , ) as u(y)=f 1 (∥y ∥)+(f·y 2 (∥y ∥)−f 1 (∥y ∥)+(f·y )y 2 (∥y ∥) −ℎ 2 (g 1 (∥y ∥)+(g·y )y 2 (∥y ∥)) + 2ℎ( ′ 1 (∥y ∥)/∥y ∥+ 2 (∥y ∥))(L×y ) + 2ℎ((g·e 1 )y 2 (∥y ∥)+(y ·e 1 )g 2 (∥y ∥) +(g·y )e 1 ′ 1 (∥y ∥)/∥y ∥+(y ·e 1 )(g·y )y ′ 2 (∥ ∥)/∥y ∥. (C.2) wherey =y−x andy =y−(−ℎ , , ),withthedipolestrengthfixedat g= 2(f·e 1 )e 1 −f androtletstrengthatL=f×e 1 . Forconvenience,wewillalsoshowonechoiceoftheregularizedfunctionsas 1 ()= 1 8 1 ( 2 + 2 ) 1/2 + 2 ( 2 + 2 ) 1/2 , 2 ()= 1 8 1 ( 2 + 2 ) 3/2 , 1 ()= 1 4 1 ( 2 + 2 ) 3/2 − 3 2 ( 2 + 2 ) 5/2 , 2 ()=− 3 4 1 ( 2 + 2 ) 5/2 . (C.3) Finally,usingabovemethodforinducedvelocitycalculationoftheBlaketensor,onecannaively integrate (2.38) by (i) finding the pairwise induced velocities inside the interacting domain, (ii) obtainafiniteincrementinrotorphaseangleaccordingtotheequation,(iii)integratewithstandard methodandachievenumericalstabilitybykeepingasufficientlysmalldiscretetimestep. However, 140 thisnumericalintegrationstepcanalso beachievedby(i)directlyupdaterotor/Stokesletlocations due to pairwise induced velocities with any standard numerical integration scheme, disregarding the trajectory constraint, (ii) enforce constraint forces by projecting updated rotors back onto the closestpointalongthetrajectory,(iii)calculatethecorrespondingnewphaseangle. Inpractice,we foundthatwhilebothmethodconvergestothesameanswerinthelimitofinfinitesimaltimesteps, thesecondprojection-basedmethodcanyieldmuchmorenumericalstabilityevenwhentimestep values are chosen to be very large (e.g.,Δ > 1) due to the nonlinear correction! Thus, given an arbitrary (isotropic) initial state of the rotors, the second method can be first used to speed up the process by quickly bring the system into its steady state wave-like behavior before shrinking time stepsformoreaccuratecomputationalresults. References 205. Ainley, J., Durkin, S., Embid, R., Boindala, P. & Cortez, R. The method of images for regularizedStokeslets.JournalofComputationalPhysics227, 4600–4616(2008). 206. Yan, W. & Shelley, M. Universal image systems for non-periodic and periodic Stokes flows aboveano-slipwall.JournalofComputationalPhysics375,263–270(2018). 141 AppendixD Activeporousmedia drivenchannelflow D.1 CFDsetup We used a mixed finite element method to solve the non-dimensional Brinkman equations by introducing vorticity as an independent variable to avoid numerical evaluation of Laplacian, and solve the fluid equations in weak form on a rectangular mesh. The mixed formulation requires the assignment of discrete vorticity to nodes of the mesh, velocities to edges, and pressure to faces. Correspondingly, bilinear, linear, and piecewise constant basis functions are used for vorticity, velocity and pressure inside each element, respectively; see [207–210]. These compatible discrete basisfunctionsallowtheincompressibilityconstraint/continuityequationtobesolvedsimultane- ously with the momentum balance equation. Note that in this formulation, the tangential velocity boundary conditions need to be enforced weakly through the vorticity equation, while the normal velocitiescanbeenforcedstronglyviadirectsubstitution. We verified the convergence of this method by first solving two relevant analytically tractable problems: flow inside a periodic channel driven by prescribed boundary flow (akin to the limit case of ℎ→ 0 and ¯ →∞) and pressure-driven Poiseuille flow (see §2.4.5). In both cases, we tested with constant non-trivial Brinkman coefficients. To ensure pointwise error remain less than 142 a b L 2 Absolute Error Grid Size 10 -3 10 -4 10 -1 10 -2 10 -7 10 -5 10 -6 10 3 10 1 10 2 10 0 10 -20 10 -10 10 -15 10 -5 10 3 10 1 10 2 Grid Size u y ω u x u y ω u x Figure D.1: Validation of the Brinkman solver. a. 2 error convergence of vorticity (black) and velocity (red and blue) for the slip boundary problem, where (, 0)= (,)= cos(2). Testedonunitsquaredomainusinggridsizeof20x20to200x200with= 1,== 1, = 1000. b. Convergence plot for Poiseuille flow with = 1, = = 1, = 16,Δ=−1. The error in isnoisyduetofloatingpointprecision(noverticalflow). (10 −6 ),theductsimulationsusemeshsizesnolargerthan 5× 10 −3 perside(e.g.,200x200for = and200x400for= 2). As an aside, since our activity term is prescribed and the pressure gradient term has only linear effects on the final flow field (Fig. 6.3a), we can obtain the flow field at any arbitrary Δ by computing a superposition of a zero adverse pressure flow field and a finite adverse pressure flow fieldwiththenecessaryweighting, i.e.,u(Δ)=u(Δ= 0)+Δ/Δ [u(Δ )−u(Δ= 0)]. D.2 Optimizationprocedures Ourgoalistofindoptimalductgeometries (,ℎ/) forfluidpumpingunderanimposedadverse pressure gradient Δ, given a material constraint ¯ ℎ. To optimize duct designs, we use the pumpingefficiency astheobjectivefunction. 143 To identify optimal duct designs over the morphospace(,ℎ/), we devise an optimization algorithm detailed in a pseudo-code format in Algorithm 1. The algorithm can be explained as follows. For given values of material constraint ¯ ℎ and adverse pressure gradient Δ, we pre-compute the flow speed over the entire morphospace(,ℎ/) and, at each(,ℎ/), we calculate theflow rate = and poweroutput, here proportionaltoJ= 3 . We thenidentify the duct design(,ℎ/) † that maximizes the output power and we record the corresponding valueJ † . Note that, by construction, this duct design also maximizes efficiency for this value of material constraint ¯ ℎ because is proportional toE = J/ ¯ ℎ and ¯ ℎ is constant. In addition to identifying the optimal duct design(,ℎ/) † , we assess the performance of all ducts, over the entiremorphospace,relativetothisoptimaldesign. Thisisdonebyassigningaperformancescore S, defined as the pumping efficiency scaled by the maximal pumping efficiency, at each point of the morphospace for the given choice of material and functional constraints( ¯ ℎ,Δ), ignoring designsthatcauseflowreversal. We repeat this process by sweeping over a range of material and functional constraints ( ¯ ℎ∈ {1, 5,..., 3000}) and (Δ ∈ {0, 10 −5 ,..., 10 7 }), returning optimal duct designs(,ℎ/) † as function of( ¯ ℎ,Δ). We also record, for the range of material constraints ¯ ℎ considered, the largestpossiblevaluesoftheperformancescoreS ∗ ,efficiency ∗ ,flowrate ∗ ,andadversepressure Δ ∗ overtheentiremorphospace(,ℎ/). D.3 Optimalductdesigns Before running Algorithm 1 for the full range of material constraints identified in Fig. 2.5B over the entire morphospace(,ℎ/), it is instructive to first examine the performance of the ciliary 144 Algorithm1Evaluatingductdesignperformance Require: listofCDmorphologicalandconstraintparameterswithCFDoutput: (,ℎ/) // lumendiametersandcilia-to-lumenratios ¯ ℎ // ciliarymaterialconstant Δ // imposedadversepressuregradient ≡(,ℎ/) // pre-computedaverageflowvelocity,storedasamatrix overtheentiremorphospaceateachcombinationof( ¯ ℎ,Δ) 1: Initializethefollowingasnullmatricesovertheentiremorphospace(,ℎ/): S ∗ ≡S ∗ (,ℎ/) // maximumperformance score E ∗ ≡E ∗ (,ℎ/) // maximum”efficiency” / ∗ ≡ ∗ (,ℎ/) // maximumflow rate Δ ∗ ≡Δ ∗ (,ℎ/) // largestΔ withoutflowreversal 2: for ¯ ℎ= 1, 5, 10,..., 3000do 3: forΔ/ −1 = 0, 10 −5 ,..., 10 7 ,... do 4: if max() < 0then 5: break //Stopifnodesignparameterscanpumpatthis( ¯ ℎ,Δ) 6: At( ¯ ℎ,Δ),initializematricesoverthemorphospace(,ℎ/): J≡J(,ℎ/);P≡P(,ℎ/);≡... ;E≡... ;S≡... J= 3 // calculate“outputpower” P= repmat(Δ/ −1 ,size(,ℎ/)) //replicateΔ = // calculateflowrate E=J/( ¯ ℎ) // calculate“efficiency” / S=0 // initializeperformancescore 7: Identifyoptimaldesignsand assignperformancescoresatgiven( ¯ ℎ,Δ): [J † ,(,ℎ/) † ]= max(J) //findparametersthatmaximize J S=J/J † //calculaterelativeperformancescore return(,ℎ/) † //Returnoptimaldesignparameters asfunctionof( ¯ ℎ,Δ). 8: Ignoredesignsthatcauseflow reversalatgiven ( ¯ ℎ,Δ): S( < 0)= 0 //assigndesignswithflowreversalanullscore P( < 0)= 0 9: Returntheextremaofthefollowingoverthemorphospace(,ℎ/): S ∗ = max(S,S ∗ ) E ∗ = max(E,E ∗ ) ∗ = max(, ∗ ) Δ ∗ = max(P,Δ ∗ ) pumps along two cross-sections of the morphospace, at ℎ/ = 0.5 and = 50 m, for a single value of the material constraint ¯ ℎ = 1000; see Fig. D.2A. We consider three values of adverse 145 FigureD.2: Performancescoreatdifferentadversepressuregradient. A. Theperformancescore of two representative slices ofℎ/ = 0.5 and = 50m (doted lines). B-D. Performance score for three representative adverse pressure gradient values. Slices whereℎ/ = 0.5 and = 50m are highlighted in black with the optimal score noted in circle when possible. Performance score for other values ofℎ/ (in increment of 0.05) and larger than 50m (up to 200m) are shown ingray linesinthe background. B Whenadverse pressuregradientis verylarge,only highℎ/ at relativelysmall cangetforwardflowforpositiveperformancescores. CWhenadversepressure gradient is lower, optimal duct diameter occurs at low, and optimal cilia-to-lumen ratio is high ℎ/. DInabsenceofanyadversepressuregradient,theoptimaldiameterisatthetestedmaximum andoptimalratioatthetestedminimum. ¯ ℎ= 1000forallpanels. pressureΔ and in each case, we identify the optimal design(,ℎ/) † , and compute the relative performance scoreS= 3 /( 3 ) † as a function of forℎ/= 0.5 (Fig.D.2B-D, left column) and as a function ofℎ/ for = 50m (Fig.D.2B-D, right column). We find that as the adverse pressure gradient increases, S increases in the direction of lower lumen diameter and higher cilia-to-lumen ratio ℎ/ as shown by the highlighted circle and solid black lines inside left and 146 Figure D.3: Optimal designs that maximize pumping efficiency. At each value of the material constraint ¯ ℎ,weidentifyafamilyofoptimalductdesigns(,ℎ/) † (whiteline)thatmaximize pumping efficiency. We superimpose as a 2D colormap values of the maximal performance score S ∗ andmaximalsustrainedpressureΔ ∗ . BrightnesscorrespondstoS ∗ andrepresentsthedegreeof optimality,whilesaturationtoredrepresentsthemaximaladversepressuregradientΔ ∗ ;colormap values are given in the far right. The dashed zone indicates locations in the morphospace where weusednumericalextrapolationinFig.4ofthemaintext. right column of Fig.D.2B-D respectively. Moreover, when adverse pressure become sufficiently large, designs with smaller cilia-to-lumen ratio can be completely absent since they could only producebackflowinthedirectionofapplied Δ,asshownbytheabsenceofthelinesthatrepresent designswithℎ/ < 0.7inFig.D.2B. Next, we run Algorithm 1 for one value of material constraint ¯ ℎ at a time. For each ¯ ℎ, we plot in Fig. D.3 the family of optimal designs(,ℎ/) † for the full range of adverse pressure gradients Δ that we consider (white line). We superimpose on the morphospace the optimal performance score S ∗ and corresponding maximal sustained pressureΔ ∗ using a 2D colormap with values of S ∗ as brightness and Δ ∗ as saturation towards pure red hue. It is evident that as the ciliary material constraint ¯ ℎ decreases, the (white) line of optimal designs shift towards 147 Figure D.4: Pumping efficiency. A. At very large adverse pressure gradient, positive efficiency valuesareachievedonlybydesignswithsmalllumendiameter,largecilia-to-lumenratio,andlarge material constraints. B. At intermediate adverse pressure gradient value, small lumen diameter, optimalpumpingefficiencyisachievedbydesignswithlargecilia-to-lumenratioandintermediate materialconstraints. C.Atzeroadversepressuregradient,largelumendiameter,optimalpumping efficiency is achieved by designs with small cilia-to-lumen ratio and relatively small material constraints. the bottom-left (smaller and more open channels), and lower values of material constraint also reduce the range of adverse pressure gradient that ciliated ducts can pump against. Note that our numerical results are limited to lumen diameter values above a minimal threshold to avoid 148 Figure D.5: Optimal performance score and pumping efficiency over the entire mor- phospace. A. Cumulative performance score S ∗ for all test material constraints ¯ ℎ = {1, 5, 10,..., 1000, 3000} overlayed on top of each other. The < 50% zones of this figure is usedasocclusionmasksinpanelBandFig.6.6B,C. numericalartifactsresultingfromextremelysmallaspectratiosbetweenthelumendiameter and theprescribedwavelength. Resultsoftheaboveoptimizationresultscanbedirectlytranslatedintohowpumpingefficiency / changes as a function of(,ℎ/), material constraint ¯ ℎ, and adverse pressureΔ. As shown in Fig.D.4, two trends are clear: (i) just like the case for performance score when material constraints are fixed to a single value, as adverse pressure gradient Δ increases (moving from C toA),lowerlumendiameter (leftcolumn)andhighercilia-to-lumenratioℎ/ (rightwardonx- axis)becomesmoreoptimalbyachievingahigherefficiency;(ii)furthermore,as Δincreases,we also see that higher material constraint constants ¯ ℎ become necessary to produce more efficient pumping. Combining results of Algorithm 1 for the full range of values of ¯ ℎ, we can obtain on the morphospace the cumulative performance score S ∗ (maximum possible score at given(,ℎ/) 149 coordinate)shownasthegraycolormapinFig.D.5A.Notethathereweusenumericalextrapolation for small values of corresponding to the region of the parameter space highlighted by dashed lines in Fig. D.3. The white lines of Fig. D.3 are overlaid next to each other as solid lines. In Fig. D.5B, we show the corresponding cumulative pumping efficiency ★ / . Here we used the 50% isoline of the cumulative performance scoreS ∗ in Fig. D.5A as a clipping mask for clarity. This same mask is also applied to Fig. 6.6B,C for the analysis ofΔ ∗ and ∗ in the main text, respectively. References 207. Bochev, P. B. & Hyman, J. M. in Compatible Spatial Discretizations 89–119 (Springer, 2006). 208. Bochev, P. & Gunzburger, M. A locally conservative mimetic least-squares finite element methodfortheStokesequationsinInternationalConferenceonLarge-ScaleScientificCom- puting(2009),637–644. 209. Kreeft, J. & Gerritsma, M. Mixed mimetic spectral element method for Stokes flow: A pointwisedivergence-freesolution.JournalofComputationalPhysics240,284–309(2013). 210. Hiemstra,R.,Toshniwal,D.,Huijsmans,R.&Gerritsma,M.I.Highordergeometricmethods with exact conservation properties. Journal of Computational Physics 257, 1444–1471 (2014). 150
Abstract (if available)
Abstract
Cilia are micron-size, hair-like protrusions present in almost all Eukaryotic cells. Their motile variety - also called Eukaryotic flagellum when appearing in small numbers - enables cells to direct fluid flow around themselves by moving and coordinating in wave-like fashions.
The origin of cilia activity arises from a complex but relatively conserved nanoscale structure composed of many semi-rigid microfilaments (microtubule doublets) that slide due to thousands of molecular motor proteins (axonemal dyneins). Notably, cilia can work individually as propellers for single-celled flagellates or animal spermatozoa, move together in one or a few pairs for algae Chlamydomonas in some ways analogous to appendages of mammals, or show up in large groups on microorganisms like Paramecium and Stentors or on epithelial tissue of mammalian brain ventricle, airway, and reproductive tracts (1 mm^2 of airway cilia tissue can have up to 1 million cilia!).
The physics behind ciliary motion have long fascinated biologists and fluid dynamicists since their discovery more than 300 years ago. Nevertheless, there remain debates on the precise mechanisms responsible for the spontaneous oscillation of individual cilium and the long-range coordination of many cilia. Perhaps more importantly, our understanding of cilia physics are often segregated by the cascading length scales inherent for any ciliated tissue. There is thus a pressing need to develop multiscale analyses that can incorporate key insights at the nano- and micro-scales into functional predictions at the tissue- and organ-level.
In this thesis, three such contributions will be expounded in detail:
(i) Using two-state molecular motor models and the continuous filament representation for the elasto-hydrodynamics of cilia, we demonstrate how boundary conditions and distribution of activities along a cilium can change the geometry and wave traveling direction of spontaneously emerging cilia oscillations.
(ii) Reducing individual cilia motion into phased oscillators, we theoretically and numerically quantify how key attributes of cilium beat (effective and recovery strokes, planarity of motion) combined with structural heterogeneity parameters at the tissue-level alter the quality and form of the emergent metachronal waves and associated flow functions.
(iii) Coarse-graining cilia activity as active porous-media, we illustrate how the morphology of a ciliated duct predicts its relative ability to perform bulk transport versus filtration, and establish universal design rules for ciliary pumps.
Frustrating as it may be, this work marks only the exciting beginning of model ciliated systems from the nanometer scale all the way to macroscopic functions. Combined with further advances in imaging techniques at different length and time scales, it is the author's hope that the numerical and theoretical techniques developed here will serve as groundwork for the resolution of the remaining mysteries about cilia, and perhaps even the development of diagnostic tools for cilia-related diseases in the near future.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Synchronization in carpets of model cilia: numerical simulations and continuum theory
PDF
Model based design of porous and patterned surfaces for passive turbulence control
PDF
Hydrodynamic stability of two fluid columns of different densities and viscosities subject to gravity
PDF
A million-plus neuron model of the hippocampal dentate gyrus: role of topography, inhibitory interneurons, and excitatory associational circuitry in determining spatio-temporal dynamics of granul...
PDF
Transitions in the collective behavior of microswimmers confined to microfluidic chips
PDF
Tracking and evolution of compressible turbulent flow structures
PDF
Pattern generation in stratified wakes
PDF
The projection immersed boundary method for compressible flow and its application in fluid-structure interaction simulations of parachute systems
PDF
Supersonic flow study of a capsule geometry using large-eddy simulation
PDF
Development and applications of a body-force propulsor model for high-fidelity CFD
PDF
Numerical study of shock-wave/turbulent boundary layer interactions on flexible and rigid panels with wall-modeled large-eddy simulations
PDF
Switching dynamical systems with Poisson and multiscale observations with applications to neural population activity
PDF
Neural network integration of multiscale and multimodal cell imaging using semantic parts
PDF
The physics of emergent membrane phenomena
PDF
Efficient inverse analysis with dynamic and stochastic reductions for large-scale models of multi-component systems
PDF
Reconstruction and estimation of shear flows using physics-based and data-driven models
PDF
Functional real-time MRI of the upper airway
PDF
Information geometry of annealing paths for inference and estimation
PDF
Modeling of bicontinuous nanoporous materials and investigation of their mechanical and hydraulic properties
PDF
Microbial interaction networks: from single cells to collective behavior
Asset Metadata
Creator
Ling, Feng
(author)
Core Title
Multiscale modeling of cilia mechanics and function
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Mechanical Engineering
Publication Date
09/28/2022
Defense Date
02/18/2022
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
active matter,axoneme,bifurcation,biological physics,cilia,cilia coordination,cilia oscillation,ciliary pumps,ciliopathies,coarse graining,collective motion,continuum theory,dynein,emergent physics,filament dynamics,flagella,fluid dynamics,fluid structure interaction,instability,microfluidic pump design,molecular motor,mucociliary clearance,multiscale modeling,OAI-PMH Harvest,phased oscillator,porous media,Stokes' flow,synchronization
Format
application/pdf
(imt)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Kanso, Eva Adnan (
committee chair
), Bermejo-Moreno, Ivan (
committee member
), Haselwandter, Christoph (
committee member
), Newton, Paul K. (
committee member
), Oberai, Assad (
committee member
)
Creator Email
fling@usc.edu,levincoolxyz@gmail.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-oUC110843497
Unique identifier
UC110843497
Document Type
Dissertation
Format
application/pdf (imt)
Rights
Ling, Feng
Type
texts
Source
20220331-usctheses-batch-918
(batch),
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
active matter
axoneme
bifurcation
biological physics
cilia
cilia coordination
cilia oscillation
ciliary pumps
ciliopathies
coarse graining
collective motion
continuum theory
dynein
emergent physics
filament dynamics
flagella
fluid dynamics
fluid structure interaction
instability
microfluidic pump design
molecular motor
mucociliary clearance
multiscale modeling
phased oscillator
porous media
Stokes' flow
synchronization