Close
The page header's logo
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
/
Flow physics and nutrients transport in microscopic living systems
(USC Thesis Other) 

Flow physics and nutrients transport in microscopic living systems

doctype icon
play button
PDF
 Download
 Share
 Open document
 Flip pages
 More
 Download a page range
 Download transcript
Copy asset link
Request this asset
Transcript (if available)
Content Flow Physics and Nutrients Transport in Microscopic Living Systems
by
Jingyi Liu
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(Mechanical Engineering)
May 2025
Copyright 2025 Jingyi Liu



Dedication
To open access and the pursuit of knowledge and discovery.
ii



Acknowledgements
I would like to express my deepest gratitude to my advisor, Professor Eva Kanso, for her invaluable guidance in conducting scientific research, approaching complex questions, and encouraging me to explore the
beauty of mathematics. Her insightful mentorship and discussions have been instrumental in shaping this
research. I am also deeply grateful to John Costello for his intellectual insights and inspiring discussions
on the beauty of biology and natural sciences.
I extend my sincere appreciation to my committee members: Professor Mitul Luhar, for his thoughtful
feedback and support from many aspects throughout this journey; Professor Ivan Bermejo-Moreno, whose
inspiring words encouraged me to pursue research in fluid dynamics — his class was my first introduction
to the subject; and Professor Christoph Haselwandter, who provided me with a deeper understanding of
physics from a theoretical perspective.
I am grateful to Dr. Hanliang Guo for his generous advice in the early stages of my research and for
many encouraging discussions along the way, and to Dr. Yi Man for mathematical insights that enriched
this thesis. Special thanks to my lab colleagues, whose camaraderie and intellectual exchanges made this
journey both enjoyable and stimulating. I am especially grateful to my peers Dr. Feng Ling, Dr. Yusheng
Jiao, Dr. Chenchen Huang, and Dr. Haotian Hang for their supportive, thought-provoking, and inspiring
discussions on science, nature, engineering, and future aspirations. Additionally, I appreciate the unique
perspectives and enriching conversations shared with Lionel Vincent, Sina Heydari, Basile Radisson, Janna
Nawroth, Anup Kanale, Tommaso Redaelli, JP Raimondi, Hao Cheng, Zitao Yu, Alyssa Chan, Kalyan Naik
iii



Banoth, and Morgan Jones—each of whom has contributed in their own way to my academic and personal
growth.
I gratefully acknowledge my undergraduate mentor, Yuyin Du, whose introduction to robotics — particularly in underwater vehicles and other fascinating projects — sparked my early interest in research
and eventually led me to pursue a deeper understanding of the physics of fluid dynamics. I also thank my
labmates Wei Zhang, Yang Yang, and Ji Li for the shared curiosity, support, and friendship that laid the
foundation for my research journey.
Finally, I extend my heartfelt appreciation to my family and friends for their unwavering support,
understanding, and belief in me. Their encouragement has been a constant source of strength throughout
this journey. My grandfather Xueqin Liu, a teacher at a University in Xi’an, embodied integrity and justice,
inspiring my pursuit of science and a true curiosity about the world. His spirit has shaped my values and
will continue to guide me as I grow into the person I aspire to be. I thank my mother, Huijuan Liu, the most
important person in my life, who raised me to be strong and independent. My beloved family members,
Xuemei Yin, Xianmin Jiang, Huimin, Liu, Yunyao Jiang, and Xingyao Liu, have always surrounded me with
love and encouragement. Finally, to Junchao Ma, who has walked beside me through the critical period of
this doctoral journey — no words are needed; thank you.
iv



Table of Contents
Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x
Chapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Oceanic microbes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 Life at low Reynolds numbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
Chapter 2: Motion in low Reynolds number: aquatic microorganisms . . . . . . . . . . . . . . . . . 4
2.1 Stokes Flow and General Solution in Spherical Coordinates . . . . . . . . . . . . . . . . . . 4
2.2 Non-ciliated sphere - sinking diatom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.3 Ciliated-sphere - ciliates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.4 Summary of analytical solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
Chapter 3: Nutrients transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.1 Spontaneous diffusion in nature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.2 Diffusion equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.3 Nutrient flux on surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
Chapter 4: Comparative study on swimming and attached ciliates . . . . . . . . . . . . . . . . . . . 18
4.1 Velocity scale and energy dissipation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.2 Advection-diffusion transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.4 Asymptotic analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
4.5 Numerical and asymptotic results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
4.6 Optimal feeding in attached ciliates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
4.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
Chapter 5: Microbial nutrients acquisition strategy: to swim or stay attached? . . . . . . . . . . . . 36
5.1 Biological observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
5.2 Mathematical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
v



5.4 Linking model prediction to biological data . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
Chapter 6: Flow Physics of Nutrient Transport Drives Functional Design of Ciliates . . . . . . . . . 49
6.1 Biological Observations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
6.2 Model modification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
6.3 Model predictions of optimal design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6.4 Linking model prediction and biological data . . . . . . . . . . . . . . . . . . . . . . . . . . 65
6.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
Chapter 7: Feeding in non-uniform concentration field . . . . . . . . . . . . . . . . . . . . . . . . . 70
7.1 Unsteady diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
7.2 Unsteady advection-diffusion transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
7.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
7.4 Unsteady transport in patchy concentration fields . . . . . . . . . . . . . . . . . . . . . . . 85
7.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
Chapter 8: Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
8.1 Beyond axisymmetric squirmers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
8.2 Sensing, Motion Control, and Learning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
Appendix A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
Notes on functions in spherical coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
A.1 Gegenbauer differential equation of degree ±
1
2
. . . . . . . . . . . . . . . . . . . . . . . . . 102
A.2 Envelope model with radial surface velocity . . . . . . . . . . . . . . . . . . . . . . . . . . 104
A.3 Bessel functions and spherical Bessel functions . . . . . . . . . . . . . . . . . . . . . . . . . 106
Appendix B . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
Numerical methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
B.1 Adjoint-based optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
B.2 Legendre spectral method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
B.2.1 Steady-state advection-diffusion equation . . . . . . . . . . . . . . . . . . . . . . . 112
B.2.2 Numerical scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
B.2.3 Unsteady advection-diffusion equation . . . . . . . . . . . . . . . . . . . . . . . . . 116
vi



List of Tables
1.1 Summary of published models and novel models . . . . . . . . . . . . . . . . . . . . . . . . 3
2.1 Comparison of Stokes flow around attached and swimming ciliate model . . . . . . . . . . 14
4.1 Summary of asymptotic analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
5.1 Survey of size a and flow measurements U in attached and swimming ciliates and sinking
diatoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.2 Summary of mathematical expressions for rigid and ciliated sphere models . . . . . . . . . 41
5.3 Summary of mathematical expressions for rigid and ciliated sphere models . . . . . . . . . 42
5.4 Summary of asymptotic results for Sh number of rigid and ciliated spheres . . . . . . . . . 44
6.1 Summary of morphological data of surveyed ciliates . . . . . . . . . . . . . . . . . . . . . . 54
6.2 Flow rates measurements for surveyed ciliates . . . . . . . . . . . . . . . . . . . . . . . . . 55
vii



List of Figures
1.1 Oceanic microbes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
2.1 Model spherical cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Modeling swimming and attached ciliates . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
4.1 Fluid flow for swimming and attached ciliates with same hydrodynamic power . . . . . . . 20
4.2 Comparison between swimming and attached mdoel ciliates . . . . . . . . . . . . . . . . . 23
4.3 Sherwood number as a function of Péclet number . . . . . . . . . . . . . . . . . . . . . . . 31
4.4 Sherwood number for hybrid surface motions . . . . . . . . . . . . . . . . . . . . . . . . . 32
4.5 Numerical optimization of surface motions that maximize feeding rates in attached ciliates 34
5.1 Phylogenetic tree . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.2 Stokeslet and Envelope models of attached and swimming ciliates . . . . . . . . . . . . . . 46
5.3 Comparison of Sh number versus Pe number for ciliates and diatom models . . . . . . . . 48
6.1 Ciliates utilize phagotrophy to ingest food particles . . . . . . . . . . . . . . . . . . . . . . 51
6.2 External and internal cost for nutrients transport for small and large cells. . . . . . . . . . 52
6.3 Mathematical modeling of ciliated cells with fractional feeding and ciliary areas . . . . . . 56
6.4 Flow and concentration fields of ciliated cells with vary feeding and ciliary areas . . . . . . 61
6.5 Optimal arrangements of cilia and feeding structures in sessile and motile cells . . . . . . . 64
6.6 Model predictions explain the diversity of cilia arrangements in sessile and swimming
ciliates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
7.1 Integration path in the complex plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
viii



7.2 Pure diffusion in concentration gradient . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
7.3 Unsteady transport in uniform and linear concentration . . . . . . . . . . . . . . . . . . . . 82
7.4 Nutrient intake of the first mode fitting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
7.5 Nutrient uptake in a concentration field with depleted upstream patch . . . . . . . . . . . 86
7.6 Nutrient uptake in a depleted environment with a rich concentration patch . . . . . . . . . 86
ix



Abstract
Nutrient acquisition strategies are fundamental to the survival and ecological roles of microorganisms. In
aquatic ecosystems, ciliated microbes have evolved specialized mechanisms to optimize nutrient uptake
by actively modifying their surrounding fluid. These organisms either swim freely or attach to surfaces,
generating feeding currents through the coordinated motion of cilia. The effectiveness of these strategies
is shaped by hydrodynamic constraints, yet the extent to which flow physics governs feeding efficiency
and ecological function remains an open question.
To understand their foraging strategies, I explore three aspects: (1) the benefits and limitations of
swimming versus attachment for nutrient acquisition, (2) the constraints imposed by flow physics and
nutrient transport on ciliate morphology, and (3) the effects of environmental variations, such as nutrient
gradients, on feeding efficiency. I solve the advection-diffusion equation governing nutrient transport in
cilia-driven flows using numerical and analytical techniques, including spherical decomposition, spectral
methods, and asymptotic analysis. To determine optimal ciliary motion and cell functional design for
maximizing nutrient uptake, I formulate and solve an inverse problem using adjoint-based optimization.
The findings reveal that, despite their contrasting behaviors, swimming and attached ciliates can
achieve hydrodynamically equivalent feeding rates, resolving a long-standing debate on the efficiency
of these two strategies. Furthermore, flow physics drive the functional design of ciliate morphology, shaping their ecological adaptations. Finally, the results suggest that ciliates may employ chemical sensing to
x



dynamically orient themselves for enhanced feeding in non-uniform nutrient environments. These findings open new directions for studying microswimmer locomotion in dynamic chemical landscapes and
investigating adaptive motion in response to environmental cues.
xi



Chapter 1
Introduction
1.1 Oceanic microbes
Oceanic microbes play an essential role in ecosystems by regulating key environmental processes, including cycling dissolved chemicals and driving the global carbon cycle [4, 6]. Their feeding activities are
fundamental to their biological fitness and ecological functions, spanning from small bacteria to larger
eukaryotes [127, 121, 50, 42, 34]. These organisms metabolize a diverse range of molecules, including
dissolved gases, organic compounds, and complex proteins, underscoring their essential role in ecosystem
dynamics[29, 17]. Understanding their feeding mechanisms, particularly among larger microbial cells with
specialized strategies, is crucial for elucidating their ecological impact.
Marine environments host an extraordinary diversity of unicellular organisms—including archaea, bacteria, and eukaryotes—which collectively constitute the majority of oceanic life. Among these, larger microorganisms display varied foraging strategies tailored to their ecological niches. For example, microbial
filter feeders such as choanoflagellates use specialized structures to optimize water filtration and prey
capture [103], while diatoms rely on sinking behavior to access nutrient-rich zones [139].
Ciliated microorganisms, shown in Fig.1.1 (purple and blue microbes), stand out for their ability to use
surface-bound cilia to generate fluid flows in highly viscous environments [36, 28, 68, 110, 125, 52]. These
ciliates display two primary behavioral modes: they either swim freely through the water column [89,
1



100 μm
1 μm
Figure 1.1: Oceanic microbes. An illustration of oceanic microbes, with focus on attached (blue) and
swimming (purple) ciliates.
90, 5, 24, 50] or attach to surfaces to generate feeding currents [126, 108, 109, 134, 143, 133, 5]. In both
cases, they actively manipulate their surrounding fluid to transport nutrients toward the cell body. This
mechanical interaction with the fluid is essential for nutrient uptake and survival, highlighting the intimate
link between biological function and physical principles.
1.2 Life at low Reynolds numbers
In the world of microscopic organisms whose typical size is in the order of 1−100µm with velocities of the
order of a few body lengths per second 10 − 1000µm · s
−1
, the inertia force is negligible, and the viscous
force is dominant. This effect can be quantified by the Reynolds number Re = ρUa/η, where U and a are
the characteristic velocity and size of the organism, ρ and η are the density and dynamic viscosity of the
fluid. Re number, representing the relative effect of inertia and viscosity of the fluid, can be in the range
of 10−4 − 10−1
for microorganisms.
2



To investigate the behavior of ciliated microbes, numerous mathematical models have been developed to study their locomotion, interactions with fluid environments, and nutrient transport mechanisms.
Table 1.1 summarizes key research problems addressed by existing models and highlights the novel contributions in this work. Despite substantial advancements, gaps remain in the study of attached ciliated
microorganisms, particularly those with partial ciliary coverage. Additionally, much of the prior research
has focused on uniform concentration fields, leaving an incomplete understanding of how ciliates operate
in non-uniform concentration environments.
This thesis addresses these gaps by presenting mathematical models for fluid dynamics and nutrient
transport. First, a model for attached ciliates is introduced, accompanied by a comparative analysis of ciliated and non-ciliated configurations. Next, the models are generalized to include cells with partial ciliary
coverage. Finally, the analysis is extended to explore spatiotemporal nutrient transport in non-uniform
concentration fields, emphasizing the dynamic interplay between ciliary motion and environmental gradients.
Table 1.1: Summary of published models and novel models presented in this work. Established studies
in previous studies ( black check marks ) and studies shown in this work (blue check marks). Notes with
’Numer.’ and ’Analy.’ represent that the solutions can be obtained numerically and analytically.
Rigid Ciliated (full cover) Ciliated (partial cover)
Sinking Stationary Swimming Attached Swimming Attached
Flow Field ✓ – ✓ ✓ ✓ ✓
Concentration ✓ ✓ ✓ ✓ ✓ ✓
(uniform) Numer. Analy. Numer. Numer. Numer. Numer.
(gradient) Numer. Analy. Numer. Numer. Numer. Numer.
Sh asymptotics ✓ – ✓ ✓ - -
3



Chapter 2
Motion in low Reynolds number: aquatic microorganisms
For incompressible fluid, the governing equation is the Navier-Stokes equations,
ρ(
∂u
∂t + (u · ∇)u) = −∇p + η∇2u + ρg,
∇ · u = 0,
(2.1)
where u is the velocity vector, p is pressure, ρ is density of the fluid, and η is the dynamic viscosity. At the
micron scale, the Reynolds number that quantifies the effects of inertial to viscous forces is Re≪ 1; that
is, the viscous forces play a dominant role, and inertia force is negligible.
2.1 Stokes Flow and General Solution in Spherical Coordinates
For zero Reynolds number, Re = 0 assumption, the governing equations for the fluid are the incompressible
Stokes equation, together with the continuity condition,
−∇p + η∇2u = 0,
∇ · u = 0.
(2.2)
where u is the velocity vector, p is pressure, ρ is density of the fluid, and η is the dynamic viscosity.
4



Stokes equations For solving fluid velocity u at zero Re limit, the incompressible Stokes equation, we
can consider following approach. Due to the continuity property of the fluid, the velocity vector u can
be expressed in terms of a vector potential Ψ as u = ∇ × Ψ [9]. Taking the curl of the Stokes equation
in Eq.(2.2), substituting u = ∇ × Ψ into the resulting equation and using the continuity condition, we can
obtain the classic result that the vector potential Ψ is governed by the bi-Laplacian ∇×∇×∇×∇×Ψ [51].
A
er
eθ
ez
ex
r
θ
ez
ex
ey
θ
φ
a
B
Figure 2.1: Model spherical cell. (A) Spherical model cell with coordinates (r, θ, ϕ), where θ ∈ [0, π] and
ϕ ∈ [0, 2π). (B) Due to axisymmetry, ϕ is an ignorable coordinate.
To solve for the fluid velocity in the fluid domain bounded internally by a spherical boundary of radius
a, it is convenient to introduce the spherical coordinates (r, θ, ϕ) and associated unit vectors (er, eθ, eϕ)
(see Fig. 2.2A). It is also convenient to introduce orthonormal frames(er, eθ, eϕ) and (ex, ey, ez) associated
with the spherical and Cartesian coordinates with ez = cos θer − sin θeθ. The fluid velocity can be
expressed in component form u ≡ (ur, uθ, uϕ). Here, we are interested only in axisymmetric flows, for
which uϕ = 0 is identically zero, and the components of the vector potential Ψ can be expressed in terms
of an axisymmetric stream function ψ (see, e.g., [9]),
Ψ =
0, 0,
ψ
r sin θ
T
. (2.3)
The non-trivial components (ur, uθ) of the fluid velocity are related to ψ(r, θ) as follows,
ur =
1
r
2 sin θ
∂ψ
∂θ , uθ = −
1
r sin θ
∂ψ
∂r .
(2.4)
5



The streamfunction ψ is governed by the biharmonic equation E2E2
(ψ) = 0 given in terms of the linear
operator E2
,
E
2 =
∂
2
∂r2
+
sin θ
r
2
∂
∂θ
1
sin θ
∂
∂θ
. (2.5)
This biharmonic equation for ψ can be solved as follows. By setting ψ1 the solution to the homogeneous
equation E2ψ1 = 0, and ψ2 the solution to the inhomogeneous equation E2ψ2 = ψ1 [51], the full solution
can be obtained as ψ = ψ1 + ψ2. This solution can be obtained analytically for arbitrary boundary
conditions in terms of r and the coordinate µ obtained via the nonlinear transformation of coordinate
µ = cos θ. Explicitly expression of ψ(r, µ) can be found in [51],
ψ(r, µ) = X∞
n=2
(anr
n + bnr
−n+1 + cnr
n+2 + dnr
−n+3)Fn(µ). (2.6)
Here, Fn(µ) = −
R
Pn−1(µ)dµ are solution functions related to the Legendre Polynomials of the first kind
Pn(µ), satisfying equation (1 − µ
2
)
d
2Fn
dµ2 + n(n − 1)Fn = 0 (see Appendix.A.1 for details), and an, bn, cn,
and dn are unknown coefficients.
By subtituting 2.6 into 2.4, we can obtain the velocity components (ur, uθ),
ur(r, µ) =
α0 + α13
1
r
3
+ α11
1
r

P1(µ) +X∞
n=2

αn+2
1
r
n+2 + αn
1
r
n

Pn(µ),
uθ(r, µ) =
−α0 + α13
1
2r
3
− α11
1
2r

V1(µ) +X∞
n=2
1
2

αn+2
n
r
n+2 + αn
n − 2
r
n

Vn(µ),
(2.7)
where α0, α11, . . . , αn are unknown coefficients, related to an, bn, cn, dn in Eq.(2.6), to be determined from
boundary conditions. The basis functions Vn(µ) are defined as
Vn(µ) = −
2
p
1 − µ2
Z
Pn(µ)dµ =
2
n(n + 1)
p
1 − µ2P
′
n
(µ), (2.8)
6



where P
′
n
(µ) = dPn(µ)
dµ . Both the Legendre polynomials Pn(µ) and basis functions Vn(µ) satisfy the
orthogonality conditions
Z +1
−1
Pn(µ)Pm(µ)dµ =
2
2n + 1
δnm,
Z +1
−1
Vn(µ)Vm(µ)dµ =
8
n(n + 1)(2n + 1) δmn. (2.9)
Analytical expressions of the pressure field The pressure is obtained by substituting Eq.(2.7) into
Stokes equation Eq.(2.2) and integrating, yielding the pressure field p(r, µ) as
p = p∞ + η
X∞
n=1
4n − 2
n + 1
αn
r
n+1Pn(µ). (2.10)
Analytical expressions of the stress field. The fluid stress tensor σ is given by σ = −pI + η(∇u +
∇Tu). For the axisymmetric flows considered here, σ admits three non-trivial stress components
σrr = −p + 2η
∂ur
∂r , σrθ = η

1
r
∂ur
∂θ + r
∂
∂r (
uθ
r
)

,
σθθ = −p + 2η

1
r
∂uθ
∂θ +
ur
r

.
(2.11)
Explicit expressions for the stress components are obtained by substituting Eq.(2.7) and Eq.(2.10) into Eq.(2.11).
Hydrodynamic force acting on the sphere The hydrodynamic force exerted by the fluid on the sphere
can be calculated by integrating the stress tensor σ over the surface S of the sphere. Due to axisymmetry,
only the force in the direction of the axis of axisymmetry, taken to be the z-axis, is non-zero,
F =
Z
S
σ · nˆdS = −4πηα11ez. (2.12)
Viscous dissipation energy The energy dissipation rate is defined as the volume integral over the
entire fluid domain V of the inner product of the velocity strain rate tensor e =
1
2
(∇u + ∇u
T) and stress
7



tensor σ, which, given proper decay at infinity, can be expressed as an integral over the surface of the
sphere by applying the divergence theorem (see, e.g., [66])
P =
Z
V
e : σdV = −
Z
S
u · (σ · nˆ)dS. (2.13)
This general solution is applicable to moving or stationary spheres with slip or no-slip boundary conditions. With the prescribed boundary conditions, we can obtain fluid solutions to the rigid spheres and
haired spheres in uniform flow, which are consistent with existing studies. In following sections, we can
apply the general solution approach to attached ciliated sphere in still fluid.
2.2 Non-ciliated sphere - sinking diatom
For modeling microorganisms with rigid skin such as diatoms, we can consider the classic problem of a
rigid sphere moving at constant velocity U in ambient fluid. Using a frame of reference attached to the
center of the sphere, the problem is equivalent to that of uniform flow past a rigid sphere. By symmetry,
the flow is axisymmetric and given by Eq. 2.7. The far-field boundary condition satisfies the uniform flow
condition, and the surface boundary condition is given by the no-slip velocity. Taking a to be the radius
of the sphere, the boundary conditions can be written as
u|r=a = 0, u|r→∞ = −Uez. (2.14)
8



By substituting boundary conditions, we can solve for the unknown coefficients and obtain α0 = 0, α13 =
U a3/2, α11 = 3U a/2 and αn+2 = 0, αn = 0,(n ≥ 2). Substituting the above coefficients back into the
general solution of velocity and pressure, we arrive at the velocity and pressure fields around the sphere
ur(r, θ) = (−1 −
a
3
2r
3
+
3a
2r
)U cos θ,
uθ(r, θ) = (1 −
a
3
4r
3
−
3a
4r
)U sin θ,
p = p∞ + η
3U a
2r
2
cos θ.
(2.15)
Finally, substituting the coefficients into general solution of force, we have the hydrodynamic force experienced by the sphere Fz = −6πηU aez.
2.3 Ciliated-sphere - ciliates
The envelope model is a classic model for ciliated micro-organisms. Cilia are a group of hair-like structure
on the surface, helping them swimming and feeding in liquid environment. Due to the surface cilia being
packed-growing and oscillating, we can consider the tips of those cilia to form a continuous deformable
surface covering the organism[75, 19]. Due to the velocity generated by cilia oscillation, the surface boundary conditions is slip boundary condition. The surface velocity is determined by ciliary motion. Consider
the cilia generated velocity can be expanded on Legendre polynomials Pn and Vn as [19]
ur|r=a =
X∞
n=0
An(t)Pn(µ), uθ|r=a =
X∞
n=1
Bn(t)Vn(µ), (2.16)
where A(t) and B(t) are time-dependent functions, relating to periodic motion of ciliary beating. Knowing
that the time scale of cilia beating is much smaller than the time scale of the fluid motion driven by cilia,
we consider the coefficients An, Bn to be constants. For ciliated microorganisms, the cilia length is much
shorter than the body length, the radial deformation can be ignored, such that only tangential velocity
9



is considered on the sphere surface. In the following model set up, we consider the surface boundary
conditions are
u|r=a =
X∞
n=1
BnVn(µ)eθ. (2.17)
Substitute the surface boundary conditions into the general solution of velocity Eq.(2.7), we have
ur|r=a = (α0 +
α13
a
3
+
α11
a
)P1 +
X∞
n=2
(
αn+2
a
n+2 +
αn
a
n
)Pn = 0,
uθ|r=a = (−α0 +
α13
2a
3
−
α11
2a
)V1 +
X∞
n=2
(
nαn+2
2a
n+2 +
(n − 2)αn
2a
n
)Vn =
X∞
n=1
BnV n.
(2.18)
By substituting the far-field boundary condition, we can solve for the constant coefficients α11, α13, αn,
(n = 0, 2, 3, 4...).
Mode 2
Mode 1
Mode 3
C
er
eθ
ez
ex
r
θ
-1 0 1
0
1
-1
/2 0
θ
µ = cos( )
2 /3 /3
Mode 2
Mode 1, motile
Mode 3
a
Mode 1, sessile
A
er
ez
ex
r
a θ
B
SwimmingAttached
Figure 2.2: Modeling swimming and attached ciliates. Spherical envelope model for (A) swimming
ciliate and (B) attached cilate. Ciliary motion is represented via a slip surface velocity. (C) First three
modes of surface velocity all at the same energy value: treadmill (mode 1), dipolar (mode 2), and tripolar
(mode 3) modes, corresponding to B1V1(µ) (B1 = 1 for attached and p
3/2 for swimming), B2V2(µ)
and B3V3(µ) with B2 =
√
3, B3 =
√
6 for both attached and swimming. Dotted lines represent lines of
symmetry of surface velocity.
Swimming ciliate To model a swimming ciliated microorganism, we can consider the swimming velocity of the squirmer is U. By using body frame of reference, the problem is equivalent to uniform flow
past a sphere with slip boundary condition, where the slip velocity is generated by cilia, and the velocity of
10



the uniform flow is the microorganism’s swimming velocity U, see Fig.2.2A. Thus, the far-field boundary
condition is set as uniform flow at infinity distance along axis ez. Thus, boundary conditions are
u|r=a =
X∞
n=1
BnVn(µ)eθ, u|r→∞ = −Uez. (2.19)
Some discussion of considering radial deformation is in appendix A.2.
Thus, similar to the rigid moving sphere problem, by substituting the far-field boundary condition, we
obtain coefficients α0 = −U,
α11 =
a
2
(3U − 2B1), α13 =
a
3
2
(−U + 2B1),
αn = −Bna
n
, αn+2 = Bna
n+2 (n ≥ 2).
(2.20)
Thus, the fluid velocity can be obtained as
ur =

− U + (−U + 2B1)
a
3
2r
3
+ (3U − 2B1)
a
2r

P1 +
X∞
n=2

a
n+2
r
n+2 −
a
n
r
n

BnPn,
uθ =

U + (−U + 2B1)
a
3
4r
3
− (3U − 2B1)
a
4r

V1 +
X∞
n=2

nan+2
2r
n+2 −
(n − 2)a
n
2r
n

BnVn.
(2.21)
The swimming velocity U is determined by surface cilia motion, thus, swimming velocity is an unknown
quantity. To solve for the flow field as well as the swimming speed, we need one more condition. Consider
this microorganism is freely swimming in the fluid, thus, force-free condition needs to be satisfied [89]
F =
Z
S
σ · ndS = 0. (2.22)
11



Substitute velocity coefficients Eq.(2.20) into general solution of the hydrodynamic force exerting on the
sphere body Eq.(2.12). To satisfy zero net force requirement, we require α11 =
a
2
(3U − 2B1) = 0. Thus,
we obtain the solution of fluid velocity and pressure as
ur(r, µ, t) =
−
2
3
+
2a
3
3r
3

B1(t)P1(µ) +X∞
n=2

a
n+2
r
n+2 −
a
n
r
n

Bn(t)Pn(µ),
uθ(r, µ, t) =
2
3
+
a
3
3r
3

B1(t)V1(µ) +X∞
n=2
1
2

nan+2
r
n+2 −
(n − 2)a
n
r
n

Bn(t)Vn(µ),
p(r, µ, t) = p∞ − η
X∞
n=2
4n − 2
n + 1
a
n
r
n+1 Bn(t)Pn(µ),
(2.23)
with the swimming speed obtained as U = 2B1/3.
Attached ciliate The attached squirmer has a stiff stalk that drags the sphere to make the sphere body
be fixed at one location. Thus, the attached squirmer is modeled as a fixed squirmer which is solely stirring
the fluid without locomotion and the stalk is infinitely long so that no wall effect on the fluid, see Fig.2.2B.
Use body frame of reference and consider ϕ-axis symmetry, we have the boundary conditions
u|r=a =
X∞
n=1
Bn(t)Vn(µ)eθ, u|r→∞ = 0. (2.24)
Substituting the far-field boundary condition, we obtain α0 = 0. By substituting α0, we can obtain coefficients as
αn = −Bna
n
, αn+2 = Bna
n+2
. (2.25)
12



Similarly, substitute the coefficients Eq.(2.25) into general solution of velocity Eq.(2.7) and pressure Eq.(2.10),
we obtain
ur(r, µ, t) = X∞
n=1

a
n+2
r
n+2 −
a
n
r
n

Bn(t)Pn(µ), (2.26)
uθ(r, µ, t) = X∞
n=1
1
2

nan+2
r
n+2 −
(n − 2)a
n
r
n

Bn(t)Vn(µ), (2.27)
p(r, µ, t) = p∞ − η
X∞
n=1
4n − 2
n + 1
a
n
r
n+1 Bn(t)Pn(µ). (2.28)
The net hydrodynamic force can be obtained by substituting the coefficients Eq.(2.25) into the general
solution of force Eq.(2.12),
Fz = 4πηaB1. (2.29)
Thus, the force generated by the stiff stalk should resolve force balance on the cilia body
F = Fstalk +
Z
S
σ · ndS = 0. (2.30)
The force generated by the stalk can be found as Fstalk = −4πηaB1ez.
2.4 Summary of analytical solutions
A summary of analytical solution of fluid around ciliated and non-ciliated model cells are listed in Table.
2.1.
13



Table 2.1: Comparison of Stokes flow around attached and swimming ciliate model. Mathematical
expressions of fluid velocity field, pressure field, hydrodynamic power, forces acting on the sphere for both
attached and swimming ciliated sphere. And the swimming speed for freely swimming ciliated sphere. All
quantities are given in dimensional form in terms of the radial distance r and angular variable µ = cos θ.
attached ciliated sphere
Fluid velocity field ur(r, µ) = P∞
n=1
a
n+2
r
n+2 −
a
n
r
n

BnPn(µ)
uθ(r, µ) = P∞
n=1
1
2

nan+2
r
n+2 −
(n − 2)a
n
r
n

BnVn(µ)
Pressure field p(r, µ) = p∞ − η
P∞
n=1
4n − 2
n + 1
a
n
r
n+1 BnPn(µ)
Energy dissipation rate P = 16πaη P∞
n=1
B2
n
n(n + 1)
Hydrodynamic force Fh = 4πηaB1ez
swimming ciliated sphere
Fluid velocity field ur(r, µ) =
−
2
3
+
2a
3
3r
3

B1P1(µ) + P∞
n=2
a
n+2
r
n+2 −
a
n
r
n

BnPn(µ)
uθ(r, µ) =
2
3
+
a
3
3r
3

B1V1(µ) + P∞
n=2
1
2

nan+2
r
n+2 −
(n − 2)a
n
r
n

BnVn(µ)
Pressure field p(r, µ) = p∞ − η
P∞
n=2
4n − 2
n + 1
a
n
r
n+1 BnPn(µ)
Energy dissipation rate P = 16πaη
1
3
B2
1 +
P∞
n=2
B2
n
n(n + 1)
Swimming speed U =
2
3
B1
non-ciliated sphere
Fluid velocity field ur(r, µ) = (−1 −
a
3
2r
3
+
3a
2r
)Uµ, uθ(r, θ) = (1 −
a
3
4r
3
−
3a
4r
)U
p
1 − µ2
Pressure field p = p∞ + η
3U a
2r
2
cos θ
Energy dissipation rate P = 6πaηU2
Drag force D = 6πµaU
14



Chapter 3
Nutrients transport
Feeding of oceanic microbes is essential for their biological fitness and ecological function [127, 121, 50,
42, 34]. The metabolic processes of microbes, from small bacteria to larger ciliates like the Stentor or
Paramecium, hinge on the absorption of particles or molecules at their surface [29, 17, 132, 134, 110, 133].
These particles or molecules vary widely depending on the organism and encompass dissolved oxygen
and other gases, lightweight molecules, complex proteins, organic compounds, and small particles and
bacteria. The typical random motion of these particles and bacteria is akin to a diffusive process at the
scale of the microorganism [84, 85, 90, 13, 14]. For simplicity, all cases will be collectively referred to as
"nutrients."
3.1 Spontaneous diffusion in nature
The energy needed for a passive small solute, such as molecule, to diffuse in water can be approximated
using the Stokes-Einstein relation; the thermal energy kBT, where kB is the Boltzman constant and T is
temperature, supports this diffusion process. For example, at T = 300K, kBT ∼ 10−21 Joules, and the
energy required to move a molecule of size b = 10−10m in water, where dynamic viscosity is η and the
diffusivity coefficient D = 10−9 m2
· s
−1
, is given by 6πηbD ∼ 10−21 Joules, thus thermal fluctuations
are sufficient to support diffusive transport.
15



3.2 Diffusion equation
Considering passive particles of nutrients with a size much smaller than that of the microorganism, the
surrounding nutrient concentration is governed by the diffusion process. Considering the concentration
of nutrients around the model cell is described as C(t, r, θ, ϕ). Here, t is time and (r, θ, ϕ) are spherical
coordinates. Assume metabolism rate of the cell is much faster than the time scale of diffusion, we can
consider concentration C of nutrients subject to zero concentration at the spherical surface [14, 17, 84, 90].
Considering diffusion of nutrients at constant diffusivity D [14] and starting from an initial concentration
field C(0, r, θ, ϕ) = Cb, the diffusion equation subject to zero concentration at the cell surface and proper
convergence to Cb as r → ∞,
∂C
∂t = D∇2C,
C|t=0 = Cb, C|r=a = 0, C|r→∞ = Cb.
(3.1)
If the background concentration is uniform C(0, r, θ, ϕ) = Cb = Co, where Co is constant, the concentration field evolves as C(t, r, θ, ϕ) = Co[1−a erfc
(r − a)/
√
4t

/r] [25]. As time t → ∞, the concentration
field is C(t → ∞) = Co[1 − a/r].
3.3 Nutrient flux on surface
To evaluate nutrient intake, we can use Fick’s law and define nutrient uptake is defined as the surface
integral of the nutrient flux over the cell surface I = −
H
n · (−D∇C)dS; here n is the outward unit
normal, dS = 2πa2
sin θdθ is the area element on the sphere, and the sign convention is chosen such that
I is positive if the cell takes up nutrients,
I(t) = Z
S
D∇C · ndS = 2πa2D
Z π
0
∂C(t)
∂r |r=a sin θdθ. (3.2)
1



where nˆ is outward normal of the sphere surface. In a uniform background concentration, the nutrient
uptake at the cell surface is given by I(t) = 4πaDCo(1 + a/√
πt). At steady state, it is I(t → ∞) =
4πaDCo.
17



Chapter 4
Comparative study on swimming and attached ciliates
Ciliated microorganisms near the base of the aquatic food chain either swim to encounter prey or attach at
a substrate and generate feeding currents to capture passing particles. Whether motile or sessile, these ciliates perform work against the surrounding fluid, creating flow fields that affect the transport of nutrients
and maintain a sufficient turnover rate of nutrients, unattainable by diffusion only [127].
Here, we can represent attached and swimming ciliates using a popular spherical model in viscous fluid
with slip surface velocity that afford analytical expressions of ciliary flows [19]. We can solve an advectiondiffusion equation for the concentration of dissolved nutrients, where the Péclet number (Pe) reflects the
ratio of diffusive to advective time scales. For freely swimming ciliates, the existing literature identifies
the "treadmill" mode of surface ciliary motion as optimal for both feeding efficiency and swimming speed
[90]. This raises two key questions: What is the optimal ciliary mode for attached ciliates, and how does it
differ from that of swimming ciliates? In this chapter, I evaluate, given a fixed amount of available energy,
the effect of surface velocities on feeding rates in attached ciliated microorganisms and it’s comparison to
swimming ciliated microorganisms ( see also [79]).
18



4.1 Velocity scale and energy dissipation
Constraint energy In order to compare attached and swimming ciliates that exert the same hydrodynamic power P on the surrounding fluid, we can introduce a characteristic velocity scale U based on the
total hydrodynamic power U =
p
P/(8πaη). To obtain non-dimensional forms of the equations and
boundary conditions, we can consider the spherical radius a = 1 and T = 1/U = 1 as the characteristic length and time scales of the problem. These considerations impose the following constraints on the
velocity coefficients Bn (Table 2.1)
attached : X∞
n=1
2B2
n
n(n + 1) = 1, swimming :
2
3
B
2
1 +
X∞
n=2
2B2
n
n(n + 1) = 1. (4.1)
Considering only the treadmill mode leads to B1 = 1 in the attached case and B1 =
p
3/2 in the swimming
case, with all other coefficients identically zero (Bn̸=1 = 0). That is, for the same hydrodynamic power
P, the attached sphere exhibits a slower surface velocity than the swimming sphere (B1 = 1 versus
B1 =
p
3/2), which can be proved using the reciprocal theorem [51]. Considering only the second mode,
we get B2 =
√
3 and Bn̸=2 = 0 in both the attached and swimming spheres, and when only the third mode
is considered, B3 =
√
6 and Bn̸=3 = 0 (Fig. 4.1). In the attached case, when a surface motion consists of
multiple modes simultaneously, the portion of energy assigned to each mode is denoted by β
2
n
, such that
Bn = βn
p
2/(n(n + 1)). For example, if the total energy budget P is equally distributed between the
first two modes, we have β
2
1 = 0.5, β2
2 = 0.5, and B1 =
√
0.5, B2 =
√
1.5.
In Fig. 4.1, it can be shown the streamlines around the swimming and attached spheres with slip surface
velocity corresponding to the treadmill mode (mode 1), dipolar mode (mode 2), and tripolar mode (mode
3). In all cases, the hydrodynamic power P/8πηa = 1 is held constant. Modes 2 and 3 produce zero net
force on the sphere, causing no net motion even in the swimming case, thus the streamlines are the same.
19



Mode 1, the treadmill mode, is the only mode that leads to motility. The associated streamlines are shown
in body-fixed frame in the swimming case and in inertial frame in the attached case.
A Mode 1, attached Mode 2Mode 1, swimming Mode 3
0
1.2
B C D Flow velocity/ U
Figure 4.1: Fluid flow for swimming and attached ciliates with same hydrodynamic power. Fluid
magnitude (colormap) and streamlines (black line) corresponding to the first three modes. First three
modes are at the same energy value, where the coefficient of mode 1 is (A) B1 =
p
3/2 for swimming (B)
B1 = 1 for attached ciliate, and that of mode 2 and mode 3 are B2 =
√
3, B3 =
√
6 for both attached and
swimming.
Reciprocal theorem Let u1 and σ1 be the fluid velocity and stress field at the surface of the attached
ciliate and let um = Uez + u1 and σm be those of the swimming ciliate; here, Uez represents the rigid
translation of the ciliate and u1 the ciliary treadmill motion, assuming both swimming and attached ciliates
exhibit the same u1. The total energy dissipation rates Ps and Pm of the attached and swimming ciliates,
respectively, are given by
Ps = −
I
u1 · σ1 · nˆdS, Pm = −
I
um · σm · nˆdS = −
I
u1 · σm · nˆdS, (4.2)
where, to obtain Pm, we can use the fact that H
ez ·σm ·nˆdS = 0 because the swimming ciliate is force-fee
in the ez direction.
Now apply the reciprocal theorem [51] to the two auxiliary problems (u1,σ1) and (um,σm),
I
u1 · σm · nˆdS =
I
um · σ1 · nˆdS. (4.3)
20



The left-hand side is equal to Pm. The right-hand side, recalling that um = Uez + u1, is equal to Ps −
U
H
ez · σ1 · nˆdS, that is,
Pm = Ps − U
I
ez · σ1 · nˆdS. (4.4)
By definition of the treadmill surface velocity u1(r = a, θ) = B1V1(cos θ)eθ and associated fluid velocity
(Table 2.1) and stress field Eq.(2.11), we can get that Pm = (16/3)πηaB2
1
and Ps = 8πηaB2
1
, as listed in
Table 2.1. We also can get that the last integral on the right-hand side is equal to U
H
ez · σ1 · nˆdS =
4πηaB1U, which for U = 2B1/3 is consistent with the results in Table 2.1. That is, for the same surface
velocity, the dissipation rate in the swimming case is lower by one-third compared to the dissipation rate
in the attached case. Alternatively, to maintain the same dissipation rate, the treadmill surface velocity in
the swimming case should be set to B1,swimming/B1,attached =
p
3/2, as done in previous text.
4.2 Advection-diffusion transport
Advection-diffusion equation To determine the effect of the advective flows generated by the ciliated
sphere on the nutrient concentration field around the sphere, we can consider the advection-diffusion
equation for the steady-state concentration C of nutrients subject to zero concentration at the spherical
surface [14, 17, 84, 90],
u · ∇C = D∆C, C(µ)|
r=a=1 = 0, C(µ)|
r→∞ = Co. (4.5)
We can normalize the concentration field by its far-field value Co at large distances away from the sphere
and consider the transformation of variable c = (Co − C)/Co [84, 90]. Writing the advection-diffusion
equation and boundary conditions Eq.(4.5) in non-dimensional form in terms of the new variable c(r, µ)
yields
Pe u · ∇c = ∆c, c(µ)|
r=a=1 = 1, c(µ)|
r→∞ = 0, (4.6)
21



where the Péclet number Pe = aU/D, which quantifies the ratio of diffusive to advective time scales.
When advection is dominant, the advective time is much smaller than the diffusive time and Pe ≫ 1,
when diffusion is dominant, the advective time is much larger and Pe ≪ 1. At Pe = 1, the two processes
are in balance. For intermediate Pe number, the concentration field is solved numerically by using Legendre
spectral method [90, 79], see details in B.2.
Sherwood number The nutrient uptake is estimated by the Sherwood number,
Sh = I
Idiffusion
= −
1
4πaDCo
I
nˆ · (−D∇C)dS = −
1
2
Z 1
−1
∇c · er|
r=a=1 dµ. (4.7)
where and the steady-state inward current due to molecular diffusion is given by Idiffusion = 4πaDCo
shown in Chapter.3.
4.3 Numerical results
By substituting the analytical solutions of the flow field u from Table 2.1 into Eq.(4.6). We can arrive at
governing equations for the concentration field c, which can be solved numerically using a spectral method
(Appendix.B.2) and analytically in the asymptotic limit of small and large Péclet numbers in the following
section.
A comparison of the steady-state concentration field around the swimming and attached spheres at Pe
= 100 and Pe = 1000 is shown in Fig.4.2. In the treadmill mode, the swimming sphere swims into regions
of higher concentration, which thins the diffusive boundary layer at its leading surface, leaving a trailing
plume or “tail” of nutrient depletion. Similar concentration fields are obtained in the attached sphere
because, although the sphere is fixed, the surface treadmill velocity generates feeding currents that bring
nutrients towards the surface of the sphere. In the dipolar and tripolar modes, feeding currents bring fresh
nutrients to the spherical surface from, respectively, two opposite and three nearly equiangular directions.
22



Concentration
Sh = 6.7 Sh = 6.0 Sh = 3.0 Sh = 7.6 0%
A
motile sessile
mode 2 mode 3
Sh = 20.8 Sh = 22.7 Sh = 9.8
mode 1 mode 1
Pe=100Pe=1000
motile sessile
Sh = 23.2
100%
B
Figure 4.2: Comparison between swimming and attached mdoel ciliates. Flow streamlines (white)
and concentration fields (colormap) at (A) Pe = 100 and (B) 1000 for the same hydrodynamic power P = 1
and distinct surface motions. In the treadmill mode, the streamlines, concentration field, and Sh number
differ between the attached and swimming spheres, but are the same in the dipolar and tripolar surface
modes.
By comparing the Sh number associated with each mode at both Pe = 100 and Pe = 1000, it can
be shown larger Pe leads to higher nutrient uptake. In the treadmill mode, evaluating Sh number for
the swimming and attached spheres led, respectively, to 7.6 and 6.7 at Pe = 100 and 23.2 and 20.8 at
Pe = 1000. Indeed, at each Pe, the swimming sphere with treadmill surface velocity produced the largest
Sh number, implying that for the same hydrodynamic power, motility maximized nutrient uptake. But the
percent difference in nutrient uptake between the swimming and attached sphere decreased with increasing Pe, from 13.4% to 11.5%. For the swimming sphere, comparing the treadmill mode and dipolar mode
(mode 2), it can be found Sh = 7.6 and 6.0 at Pe = 100 and Sh = 23.2 and 22.7 at Pe = 1000. That is, for
the swimming sphere, the difference in nutrient uptake between mode 1 and mode 2 also decreased with
increasing Pe from 26.7% to 2.2%. Interestingly, for the attached sphere, the nutrient uptake in mode 2
(Sh = 22.7) exceeded that of mode 1 (Sh = 20.8) at Pe = 1000, with a 9.1% increase.
23



4.4 Asymptotic analysis
To understand the asymptotic behavior of Sh in the limit of large Pe, I performed an asymptotic analysis
to obtain the scaling of Sh with Pe for the attached ciliated sphere. The approach is similar to that used
in [84, 90] for the treadmill mode in the swimming ciliated sphere. In this section, I will show derivations of
asymptotic solutions of the Sherwood number for a attached ciliated sphere. Particularly, I seek asymptotic
expressions associated with mode 1 and mode 2.
Large Pe limit Here, we can can start with a general expression of velocity field corresponding to the
n
th mode surface velocity at the unit sphere a = 1, for which uθ(µ)|
r=1 = BnVn(µ), where n = 1, n = 2
are considered later in this paper. The flow field corresponding to n
th mode is given by (see Tables 2.1)
ur =

1
r
n+2 −
1
r
n

BnPn(µ), uθ =

n
r
n+2 −
n − 2
r
n

Bn
2
Vn(µ), (4.8)
Take the Taylor series expansion of flow field at r = 1 and keep only the leading order terms,
ur(r, µ) = ur|
r=1 +
∂ur
∂r


r=1 (r − 1) + · · · = −2Bn(r − 1)Pn(µ) + . . . ,
uθ(r, µ) = uθ|
r=1 +
∂uθ
∂r



r=1
(r − 1) + · · · = BnVn(µ) + . . . .
(4.9)
We can define the temporary variable y = r − 1 (not to confused with the y-coordinated in the inertial
(x, y, z) space). The region y ≪ 1 represents a thin boundary layer around the spherical surface. Since the
concentration boundary layer is expected to be thinner as Pe increases, we can rescale r−1 = y = Pe−mY ,
where Y is a new variable.
24



Next, we can substitute the linearized flow field u from Eq.(4.9) into the advection-diffusion equation Eq.(4.5) and use the new variable r − 1 = Pe−mY . Keeping only the leading-order terms, we can
obtain the advection operator u · ∇,
ur
∂
∂r = −2BnPnY
∂
∂Y , uθ
1
r
∂
∂θ = −BnVn
p
1 − µ2
1
Pe−mY + 1
∂
∂µ. (4.10)
Similarly, rewriting the Laplacian operator using the new variable r − 1 = Pe−mY , we arrive at
∇2 = Pe2m ∂
2
∂Y 2
+
2Pem
Pe−mY + 1
∂
∂Y +
1
(Pe−mY + 1)2
∂
∂µ
(1 − µ
2
)
∂
∂µ
. (4.11)
The leading-order term in the Laplacian operator scales with Pe2m, while the leading-order term in the
advection operator scales with Pe0
. Matching order on both sides of the advection-diffusion equation, we
have 2m = 1. Thus, m =
1
2
, and we have r − 1 = Pe−1/2Y .
We now substitute Eq.(4.10) and Eq.(4.11), with m = 1/2 into the dimensionless advection-diffusion
equation Eq.(4.6), keeping only the leading-order terms, we arrive at
−Bn

2PnY
∂c
∂Y + Vn
p
1 − µ2
∂c
∂µ
=
∂
2
c
∂Y 2
. (4.12)
We can define a similarity variable Z = Y /g(µ) such that, by the chain rule,
∂
∂Y =
1
g
∂
∂Z ,
∂
∂µ = g
′ ∂
∂g = −
Zg′
g
∂
∂Z . (4.13)
We can substitute into Eq.(4.12) and rearrange terms to arrive at the ordinary differential equation
∂
2
c
∂Z2
+ Bn

2Png
2 −
1
n(n + 1)(1 − µ
2
)P
′
n
gg′

Z
∂c
∂Z = 0. (4.14)
25



For a similarity solution to exist, the term 2µg2 − (1 − µ
2
)gg′ needs to be equal to a constant. Setting the
value of this constant to 2, the problem becomes that of solving
d
2
c
dZ2
+ 2Z
dc
dZ = 0, Bn

2Png
2 −
1
n(n + 1)P
′
n
(1 − µ
2
)(g
2
)
′

= 2. (4.15)
For the first mode, n = 1 and P1 = µ, the solution is in the form,
c(Z) = C1erf(Z) + C2, g2
(µ) = C3 − 12µ + 4µ
3
3B1(µ2 − 1)2
.
(4.16)
Here, C1, C2 and C3 are unknown constants to be determined from the boundary conditions c(Z = 0) = 1,
c(Z → ∞) = 0, and from the condition that at µ = 1, the concentration field and the function g(µ) must
be bounded [84]. Put together, we can get that C3 = 8, C2 = 1, and C1 = −1. Thus, the asymptotic
solution of the concentration field in the limit of large Péclet, Pe ≫ 1, is given by
c(Z) = 1 − erf(Z) = erfc(Z),
(4.17)
where
Z = Pe
1
2
(r − 1)
g(µ)
, g(µ) = s
8 − 12µ + 4µ3
3B1(1 − µ2)
2
. (4.18)
The Sherwood number at large Pe is given by,
Shmode 1 = −
1
2
Z 1
−1
∂c
∂r




r=1
dµ =
1
2
Pe
1
2
Z 1
−1
2
√
π
e
−Z2 1
g




r=1
dµ
=
1
√
π
Pe
1
2
Z 1
−1
s
3B1(1 − µ2)
2
8 − 12µ + 4µ3
dµ =
r
4B1
3π
Pe
1
2 ,
(4.19)
26



Similarly, for the second mode, n = 1 and P2 =
1
2
(3µ
2 − 1), we can obtain the solution
Z = Pe
1
2
(r − 1)
g(µ)
, g(µ) = s
1 + µ4 − 2µ2
B2(1 − µ2)
2µ2
. (4.20)
The Sherwood number at large Pe is
Shmode 2 =
r
B2
π
Pe
1
2 ,
(4.21)
Small Pe Limit At small Péclet number, we can expand the concentration field as
c = Pe0
c0 + Pe1
c1 + Pe2
c2 + . . . (4.22)
By substituting the expanded concentration into the dimensionless advection-diffusion equation Eq.(4.6),
we can arrive at the following system of equations, to be solved at each order in Pe,
Order 0 : 0 = ∇2
c0, c0|
r=1 = 1, c0|
r→∞ = 0.
Order 1 : u · ∇c0 = ∇2
c1, c1|
r=1 = 0, c1|
r→∞ = 0.
Order 2 : u · ∇c1 = ∇2
c2, c2|
r=1 = 0, c2|
r→∞ = 0.
(4.23)
At the leading order, the solution is simply c0 = 1/r. To find the solution at higher orders, we can substitute
the velocity field Eq.(4.8) into the higher order equations in Eq.(4.23). At order Pe1
, we can get
−Bn

1
r
n+2 −
1
r
n

1
r
2
Pn(µ) = ∇2
c1, c1(r = 1) = 0. (4.24)
27



Recall that the Legendre polynomials satisfy the Legendre differential equation d
dµ[(1 − µ
2
)P
′
m] + m(m +
1)Pm = 0. Expanding c1 in terms of Legendre polynomials c1(r, µ) = P∞
m=0 Rm(r)Pm(µ), and substituting back into the above equation, we can get
−Bn

1
r
n+2 −
1
r
n

Pn(µ) = X∞
m=0

d
dr (r
2
dRm
dr ) − m(m + 1)Rm

Pm(µ). (4.25)
By equating the Legendre polynomials on both sides in the above equation, we find that only the term
Rn(r) survives,
−Bn

1
r
n+2 −
1
r
n

= r
2R
′′
1n + 2rR′
1n − 6R1n. (4.26)
For the first mode, solving the above ordinary differential equations with n = 1, taking into consideration
the boundary conditions, we find the solution for c1,
c1(r, µ) =
3
4
r
−2 −
1
4
r
−3 −
1
2
r
−1

B1µ. (4.27)
Repeating the same procedure at order Pe2
, we find that
∇2
c2 = B
2
1

4 − 9r + 2r
2 + 3r
3
12r
7
P0(µ) −
(r − 1)(5 − 4r − 9r
2 + 6r
3
)
8r
7
2
3
P2(µ)

. (4.28)
We can expand c2 =
P∞
m=0 Rm(r)Pm(µ) in terms of a Legendre polynomial expansion with unknown
Rm(r) and substitute back in the above equation; we get that only the terms R0(r) and R2(r) survive,
such that the general solution for c2(r, µ) is given by
c2(r, µ) = R0(r)P0(µ) + R2(r)P2(µ). (4.29)
28



One can readily verify that the ordinary differential equation governing R0(r) is given by
r
2R
′′
0 + 2rR′
0 = B
2
1
4 − 9r + 2r
2 + 3r
3
12r
5
. (4.30)
The solution to this equation, taking into account the boundary conditions, is of the form
R0(r) = B
2
1

−
77
720r
+
1
60r
5
−
1
16r
4
+
1
36r
3
+
1
8r
2

. (4.31)
It turns out that, for computing the Sherwood number below, only R0(r) is needed, as shown in
−
1
2
Z 1
−1
∂c2
∂r




r=1
dµ = −
1
2
X∞
m=0
dRm
dr




r=1
Z 1
−1
Pm(µ)dµ =
dRm
dr




r=1
δm0 (4.32)
with implementing the property of Legendre polynomials R 1
−1
Pm(µ)dµ = 2δm0. We thus need not calculate R2(r). For the first mode, Sh in the small Pe limit is given by
Shmode 1 = −
1
2
Z 1
−1
∂c0
∂r




r=1
dµ −
1
2
Z 1
−1
∂c1
∂r




r=1
dµ
Pe −
1
2
Z 1
−1
∂c2
∂r




r=1
dµ
Pe2
= 1 +
43B2
1
720
Pe2
.
(4.33)
Similarly, we can compute Sherwood number for surface velocity containing only the second mode by
replacing velocity field with n = 2,
Shmode 2 = 1 +
41B2
2
25200
Pe2
. (4.34)
4.5 Numerical and asymptotic results
Summary of asymptotics In Table 4.1, I summarize the results of the above asymptotic analysis for
mode 1 and mode 2 of the attached sphere. Note that the asymptotic analysis of mode 2 applies equally to
29



Table 4.1: Summary of asymptotic analysis Asymptotic expressions for Sh as a function of Pe for attached and swimming ciliated sphere model. The velocity coefficients associated with each mode are
chosen satisfying same constraint hydrodynamic power.
Large Pe limit Small Pe limit
attached swimming attached swimming
Mode 1 Sh =
2
√
3π
Pe
1
2 Sh =
2
√
3π

3
2
1
4
Pe
1
2 , Sh = 1 +
43
720
Pe2 Sh = 1 + 1
3
r
3
2
Pe,
(present) [84, 90] (present) [84, 90]
Mode 2 Sh =
1
√
π
3
1
4Pe
1
2 (present) Sh = 1 +
41
8400
Pe2
(present)
the attached and swimming spheres. For comparison, I also include in this table the asymptotic results of
[84, 90] for the swimming ciliated sphere. These asymptotic results are superimposed onto the numerical
computations in the two insets in Fig.4.3A and B.
At small Pe (Pe ≪ 1), in the treadmill mode, Sh scales with Pe2
and Pe1
, respectively, for the attached
and swimming spheres and, in the dipolar mode, Sh scales with Pe2
. That is, at Pe ≪ 1, the treadmill mode
outperforms the dipolar mode in the swimming case, and the swimming sphere outperforms the attached
sphere in the treadmill mode.
At large Pe (Pe ≫ 1), in the treadmill mode, Sh scales with √
Pe in both the attached and swimming
ciliated spheres. That is, there is no distinction in the scaling of Sh with Pe between the swimming and
attached spheres. Interestingly, it can be found that in the dipolar mode, Sh also scales with √
Pe, indicating
no distinction between the treadmill and dipolar modes. The constant coefficients in the asymptotic scaling
differ slightly: Sh ≈ 0.65√
Pe (treadmill, attached) and Sh ≈ 0.72√
Pe (treadmill, swimming), and Sh ≈
0.74√
Pe (dipolar, both). Thus, in the swimming sphere, the treadmill and dipolar modes perform nearly
similarly in the large Pe limit, while in the attached sphere, the dipolar mode outperforms the treadmill
mode by a distinguishable difference (by about 10%).
30



Comparison of results To probe the trends in Sh number over a larger range of Pe values, I computed,
for the same hydrodynamic power, the steady state concentration field for the attached and swimming
sphere, and for mode 1, 2, and 3, for Pe ∈ [0, 1000]. In Fig.4.3A, the results are Sh for the attached sphere.
The solid lines in color blue, purple and grey represent mode 1, mode 2, and mode 3, respectively. Mode 1
exhibits the best feeding performance (highest Sh) for Pe ≤ 284. Mode 2 exceeds mode 1 after this critical
Pe value. At Pe = 1000, Sh number of mode 2 is about 10% higher than that of mode 1. Numerical results
for the swimming sphere are shown in Fig.4.3B. Mode 1 outperforms modes 2 and 3 for the entire range
of Pe values, but the difference between mode 1 and mode 2 seems to decrease and the two modes seem to
approach each other asymptotically at larger Pe.
Sherwood number, Sh
B
Small Pe
2
1
=
1
2
1
2
10 0
10-4
10-2
mode1
mode2
mode3
A
5
10
15
20
25
1
5
10
15
20
25
1
Péclet number, Pe
0
mode1
mode2
mode3
Péclet number, Pe
0 500 1000
100
10-1 101
2
100 2
10-4
10-2
Sh - 1
10 0
10 -1 10 1
Small Pe
1
2
1
2
Sessile Motile
15
20
25
500 1000
15
20
25
103
500
Large Pe Large Pe Sherwood number, Sh
Sh
Pe
Sh - 1 Sh
slope
Pe
500 1000
Sessile Motile
Figure 4.3: Sherwood numebr as a function of Péclet number. (A) attached ciliate model and (B)
swimming ciliate model for the same hydrodynamic power P = 1. Solid lines are numerical calculations
for mode 1 (blue), mode 2 (purple), and mode 3 (grey). Dashed lines and scaling laws in the limit of large
and small Pe are obtained from asymptotic analysis for mode 1 (blue) and mode 2 (purple).
4.6 Optimal feeding in attached ciliates
By keeping the total hydrodynamic power constant, we can investigate numerically how Sh number varies
when multiple surface modes coexist in both open-loop and closed-loop approaches.
Open-loop For an open-loop approach considering two modes, we can distribute the total hydrodynamic power into the first two modes, assigning an energy fraction β
2
1
to mode 1 and the remaining
31



1
2
20
25
15
Pe = 1000Pe = 10
A
0
5
10
15
20
0
5
10
800
15
400
20
0
25
0
1
0.5
Sherwood number, Sh
Pe
0.50 10.50 1
-30
-10
0
10
0.50 1
B C Pe = 1000
↓20%
↓3%
↓2%
.9
10%
.9 .1
↓
↓
↓
↓
10% 10% ↓
↓
more
robust
less
robust
β
2
1 β
2
1 β
2
1
∂Sh
∂β2
1
Sh Sh
β
2
1
Figure 4.4: Sherwood number for hybrid surface motions. (A) Hybrid surface motions with two
fundamental modes, treadmill and dipolar, and constraint hydrodynamic power P/8πηa = β
2
1 + β
2
2 = 1,
where β
2
1
represents the portion of the energy assigned to mode 1. Sh versus Pe and β
2
1
as β
2
1
vary from 0
to 1. Close up at Sh number versus β
2
1
at (B) Pe = 10 and (C) Pe = 1000 (left), and partial derivative of Sh
with respect to β
2
1
(right).
fraction 1 − β
2
1
to mode 2. I varied β
2
1
from 1 (all hydrodynamic power assigned to mode 1) to 0 (all hydrodynamic power assigned to mode 2) at fixed intervals ∆β
2
1 = 0.05. I also varied Pe from 0 to 1000 at
∆Pe = 0.1. For each combination of (Pe, β2
1
), I computed the steady-state concentration field and calculated the resulting Sh number. I found that at small Pe, Sh increased monotonically as β
2
1
varied from 0 to
1, indicating that mode 1 is optimal. At larger Pe, a new local maximum appeared at β
2
1 = 0 (mode 2). This
change is evident when comparing Fig. 4.4B and Fig. 4.4C, which illustrate Sh as a function of β
2
1
at Pe = 10
and Pe = 1000, respectively. At Pe = 10, the maximal Sh at β
2
1 = 1 is a global optimum. At Pe = 1000, two
local optima in Sh number are obtained at β
2
1 = 1 and β
2
1 = 0, with Sh|β
2
1=0 > Sh|β
2
1=1, indicating that
mode 2 is a global optimum. Interestingly, when calculating the sensitivity ∂Sh/∂β2
1
of these maxima to
variations in β
2
1
, I found that the maximum at β
2
1 = 0 (mode 2) is more sensitive to variations in β
2
1
, with

∂Sh/∂β2
1


β
2
1=0 ≫

∂Sh/∂β2
1


β
2
1=1. That is, a small variation in β
2
1
leads to a larger drop in Sh at β
2
1 = 0
(mode 2), while the same variation in β
2
1
leads to a small drop in Sh at β
2
1 = 1 (mode 1). In Fig.4.4C, it is
shown that a 10% variation of β
2
1
leads to a 3% drop from the optimal value at mode 1 and 20% drop from
optimal value at mode 2. Similar results can be obtained by considering three modes coexist [79].
Closed-loop The optima in Fig. 4.4 are identified in an open-loop search over the parameter spaces
(β
2
1
,Pe) and (β
2
1
, β2
2
,Pe). Such open-loop search becomes unfeasible when considering higher-order
32



surface modes. A closed-loop optimization algorithm that seeks surface velocities that optimize Sh is
needed. Here, I adapted the adjoint optimization method with the gradient ascent algorithm shown in
Appendix.B.1 [90].
In Fig. 4.5, I considered initial surface velocities with 10 modes, satisfying the constraint on the total
hydrodynamic power PN
n=1 β
2
n = 1. I used the closed loop optimization algorithm to identify optimal
surface velocities that maximize Sh. In Fig. 4.5A, I show the optimization results at Pe = 10 and two distinct
initial conditions (grey line). The optimization algorithm converges to an optimal solution (black line)
that is close to mode 1 (superimposed in blue). The energy distribution among all modes as a function of
iteration steps shows that, while the initial energy was distributed among multiple modes, in the converged
solution, energy is mostly assigned to mode 1. Indeed, comparing Sh number (black marker ’∗’) to Sh
number of mode 1 (blue line) shows that the optimal Sh converges to that of mode 1. Flow streamlines and
concentration fields at these optima are shown in the bottom row of Fig. 4.5A.
In Fig. 4.5B, I show the optimization results at Pe = 1000 and the same two initial conditions (grey line)
considered in Fig. 4.5A. Here, unlike in Fig. 4.5A, the optimization algorithm converges to two different
optimal solutions depending on initial conditions: one optimal (black line) is close to mode 1 (superimposed
in blue) and the other optimal is close to mode 2 (superimposed in purple), as reflected in the energy
distribution (second row) and in Sh values (third row). Flow streamlines and concentration fields at these
two distinct optima are shown in the bottom row of Fig. 4.5B. The second optimal, the one corresponding
to mode 2, exhibits higher Sh value. These results are consistent with the open loop analysis presented in
Figs. 4.3 and 4.4.
4.7 Summary
For swimming ciliates, it is found that assigning all energy into a treadmill surface velocity that induces
swimming is an optimal strategy for maximizing feeding rate but only in a finite Pe range. This is in
33



10
20
10
20 Sh
Sherwood number, Sh Energy
1
2
optimal
initial
optimal
mode 1
-1 10
3
Sh
Shmode 1
0
1
2
0
1
Iteration
0 3 1 4 2
Iteration
0 21 3 4mode 1
5
Pe = 10 Pe = 1000 B
Sh = 2.2 Sh = 2.2 Sh = 20.3
Sh
Shmode 2
mode 2
mode 2
Iteration
0 1 4 32
Sh
Shmode 1
mode 1
1
2
3
Surface velocity u
mode 1
0
1
0
1
0
1
initial
optimal
mode 1
-1 10
initial
optimal
mode 1
-1 10 -1 10
initial
optimal
A
Iteration
0 321
Shmode 1
0
1
-1
0
1
-1
0
1
2
Sh = 21.9
Figure 4.5: Numerical optimization of surface motions that maximize feeding rates in attached
ciliates. (A) Optimization results at Pe = 10 for two different initial guesses. (B) Optimization results
at Pe = 1000 for the same two different initial guesses. Top row shows initial surface velocity (grey),
optimal surface motion (black), and mode 1 (blue). Second row shows the energy distribution among
ten different velocity modes at each iteration of the numerical optimization process. Third row shows
Sherwood number (black) at each iteration, Sh number for only mode 1 Shmode1 (blue), and the difference
in Sh number (Shmode1-Sh) (grey). Last row shows fluid and concentration fields under optimal surface
conditions. In panel B, at Pe = 1000, the numerical optimization algorithm converges to one of two distinct
solutions, depending on initial conditions, that are either close to mode 1 (blue) or mode 2 (purple).
34



contrast to the findings in [90] that claimed optimality of the treadmill mode for all Pe. In the limit of large
Pe, there is no distinction in nutrient uptake between the treadmill mode and the symmetric dipolar mode
that applies zero net force on the ciliated body inducing no swimming motion.
For attached ciliates, is is found that the treadmill mode achieves optimal feeding rate at relatively low
Pe values below a critical value Pecr ≈ 280, while the dipolar mode becomes optimal for Pe exceeding
this threshold. The asymptotic analysis supports that, in the large Pe limit, the dipolar surface mode
outperforms other modes in terms of feeding rate. However, the sensitivity analysis shows that although
at large Pe the treadmill mode leads to lower Sherwood values, it is more robust to perturbations in the
surface velocity. The dipolar mode leads to higher Sherwood values but it is significantly more sensitive to
perturbations, with feeding efficiency dropping rapidly even with small perturbations in surface motion.
This finding challenges previous assumptions that motility inherently improves feeding rate in ciliates [90, 5].The optimal cilia-driven surface velocity for maximizing feeding rate varies significantly depending on the Péclet number, with distinct advantages observed for both swimming and attached ciliates
under different Pe numbers. This study enriches the understanding of the complexity of feeding strategies
in ciliated microorganisms and highlights the importance of considering various environmental conditions
when evaluating the ecological roles and evolutionary adaptations of these microbes.
35



Chapter 5
Microbial nutrients acquisition strategy: to swim or stay attached?
Nutrient uptake is a fundamental challenge for all microbes, driven by the patchy and dynamic nature
of nutrient distributions in marine environments. Among small, single-celled protists—particularly ciliates—two main strategies have evolved: actively swimming to seek out nutrient-rich regions, or remaining attached to a surface while generating feeding currents to draw in nutrients. Both strategies—“swim”
and “stay”—have independently emerged and diversified across many branches of the microbial tree of
life 5.1. This widespread occurrence raises a paradox: if one strategy were truly optimal, we might expect it to be consistently favored and reflected across evolutionary lineages. Yet, both behaviors persist in
modern species, suggesting that optimality is context-dependent and shaped by ecological and physical
constraints.
In this chapter, we investigate this evolutionary trade-off through a comparative study of motile and
sessile ciliates. By modeling their fluid-driven nutrient uptake under varying environmental conditions,
we examine whether one strategy offers a universal advantage—or whether the coexistence of both points
to a deeper adaptive balance shaped by the physics of the microscopic world.
In this chapter, I additionally include a comparison with diatoms to distinguish the effects of relative body motion independent of cilia-driven feeding currents. Furthermore, a combination of a survey
36



of existing biological observations with model predictions will be presented to provide a study of these
evolutionary adaptations (see also [78]).
5.1 Biological observations
Tree of life for microbes To observe the biological data, a survey on the morphology, flows, and phylogenetic lineage of ciliates and diatoms is presented in Fig. 5.1, [65, 106, 53, 48]. Attached ciliates, such as
the Stentor [134], Opercularia [125, 143], and Vorticella [126, 104, 133, 98, 109], are characterized by a ciliary
crown, where the motion of beating cilia entrains fluid toward the cell. The cell body and ciliary crown are
positioned away from the surface they live upon, usually with a stalk, to minimize the effect of that surface
on slowing down the cilia-driven microcurrents [73, 21, 108, 109]. Further, to avoid generating recirculating microcurrents and reduce reprocessing of depleted water [125, 133, 110], attached ciliates actively
regulate their orientation to feed at an angle relative to the substratum [109]. At optimal inclination, the
effective cilia-generated force is nearly parallel to the bounding substrate and creates quasi-unidirectional
flows that drive nutrients and particles past the cell feeding apparatus [109, 135, 67]. In swimming ciliates,
such as the Paramecium and Volvox, the surface of the organism is often entirely covered with cilia that
beat in a coordinated manner and power the organism to swim through the surrounding fluid [24, 23, 71,
76]. Diatoms lack motility apparatus and sink by regulating their buoyancy [63, 93, 43] (Fig. 5.1B).
Kinematic measurements Empirical flow measurements around attached [108, 107, 143, 52, 98, 135]
and swimming [36, 32] ciliates are sparse. Here, considering morphometric and flow data from ten species
of attached ciliates [133, 126, 125, 98, 109, 134], ten species of swimming ciliates [24, 23, 76], and seven
species of diatoms [93, 43], a summary of the ranges of sizes and characteristic speeds are shown in Table 5.1. (detailed measurements are shown in [78]). Size is represented by the volume-equivalent spherical
radius a. The characteristic speeds U for attached ciliates are based on the maximal flow speeds measured
37



Urocentrum Paramecium Opercularia VorticellaEpistylisZoothamniun
Eukaryotes
Heterotrichea
Sessilida
Stentor
Peniculida
Spirotrichea
Opalinata
Viridiplantae SAR
Volvox
Tintinnida
Class PhylumOrderFamilyGenus Clade
Oligohymenophorea
Hymenostomatida
Chlorophyta Ciliophora Ochrophyta
Chlorophyceae Bigyra
Alveolata Stramenopiles
Chlamydomonadales Heterotrichida Sporadotrichida
Volvocaceae Stentoridae Oxytrichidae Codonellidae Tintinnidae Urocentridae Parameciidae Tetrahymenidae Zoothamniidae Operculariidae Vorticellidae Opalinidae
Bacillariophyceae
(=Diatomeae)
Stylonychia Tetrahymena Eutintinnus Tintinnopsis Opalina Diatom Swimming ciliate
Sessile ciliate
Sinking Diatom
A
Uswim g
F cilia
Z
X
Y
B
-Uswim
Lab/Body frame Lab frame Body frame Body frame Lab frame
Epistylididae
Usink
-Usink
Figure 5.1: (A) Phylogenetic tree showing microorganisms known to feature cilia that generate feeding
currents in either attached (blue) or free swimming (purple) states. The class of diatoms – non-swimming
cells that sink when experiencing nutrient limitation – is shown for comparison. (B) Flow fields around a
attached ciliate, swimming ciliate, and sinking diatom, in lab and body frame of references. Streamlines
are shown in blue in lab frame (X, Y, Z).
near the ciliary crown. For swimming ciliates and sinking diatoms, we collected measured swimming and
sinking speeds, which, given the no-slip boundary condition in this viscous regime [112, 71], also represent
flow speeds near the surface of these microorganisms.
Microcurrents and Péclet number The microcurrents generated by these ciliates improve solute transport to and from the surface of the microorganism. For a characteristic microcurrent of speed U =
100 µm · s
−1
, small molecules and particles would be transported over a characteristic distance a = 100
µm in approximately a/U = 1 s. In contrast, the same substance transported by diffusion alone takes a
considerably longer time to traverse the same distance. For example, diffusive transport of oxygen and
small molecules, with diffusivities that are in the order of D = 10−9 m2
· s
−1
, takes about a
2/D = 10 s,
while live and dead bacterial particles with respective diffusivity D = 4×10−10 m2
·s
−1
and D = 2×10−13
m2
· s
−1
[13] take about a
2/D = 25 s and 10000 s, respectively. The ratio of diffusive a
2/D to advective a/U timescales defines the Péclet number, Pe = aU/D. For Pe ≪ 1, mass transport is controlled by
38



Table 5.1: Survey of size a and flow measurements U in attached and swimming ciliates and
sinking diatoms. Size a is calculated using the volume-equivalent spherical radius. Corresponding ranges
of Pe numbers are based on the diffusivity of oxygen, D = 10−9 m2
·s
−1
, live bacteria, D = 4×10−10
,
m2
·s
−1
), and dead bacteria D = 2×10−13 m2
·s
−1
.
Empirical measurements Péclet number, Pe
microorganism characteristic live dead
size speed oxygen bacteria bacteria
a (µm) U (µm·s
−1
) diffusivity diffusivity diffusivity
Attached ciliates 15–60 50–2500 1–80 2–210 (5–400)×103
Swimming ciliates 15–180 50–3200 3–160 8–390 (17–800)×103
Sinking diatoms 10–120 40–210 0.4–23 – –
molecular diffusion. For the microorganisms that we surveyed, we obtain Pe ranging from nearly 0 to as
large as 103 − 105 depending on the nutrient diffusivity (Table 5.1). This dimensional analysis suggests
that the flows generated by the microorganisms substantially enhance the transport of solutes to and from
their surface, and while it clearly shows that diatoms typically occupy a smaller range of Pe numbers, this
analysis does not reveal a clear distinction between attached and swimming ciliates. To further explore
such distinction, and to assess whether ciliates are disadvantaged by flow physics in their attached state
compared to their swimming state as suggested in [5], cross comparison among mathematical models are
shown in this chapter and in [78].
5.2 Mathematical models
To quantify and compare nutrient uptake across microorganisms, the cell body is approximated as a sphere
of radius a, as typically done in modeling attached and swimming ciliates [19, 84, 90, 5] and sinking
diatoms [97, 63, 62].
Fluid model: envelope model The fluid velocity u around the sphere is governed by the incompressible
Stokes equations, −∇p + η∇2u = 0 and ∇ · u = 0, where p is the pressure field and η is viscosity. The
Stokes equations can be solved for two models of cilia activity: cilia represented as a Stokeslet force Fcilia
39



placed at a distance L and pointing towards the center of the sphere and no-slip velocity at the spherical
surface [18, 66, 140, 5] (Fig. 5.2A), and densely-packed cilia defining an envelope model with a slip velocity
at the spherical surface where all cilia exert tangential forces pointing from one end of the sphere to the
opposite end [19, 89, 90] (Fig. 5.2B). For envelope model, consider the "treadmill" slip velocity at the surface
of the sphere for both attached and swimming ciliate
uθ|
r=a = B
p
1 − µ2, ur|
r=a = 0, (5.1)
where B is a constant parameter that defines a velocity scale B = U = 1. In Table 5.2, the mathematical
expressions are listed for the boundary conditions, fluid velocity field, pressure field, forces acting on the
sphere, hydrodynamic power, and swimming speed for freely moving ciliated sphere. Results are consistent
with the results of [89] for a swimming sphere.
Fluid model: Stokeslet model Following [5, 30], we can represent the ciliary activity as a point force
F = −Fciliaez located at Lez, at a distance L from the center of the sphere. For the swimming sphere, the
flow field is obtained by superposition of two solutions of the Stokes equation, to obtain a solution with
no slip boundary condition at the sphere’s surface: u = u1 + u2, where u1 = Go · F and Go is Green’s
tensor due to a point force near a rigid sphere [18, 66] and u2 is the solution due to uniform flow −Uez
pass a rigid sphere, expressed in the sphere’s frame of reference [66, 5].
In the Stokes regime, force balance on the swimming sphere is given by
0 = T + K + D, (5.2)
where T = −F is the thrust force generated by the flagella (Stokeslet) in the swimming direction ez, D is
the drag force due to a moving sphere in fluid with speed U, and K is the hydrodynamic force acting on
40



Table 5.2: Summary of mathematical expressions for rigid and ciliated sphere models. Envelope
model subject to treadmill slip velocity. Mathematical expressions of boundary conditions, fluid velocity
field, pressure field, forces acting on the sphere, hydrodynamic power, and speed for each representative
case: attached ciliated sphere, freely swimming ciliated sphere, and sinking (non-ciliated) sphere. All
quantities are given in dimensional form in terms of the radial distance r and angular variable µ = cos θ.
B.C. at the surface of the sphere B.C. at infinity
Attached ciliated sphere uθ|
r=a = B
p
1 − µ2, ur|
r=a = 0 u|
r→∞ = 0
Swimming ciliated sphere uθ|
r=a = B
p
1 − µ2, ur|
r=a = 0 u|
r→∞ = −Uez
Sinking (non-ciliated) sphere u|
r=a = 0 u|
r→∞ = −Uey
Fluid velocity field
Attached ur(r, µ) =
a
3
r
3
−
a
r

Bµ, uθ(r, µ) = 1
2

a
3
r
3
+
a
r

B
p
1 − µ2
Swimming ur(r, µ) =
−
2
3
+
2a
3
3r
3

Bµ, uθ(r, µ) =
2
3
+
a
3
3r
3

B
p
1 − µ2,
Sinking ur(r, µ) =
−1 +
3a
2r
−
a
3
2r
3

Uµ, uθ(r, µ) =
1 −
3a
4r
−
a
3
4r
3

U
p
1 − µ2
Fluid velocity field in lab frame Far-field signature
Attached same as above Force monopole (u ∼ 1/r)
(Stokeslet)
Swimming ur(r, µ) = 2a
3
3r
3
Bµ, uθ(r, µ) = a
3
3r
3
B
p
1 − µ2 Potential dipole (u ∼ 1/r3
)
Sinking ur(r, µ) =
3a
2r
−
a
3
2r
3

Uµ, Force monopole (u ∼ 1/r)
uθ(r, µ) =
−
3a
4r
−
a
3
4r
3

U
p
1 − µ2 (Stokeslet)
Fluid pressure field Forces on sphere
Attached p(r, µ) = p∞ − η
a
r
2
Bµ F = 4πηaBez
Swimming p(r, µ) = p∞ − η
a
r
2
Bµ F = 4πηaBez, D = −6πηaUez
Sinking p(r, µ) = p∞ + η
3a
2r
2
Uµ W = −
4
3
πa3
(δρ)gey, D = 6πηaUey
41



Table 5.3: (table extension).
Hydrodynamic Power Swimming Speed
Attached P = 8πaηB2 U = 0
Swimming P = 8πaη
2
3
B2 U =
2
3
B
Sinking P = 6πaηU2 U =
2ga2
(δρ)
9η
the sphere by the flow generated by the point force F. Evaluating these terms leads to a scalar equation
in the ez-direction that can be solved to obtain the swimming speed
U =
Fcilia
6πµa
1 −
3a
2L
+
a
3
2L3

. (5.3)
For a attached sphere, The force balance on the sphere is given by:
0 = Ftether + T + K, (5.4)
where T = −F and K are defined as above. This equation can be solved to obtain the force Ftether
provided by a tether to prevent the sphere from moving in the ez direction. Namely,
Ftether =

1 −
3a
2L
+
a
3
2L3

Fciliaez, (5.5)
To non-dimensional these solutions in Eq.(5.3) and Eq.(5.5), we chose a as the characteristic length
scale and U = Fcilia/(8πµa) as the characteristic velocity scale. Results are consistent with the results
of [5].
Nutrients transport model To systematically compare the nutrient transport of model organisms, we
can consider a uniform concentration background C|r→∞ = Co. For non-dimensionalization, we can
scale the concentration using c = (C∞ − C)/C∞ assuming zero concentration at an infinite distance and
42



uniform concentration on the sphere’s surface. By choosing characteristic velocity scale U and length scale
a, we obtain the dimensionless governing equation and boundary conditions Eq.(4.6). The concentration
field can be solved numerically by using the Legendre spectral method, see Appendix. B.2.
Sherwood number and clearance rate To assess the effects of these cilia-generated flows on the transport of nutrients to the cell surface, we can use two common metrics of feeding.
For diffusion coupled with advective microcurrents, the nutrient uptake is quantified by the Sherwood
number Sh, which is equivalent to a dimensionless nutrient uptake, where I is scaled by Idiffusion. We can
apply both metric to each of the Stokeslet and envelope models.
The clearance rate is calculated by integrating the z-component of the flow velocity uz past an annular
disk around the equator of the sphere. Normalizing the result by a flux πR2U of a uniform flow U through
a disk of area πR2
, we get the clearance rate (with minus sign to indicate positive clearance value)
Q = −
2
R2U
Z R
a
uz|z=0rdr. (5.6)
When applied to the Stokeslet model, this clearance rate is consistent with that in [5]. By applying it to the
fluid solution in the envelope model, we can compare feeding rates across the two models ciliary models
which reflect different distributions of the ciliary force – concentrated in the Stokeslet model versus fully
distributed in the envelope model.
Asymptotic results For analyzing the advection-diffusion equation under extreme Péclet number Pe ≪
1 and Pe ≫ 1, we can use the approach employed in [3, 2] for a rigid sphere and the same approach in
Chapter.4 for the spherical envelope model. The results are summarized in Table.5.4.
43



Table 5.4: Summary of asymptotic results for Sh number of rigid and ciliated spheres. Expressions
for Sh as a function of Pe for attached and swimming ciliated sphere model, compared to a sinking sphere.
Large Pe limit Small Pe limit
Sherwood
number
Reference Sherwood
number
Reference
attached Sh = √
2
3π
Pe
1
2 present study Sh = 1 + 43
720Pe2
present study
Swimming Sh = √
2
3π
Pe
1
2 [84], [90] Sh = 1 + 1
3
Pe [84], [90]
Sinking Sh = 0.55Pe
1
3 [3, 2, 50] Sh = 1 + 1
3
Pe [3, 50]
5.3 Numerical results
In Fig. 5.2, the Stokeslet model is solved by considering the Stokeslet force Fcilia placed distance L = 2
(Fig. 5.2A), and envelope model with a slip velocity B = 1, U = 2/3 (Fig. 5.2B). In Fig. 5.2C,D, flow
streamlines (white) and concentration fields (colormap at Pe = 100) are shown in the Stokeslet and envelope models. In the attached sphere, ciliary flows drive fresh nutrient concentration from the far-field
towards the ciliated surface. These fresh nutrients thin the concentration boundary layer at the leading
surface of the sphere, where typically the feeding apparatus is found in sessile ciliates, with a trailing plume
or "tail" of nutrient depletion. Similar concentration fields are obtained in the swimming case, albeit with
narrower trailing plume.
In Fig. 5.2E, it is shown the clearance rate Q in the context of the Stokeslet model as a function of
the ciliary force location L/a for a small annular encounter zone of radius R = 1.1a extending away
from the cell surface. Swimming is always more beneficial. However, the increase in clearance rate due to
swimming is less than 10%. This is in contrast to the several fold advantage obtained in [5] for L = 4a and
R = 10a [78]. By employing the same metric Q in the envelope model, it can be found that motility is also
more advantageous, albeit at less than 5% benefit (Fig. 5.2H). The choice of the size of the encounter zone
are considered as following. Nutrient encounter and feeding in ciliates occur near the leading edge of the
44



ciliary band [44, 130, 57, 58]. Cilia are typically of the order of 10 microns in length, and the cell body of
a ciliate is typically in the range of 10-1000 microns. The choice R = 1.1a indicates encounter within an
annular protrusion extending 10% beyond the body radius because it falls within the biological range and
because a larger encounter zone would induce additional drag on the body that needs to be accounted for
in the model.
In Fig. 5.2F and G, the Sh number based on the Stokelet and envelope models are reported respectively.
In the Stokeslet model (Fig. 5.2F), sessile spheres do better when the cilia force is close to the cell surface
(L−a)/a ⪅ 1.25. In the envelope model (Fig. 5.2G), motile spheres do slightly better for all Pe⪅ 100. The
difference ∆Sh between the sessile and motile spheres favors, by less than 20%, the sessile strategy in the
Stokeslet model and the swimming strategy in the envelope model (Fig. 5.2H).
Comparing Sh between the Stokeslet and envelope models (Fig. 5.2C and D), it can be shown that,
at Pe = 100, Sh = 2.7 (sessile) and 2.6 (motile) in the Stokeslet model compared to Sh = 6.7 (sessile) and
6.9 (motile) in the envelope model. This is over a two-fold enhancement in nutrient uptake at the same
swimming speed U = 2/3 simply by distributing the ciliary force over the entire surface of the cell!
Indeed, this improvement occurs because the ciliary motion in the envelope model significantly thins the
concentration boundary layer along the entire cell surface as opposed to only near where the cilia force is
concentrated in the Stokeslet model.
5.4 Linking model prediction to biological data
Consider a wider range of Pe number Pe ∈ [0, 1000] for the attached and swimming spheres, and a sinking
sphere, the numerical predictions are shown in Fig. 5.3A ( solid lines, log-log scale). It is shown that at small
Pe, swimming is more advantageous than attachment; in fact, any motion, even sinking, is better than no
motion at all [127]. However, at larger Pe, there is no distinction in Sh number between the attached and
swimming sphere.
45



H
Stokeslet
Sessile
Envelope
Motile
Motile
Sessile
1
2
3
Sherwood
number, Sh
4
40
1.1 2 3
%
Clearance
rate, Q /πR U2
D
Sh = 2.7 Sh = 2.6
Sh = 6.7 Sh = 6.9
Peclet number, Pe
0 50 100
A
a
L
Encounter
R zone
0
20
B
EC
a
F
1
3
5
7
Sherwood
number, Sh
ΔSh
Force location, L/a
G
ΔQ
πR
2U
%
Envelope
0
40
20
20
0
0.4
0.2
0.2
Envelope
L/a=1.1
L/a=2
L/a=1.1
L/a=2
4
0% 50% 100%
Concentration
level
Fluid flow
Figure 5.2: Stokeslet and Envelope models of attached and swimming ciliates. (A) Stokeslet model
where ciliary activity is represented by a Stokeslet force Fcilia is located at a distance (L−a)/a outside the
spherical cell surface with no-slip surface velocity. (B)Envelope model where cilia activity is distributed
over the entire cell surface with slip surface velocity. (C,D) Fluid streamlines (white) and nutrient concentration fields (colormap) in the attached and swimming cases. Here, L/a = 2, a = 1 and Fcilia is chosen to
generate a swimming speed U = 2/3 in the swimming case to ensure consistency with the envelope model.
(E) Nutrient uptake in attached and swimming Stokeslet-sphere model based on calculation of clearance
rate Q of a fluid volume passing through an annular disk of radius R/a = 1.1 and Sherwood number Sh.
In the latter, Pe is 100. (F) Nutrient uptake in attached and swimming envelope model based on calculation
of Sherwood number Sh as a function of Pe. (G) Difference in clearance rate ∆Q = Qswimming −Qattached
and Sherwood number ∆Sh = ∆I/Idiffusion = Shswimming − Shattached in the Stokeslet-sphere model for
L/a = 1.1 and L/a = 2 and in the envelope model. In both metrics, the difference is less than 20%: ∆Q
is less than 20% the advective flux πR2U and ∆I is less than 20% of the corresponding diffusive uptake
Idiffusion = 4πRDC∞. Shaded grey area denotes when the attached strategy is advantageous.
46



Next, by using the Pe numbers obtained from experimental measurements of attached (blue ⃝) and
swimming (purple △) ciliates and sinking diatoms (green □), we can compute the corresponding values
of Sh number (Fig. 5.3A). Sinking diatoms are characterized by smaller values of Sh number, whereas with
increasing Pe, the Sh values of attached ciliates span similar ranges as those of swimming ciliates.
In Fig. 5.3B and C, with superimposed asymptotic results, it can be shown that, at small Pe ≪ 1, the
Sh numbers for swimming and sinking spheres scale similarly with Pe (Sh ∼ Pe), whereas Sh scales worse
(Sh ∼ Pe2
) for the attached sphere. The thorough literature survey indicates no data points for attached
microorganisms in this limit. At large Pe ≫ 1, the Sh numbers of the attached and swimming spheres
scale similarly with Pe (Sh ∼ Pe
1
2 ), whereas the sinking sphere scales worse (Sh ∼ Pe
1
3 ). Similar scaling is
found in swimming ciliate models [84, 123, 90]. These results confirm that, hydrodynamically, attached and
swimming ciliates are equivalent in the limit of large Pe. When cilia generate strong feeding currents that
drive nutrients and particulates toward the cell body, attached microorganisms can be equally competitive
with swimming microorganisms that swim to feed.
5.5 Summary
This study provides a new perspective on evaluating the role of flow physics in the feeding strategies of
microorganisms. Prior methods in support of either the swimming or attached strategies as optimal drew
general conclusions from focused analyses. Support for swimming came principally from flow-based models of idealized organisms propelling themselves through water [89, 68, 102, 5]. Support for maximum feeding by attached protists came from empirical measurements of prey removal by swimming versus attached
individuals [61, 28]. The presented approach shows that, while feeding rates may vary between organisms
and mathematical models, given a cellular (ciliary) machinery that allows microorganisms to manipulate
the surrounding fluid and generate flows, flow physics itself imposes no constraint on what is achievable by
the swimming versus attached strategies – they can both be equally competitive in transporting nutrients
47



100
102
10-3
10-2
10-1
100
101
101
Sherwood number, (Sh - 1)
Péclet number, Pe
Sh
Péclet number, Pe
Sh - 1
10-1 100
10-3
10-2
10-1
= 2
= 1
1=
Small Péclet limit
A B
Swimming
Sessile
Sinking
oxygen
Live
bacteria
C
103
10-1
Pe = 1
Coscinodiscus.
Porosira.
Palmerina.
Paramecium.
Stylonychia.
Urocentrum.
Tillina.
Tetrahymena.
Opalina.
Vorticella.
Epistylis.
Opercularia.
Pseudovorticella.
Stentor. 101
102
103
100
101
=
1
3
=
1
2
1
2
Large Péclet limit
Figure 5.3: Comparison of Sh number versus Pe number for ciliates and diatom models. for the
sinking (green) diatom and the swimming (purple) and attached (blue) ciliates based on the envelope model.
(A) Shifted Sherwood number (Sh - 1) versus Péclet number in the logarithmic scale for a range of Pe from
0 to 1000. Pe numbers associated with experimental observations of diatoms (square), swimming ciliates
(triangle), and attached ciliates (circle) are superimposed. Corresponding Sh numbers are calculated based
on the mathematical model. Empty symbols are for oxygen diffusivity D = 1 × 10−9m2
· s
−1
and the
solid symbols correspond to the diffusivity D = 4 × 10−10m2
·s
−1 of live bacteria [13]. (B-C) Asymptotic
analysis (dashed lines) of Sherwood number in the large Péclet limit (B) and small Péclet limit (C).
and wastes to and from the cell surface in the large Pe limit where nutrient advection is dominant. These
findings suggest that the choice of feeding strategy was likely influenced by evolutionary, ecological or behavioral variables other than flow physics, such as metabolic or sensory requirements [16, 94, 20], predator
avoidance [111], symbiotic relations [62], and nutrient availability or environmental turbulence [111, 72,
70].
48



Chapter 6
Flow Physics of Nutrient Transport Drives Functional Design of Ciliates
Nutrient transport is essential for cellular function and growth. Early life forms relied solely on diffusion
to absorb dissolved nutrients, a constraint that limited cell size and complexity [14, 123]. The evolution
of phagotrophy—the ability to ingest particulate food—marked a major transition in life’s history [26, 74].
This transition coincided with the development of both internal transport systems and external structures
specialized for food capture [95, 22, 35, 22]. Among modern single-celled eukaryotes, ciliates exemplify
this integration, using surface cilia to drive fluid flow and specialized organelles to process ingested particles [45].
Despite the shared molecular machinery underlying both intracellular and extracellular transport, the
organization of cilia varies widely across species. Swimming ciliates typically exhibit globally distributed
cilia that propel them through the water, while sessile ciliates arrange ciliary bands around the oral apparatus to create localized feeding currents. This structural divergence raises a key question: how does cilia
arrangement affect feeding performance, and are certain configurations functionally optimal depending
on the organism’s lifestyle?
In this chapter, we address this question through a systematic modeling study of ciliate morphologies
and flow-driven nutrient uptake. By quantifying the relationship between form, flow, and function, we
49



investigate how physical principles constrain and shape the design of these microorganisms (see also [77]).
6.1 Biological Observations
Phagotrophy Phagotrophy was accompanied by development of specialized mechanisms to both the
cell exterior - for food particle collection - and to the cell interior - for material transport within the
cell [74, 95]. These external and internal features fundamentally distinguished single-celled eukaryotes
from earlier archaea and bacteria [35], allowing them to attain much larger size and structural complexity.
In turn, these changes led to the emergence of diverse unicellular predation strategies [74], the evolution
of multicellularity [22], and the development of the modern food webs that determine Earth’s habitability [49, 141, 62]. However, despite its key role in the evolution of life, the link between intracellular and
extracellular transport in single-celled eukaryotes remains largely unexplored.
Intracellular and Extracellular Transport in Ciliates Ciliates are important aquatic protistan grazers, characterized by the presence of hair-like organelles called cilia. Cilia protrude from the cell surface
and act as microactuators that move fluid across the cell surface to enable the cell to swim or generate feeding currents (Fig. 6.1) [45, 62]. The ciliary activity directs fluid towards the cell "mouth" [116] or cytostome,
where food particles are ingested and subsequently enclosed in membrane-bound food vacuoles [100]. The
vacuoles are then transported along complex microtubular highways, powered by molecular motors, to
various parts of the cell [95], at rates exceeding diffusive transport over larger distances. Transport destinations include the mitochondria where energy generation occurs. A portion of this energy is, in turn,
allocated to power the cilia. Interestingly, while intracellular transport in eukaryotes is a highly organized
process that is conserved across cell types and evolutionary lineages [95], the arrangements of cilia in cells
that swim are markedly different from those in cells that attach and generate feeding currents. The latter
50



100 μm
cytosome
food
vacuoles
A B food
particulates
ciliary
bands
ciliary
area food
vacuoles
cytosome
Figure 6.1: Ciliates utilize phagotrophy to ingest food particles. An illustration of A. swimming
Paramecium and B. attached Stentor, where cilia beat in concert to direct food particles and prey to the cell
feeding apparatus. The Paramecium cell surface is densely covered by cilia, while the Stentor generates
feeding currents using ciliary bands that are localized around the feeding apparatus.
exhibit ciliary bands that are localized around the feeding apparatus (Fig. 6.1B). This raises the following
question: do optimal cilia arrangements vary according to the life strategy of the cell, swimming versus
attached?
Energy cost of enhancing feeding by moving fluid Enhancing nutrient transport to the cell surface
through fluid movement – whether by swimming or attaching and creating feeding currents – depends on
the Péclet number Pe = aU/D. Transport improves when the fluid motion over a length a occurs on a time
scale a/U that is significantly shorter than the diffusion time scale a
2/D, where D is nutrient diffusivity.
According to Pe, this benefit increases linearly with cell size, but it actually depends more strongly on cell
size than suggested by Pe alone. This becomes apparent when considering the cost of energy dissipated in
the viscous fluid, an energy that the cell itself must supply [14]. The rate of dissipation energy P per cell
volume V required to improve the cell’s nutrient uptake by fluid movement scales with ηD2/a4
, while that
required to propel the same cell at speed U scales with ηU2/a2
(see Fig.6.2), where η is the fluid viscosity.
The ratio defines a relative cost C ∼ Pe−2 of improving nutrient uptake through fluid movement: the cost
decreases as an inverse square law with increasing Pe and thus cell size a. Compared to a bacterial cell
of size a = 1 µm swimming at speed U = 10a, for a ciliate of size a = 100 µm, the cost of improving
nutrient transport by moving fluid reduces by a factor of 104
.
51



C
Average time, [s]
1
2
100
101
102
103
10-2
100
102
104
Cell radius, a [µm]
Kinesin
Myosin
Dynein
Diffusivity, [µm2/s]
0.4 - 9.8 Molecular
diffusion
649 - 711
836
Average speed, [nm/s]
Kinesin
Myosin
Dynein
241 - 381
422
Average time, [s]
0 1000
0
1
2
x104
500
3D Brownian
motion
Molecular motor
transport
Cell radius, a [µm]
B
E
improving nutrient uptake by
moving fluid
propelling cell at constant
speed U
Power per cell volume,
P/V [picoWatt/μm3]
P/V ~ ηD
2
a
-4
P/V ~ ηU
2
a
-2
10-10
10-6
10-2
10-14 10
-8
10
-4
10
0
10
4
C ~ Pe-2
Cost, C
U = 10a
10 x
-4
10 x
A
U = 10a
100
101
102
103
100
101
102
103
D
Figure 6.2: External and internal cost for nutrients transport for small and large cells. External(Top
row): Estimates of viscous dissipation power per unit volume needed for a cell of radius a to double its
nutrient uptake, compared to diffusion only, through either swimming or through stirring the surrounding
fluid (purple lines) [14]. Estimates of viscous dissipation power needed to propel the same cell at speed U =
10a (dashed line), with grey shaded areas highlighting the range of speed from U = a to U = 40a. (A) For
small cell a ≤ 5µm, the power needed for improving nutrient uptake through fluid movement costs more
than simply propelling the cell; for larger cells, a between 100 and 1000 µm, the power requirements for
improving nutrient uptake through fluid movement are smaller than those of moving the cell. (B)The cost
of improving nutrient uptake through fluid movement decreases as an inverse square law with increasing
cell size, making it an efficient strategy for larger cells. Internal(bottom row): (C,D) Average time for a
microscopic particle to move a distance a equal to the cell radius through diffusion in the 3D cytoplasm
(grey region) and through protein-based cytoskeleton transport (pink lines) shown in (C) regular and (D)
logarithmic scales. The logarithmic scale reveals the standard scaling of average time with distance square
in diffusive transport and with distance in advective transport. (E) Ranges of parameter values of diffusivity
coefficient [101] and molecular motor speeds [88, 136, 8, 120] are obtained from a literature survey.
Essentially, as cell size increases from a few microns to several tens of microns, the energetic cost of
moving fluid for enhancing the cell nutrient uptake, via ciliary-driven feeding current, decreases dramatically, making it an energetically feasible tactic for both swimming and attached cells (see Fig.6.2). From an
energetic perspective, this justifies the investment that ciliates made in maintaining active cilia that drive
fluid across the cell surface. Meanwhile, phagotrophy and orchestrated intracellular material transport
52



favor designation of a specific portion of the cell surface as a ’mouth’ for efficient food intake. We thus
asked: what optimal arrangements of cilia and feeding areas maximize the cell nutrient uptake?
Morphology data I collected morphological data from a total of 31 ciliates. For each ciliate, the spherical
equivalent radius and the cell surface fractions allocated to feeding Af and cilia coverage Ac are estimated
from multiple images. In addition to directly measuring Ac from images, we also approximated its value by
considering that the bandwidth of the cilia ring is 10 − 20µm, which is twice the length of a single cilium
(5−10µm). A summary of the mean and standard deviation of these measurements is shown in Table.6.1.
Meanwhile, wherever available, I collected measurements of maximal flow speeds in sessile ciliates and
average swimming speeds in motile ciliates, and computed the associated Pe values using the average value
of the equivalent sphere radius, along with diffusivity coefficients for Oxygen and dead bacteria, shown in
Table.6.2.
6.2 Model modification
To model the various ’mouth’ and cilia coverages, we can extend the established model in previous chapters
and extend as following settings. With prohibiting overlap between cilia and mouth, we can define a
cilia covered a fraction Ac of the cell surface and a food intake surface fraction Af
, where, by design,
Ac + Af ≤ 100% (Fig. 6.3 B). And for systematic study, we can keep energy cost P and cell volume
V constant throughout. This spherical model bridges the two archetypes of motile and sessile cells in a
concise way and enables efficient exploration of the entire cell design space spanned by (Ac, Af).
Describing surface velocity induced by ciliary bands. To describe ciliary activity over a ciliary band
covering a fraction of the spherical surface, specified by r = a = 1 and µ ∈ [µ2, µ1], where µ2 and µ1 are
53



Table 6.1: Summary of morphological data of surveyed ciliates. Mean ⟨·⟩ and standard deviation △·
of the ciliates’ size, feeding, and cilia area percentages.
Equivalent sphere radius feeding fraction Cilia fraction
⟨a⟩(µm) △a(µm) ⟨Af⟩ △Af ⟨Ac⟩ △Ac
Paramecium 34.36 7.57 0.04 0.01 0.96 0.01
Tetrahymena 17.82 3.72 0.07 0.05 0.93 0.05
Strombidium 21.44 6.35 0.17 0.02 0.26 0.09
Didinium 65.32 4.5 0.16 0.01 0.19 0.06
Cinetochilum 11.76 1.72 0.07 0.020 0.93 -
Merseres 24.97 0.71 0.13 0.007 0.10 -
Halteria 14.42 0.42 0.19 0.064 0.15 -
Stylonychia 27.51 - 0.15 - 0.06 -
Zosterodasys 48.43 - 0.02 - 0.98 -
Chilodontopsis 29.92 - 0.03 - 0.97 -
Spirostrombidium 21.61 0.22 0.15 0.006 0.26 0.09
Rimostrombidium 24.55 0.42 0.03 0.001 0.18 -
Trichototaxis 58.86 - 0.08 - 0.92 -
Uncinata 52.48 - 0.05 - 0.95 -
Lembadion 28.87 - 0.19 - 0.81 -
Blepharisma 15.14 - 0.08 - 0.92 -
Euplotes 32.37 11.13 0.11 0.011 0.24 0.08
Cyclidium 5.24 0.00 0.29 - 0.71 0.00
Stentor 149.49 111 0.24 0.05 0.08 0.03
Vorticella 29.09 5.2 0.15 0.01 0.18 0.12
Opercularia 21.53 1.07 0.11 0.01 0.17 0.13
Pseudovorticella 19.1 7.15 0.42 0.04 0.20 0.06
Epistylis 17.83 2.33 0.11 0.018 0.13 0.04
Rhabdostyla 15.53 1.73 0.15 0.017 0.18 0.06
Baikalonis 10.11 0.66 0.14 - 0.17 0.06
Zoothamnium 20.09 0.43 0.15 - 0.18 0.06
Paravorticella 8.56 - 0.45 - 0.27 0.09
Scyphidia 26.08 0.77 0.17 0.002 0.10 0.03
Opisthonecta 21.34 0.17 0.07 0.011 0.18 0.06
Pyxicola 19.47 1.47 0.06 0.009 0.19 0.06
Carchesium 15.72 - 0.15 - 0.18 0.06
54



Table 6.2: Flow rates measurements for surveyed ciliates. Flow speed U and clearance rates Q, where
the flow speed U for sessile ciliates is measured near the ciliary band in the absence of background flow,
while for motile ciliates, it corresponds to their swimming speed. The Péclet number is defined as Pe
= aU/D, where D1 = 1×10−9m2/s is the diffusivity of the oxygen molecule, and D2 = 4×10−10m2/s
is that of dead bacteria [12]. The clearance rate is the volume of water cleared per hour in unit body
volume per hour, obtained by measuring food particles ingestion rate and derived by dividing the ingestion
rate (number of particles ingested per unit time) by the average food particle concentration during the
experiment [38].
Flow speed Ref clearance Ref Pe Pe
rate (Oxygen) (dead bacteria)
U(µm/s) Q(V/h)
⟨a⟩U
D1
⟨a⟩U
D2
Paramecium 350- 2500 [24, 41] 2 × 104
[38] 12.1 - 87.0 30.2 - 217.4
Tetrahymena 180 - 500 [23, 138] 103
[38] 3.2 - 8.9 8.0 - 22.3
Strombidium 900 - 2800 [59] - - 19.3 - 60.0 48.2 - 150.1
Didinium 0.4 - 1500 [15] - - 0 - 98.0 0.1 - 244.9
Rimostrombidium 160 - 200 [15] - - 0 - 98.0 0.1 - 244.9
Stentor 160 - 360 [134] 104
[38] 23.9 - 53.8 59.8 - 134.5
Vorticella 320 - 2500 [126, 98, 109] - - 9.3 - 72.7 23.3 - 181.8
Opercularia 25 - 47 [143] - - 0.5 - 1.0 1.3 - 2.5
Pseudovorticella 100 - 500 [62] - - 1.9 - 9.5 4.8 - 23.9
Stylonychia - - 103
[38] - -
Lembadion - - 105
[38] - -
Blepharisma - - 104
[38] - -
Euplotes - - 2 × 104
[38] - -
Cyclidium - - 6 × 103
[38] - -
55



ex
ez
eθ
θ
a
er
cilia
coverage
cell
surface
feeding
apparatus
µ1
= cos(θ1
)
µ2
= cos(θ2
)
A
0
100
0 100
0 100
Occupied fraction
of cell surface, %
Feeding apparatus Af, %
Cilia coverage Ac
, %
50
50
B
Figure 6.3: Mathematical modeling of ciliated cells with fractional feeding and ciliary areas. (A)
Ciliated cell modeled as a sphere with a fraction of the spherical surface allocated for feeding Af and an
adjacent, non-overlapping fraction allocated for cilia coverage Ac. (B) Cell design space (Ac, Af) describes
the fractions of the cell surface dedicated, respectively, to the mouth and cilia coverage.
limit values shown in Fig.6.3.A, we can write u|r=a = αg(µ)eθ, where α is a constant and g(µ) is given
by
g(µ) =



sin
π
µ1 − µ
µ1 − µ2

, µ ∈ [µ2, µ1].
0, µ ∈ [−1, µ2) ∪ (µ1, 1].
(6.1)
By projecting αg(µ) on the velocity basis Vn(µ), we can obtain
Bn =
n(n + 1)(2n + 1)
8
Z 1
−1
αg(µ)Vn(µ)dµ. (6.2)
In practice, we consider a total of N velocity modes, n = 1, . . . , N, where N is chosen such that the error
is less than 10−4
. The constant α is determined by imposing the constraint of constant energy dissipation
rate
P = constant = 8πηa



P∞
n=1
2B2
n
n(n + 1), sessile cell,

2
3
B2
1 +
P∞
n=2
2B2
n
n(n + 1)
, motile cell.
(6.3)
Advection-diffusion equation To analyze the nutrient concentration around the sphere, I considered
the steady-state advection-diffusion equation for a concentration field C of nutrients with diffusivity coefficient D. We write the solution in non-dimensional form, choosing a and U =
p
P/(8πηa) as the
56



characteristic length and velocity scales, which we normalize to unity. The concentration field is normalized by the far-field background value C∞, with the dimensionless concentration c = (C∞ −C)/C∞. The
governing equation and boundary conditions are
Pe u · ∇c = ∆c, c|
r→∞ = 0,



c = 0, surface fraction designated as ’mouth’,
∇rc = 0, remaining cell surface.
(6.4)
Here, the Péclet number Pe = aU/D is defined based on the characteristic length scale a and the characteristic speed U.
Surface fraction designated as ’mouth’. Let the absorbing fraction of the surface be defined as
Af =
1
4πa2
Z 2π
0
Z θ1
0
a
2
sin θdθdϕ =
1
2
Z 1
cos θ1
d(cos θ) = 1 − µ1
2
, (6.5)
where the ’mouth’ area always starts from the pole of the spherical cell at θ = 0 and ends at θ1 (µ1 = 1),
shown in Fig.6.3. A. The adjacent ciliary band coverage near the ’mouth’ is in the range µ ∈ [µ1, µ2].
Feeding flux To estimate the feeding rate, we calculate the flux Q through a disk placed at a distance
0.1a away from the cell ’mouth’ and of radius Rd = a
p
(1 − µ
2
1
) (with rdrdϕ being an area element)
Q =
Z 2π
0
Z Rd
0
uz|
z=1.1a
rdrdϕ. (6.6)
Sherwood number By using Fick’s law, the nutrient flux at the sphere surface is equal to the area
integral of the concentration flux over the spherical surface,
I = −
Z Z
’mouth’
eˆr · (−D∇rC)dS =
Z θ2
0
D
∂C
∂r 2πa2
sin θdθ. (6.7)
57



In the case of pure diffusion, the inward flux Idiffusion due to molecular diffusion can be obtained by numerically solving the diffusion equation with the above chosen boundary conditions. Accounting for both
advective and diffusive transport, the Sherwood number Sh= I/Idiffusion is defined as dimensionless nutrient flux through the sphere surface.
Forward numerical approach The nutrient concentration field c is obtained by numerically solving the
governing equation Eq.(6.4) using the a finite difference method. To reach the far-field boundary condition,
we used a non-uniform radial mesh such that the grid is denser near the sphere surface and more sparse
in the far field. We chose to discretize θ on a uniform mesh. See [89, 90, 78, 79] for more details.
Optimization approach To search for optimal surface motions that maximize feeding in the ciliated
sphere model with partial oral fraction, I considered an optimization method based on variational analysis
and steepest descent [90, 79], see Appendix B.1. The problem consists of a PDE-constrained optimization
problem, where the goal is to find an optimal range of cilia coverage from µ1 to µ2 that maximizes the
nutrients flux, subject to the concentration field c satisfying the advection-diffusion equation and surface
velocity satisfying the constant energy constraint. The value of µ1 is dictated by the end of the ’mouth’
area to satisfy the adjacency and no overlap condition between ’mouth’ and cilia. The goal is thus to find
the optimal cilia endpoint µ2,
max
µ2,c
I(µ2, c), subject to L[c] = 0 and X
N
n=1
β
2
n = 1. (6.8)
Here, the linear operator L = Peu · ∇ − ∇2
is that of the advection-diffusion equation along with the
corresponding boundary conditions. The coefficients βn = Bn
p
2/(n(n + 1)) are weighted coefficients,
satisfying the constraint on the energy dissipation rate, P
n
β
2
n = 1. Note that for numerical solution, I
consider a finite number N of modes.
58



The goal is seeking variation of nutrient flux δI due to variation of the endpoint δµ2.
δI
δµ2
= −2πPe Z ∞
1
Z 1
−1
c

δu
δµ2
· ∇h

r
2
drdµ, (6.9)
where h(r, µ) is the adjoint function associated with the advection-diffusion equation, satisfying
Peu · ∇h + ∇2h = 0, h(r = 1) = 1, h(r → ∞) = 0, (6.10)
and the consistency equation
Z
S
∇δc · ndS = −Pe Z
V
cδu · ∇h dV. (6.11)
From functional analysis, we have
δu
δµ2
=
X
N
n=1
δur
δβn
δβn
δµ2

er +

δuθ
δβn
δβn
δµ2

eθ

,
=
X
N
n=1
δur
δBn
δBn
δµ2

er +

δuθ
δBn
δBn
δµ2

eθ

.
(6.12)
Using the Leibniz integral rule and recalling Eq.(6.1) and Eq.(6.2), the variation in velocity coefficients Bn
generated by the variation in ending point µ2 are given by
δBn
δµ2
=
n(n + 1)(2n + 1)
8

−g|µ=µ2 Vn(µ) + Z µ1
µ2
∂g
∂µ2
Vn(µ)dµ
=
n(n + 1)(2n + 1)
8
Z µ1
µ2
π(µ1 − µ)
(µ1 − µ2)
2
cos
π(µ1 − µ)
µ1 − µ2

Vn(µ)dµ.
(6.13)
By substituting Eq.(6.13) into Eq.(6.12), we obtain the derivative of the velocity field with respect to µ2.
At each iteration for updating µ2, we solve for the adjoint function h(r, µ) numerically by solving the
associated adjoint equation along with boundary conditions Eq.(6.10). By substituting h together with
59



Eq.(6.12) into Eq.(6.9), we obtain the expression for the gradient δI
δµ2
, which we can use to update the value
of I to identify µ2 that maximizes I. The step size s can be chosen dynamically at each iteration step j for
updating µ2,
µ
(j+1)
2 = µ
(j)
2 + s

δI
δµ2
(j)
, (6.14)
where s = (µ
(j)
2 − µ
(j−1)
2
)/(
δI
δµ2
(j)
−
δI
δµ2
(j−1)). This choice is based on the Barzilai-Borwein method to
accelerate the computation time.
6.3 Model predictions of optimal design
With two introduced measures of feeding performance: the feeding flux Q defined as the integrated fluid
velocity over a cross-sectional area leading to the cell’s designated mouth [68, 118, 117, 108, 109], and the
Sherwood number Sh, defined as the ratio of total nutrient flux directed to the mouth region divided by
the corresponding diffusive flux [63, 84, 123, 90, 91, 78, 79, 62], we seek for optimal configuration in design
space.
Case study: 10%, 30%, 50% ’mouth’ fraction Fixing Af = [10%, 30%, 50%], the optimal cilia arrangement results are computed with open-loop and closed-loop methods, shown in Fig.6.4. It is found
that, for Af = [30%, 50%], optimal feeding in sessile cells occurs when cilia are arranged in a narrow band
around the feeding apparatus, covering only a narrow fraction Ac of the cell surface, while in motile cells,
feeding is optimal when cilia covered the remaining cell surface Ac = 100 − Af
(also see in Fig. 6.5A,B,
open yellow markers). For Af = 10%, the optimal feeding occurs when cilia are arranged in a narrow band
around the feeding apparatus for both sessile and motile. These results aligned with calculations based on
a closed-loop optimization approach for maximizing Sh.
60



0
1
2
Surface velocity u
-1 10
0
1
2
0
1
2
-1 10
-1 10
initial
guess
Af
= 50%
Af
= 30%
Af
= 10%
sessile motile
sessile motile
sessile motile
µ = cos(θ)
motile
sessile
Surface velocity u Surface velocity u
Feeding flux, Q/πa
Sherwood number, Sh
2U
0
0.1
2
2.4
Af
= 10% 2.8
Feeding flux, Q/πa
Sherwood number, Sh
2U
0
0.1
2
2.4
Af
= 30% 2.8
Feeding flux, Q/πa
Sherwood number, Sh
2U
0
0.1
0
Cilia coverage Ac, %
45 90
2
2.4
Af
= 50% 2.8
0
Cilia coverage Ac, %
45 90
Af
= 10%
Af
= 30%
Af
= 50%
B C
Sessile
Motile
D
1
optimal
0
A
sessile motile
1
0
(Af = 30%)
C/C∞
u/U
(lab frame) (lab frame) (cell frame) (cell frame)
sessile motile
Figure 6.4: Flow and concentration fields of ciliated cells with vary feeding and ciliary areas.
(A) Flow field (lab frame) and concentration field (cell frame) corresponding to Af = 30% and Ac =
21.5% (sessile) and Af = 30% and Ac = 59.6% (motile), with Sh= 2.57 and Sh= 2.76, respectively.
(B) Feeding flux Q/πa2U and (C) Sherwood number Sh as a function of cilia fraction Ac for feeding
fractions Af = [10%, 30%, 50%]. Blue circles and purple squares represent sessile and motile ciliates,
respectively. Filled yellow markers indicate optimal values. (D) Closed-loop optimization to compute
optimal Ac that maximizes Sh; in fact, we compute optimal µ2 given µ1 corresponding to "mouth" size
Af = [10%, 30%, 50%]. All calculations are for Pe = 100.
Optimal configuration in design space Next, the feeding metrics Q and Sh can be computed for the
entire range of Af and Ac from 5% to 95% (Fig. 6.5A-D) and identified for each Af
, the cilia coverage Ac
that locally maximized nutrient uptake (Fig. 6.5, solid yellow lines). The value of Ac that maximizes Q (or
Sh) at a given choice of feeding fraction Af
is a local optimum in that it is tailored to that specific Af
. We
discovered that, regardless of the choice of performance metric, Q or Sh, the local optima of Ac depended
61



on cell motility: for Af ⪅ 20%, the local minima overlapped for sessile and motile cells, with ciliary bands
around the feeding fraction optimized feeding efficiency. At around Af ≈ 20%, the local optima in motile
cells bifurcated into a new branch where cilia covered the remaining cell surface not dedicated to feeding
(Fig. 6.5A,C), whereas in sessile cells, no such bifurcation occurred. The branch of local optima continued
smoothly (Fig. 6.5B,D). The combination (Ac, Af) that maximized Q (or Sh) over the entire design space is
a global optimum (Fig. 6.5, solid yellow markers). In sessile cells, the global optimum corresponds to ciliary
bands, while in motile cells, it features cilia spread over the entire non-oral cell surface.
A closer examination of the design space (Ac, Af) shows that at the global optima, the surface fraction
Af allocated to the feeding apparatus is less than 50% in both sessile and motile cells (Fig. 6.5, solid yellow
markers [77]). That is, unlike in diffusive transport, where maximal intake is achieved when receptors
are randomly distributed over the entire cell surface [14], in flow-driven transport, the designation of a
specific portion of the cell surface as a ‘mouth’ not only facilitates the intracellular processing of ingested
food particles, but also optimizes external transport to the cell surface. The Q- and Sh-landscapes in both
sessile and swimming cells are broad peaked. This feature implies that nutrient uptake is robust to small
variations in cell design and insensitive to perturbations in the specific feeding and ciliary arrangements.
To underscore the significance of this feature, we calculated the contour lines corresponding to Q =
104 bodyvolume/h and highlighted the regions of the cell morphospace (Ac, Af), where Q exceeding 104
bodyvolume/h. We choose this Q value because empirical observations in ciliates range from Q = 103
to
105 bodyvolume/h [38]. To map our dimensionless model results to units of bodyvolume per hour, we used
a maximal speed of 50 bodylength/s, consistent with empirical measurements in ciliates (Supplementary
Table S.5, Fig. 6.5E, inset) [27]. Our model predicts clearance rates higher than Q = 104 bodyvolume/h
over broad regions of (Ac, Af) space (Fig. 6.5A,B inset). Likewise, for Sh number, over twofold increase in
nutrient uptake is attainable nearly everywhere in the cell morphospace (Fig. 6.5C,D inset).
62



Multi-objective optimal design Because ciliary activity in motile cells serves to also power the cell
motion, we explored how the arrangements of feeding and cilia coverage affect the cell swimming speed U.
Not surprisingly, U is maximum at maximal cilia coverage (Fig. 6.5E). In sessile cells, ciliary activity exerts
a force on the surrounded fluid that should be balanced by an attachment force, sustained by a specialized
stalk [119, 79, 78]; we thus computed the attachment force over the entire cell morphospace and found
that this force is minimized at minimal cilia coverage (Fig. 6.5F). That is, in motile cells, optimal designs
that maximize Q and Sh by distributing cilia across the entire surface fraction unoccupied by the feeding
apparatus also improve the propulsive force and swimming speed U of the cell (Fig. 6.5E), while in sessile
cells, optimal designs that maximal Q and Sh by concentrating ciliary activity in a band around the cell’s
feeding area also reduce the force needed to keep the cell attached (Fig. 6.5F).
To rigorously account for these multiple objectives, we defined a multi-objective optimization problem
using the weighted-sum method, where the objective function J depended on cell motility: in motile
cells, J = αSh + (1 − α)U aimed to maximize both feeding rate Sh and swimming speed U, whereas
in sessile cells, J = αSh − (1 − α)F aimed to maximize feeding while minimizing the attachment force
F. The parameter α describes the weight given to nutrient uptake; at α = 1, J account for only feeding,
recovering the optimization problem in Fig. 6.5C,D. At α = 0, feeding does not enter in the optimization.
As we varied α from 0 to 1, the optimal solutions describe a Pareto front [86]. In motile cells, the Pareto
front starts at maximal cilia coverage at α = 0, which maximizes swimming speed, and ends at maximal
nutrient uptake at α = 1, thus favoring maximal cilia coverage over the entire cell surface not allocated to
feeding (Fig. 6.6A). In sessile cells, the Pareto front jumps abruptly from minimal cilia coverage at α = 0,
when feeding is discounted, to the branch of optimal solutions characterized by a localized ciliary band
surrounding the feeding area, ending with maximal nutrient uptake at α = 1 (Fig. 6.6B). That is, the
Pareto front in sessile cells follows closely the Sh-optimal branch, indicating no conflict between the two
objectives of maximizing uptake and minimizing attachment force.
63



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
1.8
Cilia coverage Ac, %
50
Sh
2.8
D
Sessile Motile
EC
F
.34
0
1.5
Sh
1
0
.18 2.6
0
a2U
0 100 0 100 50
F
4̟ηa
0
50
100
Feeding apparatus Af, %
.7
0
0.2
0.4
0.6
0.8
0.6
Swimming speed
Attachment force
0.2
0
50
100
Feeding apparatus Af, %
Sh>2
Sh>2 Q > 104
B
A
Cilia coverage Ac, %
50
Q > 104
0 100
Cilia coverage Ac, %
bodyvolume
hour
bodyvolume
hour
0.4
Q
a2U
Q
U
0 100 50 0 100 50 0 100 50
0
50
100
0
50
100
0
50
100
0
50
100 bodylength
10
30
second
Figure 6.5: Optimal arrangements of cilia and feeding structures in sessile and motile cells. A,
B. Feeding flux Q and C,D. Sherwood number Sh as colormaps over the morphospace (Ac, Af) of motile
(top) and sessile (bottom) cells; darker colors indicate larger values, white contour lines indicate designs
that produce same Q or Sh values, solid yellow lines indicate local optima at each Af
; filled circles and
squares indicate global optima of sessile and motile designs over the entire design space. Insets in A and
B: contour lines of Q = 104
in units of bodyvolume per hour for motile (purple) and sessile (blue); Insets
in C and D: contour lines of Sh = 2 for motile (purple) and sessile (blue) shows that Sh is more than
doubled nearly everywhere in the morphospace. E. Cell swimming speed U increases with increasing cilia
coverage. Inset: swimming speed in units of bodylength per second, highlighting contour lines of 10 and
30 bodylength per second, with maximal value of 50 bodylength per second in the rightmost corner. F.
Force required to keep the cell attached increases with cilia coverage. All calculations are done at constant
hydrodynamic power P, cell volume V, and Pe = 100 chosen based on typical particle diffusivity and flow
speed measurements in ciliates.
64



6.4 Linking model prediction and biological data
How do these model predictions help elucidate the diversity in cilia arrangements observed in sessile and
swimming ciliates? To answer this question, we surveyed 31 genera of ciliates within the phylum Ciliophora, encompassing the classes Oligohymenophorea, Heterotrichea, Litostomatea, Spirotrichea, and
Phyllopharyngea listed in Fig. 6.6 and Supplementary Table S.3. Ciliates that attach and generate feeding currents, such as Stentor [121, 134, 124, 113], Vorticella [126, 98, 109], Opercularia [143, 131, 114] and
Pseudovorticella [62, 60, 47], feature ciliary bands localized around the feeding apparatus. According to
our model, this design maximizes feeding efficiency while minimizing the required attachment force. Of
the swimming ciliates, Didinium [55, 11] and Strombidium [96, 87, 59] feature localized ciliary bands while
the entire non-oral cell surface of Paramecium [24, 41] and Tetrahymena [23, 40, 138, 10] is covered by
cilia. Mapping these empirical observations onto our morphospace (Ac, Af), we found that sessile ciliates
(blue and green circles) grouped into a single cluster, characterized by localized ciliary bands surrounding
the feeding apparatus, whereas motile ciliates (red and purple circles) formed two distinct clusters, one
characterized by localized ciliary bands and the other characterized by cilia covering the entire cell surface (Fig. 6.6A,B). The motile and sessile ciliates that feature localized ciliary bands followed closely the
branches of local optima of Q and Sh predicted by our model, which in sessile cells also coincide with the
Pareto front and lead to the global optimum of nutrient uptake. The motile ciliates, such as the Paramecium and Tetrahymena, closely followed the predictions of our Pareto front for motile cells. These results
indicate that cilia allocation is a balance between swimming and feeding.
Theoretically, for the same energetic cost, flow physics permits larger increase in nutrient uptake for
the swimming cells that follow the Pareto front but at the expense of reducing their swimming speed.
Larger increase in both nutrient uptake and swimming speed are possible for the swimming cells that
exhibit localized ciliary bands around their feeding apparatus, yet at the cost of having to build more cilia.
Such improvements are possible at no additional energetic cost for powering the cilia. This is evident
65



Pareto
front
α = 1
B α = 0
A
0
50
100 Feeding apparatus Af, %
Paramecium
Tetrahymena
Strombidium
Didinium
Cinetochilum
Meseres
Halteria
Stylonychia
Zosterodasys
Chilodontopsis
Spirostrombidium
Rimostrombidium
Trichototaxis
Uncinata
Lembadion
Blepharisma
Euplotes
Pseudovorticella
Opercularia
Vorticella
Stentor
Epistylis
Rhabdostyla
Baikalonis
Zoothamnium
Paravorticella
Scyphidia
Opisthonecta
Pyxicola
50 Carchesium
Cilia coverage Ac, %
0 100
Motile
Sessile
sessile
Sh-optimal: motile
sessile
Q-optimal: motile
D
C
2
2
0
Q = 104
Q = 104
Cilia coverage Ac, %
0 100 50
Surveyed ciliates
0
α = 0
α = 1
Pareto
front
0
50
100 Feeding apparatus Af, %
0
50
100 Feeding apparatus Af, %
0
50
100 Feeding apparatus Af, %
Pareto fronts:
sessile
motile
Sh-optimal
Sh-optimal
motile sessile
Figure 6.6: Optimal designs explain the diversity of cilia arrangements in sessile and swimming
ciliates. A. and B. Solid yellow lines are local optima at each Af
from Fig.??. Pareto fronts that maximize
J = (αSh + (1 − α)U) for swimming cells (purple solid line) and J = (αSh − (1 − α)F) for sessile
cells (blue solid line). Here, we normalized Sh, U and F to each lie in the range from 0 to 1, and for each
value of α ∈ [0, 1], we computed the corresponding optimal cell design (Ac, Af ) to obtain the Pareto front
as α increases from 0 to 1. We surveyed 31 motile (purple and red) and sessile (blue and green) ciliates,
listed by their genus (Supplementary Tables S.3 and S.4). Empirical measurements are mapped onto the
design space (Ac, Af): sessile ciliates cluster in one region of the morphospace, characterized by a ciliary
crown surrounding the feeding apparatus, while motile ciliates exhibit two types of cell morphologies,
one similar to the sessile ciliates and one that maximizes cilia coverage over the entire cell surface. C. All
surveyed ciliates lie in a region of the morphospace where the model predicts a feeding rate Q of over 104
bodyvolume per hour, consistent with experimental measurements (Supplementary Table S.5) (contour
lines of Q and optimal designs taken from Fig. 6.5A,B). D. Likewise, all surveyed ciliates correspond to
over twofold increase in nutrient uptake compared to diffusion alone at Pe = 100; (contour lines of Sh and
optimal designs taken from Fig. 6.5C,D). Superimposing optimal solutions, including the Pareto fronts,
show remarkable consistency between optimal model predictions and surveyed ciliates.
66



because our modeling results are all obtained at the same hydrodynamic power P and cell volume V. Why
don’t all swimming cells then cluster at the Pareto front that optimizes feeding and swimming?
To explore this question, we superimposed all surveyed ciliates, both motile and sessile, onto our model
predictions of nutrient flux Q and Sherwood number Sh (Fig. 6.6C,D). We found that all surveyed ciliates,
whether motile or sessile, with ciliary bands or covered by cilia, occupied regions of the morphospace
where nutrient uptake was more than 104 bodyvolume per hour, consistent with empirical observations;
in these regions, nutrient uptake was more than doubled compared to diffusion only. That is, swimming
and sessile ciliates achieve similar levels of nutrient uptake through different body plans, and beyond a
threshold of doubling nutrient uptake, further increases in Sh number do not seem to be a dominant selective force in cell design. Because of the coevolution of extracellular and intracellular transport mechanisms
in eukaryotic cells, nutrient acquisition does not depend solely upon external transport. Intracellular processing of transported nutrients may have played a critical role in cell adaptations. Such considerations
support the hypothesis that, unlike cells that rely solely on diffusion, the limiting factor for nutrient uptake in unicellular eukaryotes may not be the external transport of nutrients to the cell surface, but rather
the cell’s capacity to process these transported nutrients [56]. This is far from being the only plausible explanation, evolutionary processes other than biophysical optimality may have limited further increase in nutrient transport. This is particularly relevant for organisms like protists, which typically have
small effective population sizes, making them more influenced by evolutionary processes driven by genetic drift [81, 137, 54]. Another factor could be the energetic costs associated with the formation and
maintenance of the ciliary structures. Since the energetic cost of structural investments scale with cell
size, evolutionary adaptations must balance the benefits of increased nutrient uptake against the costs of
building and maintaining additional ciliary structures [82]. Or simply, a twofold increase in nutrient acquisition over diffusion might be sufficient for survival, and further increases might not confer a practical
selective advantage.
67



6.5 Summary
We established a quantitative map linking nutrient uptake in ciliates to the cell’s design or surface allocation to feeding structures and cilia activity. The evolutionary transition from small prokaryotic cells,
which relied on diffusion for nutrient transport [14, 141], to larger eukaryotic cells, capable of consuming
particulates and smaller cells via phagotrophy [74], entailed increases in cell volume of up to eight orders
of magnitude [115]. Sustaining these larger cells required expanding the mechanisms for nutrient transport to the cell surface beyond simple diffusion [74, 123]. We showed that the energetic cost of moving
fluid for enhancing nutrient uptake, for example by ciliary activity, decreases dramatically in larger cells,
as an inverse square law in Péclet number, making it an energetically viable strategy for attached and
swimming unicellular eukaryotes. The designation of a small surface area as a ’mouth’ enhances nutrient
transport to the cell surface, thanks to the cell’s ability to actively direct fluid to the designated mouth.
Importantly, the diversity of cilia arrangements relative to the oral area in ciliates is not accidental but,
instead, is shaped by flow physics. Sessile ciliates maximize prey encounter while also minimizing the
force necessary to keep the ciliate attached to the substrate by concentrating their cilia in a band around
their oral area. Swimming ciliates display more diverse designs: some have localized ciliary bands around
their oral regions, while others are entirely covered with cilia. At first glance, this diversity is puzzling
because of the view that optimal feeding corresponds to optimal swimming [90, 5, 78] and localized ciliary
bands are less effective for swimming. But the flow physics of nutrient transport at the microscopic scale
is more nuanced.
Our study uncovered an important, previously unrecognized, connection between the cell surface area
allocated to feeding and optimal cilia coverage, and showed a bifurcation in optimal designs of swimming
ciliates from localized ciliary bands to cells covered in cilia. These bifurcated designs achieve similar
benefits in nutrient uptake. The implications are important for swimming - ciliates with localized ciliary
bands may not reach the same swimming speeds as those covered in cilia, but both can achieve similar
68



levels of nutrient uptake. Indeed, surveyed sessile and swimming ciliates lay within the same range of
nutrient uptake, despite the theoretical potential for further increases of nutrient uptake in swimming
ciliates, at no additional energetic cost. These findings indicate that, unlike smaller prokaryotes [14], the
limiting factor for unicellular eukaryotes is not the energy required for generating feeding currents. The
limiting factor could be the cell’s feeding structures and their capacity to process transported nutrients [56]
or the investment in building and maintaining more cilia [82]. Alternatively, evolutionary processes other
than biophysical optimality might have been relevant [81, 137, 54], or simply that a twofold increase in
nutrient uptake is "good enough" for survival.
Adaptations to the surrounding fluids for nutritional transport by individual cells capable of phagotrophy, as we note here for ciliates, has shaped the lifestyles of these protists as free-living unicellular predators consuming other organisms [74] and enabled their proliferation as one of the most influential planktonic protistan lineages. These adaptations also formed the building blocks that led to more complex
collections of cells comprising colonies [121, 123, 127, 118]. Cooperation means sharing of resources
among individuals with different abilities [121], and is facilitated when having an abundance of these
resources [105, 83]. Our findings, which strongly suggest that ciliates have the energetic capacity to transport more nutrient-laden fluids than they are able to individually process, offer a fresh perspective and
opportunities for future research on the origins of cooperation and teamwork between unicellular protists [121, 118, 92, 123, 127, 62].
69



Chapter 7
Feeding in non-uniform concentration field
Marine microbe often inhabit environments with patchy or non-uniform nutrient distribution [39, 80, 33,
64, 129, 7, 128]. Motile microbes navigate toward nutrient-dense regions through chemotaxis, actively
seeking out localized nutrient hotspots. Sessile ciliates attach to substrates and generate feeding currents
to capture passing particulates and dissolved nutrients. Optimal ciliary activity that maximizes nutrient
flux at the cell surface while minimizing the rate of hydrodynamic energy dissipation is well characterized
in uniform nutrient fields [90, 91, 78, 79]. However, it is unclear how ciliary motion should change when
nutrients are non-uniform or patchy. To address this question, in this chapter, I will present the analysis
using a mathematical model that couples cilia-driven flows with unsteady nutrient concentration.
7.1 Unsteady diffusion
We can describe the concentration C(t, r, θ, ϕ) of nutrients around a spherical cell of radius a. Here, t
is time and (r, θ, ϕ) are spherical coordinates. Considering diffusion of nutrients at constant diffusivity
70



D [14] and starting from an initial concentration field C(0, r, θ, ϕ) = Cb, the unsteady diffusion equation
subject to zero concentration at the cell surface and proper convergence to Cb as r → ∞ is given by,
∂C
∂t = D∇2C,
C|t=0 = Cb, C|r=a = 0, C|r→∞ = Cb (or ∇C|r→∞ = ∇Cb).
(7.1)
It is shown in Chapter.3, if the background concentration is uniform Cb = Co, where Co is constant, the
concentration field evolves as C(t, r, θ, ϕ) = Co[1 − a erfc
(r − a)/
√
4t

/r]. The nutrient uptake at the
cell surface is given by I(t) = 4πaDCo(1 + a/√
πt). At steady state, it is I(t → ∞) = 4πaDCo.
To address a linear concentration background, consider a background concentration Cb = Co(1±z/L)
that changes linearly along the z-axis from the base level Co; here, 1/L represents the constant slope of
the concentration gradient, taken such that the length scale L ≫ a is larger than the cell radius a. As
L → ∞, we can recover the uniform concentration field.
Non-dimensional equations Consider the length scale a, time scale a
2/D, and concentration scale
Co. Introduce the dimensionless concentration c = C/Co, and redefine time t := t/(a
2/D) and all spatial
coordinates (·) := (·)/a. Letting ϵ = a/L be the ratio between the cell size and concentration gradient,
the dimensionless counterparts to the unsteady diffusion equation and boundary conditions in Eq.(7.1) are
given by
∂c
∂t = ∇2
c, c(t = 0) = 1 + ϵz, c(r = 1) = 0, ∇c(r → ∞) = ϵez. (7.2)
The dimensionless counterpart to the nutrient uptake at the cell surface is written as
J(t) = I(t)
4πaDCo
=
1
2
Z π
0
∂c(t)
∂r |r=1 sin θdθ. (7.3)
7



In a uniform background concentration [63, 84, 90, 91], the dimensionless nutrient uptake at the cell surface
is given by J(t) = (1 + 1/
√
4πt). At steady state, it is J = 1.
Analytical approach To obtain the time evolution of the concentration field, I solved Eq.(7.2) analytically for c(t, r, θ); the ϕ coordinate is ignorable because of axisymmetry about the z-axis (Fig. 7.2A).
Specifically, we introduced the nonlinear transformation µ = cos θ and used the Laplace transform to get
the unsteady solution.
We consider the concentration field c
′ = c−(1+ϵz). Substitute back into Eq.(7.2), using the nonlinear
transformation µ = cos θ, the concentration field c
′
satisfies
∂c′
∂t = ∇2
c
′
,



c(t = 0, r, µ) = 0,
c(t, r = 1, µ) = −1 − ϵµ,
∇c(t, r → ∞, µ) = 0,
(7.4)
Take the Laplace transform of c
′
to be c¯(s, r, µ) := L{c
′}(s, r, µ) = R ∞
0
c
′
(t, r, µ)e
−stdt; apply the Laplace
transform to the governing equation and boundary conditions to get
sc¯ =
∂
2
c¯
∂r2
+
2
r
∂c¯
∂r +
1
r
2
∂
∂µ
h
(1 − µ
2
)
∂c¯
∂µ
i
, (7.5)



c¯(s, r = 1, µ) = −
1
s
(1 + ϵµ),
c¯(s, r → ∞, µ) = 0,
(7.6)
Use separation of variables c¯ = Φ(µ)R(s, r) and substitute into Eq.(7.5),
∂
2R
∂r2
Φ + 2
r
∂R
∂r Φ + 1
r
2
∂
∂µ
h
(1 − µ
2
)
∂Φ
∂µ
i
R − sΦR = 0. (7.7)
72



Multiply by r
2
Φ(µ)R(s, r)
, take all terms that depend on µ to one side, and define the integer number n ≥ 0,
r
2
d
2R
dr2
1
R
+ 2r
dR
dr
1
R
− sr2 = −
d
dµ
h
(1 − µ
2
)
dΦ
dµ
i
1
Φ
= n(n + 1). (7.8)
We arrive at two differential equations governing R(s, r) and Φ(µ) separately
r
2
d
2R
dr2
+ 2r
dR
dr −

sr2 + n(n + 1)
R = 0, (7.9)
d
dµ
h
(1 − µ
2
)
dΦ
dµ
i
+ n(n + 1)Φ = 0. (7.10)
The function Φ(µ) satisfies the Legendre differential equation, whose solution is given by the Legendre
polynomials Pn(µ); the function R(s, r) ≡ R(
√
sr) satisfies the modified spherical Bessel equation with
√
sr as the independent variable, whose solution is given by the n
th modified spherical Bessel functions of
√
sr. The terms related to modified Bessel function of the first kind in(·) are dropped to ensure convergence
of the solution in the three-dimensional domain of interest is infinite, which extends to infinity and is
bounded internally by a sphere. Putting the two solutions together, we get
c¯(s, r, µ) = X∞
n=0
Ankn(
√
s r)Pn(µ), (7.11)
where An are constant coefficients, Pn(µ) the n
th Legendre polynomial, kn(·) the n
th modified spherical
Bessel function of the second kind.
Substitute the general solution in Eq.(7.11) back into the boundary conditions in Eq.(7.6),
−
1
s
= A0k0(
√
s), −
ϵ
s
= A1k1(
√
s), An = 0 for n ≥ 2. (7.12)
7



The first two modified Bessel functions of the second kind are
k0(z) = e
−z
z
, k1(z) = e
−z
(z + 1)
z
2
. (7.13)
Solve for the unknown coefficients in Eq.(7.12) and substitute back into 7.11,
c¯(s, r, µ) = −
k0(
√
s r)
sk0(
√
s)
P0(µ) −
ϵk1(
√
s r)
sk1(
√
s)
P1(µ). (7.14)
Use the Inversion Theorem to obtain the concentration field
c
′
(t, r, µ) = −
1
2πi Z γ+i∞
γ−i∞

k0(
√
s r)
sk0(
√
s)
P0(µ) + k1(
√
s r)
sk1(
√
s)
ϵP1(µ)

e
stds. (7.15)
Re
Im
γ + i∞
γ - i∞
C D
EF
A
B
Figure 7.1: Integration path in the complex plane. Because the origin is a singular point of the integrand in Eq.(7.15), the path of integration from γ − i∞ to γ + i∞ along the grey line AB is replaced with
integration along the half circle in left plane [25]; specifically, the integral path is changed from the grey
line AB to the blue line A-CDEF-B. Lines CD and EF are very close to the real axis, and the singular point
at origin is avoided by the small circle DE.
74



The integration path from γ − i∞ to γ + i∞ can be replaced by the integral along a half circle with
infinite radius in the left half domain, along paths ACDEFB shown in Fig.7.1, [25]. Note that integral along
the small circle DE of any function f(z) is (2πi Ress=0f(z)), where Res is the residual. Using Eq.(7.13),
the residuals of the two terms in the integral Eq.(7.15) are given by
Ress=0
e
st k0(
√
s r)
sk0(
√
s)

= lim
s→0

e
st k0(
√
s r)
k0(
√
s)

=
1
r
,
Ress=0
e
st k1(
√
s r)
sk1(
√
s)

= lim
s→0

e
st k1(
√
s r)
k1(
√
s)

=
1
r
2
.
(7.16)
Substituting the residuals Eq.(7.16) into the integration Eq.(7.15), we get
c
′ = −
1
2πi
I P0(µ) + II ϵP1(µ) + 2πi
1
r
P0(µ) + 1
r
2
ϵP1(µ)
,
= −
1
2πi
I P0(µ) + II ϵP1(µ)

−
1
r
− ϵµ
1
r
2
,
(7.17)
where
I = Z
CD+EF
k0(
√
s r)
sk0(
√
s)
ds, II = Z
CD+EF
k1(
√
s r)
sk1(
√
s)
ds. (7.18)
Considering s = u
2
e
iπ for the integration in the upper left plane [25], we have the integral along line EF
of the first term
IEF =
Z ∞
0
e
−u
2
t
k0(rueiπ/2
)
k0(ueiπ/2)
2du
u
. (7.19)
The relation of the modified spherical Bessel functions and the spherical Hankel functions are given by [1]
kn(x) =



−
1
2
πinh
(1)
n (ix), −π < arg x ≤
π
2
,
−
1
2
πi−nh
(2)
n (−ix), −
π
2 < arg x ≤ π,
(7.20)
75



where h
(1)
n and h
(2)
n are spherical Hankel functions of the first and second kind; their relations to spherical
Bessel functions jn and yn are given by
h
(1)
n
(x) = jn(x) + iyn(x),
h
(2)
n
(x) = jn(x) − iyn(x).
(7.21)
As a side note, the spherical Bessel and Hankel functions are related to the Hankel and Bessel functions
Jn, Yn, Kn, In, H(1)
n , H(2)
n via jn(z) = p π
2z
Jn+1/2
(z), yn(z) = p π
2z
Yn+1/2
(z), h
(1)
n (z) = p π
2z H
(1)
n+1/2
(z),
h
(2)
n (z) = p π
2z H
(2)
n+1/2
(z) (see more in Append.A.3). Using the relations in Eq.(7.20) and Eq.(7.21), the integration in Eq.(7.19) can be written as
IEF =
Z ∞
0
e
−u
2
t h
(2)
0
(ur)
h
(2)
0
(u)
2du
u
=
Z ∞
0
e
−u
2
t
j0(ur) − iy0(ur)
j0(u) − iy0(u)
2du
u
. (7.22)
For the integration along CD in the lower left plane, we set s = u
2
e
−iπ. By similar procedure, we obtain
ICD =
Z 0
∞
e
−u
2
t
k0(ure−iπ/2
)
k0(ue−iπ/2)
2du
u
= −
Z ∞
0
e
−u
2
t
j0(ur) + iy0(ur)
j0(u) + iy0(u)
2du
u
. (7.23)
Put together, the first integral in Eq.(7.18) is given by
I = IEF + ICD = 4i
Z ∞
0
e
−u
2
t
j0(ur)y0(u) − y0(ur)j0(u)
j
2
0
(u) + y
2
0
(u)
du
u
. (7.24)
Similarly, the second integral in Eq.(7.18) can be obtained
II = 4i
Z ∞
0
e
−u
2
t
j1(ur)y1(u) − y1(ur)j1(u)
j
2
1
(u) + y
2
1
(u)
du
u
. (7.25)
76



Putting all the pieces together, we get the solution of the concentration field
c
′ = −
1
r
−
ϵµ
r
2
−
2
π
Z ∞
0
e
−u
2
t

j0(ur)y0(u) − y0(ur)j0(u)
j
2
0
(u) + y
2
0
(u)
P0(µ)+
j1(ur)y1(u) − y1(ur)j1(u)
j
2
1
(u) + y
2
1
(u)
ϵP1(µ)

du
u
.
(7.26)
The actual concentration field c = c
′ + 1 + ϵz is given by
c = 1 −
1
r
+ ϵµ(r −
1
r
2
) −
2
π
Z ∞
0
e
−u
2
t

j0(ur)y0(u) − y0(ur)j0(u)
j
2
0
(u) + y
2
0
(u)
P0(µ)+
j1(ur)y1(u) − y1(ur)j1(u)
j
2
1
(u) + y
2
1
(u)
ϵP1(µ)

du
u
.
(7.27)
Alternatively, by rearranging terms, the solution can be written as
c(t, r,µ) = 1 −
1
r
erfc
r − 1
√
4t

+ ϵµ
r −
1
r
2

+ ϵµ
−
2
πr2
Z ∞
0
e
−u
2
t u(r − 1) cos (u(1 − r)) + (1 + ru2
) sin (u(1 − r))
u(1 + u
2)
du
.
(7.28)
Rearrange the terms, we can express the final analytical solution as
c(t, r, µ) =1 −
1
r
erfc
r − 1
√
4t

+ ϵµ
r −
1
r
2

−
2ϵµ
πr2
Z ∞
0
e
−u
2
t u(r − 1) cos (u(1 − r)) + (1 + ru2
) sin (u(1 − r))
u(1 + u
2)
du.
(7.29)
Here, u is a placeholder variable for integration. To obtain the steady-state solution c(r, µ), we take the
limit t → ∞,
c(r, µ) = 1 −
1
r
+ ϵµ
r −
1
r
2

. (7.30)
For ϵ = 0, the unsteady and state-state solutions are consistent with the classical solutions of the diffusion
equation in uniform concentration [25].
77



Numerical approach Next, we solved for the unsteady diffusion equation numerically using the Legendre spectral method in θ and employing a finite difference discretization in r with a stretched mesh [90,
79, 78]. For time stepping, we evolved the unsteady diffusion equation using the Backward Euler method
(see details in B.2).
A B
C
t = 5 t = 100
.8
.6
.6
.8
.6
.8
.6
.8
1
2
.8
.6
.6
.8
1.2
time
0 100
steady state
uniform
gradient
.8 1.2
t = 0
50
ez
ex
ey
φ a θ
Nutrient intake, I/4πaD
Co
-20a 0 20a
concentration,
1
2
0
uniform
gradient
D
z coordinate
2a
uniform / gradient
C/Co
Cb = Co
(1 + z/L)
Cb = Co
0
1.2
C/Co
Figure 7.2: Pure diffusion in concentration gradient. A. Mathematical model of sessile ciliated cell as
a sphere of radius a, with spherical (r, θ, ϕ) and Cartesian (x, y, z) coordinates. B. Concentration field of
dissolved nutrients following a diffusion process in an initially (top row) uniform concentration (ϵ = 0)
and (bottom row) gradient (ϵ = 0.05). Snapshots shown at time instances t = [0, 5, 100] and at steadystate. C. Nutrient intake J = I/4πaDCo in pure diffusion in uniform and linear gradients. D. Steady-state
solution of concentration c(x = 0, y = 0, z) as a function of the z coordinate.
Summary In Fig. 7.2B, we show the time evolution of the concentration field in a uniform background
concentration ϵ = 0 (top row) and in a linear concentration gradient ϵ = 0.05 (bottom row). By choice, the
linearly changing background concentration Cb is equal to Co at the cell center z = 0; it increases linearly
for positive z and decreases by the same amount for negative z. As a result, the concentration field is
asymmetric along the z-axis, compared to the uniform background concentration. But, by construction,
this asymmetry leads to the same time evolution of the nutrient uptake J(t)(Fig. 7.2C) because the nutrient
78



flux on the positive and negative half-spheres balance each other. Indeed, from the property of Legendre
polynomials R 1
−1
Pn(µ)dµ = 2δn0, the terms in Eq.(7.29) associated with concentration gradient ϵ do not
contribute to the integral in Eq.(7.3). That is, although the spatial concentration fields differ between
the uniform and linear gradient background (Fig. 7.2B,D), by construction, they both give rise to the same
nutrient uptake at the cell surface (Fig. 7.2C). The difference between the two concentration fields is evident
when examining their variation in the z-direction at steady state (Fig. 7.2D).
7.2 Unsteady advection-diffusion transport
Model set up Coupled with the ciliate’s induced flow field, the nutrient concentration field is governed
by an advection-diffusion equation. With focusing on transient nutrient concentration evolution, we consider far-field boundary condition preserved within the considered time. The governing equation, initial
condition and boundary conditions are
∂C
∂t + u · ∇C = D∇2C,
C|t=0 = Cb, C|r=a = 0, ∇C|r→∞ = ϵez.
(7.31)
Non-dimensional equations Consider the length scale a, concentration scale Co, and velocity scale U
based on the energy dissipation rate U =
p
P/(8πaη) [79]. Here, we have two time scales: the advection
time scale τadvection = a/U and the diffusion time scale τdiffusion = a
2/D. Choosing the advection time
scale τadvection, we obtain the dimensionless transport equation, together with the initial and boundary
conditions,
∂c
∂t + u · ∇c =
1
Pe
∇2
c,
c|t=0 = 1+ϵz, c|r=1 = 0, ∇c|r→∞ = ϵez,
(7.32)
79



where the Péclet number Pe = aU/D is the ratio between the advection and diffusion rates.
Decomposition of the concentration field To simplify the solution to the advection-diffusion equation in Eq.(7.32), we introduce the concentration field c
′ = c − (1 + ϵz), which represents the disturbance
to the background concentration due to the presence of the cell. This field satisfies
∂c′
∂t + u · ∇c
′ + ϵuz =
1
Pe
∇2
c
′
,
c
′
|t=0 = 0, c′
|r=1 = − (1 + ϵµ), ∇c
′
|r→∞ = 0,
(7.33)
where uz = u · ez is the fluid velocity component along ez direction. By linearity of the advectiondiffusion equation, we further decompose the concentration field c
′
into two parts c
′ = ch + ϵc˜, where ch
is the solution associated with homogeneous ambient concentration background, and c˜ is the correction
due to the concentration gradient. Substituting back into Eq.(7.33), we get two sets of differential equations
and boundary conditions. The equations governing ch(t, r, θ) are given by
∂ch
∂t + u · ∇ch =
1
Pe
∇2
ch,
ch|t=0 = 0, ch|r=1 = −1, ∇ch|r→∞ = 0,
(7.34)
The equations governing c˜(t, r, θ) are given by
∂c˜
∂t + u · ∇c˜+ uz =
1
Pe
∇2
c, ˜
c˜|t=0 = 0, c˜r=1 = −µ, ∇c˜r→∞ = 0.
(7.35)
That is, the concentration c˜(t, r, µ) follows an advection-diffusion process with a distributed ’source’ term
represented by uz(r, µ). This source term describes a supply of nutrients by the component uz(r, µ) of
cilia-induced flow that is aligned with the background concentration gradient in the ez direction. The
80



strength of this source term decays as 1/r, indicating that the cell-induced effect decays inversely with
distance from the cell.
The dimensionless nutrient intake is given as the sum of nutrient intake Jh in a uniform concentration
field and nutrient intake J˜ in a concentration gradient field, calibrated to have zero background concentration at the cell center.
J(t) = Jh(t) + ϵJ˜(t) = 1
2
Z 1
−1
∂ch
∂r



r=1
dµ + ϵ
1
2
Z 1
−1
∂c˜
∂r



r=1
dµ. (7.36)
7.3 Numerical results
Unsteady transport in uniform field Consider a sessile ciliated sphere generating feeding currents in
a uniform concentration with ϵ = 0. The steady-state limit of this problem has been studied [112, 63, 84,
90, 62, 79, 78]. In this section, we will focus on the unsteady dynamics of nutrient transport. By definition,
the concentration field is c(t, r, µ) = 1 + ch(t, r, µ). We can numerical solved equation Eq.(7.34) for
ch(t, r, µ). In Fig. 7.3A, it is shown that the streamlines (blue) of flow fields generated by the ciliated sphere
with surface velocity corresponding to one of the first four modes: ’treadmill’ (mode 1), dipolar (mode 2),
tripolar (mode 3), and quadripolar (mode 4). In all cases, the total hydrodynamic power P/8πηa = 1
are constrained as a constant. Only the treadmill mode (n = 1) contributes a net hydrodynamic force
that needs to be balanced by an attachment force. By symmetry, the subsequent three higher-order modes
produce zero force. Notice that odd modes (n = 1, 3) produce flows with lateral (left-right) symmetry
while even modes (n = 2, 4) result in flows with both lateral and front-back symmetry.
At Pe = 10, the snapshots of the concentration field c(t, r, µ) for each mode at time t = 20 are shown
in Fig. 7.3B. The concentration fields reflect the symmetries of the underlying flows. The time evolution
of the nutrient intake J(t) corresponds to each of these four modes is shown in Fig. 7.3E. Compared to
higher-order modes, the treadmill mode, at the same hydrodynamic power, provides the largest nutrient
81



B
A
.8 .6
mode 1 mode 2 mode 3 mode 4
C
c˜
-4
4
-1.5
1.5
-1
1
-1
1
1
.1
2
2
-.1 .1 .1 -.1
-1
perturbed field, ˜c
D
0
1.5
.8 1.2
0
10
|| ∆
c˜||
concentration field,
C/Co
2a
2a
2a
2a
flow streamlines concentration field,
C/Co
0
1
2a
.8
.6
C/Co
C/Co
10 20 0
40
20
mode 1
mode 2 ~t
~t0.5
~t0.3
0
~t0.1
E H Nutrient intake, J(t)
mode 1
mode 2, 3, 4
mode 4
mode 3
mode 2
mode 1
3
2
1
100 200
0
.5
1
0
time, t
||J(t) - Jsteady ||
Jsteady
100 2000
time, t
100 200
0
20
40
mode 1
200
0
0.01
0
mode 3
0
Nutrient intake, (t) ˜J
~ t0.5
F
uniform uniform G gradient gradient R(t)
time, t time, t
Figure 7.3: Unsteady transport in uniform and linear concentration. A. Flow streamlines around a
sessile ciliated cell generated by the first four surface velocity modes. Concentration field c in B. uniform
and D. linear gradient, where c = ch + ϵc˜+ (1 + ϵz). C. snapshots of c˜ with insets showing |∇c˜| around
the cell. E. Time evolution of nutrient intake Jh(t) in uniform concentration. In uniform background
concentration, E. mode 1 exhibits the largest nutrient intake and F. fastest convergence to steady-state. In
linear gradients, G. the component J˜(t) is identically zero for even modes and largest for the first mode.
H. the distance R(t) where c˜ begins to reach 10% grows with time.
82



intake to the cell, consistent with findings in [90, 79] at steady-state. Importantly, depending on the surface
mode, the rate of change of nutrient intake approaches steady-state at different rates: the treadmill mode
approaches steady solution at a faster rate compared to higher-order modes (Fig. 7.3F).
Unsteady transport in linear field Consider now the sessile ciliated sphere in a linear concentration
gradient with ϵ = 0.05. Given the numerical solution ch(t, r, µ) of Eq.(7.34), we can separately solve
Eq.(7.35) for c˜(t, r, µ). By substituting the numerical solutions into c = ch+ϵc˜+(1+ϵz), we can recover the
full concentration field. The concentration field c˜ is shown in Fig. 7.3C for all four surface velocity modes.
The insets show a zoom in on the concentration gradient |∇c˜| around the cell, with dashed contour lines
|∇c˜| = 2 shown to highlight the fact that these contour lines are front-back asymmetric for odd modes
and front-back symmetric for even modes – with opposite sign of the concentration gradient. Fig. 7.3G
shows the corresponding nutrient intake J˜(t) for each mode. As c˜(t, r, µ) changes in time, the nutrient
intake J˜(t) associated with the even modes is identically zero for all time because of the front-back flow
symmetries of these modes; J˜(t) associated with the first mode increases with the square root of time, but
J˜(t) associated with the third mode converges to a constant as time increases.
To emphasize the role of the ’forcing term’ uz(r, µ) due to cilia-induced flows in driving nutrients to
the cell surface in Eq.(7.35), I introduced a distance R(t)that captures how, starting from an identically zero
initial value c˜(t = 0, r, µ) = 0, uz(r, µ) creates a region of high concentration around the cell. Specifically,
I defined R(t) such that {c <˜ 0.1 for r > R}; the distance R(t) marks the region within which c˜increases
from zero up to a percentage of the normalized baseline C/Co = 1, beyond R(t), this concentration drops
to below 10%. This distance can be used for approximating the growth domain of the concentration field
c˜(t, r, µ). The time evolution of this distance R(t) is shown in Fig. 7.3H; R scales with time as follows:
for mode 1, R ∼ t grows linearly in time; for mode 2, R ∼
√
t grows proportionally to the square root
of time, for modes 3 and 4, R ∼ t
0.3
and R ∼ t
0.1
, respectively, indicating, as expected, that R decreases
with higher-order modes. Mode 1 maximizes transport to the cell surface.
83



A
100
6
12
0
2000
time
Fiiting coefficients
~t0.5
0
20
40
0
numerics
500
0
20
40
2500
Nutrient intake, J(t)
time
estimation
Pe = 500
Pe = 200
Pe = 100
B
10
J˜(t) = a(t)Peb(t)
5
Pe = 500
Pe = 200
Pe = 100 0.5
0.55 b(t) a(t)
Nutrient intake, J(t)
C
time
Figure 7.4: Nutrient intake of the first mode fitting. A. Coefficients of curve fitting of nutrient intake
J˜(t) corresponding to mode 1 versus Pe at every instant time. B. Nutrient intake comparison between
numerical computation J(t) = Jh(t) + ϵJ˜(t) and estimation J(t) = Jh|t→∞ + 0.94ϵ
√
t Pe under Pe
= [100, 200, 500] with time range t ∈ [0, 500]. C. Zoom-in version for nutrient intake comparison for time
range t ∈ [0, 10]. All calculations are with gradient ϵ = 0.05.
Prediction of nutrient uptake Next, with computed the concentration field c(t, r, µ) and nutrient
intake J(t, r, µ) for mode 1 for a range of Péclet number Pe ∈ [1, 100] and a range of time t ∈ [0, 200],
we can fit the component J˜ of the nutrient intake using a power law model J˜ = a(t)Peb(t)
. The fitting
coefficients a(t) and b(t) are shown in Fig. 7.4A, where at all time the R-squared values are larger than 0.99:
b(t) converges to 0.54 after time t > 50 and a(t) scales with t
0.5
at all time, giving rise to J˜ ∼ t
0.5Pe0.54
.
Thus, we can simplify the scaling law J˜(t) = c(
√
tPe)d
and found that J˜(t) = 0.94(t Pe)0.55, with an
R-squared value larger than 0.99, which can be approximated further to J˜(t) = 0.94√
t Pe. Substituting
back into Eq.(7.36), we got an approximate expression for the optimal nutrient intake J(t) associated with
mode 1 in linear concentration gradients,
J(t) = Jh(t) + 0.94ϵ
√
t Pe. (7.37)
In Fig. 7.4B, it is tested that the ability of this power law to predict nutrient uptake beyond the ranges
of Péclet number Pe ∈ [1, 100] and time t ∈ [0, 200] which were used to obtain the fitted coefficient.
In particular, I compared the model prediction in Eq.(7.37) to results from full numerical simulations for
Pe = 100, 200, and 500, for a time range extending up to t = 500. The agreement between the model
84



in Eq.(7.37) and the numerical simulations beyond a short time range is excellent. To interrogate the initial
mismatch at short time scales, see the zoom in version in Fig. 7.4C; again, beyond an initial time window
not exceeding t = 5, the model exhibits remarkable agreement with simulations.
7.4 Unsteady transport in patchy concentration fields
In this section, I will show two case studies that represent the feeding of sessile ciliate in a patchy concentration background. First, consider an initial concentration field c(0, r, µ) = 1 − Rb/|r − r0,b| that
is a steady-state solution of the diffusion equation outside an empty blob of size Rb centered at r0,b. In
Fig. 7.5A, it can be shown this initial concentration field for Rb = 1 and r0,b = 3ez, representing a depleted
concentration patch upstream of the mode ciliate. Using this initial concentration field, we can directly
solve the unsteady state advection-diffusion equation 7.31 numerically at Pe = 10 for mode 1 and mode 2
at the same constant energy dissipation rate, as done before. The time evolution of the concentration field
associated with each mode is shown in Fig. 7.5B, C for mode 1 and mode 2, respectively. The corresponding
nutrient intake J(t) is shown in Fig. 7.5D. These results show that mode 2 outperforms mode 1 at an initial
transient time, but mode 1 takes over as the optimal mode in the long term.
Next, we can consider an initial concentration field, where the entire domain is depleted except for a
nutrient-rich blob (shown in black in Fig. 7.6A) placed on one side of the ciliated sphere. The time evolution
of this nutrient-rich patch due to cilia-induced flows is shown in Fig. 7.6B, C for mode 1 and mode 2, and
the nutrient intake versus time is shown in Fig. 7.6D. In contrast to the depleted nutrient patch study, mode
1 is optimal during the initial transient phase, while mode 2 becomes more advantageous over the long
term.
Taken together, these results show that in a linear concentration gradient with potentially infinite nutrient supply, mode 1 is optimal when the direction of cilia-induced velocity is aligned with the direction of
85



A B time = 20
10 20
0
3
0
time = 2 time = 10
nutrient intake, J(t)
time = 0
depleted patch
ahead
mode 1mode 2
D C
time
1
2
.5
mode 1
mode 2
2a
.7 .7
.5
0
1
C/Co
Figure 7.5: Nutrient uptake in a concentration field with depleted upstream patch. A. The initial
depleted patchy concentration upstream of the model cell. B-C. The concentration field snapshots at time
t = [0, 2, 10, 20] corresponding to surface velocity with B. mode 1 and C. mode 2. D. The nutrient intake
J(t) versus time corresponds with input surface velocity with only mode 1 (black) and mode 2 (grey). All
case studies are under Pe = 10.
concentration,
B
0
1
10 20 0
0.4
0
mode 1
mode 2
Nutrient intake, J(t) concentration patch ahead
mode 1mode 2
.7
.5
time
.01
0.2
C/Co
A time = 0 time = 2 time = 10 time = 20
CD
Figure 7.6: Nutrient uptake in a depleted environment with a rich concentration patch. A. The
initial field contains a concentrated blob upstream of the model cell. B-C. The concentration snapshots at
time t = [0, 2, 10, 20] corresponding to surface velocity with B. mode 1 and C. mode 2. D. Nutrient intake
versus time corresponds to surface velocity with mode 1 (black) and mode 2(grey). All case studies are
under Pe = 10.
86



the nutrient concentration gradient. However, in nutrient-depleted domains with sparse nutrient availability, mode 1 may become suboptimal. An adaptive adjustment of beating modes needs to be considered for
optimizing nutrient intake in non-uniform fields. These findings underscore that in dynamically changing
environments, there is no unique strategy or surface activity that maximizes uptake at all time.
7.5 Summary
In still water, in the absence of feeding currents, the problem is governed by nutrient diffusion only. In
this diffusive limit, starting from an initial linear concentration field, an analytical solution for the unsteady concentration field around a spherical cell is derived. By symmetry, the diffusion process in such
linear gradients provides no benefits over what is achievable by a sessile cell in a uniform background
concentration.
Given the ability of the ciliated cell to generate feeding currents and considering ciliary surface activities that are axisymmetric about the direction of concentration gradient, the envelope model is considered
for sessile ciliate. To properly compare feeding efficiencies across different modes, I fixed the hydrodynamics dissipation power to be constant. I then formulated the feeding problem as an advection-diffusion
process in linear concentration gradients. By decomposing the concentration field into two fields governed
by a homogeneous and a "forced" advection-diffusion process, I solved for each field separately using a
spectral method. It can be found that, in linear concentration gradients, only the odd modes contribute to
nutrient uptake. Odd modes generate asymmetric flows along the direction of concentration gradient that
increase nutrient intake; even modes generate symmetric flows and thus provide no advantage beyond
what is achievable in a uniform concentration field. Moreover, only the treadmill surface activity (mode
1) contributes significantly to sustained nutrient intake, while higher-order modes have only a minor impact on nutrient intake. Focusing on the treadmill mode and analyzing how nutrient uptake scales with
87



Péclet and time, I obtained an analytic expression J(t) = Jh|t→∞ + 0.94ϵ
√
t Pe, where ϵ is the concentration gradient. This approximate model accurately predicts the optimal rate of nutrient intake in the long
time regime and large Péclet for a sessile ciliate in linear concentration gradients. Lastly, by considering
both depleted and nutrient-rich patches upstream of the cell, it can be found that shifting between ciliary
beating modes optimizes cumulative nutrient intake beyond what is achievable by a single mode.
These findings demonstrate that directional cilia beating aligned with the concentration gradient is
essential for maximizing feeding performance. As such, they offer new perspectives on the mechanisms
by which sessile ciliates could regulate their orientation in response to concentration gradients, say relative
to the substrate to which they are attached [109, 107], to maximize feeding flux. Extending the approach to
motile ciliates in concentration gradients would open up exciting opportunities to study chemotaxis [46,
69, 99, 142, 37]. More generally, these findings are important in predicting the ecological impact of sessile
ciliates on both local and global nutrient distributions within microbial communities [31, 122].
88



Chapter 8
Future work
This thesis has focused on modeling fluid dynamics and nutrient transport in both swimming and sessile
ciliated microorganisms. Building on this work, several promising directions emerge for future research.
These include extending current models to non-axisymmetric squirmers, exploring motion strategies in
dynamic chemical environments, and incorporating adaptive learning mechanisms. In addition to advancing our understanding of biological systems, these investigations hold potential applications in the
development of bio-inspired micro-robots for biomedical purposes.
8.1 Beyond axisymmetric squirmers
To capture a broader range of biological and artificial swimmer behaviors, future research can relax the
assumption of axisymmetry in squirmer models. This extension would enable the investigation of optimal
stroke patterns not only for forward locomotion but also for turning, maneuvering, and responding to
environmental stimuli—such as rapid escape responses triggered by flow or chemical cues. These models
will incorporate asymmetric surface motions and swimmer-environment interactions to better understand
how locomotion strategies evolve in complex settings. Insights gained from this work could inform the
design of robotic swimmers that operate efficiently in unstructured or unpredictable environments.
89



8.2 Sensing, Motion Control, and Learning
Another exciting direction is the incorporation of adaptive control strategies into micro-swimmer models. Moving beyond deterministic simulations, reinforcement learning (RL) frameworks can be employed
to discover optimal behaviors for nutrient acquisition in spatially and temporally varying chemical landscapes. By embedding local sensing mechanisms, it becomes possible to simulate how micro-agents dynamically adjust their motion to maximize resource uptake. These studies can reveal principles of biological adaptation while also guiding the development of intelligent micro-robots capable of real-time
decision-making. Such bio-inspired control, integrated with advancements in micro-robotic actuation, has
promising applications in targeted drug delivery and other biomedical tasks requiring precise navigation
in complex environments.
90



Bibliography
[1] Milton Abramowitz and Irene A Stegun. Handbook of mathematical functions with formulas,
graphs, and mathematical tables. Vol. 55. US Government printing office, 1968.
[2] Andreas Acrivos and JD Goddard. “Asymptotic expansions for laminar forced-convection heat
and mass transfer”. In: Journal of Fluid Mechanics 23.2 (1965), pp. 273–291.
[3] Andreas Acrivos and Thomas D Taylor. “Heat and mass transfer from single spheres in Stokes
flow”. In: The Physics of Fluids 5.4 (1962), pp. 387–394.
[4] TD Ainsworth, AJ Fordyce, and EF Camp. “The other microeukaryotes of the coral reef
microbiome”. In: Trends in Microbiology 25.12 (2017), pp. 980–991.
[5] Anders Andersen and Thomas Kiørboe. “The effect of tethering on the clearance rate of
suspension-feeding plankton”. In: Proceedings of the National Academy of Sciences 117.48 (2020),
pp. 30101–30103.
[6] Frank O Aylward, John M Eppley, Jason M Smith, Francisco P Chavez, Christopher A Scholin, and
Edward F DeLong. “Microbial community transcriptional networks are conserved in three
domains at ocean basin scales”. In: Proceedings of the National Academy of Sciences 112.17 (2015),
pp. 5443–5448.
[7] Roman Babko, Tetiana Kuzmina, Yaroslav Danko, Joanna Szulżyk-Cieplak, and Grzegorz Łagód.
“Oxygen gradients and structure of the ciliate assemblages in floodplain lake”. In: Water 12.8
(2020), p. 2084.
[8] Štefan Bálint, Ione Verdeny Vilanova, Ángel Sandoval Álvarez, and Melike Lakadamyali.
“Correlative live-cell and superresolution microscopy reveals cargo transport dynamics at
microtubule intersections”. In: Proceedings of the National Academy of Sciences 110.9 (2013),
pp. 3375–3380.
[9] Cx K Batchelor and GK Batchelor. An introduction to fluid dynamics. Cambridge university press,
1967.
[10] Brian A Bayless, Francesca M Navarro, and Mark Winey. “Motile cilia: innovation and insight
from ciliate model organisms”. In: Frontiers in cell and developmental biology 7 (2019), p. 265.
91



[11] C DALE Beers. “Diet in relation to depression and recovery in the ciliate Didinium nasutum.” In:
(1933).
[12] Howard C Berg. Random walks in biology. Princeton University Press, 1993.
[13] Howard C Berg. “Random walks in biology”. In: Random Walks in Biology. Princeton University
Press, 2018.
[14] Howard C Berg and Edward M Purcell. “Physics of chemoreception”. In: Biophysical journal 20.2
(1977), pp. 193–219.
[15] Oliver S Beveridge, Owen L Petchey, and Stuart Humphries. “Mechanisms of
temperature-dependent swimming: the importance of physics, physiology and body size in
determining protist swimming speed”. In: Journal of Experimental Biology 213.24 (2010),
pp. 4223–4231.
[16] Luis Alberto Bezares-Calderón, Jürgen Berger, and Gáspár Jékely. “Diversity of cilia-based
mechanosensory systems and their functions in marine animal behaviour”. In: Philosophical
Transactions of the Royal Society B 375.1792 (2020), p. 20190376.
[17] William Bialek. Biophysics: searching for principles. Princeton University Press, 2012.
[18] John R Blake. “A note on the image system for a stokeslet in a no-slip boundary”. In:
Mathematical Proceedings of the Cambridge Philosophical Society. Vol. 70. 2. Cambridge University
Press. 1971, pp. 303–310.
[19] John R Blake. “A spherical envelope approach to ciliary propulsion”. In: Journal of Fluid
Mechanics 46.1 (1971), pp. 199–208.
[20] Robert A Bloodgood. “Sensory reception is an attribute of both primary cilia and motile cilia”. In:
Journal of cell science 123.4 (2010), pp. 505–509.
[21] Jens Boenigk and Hartmut Arndt. “Bacterivory by heterotrophic flagellates: community structure
and feeding strategies”. In: Antonie Van Leeuwenhoek 81 (2002), pp. 465–480.
[22] Martin E Boraas, Dianne B Seale, and Joseph E Boxhorn. “Phagotrophy by a flagellate selects for
colonial prey: a possible origin of multicellularity”. In: Evolutionary Ecology 12 (1998),
pp. 153–164.
[23] Christopher Brennen and Howard Winet. “Fluid mechanics of propulsion by cilia and flagella”. In:
Annual Review of Fluid Mechanics 9 (1977), pp. 339–398.
[24] WE Bullington. “A further study of spiraling in the ciliate Paramecium, with a note on
morphology and taxonomy”. In: Journal of Experimental Zoology 56.4 (1930), pp. 423–449.
[25] Horatio Scott Carslaw and John Conrad Jaeger. Conduction of heat in solids. Tech. rep. Clarendon
Press, 1959.
92



[26] Thomas Cavalier-Smith. “The phagotrophic origin of eukaryotes and phylogenetic classification
of Protozoa.” In: International journal of systematic and evolutionary microbiology 52.2 (2002),
pp. 297–354.
[27] Ray Chang and Manu Prakash. “Biophysical limits of ultrafast cellular motility”. In: bioRxiv
(2024), pp. 2024–08.
[28] Karen K Christensen-Dalsgaard and Tom Fenchel. “Increased filtration efficiency of attached
compared to free-swimming flagellates”. In: Aquatic microbial ecology 33.1 (2003), pp. 77–86.
[29] Horst Werner Doelle. Bacterial metabolism. Academic Press, 2014.
[30] Julia Dölger, Lasse Tor Nielsen, Thomas Kiørboe, and Anders Andersen. “Swimming and feeding
of mixotrophic biflagellates”. In: Scientific Reports 7.1 (2017), pp. 1–10.
[31] Andrew Dopheide, Gavin Lear, Rebecca Stott, and Gillian Lewis. “Relative diversity and
community structure of ciliates in stream biofilms according to molecular and microscopy
methods”. In: Applied and Environmental Microbiology 75.16 (2009), pp. 5261–5272.
[32] Knut Drescher, Raymond E Goldstein, Nicolas Michel, Marco Polin, and Idan Tuval. “Direct
measurement of the flow field around swimming microorganisms”. In: Physical Review Letters
105.16 (2010), p. 168101.
[33] William M Durham, Eric Climent, Michael Barry, Filippo De Lillo, Guido Boffetta,
Massimo Cencini, and Roman Stocker. “Turbulence drives microscale patches of motile
phytoplankton”. In: Nature communications 4.1 (2013), p. 2148.
[34] Jens Elgeti, Roland G Winkler, and Gerhard Gompper. “Physics of microswimmers—single
particle motion and collective behavior: a review”. In: Reports on progress in physics 78.5 (2015),
p. 056601.
[35] Laura Eme, Anja Spang, Jonathan Lombard, Courtney W Stairs, and Thijs JG Ettema. “Archaea
and the origin of eukaryotes”. In: Nature Reviews Microbiology 15.12 (2017), pp. 711–723.
[36] Richard B Emlet. “Flow fields around ciliated larvae: Effects of natural and artificial tethers.” In:
Marine ecology progress series. Oldendorf 63.2 (1990), pp. 211–225.
[37] Wen-Zhen Fang, Tongzhao Xiong, On Shun Pak, and Lailai Zhu. “Data-driven intelligent
manipulation of particles in microfluidics”. In: Advanced Science 10.5 (2023), p. 2205382.
[38] Tom Fenchel. “Suspension feeding in ciliated protozoa: feeding rates and their ecological
significance”. In: Microbial Ecology 6 (1980), pp. 13–25.
[39] Tom Fenchel and Nicholas Blackburn. “Motile chemosensory behaviour of phagotrophic protists:
mechanisms for and efficiency in congregating at food patches”. In: Protist 150.3 (1999),
pp. 325–336.
[40] Joseph Frankel. “Cell biology of Tetrahymena thermophila”. In: Methods in cell biology 62 (1999),
pp. 27–125.
93



[41] Anette Funfak, Cathy Fisch, Hatem T Abdel Motaal, Julien Diener, Laurent Combettes,
Charles N Baroud, and Pascale Dupuis-Williams. “Paramecium swimming and ciliary beating
patterns: a study on four RNA interference mutations”. In: Integrative Biology 7.1 (2015),
pp. 90–100.
[42] Josep M Gasol and David L Kirchman. Microbial ecology of the oceans. John Wiley & Sons, 2018.
[43] Brad J Gemmell, Genesok Oh, Edward J Buskey, and Tracy A Villareal. “Dynamic sinking
behaviour in marine phytoplankton: rapid changes in buoyancy may aid in nutrient uptake”. In:
Proceedings of the Royal Society B: Biological Sciences 283.1840 (2016), p. 20161126.
[44] THJ Gilmour. “Ciliation and function of the food-collecting and waste-rejecting organs of
lophophorates”. In: Canadian Journal of Zoology 56.10 (1978), pp. 2142–2155.
[45] William Gilpin, Matthew Storm Bull, and Manu Prakash. “The multiscale physics of cilia and
flagella”. In: Nature Reviews Physics 2.2 (2020), pp. 74–88.
[46] Raymond E Goldstein. “Green algae as model organisms for biological fluid dynamics”. In: Annual
review of fluid mechanics 47.1 (2015), pp. 343–375.
[47] Fernando Gómez, Lu Wang, and Senjie Lin. “Morphology and molecular phylogeny of peritrich
ciliate epibionts on pelagic diatoms: Vorticella oceanica and Pseudovorticella coscinodisci sp.
nov.(Ciliophora, Peritrichia)”. In: Protist 169.2 (2018), pp. 268–279.
[48] Jean-David Grattepanche, Laura M Walker, Brittany M Ott, Daniela L Paim Pinto,
Charles F Delwiche, Christopher E Lane, and Laura A Katz. “Microbial diversity in the eukaryotic
SAR clade: Illuminating the darkness between morphology and molecular data”. In: BioEssays
40.4 (2018), p. 1700198.
[49] Richard K Grosberg and Richard R Strathmann. “The evolution of multicellularity: a minor major
transition?” In: Annu. Rev. Ecol. Evol. Syst. 38.1 (2007), pp. 621–654.
[50] Jeffrey S Guasto, Roberto Rusconi, and Roman Stocker. “Fluid mechanics of planktonic
microorganisms”. In: Annual Review of Fluid Mechanics 44 (2012), pp. 373–400.
[51] John Happel and Howard Brenner. Low Reynolds number hydrodynamics: with special applications
to particulate media. Vol. 1. Springer Science & Business Media, 2012.
[52] Christoph Hartmann, Özlem Özmutlu, Hannes Petermeier, Johannes Fried, and Antonio Delgado.
“Analysis of the flow field induced by the sessile peritrichous ciliate Opercularia asymmetrica”.
In: Journal of biomechanics 40.1 (2007), pp. 137–148.
[53] Cody E Hinchliff, Stephen A Smith, James F Allman, J Gordon Burleigh, Ruchi Chaudhary,
Lyndon M Coghill, Keith A Crandall, Jiabin Deng, Bryan T Drew, Romina Gazis, et al. “Synthesis
of phylogeny and taxonomy into a comprehensive tree of life”. In: Proceedings of the National
Academy of Sciences 112.41 (2015), pp. 12764–12769.
[54] Alexander T Ho and Laurence D Hurst. “Effective population size predicts local rates but not local
mitigation of read-through errors”. In: Molecular Biology and Evolution 38.1 (2021), pp. 244–262.
94



[55] Portia A Holt and George B Chapman. “The fine structure of the cyst wall of the ciliated
protozoon Didinium nasutum 1”. In: The Journal of Protozoology 18.4 (1971), pp. 604–614.
[56] Masaki Ishida, Richard D Allen, and Agnes K Fok. “Phagosome formation in Paramecium: roles of
somatic and oral cilia and of solid particles as revealed by video microscopy”. In: Journal of
Eukaryotic Microbiology 48.6 (2001), pp. 640–646.
[57] Houshuo Jiang and Edward J Buskey. “Relating ciliary propulsion morphology and flow to
particle acquisition in marine planktonic ciliates I: the tintinnid ciliate Amphorides
quadrilineata”. In: Journal of Plankton Research (2024), fbae012.
[58] Houshuo Jiang and Edward J Buskey. “Relating ciliary propulsion morphology and flow to
particle acquisition in marine planktonic ciliates II: the oligotrich ciliate Strombidium capitatum”.
In: Journal of Plankton Research (2024), fbae011.
[59] Houshuo Jiang and Edward J Buskey. “Relating ciliary propulsion morphology and flow to
particle acquisition in marine planktonic ciliates II: the oligotrich ciliate Strombidium capitatum”.
In: Journal of Plankton Research (2024), fbae011.
[60] Mengmeng Jiang, Tao Hu, Zhaoyi Wang, Ziyao Liang, Jiqiu Li, and Xiaofeng Lin. “Morphology
and phylogeny of three Pseudovorticella species (Ciliophora: Peritrichia) from brackish waters of
China”. In: Journal of Eukaryotic Microbiology 66.6 (2019), pp. 869–881.
[61] Per R Jonsson, Mona Johansson, and Richard W Pierce. “Attachment to suspended particles may
improve foraging and reduce predation risk for tintinnid ciliates”. In: Limnology and
oceanograpghy 49.6 (2004), pp. 1907–1914.
[62] Eva A Kanso, Rubens M Lopes, J Rudi Strickler, John O Dabiri, and John H Costello. “Teamwork
in the viscous oceanic microscale”. In: Proceedings of the National Academy of Sciences 118.29
(2021), e2018193118.
[63] L Karp-Boss, E Boss, PA Jumars, et al. “Nutrient fluxes to planktonic osmotrophs in the presence
of fluid motion”. In: Oceanography and marine biology 34 (1996), pp. 71–108.
[64] Johannes M Keegstra, Francesco Carrara, and Roman Stocker. “The ecological roles of bacterial
chemotaxis”. In: Nature Reviews Microbiology 20.8 (2022), pp. 491–504.
[65] Patrick J Keeling, Gertraud Burger, Dion G Durnford, B Franz Lang, Robert W Lee,
Ronald E Pearlman, Andrew J Roger, and Michael W Gray. “The tree of eukaryotes”. In: Trends in
ecology & evolution 20.12 (2005), pp. 670–676.
[66] S Kim and SJ Karrila. Microhydrodynamics Butterworth. 1991.
[67] Thomas Kiørboe. “Predation in a Microbial World: Mechanisms and Trade-Offs of Flagellate
Foraging”. In: Annual Review of Marine Science 16 (2023).
[68] Julius B Kirkegaard and Raymond E Goldstein. “Filter-feeding, near-field flows, and the
morphologies of colonial choanoflagellates”. In: Physical Review E 94.5 (2016), p. 052401.
95



[69] Julius B Kirkegaard and Raymond E Goldstein. “The role of tumbling frequency and persistence
in optimal run-and-tumble chemotaxis”. In: IMA Journal of Applied Mathematics 83.4 (2018),
pp. 700–719.
[70] Alexandra Klimenko, Yury Matushkin, Nikolay Kolchanov, and Sergey Lashin. “Leave or stay:
Simulating motility and fitness of microorganisms in dynamic aquatic ecosystems”. In: Biology
10.10 (2021), p. 1019.
[71] Eric Lauga and Thomas R Powers. “The hydrodynamics of swimming microorganisms”. In:
Reports on Progress in Physics 72.9 (2009), p. 096601.
[72] Federico M Lauro, Diane McDougald, Torsten Thomas, Timothy J Williams, Suhelen Egan,
Scott Rice, Matthew Z DeMaere, Lily Ting, Haluk Ertan, Justin Johnson, et al. “The genomic basis
of trophic strategy in marine bacteria”. In: Proceedings of the National Academy of Sciences 106.37
(2009), pp. 15527–15533.
[73] Johanna Laybourn. “Energy budgets for Stentor coeruleus Ehrenberg (Ciliophora)”. In: Oecologia
22 (1976), pp. 431–437.
[74] Brian S Leander. “Predatory protists”. In: Current Biology 30.10 (2020), R510–R516.
[75] MJ Lighthill. “On the squirming motion of nearly spherical deformable bodies through liquids at
very small Reynolds numbers”. In: Communications on Pure and Applied Mathematics 5.2 (1952),
pp. 109–118.
[76] Maciej Lisicki, Marcos F Velho Rodrigues, Raymond E Goldstein, and Eric Lauga. “Swimming
eukaryotic microorganisms exhibit a universal speed distribution”. In: Elife 8 (2019), e44907.
[77] Jingyi Liu, John H Costello, and Eva Kanso. “Flow physics of nutrient transport drives functional
design of ciliates”. In: (2025). Manuscript accepted for publication.
[78] Jingyi Liu, Yi Man, John H Costello, and Eva Kanso. “Feeding rates in sessile versus motile ciliates
are hydrodynamically equivalent”. In: eLife 13 (2024).
[79] Jingyi Liu, Yi Man, John H Costello, and Eva Kanso. “Optimal feeding in swimming and attached
ciliates”. In: Journal of Fluid Mechanics 1003 (2025), A26.
[80] Richard A Long and Farooq Azam. “Microscale patchiness of bacterioplankton assemblage
richness in seawater”. In: Aquatic Microbial Ecology 26.2 (2001), pp. 103–113.
[81] Michael Lynch and Wilfried Gabriel. “Mutation load and the survival of small populations”. In:
Evolution 44.7 (1990), pp. 1725–1737.
[82] Michael Lynch, Paul E Schavemaker, Timothy J Licknack, Yue Hao, and Arianna Pezzano.
“Evolutionary bioenergetics of ciliates”. In: Journal of Eukaryotic Microbiology 69.5 (2022), e12934.
[83] DW Macdonald and DDP Johnson. “Patchwork planet: the resource dispersion hypothesis,
society, and the ecology of life”. In: Journal of Zoology 295.2 (2015), pp. 75–107.
96



[84] Vanesa Magar, Tomonobu Goto, and Timothy J Pedley. “Nutrient Uptake by a Self-Propelled
Steady Squirmer”. In: The Quarterly Journal of Mechanics and Applied Mathematics 56.1 (2003),
pp. 65–91.
[85] Vanesa Magar and Timothy J Pedley. “Average nutrient uptake by a self-propelled unsteady
squirmer”. In: Journal of fluid mechanics 539 (2005), pp. 93–112.
[86] R Timothy Marler and Jasbir S Arora. “Survey of multi-objective optimization methods for
engineering”. In: Structural and multidisciplinary optimization 26 (2004), pp. 369–395.
[87] Maira Maselli, Andreas Altenburger, Diane K Stoecker, and Per Juel Hansen. “Ecophysiological
traits of mixotrophic Strombidium spp”. In: Journal of Plankton Research 42.5 (2020), pp. 485–496.
[88] Amit D Mehta, Ronald S Rock, Matthias Rief, James A Spudich, Mark S Mooseker, and
Richard E Cheney. “Myosin-V is a processive actin-based motor”. In: Nature 400.6744 (1999),
pp. 590–593.
[89] Sébastien Michelin and Eric Lauga. “Efficiency optimization and symmetry-breaking in a model
of ciliary locomotion”. In: Physics of fluids 22.11 (2010), p. 111901.
[90] Sébastien Michelin and Eric Lauga. “Optimal feeding is optimal swimming for all Péclet
numbers”. In: Physics of Fluids 23.10 (2011), p. 101901.
[91] Sébastien Michelin and Eric Lauga. “Unsteady feeding and optimal strokes of model ciliates”. In:
Journal of Fluid Mechanics 715 (2013), pp. 1–31.
[92] Richard E Michod and Denis Roze. “Cooperation and conflict in the evolution of multicellularity”.
In: Heredity 86.1 (2001), pp. 1–7.
[93] Kevin A Miklasz and Mark W Denny. “Diatom sinkings speeds: Improved predictions and insight
from a modified Stokes’ law”. In: Limnology and Oceanography 55.6 (2010), pp. 2513–2525.
[94] David R Mitchell. “The evolution of eukaryotic cilia and flagella as motile and sensory
organelles”. In: Eukaryotic membranes and cytoskeleton: Origins and evolution (2007), pp. 130–140.
[95] Saurabh S Mogre, Aidan I Brown, and Elena F Koslover. “Getting around the cell: physical
transport in the intracellular world”. In: Physical Biology 17.6 (2020), p. 061003.
[96] David JS Montagnes. “Growth responses of planktonic ciliates in the genera Strobilidium and
Strombidium”. In: Marine Ecology Progress Series 130 (1996), pp. 241–254.
[97] W. J. de Munk and Gordon A. Riley. “Absorption of nutrients by aquatic plants”. In: (1952).
[98] Moeto Nagai, Masamichi Oishi, Marie Oshima, Hiroshi Asai, and Hiroyuki Fujita.
“Three-dimensional two-component velocity measurement of the flow field induced by the
Vorticella picta microorganism using a confocal microparticle image velocimetry technique”. In:
Biomicrofluidics 3.1 (2009), p. 014105.
97



[99] Gabela Nelson, Alexis Strain, Atsuko Isu, Alireza Rahnama, Ken-ichi Wakabayashi,
Adam T Melvin, and Naohiro Kato. “Cells collectively migrate during ammonium chemotaxis in
Chlamydomonas reinhardtii”. In: Scientific Reports 13.1 (2023), p. 10781.
[100] Philip Charles Nelson and Philip Nelson. Biological physics. WH Freeman New York, 2004.
[101] Anja Nenninger, Giulia Mastroianni, and Conrad W Mullineaux. “Size dependence of protein
diffusion in the cytoplasm of Escherichia coli”. In: Journal of bacteriology 192.18 (2010),
pp. 4535–4540.
[102] Hoa Nguyen, MAR Koehl, Christian Oakes, Greg Bustamante, and Lisa Fauci. “Effects of cell
morphology and attachment to a surface on the hydrodynamic performance of unicellular
choanoflagellates”. In: Journal of the Royal Society Interface 16.150 (2019), p. 20180736.
[103] Lasse Tor Nielsen, Seyed Saeed Asadzadeh, Julia Dölger, Jens H Walther, Thomas Kiørboe, and
Anders Andersen. “Hydrodynamics of microbial filter feeding”. In: Proceedings of the National
Academy of Sciences 114.35 (2017), pp. 9373–9378.
[104] Lowell E Noland and Harold Eugene Finley. “Studies on the taxonomy of the genus Vorticella”. In:
Transactions of the American Microscopical Society 50.2 (1931), pp. 81–123.
[105] Martin A Nowak. “Five rules for the evolution of cooperation”. In: science 314.5805 (2006),
pp. 1560–1563.
[106] Laura Wegener Parfrey, Daniel JG Lahr, and Laura A Katz. “The dynamic nature of eukaryotic
genomes”. In: Molecular biology and evolution 25.4 (2008), pp. 787–794.
[107] Rachel E Pepper, Emily E Riley, Matthieu Baron, Thomas Hurot, Lasse Tor Nielsen, MAR Koehl,
Thomas Kiørboe, and Anders Andersen. “The effect of external flow on the feeding currents of
sessile microorganisms”. In: Journal of the Royal Society Interface 18.175 (2021), p. 20200953.
[108] Rachel E Pepper, Marcus Roper, Sangjin Ryu, Paul Matsudaira, and Howard A Stone. “Nearby
boundaries create eddies near microscopic filter feeders”. In: Journal of The Royal Society Interface
7.46 (2010), pp. 851–862.
[109] Rachel E Pepper, Marcus Roper, Sangjin Ryu, Nobuyoshi Matsumoto, Moeto Nagai, and
Howard A Stone. “A new angle on microscopic suspension feeders near boundaries”. In:
Biophysical journal 105.8 (2013), pp. 1796–1804.
[110] Michala E Pettitt, Belinda AA Orme, John R Blake, and Barry SC Leadbeater. “The hydrodynamics
of filter feeding in choanoflagellates”. In: European Journal of Protistology 38.4 (2002), pp. 313–332.
[111] Richard W Pierce and Jefferson T Turner. “Ecology of planktonic ciliates in marine food webs”.
In: Rev. Aquat. Sci 6.2 (1992), pp. 139–181.
[112] Edward M Purcell. “Life at low Reynolds number”. In: American journal of physics 45.1 (1977),
pp. 3–11.
98



[113] Deepa Rajan, Tatyana Makushok, Asa Kalish, Lilibeth Acuna, Alex Bonville,
Kathya Correa Almanza, Brenda Garibay, Eric Tang, Megan Voss, Athena Lin, et al. “Single-cell
analysis of habituation in Stentor coeruleus”. In: Current Biology 33.2 (2023), pp. 241–251.
[114] Mireya Ramírez-Ballesteros, Alfonso Lugo-Vázquez, and Rosaura Mayén-Estrada. “Some
symbiotic, ecological, cytological and geographical distribution aspects of Opercularia articulata
(Operculariidae: Peritrichia) including new data from Mexico”. In: Biologia 76.5 (2021),
pp. 1501–1508.
[115] John A Raven. “Allometry and stoichiometry of unicellular, colonial and multicellular
phytoplankton”. In: New phytologist 181 (2008), pp. 295–309.
[116] Mads Rode, Thomas Kiørboe, and Anders Andersen. “Feeding flow and membranelle filtration in
ciliates”. In: Physical Review Fluids 7.2 (2022), p. 023102.
[117] Mads Rode, Giulia Meucci, Kristian Seegert, Thomas Kiørboe, and Anders Andersen. “Effects of
surface proximity and force orientation on the feeding flows of microorganisms on solid
surfaces”. In: Physical Review Fluids 5.12 (2020), p. 123104.
[118] Marcus Roper, Mark J Dayel, Rachel E Pepper, and MAR Koehl. “Cooperatively generated
stresslet flows supply fresh fluid to multicellular Choanoflagellate colonies”. In: Physical review
letters 110.22 (2013), p. 228104.
[119] Sangjin Ryu, Rachel E Pepper, Moeto Nagai, and Danielle C France. “Vorticella: a protozoan for
bio-inspired engineering”. In: Micromachines 8.1 (2016), p. 4.
[120] Mark J Schnitzer and Steven M Block. “Kinesin hydrolyses one ATP per 8-nm step”. In: Nature
388.6640 (1997), pp. 386–390.
[121] Shashank Shekhar, Hanliang Guo, Sean P Colin, Wallace Marshall, Eva Kanso, and
John H Costello. “Cooperative hydrodynamics accompany multicellular-like colonial
organization in the unicellular ciliate Stentor”. In: bioRxiv (2023).
[122] Evelyn B Sherr and Barry F Sherr. “Significance of predation by protists in aquatic microbial food
webs”. In: Antonie Van Leeuwenhoek 81 (2002), pp. 293–308.
[123] Martin B Short, Cristian A Solari, Sujoy Ganguly, Thomas R Powers, John O Kessler, and
Raymond E Goldstein. “Flows driven by flagella of multicellular organisms enhance long-range
molecular transport”. In: Proceedings of the National Academy of Sciences 103.22 (2006),
pp. 8315–8319.
[124] Mark M Slabodnick and Wallace F Marshall. “Stentor coeruleus”. In: Current Biology 24.17 (2014),
R783–R784.
[125] Vladimír Sládecek. “Indicator value of the genus Opercularia (Ciliata)”. In: Hydrobiologia 79.3
(1981), pp. 229–232.
[126] MA Sleigh and D Barlow. “Collection of food by Vorticella”. In: Transactions of the American
Microscopical Society (1976), pp. 482–486.
99



[127] Cristian A Solari, Sujoy Ganguly, John O Kessler, Richard E Michod, and Raymond E Goldstein.
“Multicellularity and the functional interdependence of motility and molecular transport”. In:
Proceedings of the National Academy of Sciences 103.5 (2006), pp. 1353–1358.
[128] Lucas J Stal and Mariana Silva Cretoiu. The marine microbiome. Springer, 2022.
[129] Roman Stocker. “Marine microbes see a sea of gradients”. In: science 338.6107 (2012), pp. 628–633.
[130] Jean-Baptiste Thomazo, Benjamin Le Révérend, Lea-Laetitia Pontani, Alexis M Prevost, and
Elie Wandersman. “A bending fluctuation-based mechanism for particle detection by ciliated
structures”. In: Proceedings of the National Academy of Sciences 118.31 (2021), e2020402118.
[131] Laura RP Utz, Taiz LL Simao, Lucia SL Safi, and Eduardo Eizirik. “Expanded phylogenetic
representation of genera Opercularia and Epistylis sheds light on the evolution and higher-level
taxonomy of peritrich ciliates (Ciliophora: Peritrichia)”. In: Journal of Eukaryotic Microbiology
57.5 (2010), pp. 415–420.
[132] Franco Verni and Paolo Gualtieri. “Feeding behaviour in ciliated protists”. In: Micron 28.6 (1997),
pp. 487–504.
[133] Kay Vopel, Christian H Reick, Günter Arlt, Martina Pöhn, and Jörg A Ott. “Flow
microenvironment of two marine peritrich ciliates with ectobiotic chemoautotrophic bacteria”.
In: Aquatic Microbial Ecology 29.1 (2002), pp. 19–28.
[134] Kirsty Y Wan, Sylvia K Hürlimann, Aidan M Fenix, Rebecca M McGillivary, Tatyana Makushok,
Evan Burns, Janet Y Sheung, and Wallace F Marshall. “Reorganization of complex ciliary flows
around regenerating Stentor coeruleus”. In: Philosophical Transactions of the Royal Society B
375.1792 (2020), p. 20190167.
[135] Hava Wandel and Roi Holzman. “Modulation of Cilia Beat Kinematics Is a Key Determinant of
Encounter Rate and Selectivity in Tintinnid Ciliates”. In: Frontiers in Marine Science 9 (2022),
p. 845903.
[136] Zhaohui Wang and Michael P Sheetz. “The C-terminus of tubulin increases cytoplasmic dynein
and kinesin processivity”. In: Biophysical journal 78.4 (2000), pp. 1955–1964.
[137] Phillip C Watts, Nina Lundholm, Sofia Ribeiro, and Marianne Ellegaard. “A century-long genetic
record reveals that protist effective population sizes are comparable to those of macroscopic
species”. In: Biology Letters 9.6 (2013), p. 20130849.
[138] Christopher R Wood, Robert Hard, and Todd M Hennessey. “Targeted gene disruption of dynein
heavy chain 7 of Tetrahymena thermophila results in altered ciliary waveform and reduced swim
speed”. In: Journal of cell science 120.17 (2007), pp. 3075–3085.
[139] Alexandra Z Worden, Michael J Follows, Stephen J Giovannoni, Susanne Wilken,
Amy E Zimmerman, and Patrick J Keeling. “Rethinking the marine carbon cycle: factoring in the
multifarious lifestyles of microbes”. In: Science 347.6223 (2015), p. 1257594.
100



[140] Jacek K Wróbel, Ricardo Cortez, Douglas Varela, and Lisa Fauci. “Regularized image system for
Stokes flow outside a solid sphere”. In: Journal of Computational Physics 317 (2016), pp. 165–184.
[141] Jonathan P Zehr, Joshua S Weitz, and Ian Joint. “How microbes survive in the open ocean”. In:
Science 357.6352 (2017), pp. 646–647.
[142] Guangpu Zhu, Wen-Zhen Fang, and Lailai Zhu. “Optimizing low-Reynolds-number predation via
optimal control and reinforcement learning”. In: Journal of Fluid Mechanics 944 (2022), A3.
[143] Bogumila E Zima-Kulisiewicz and Antonio Delgado. “Synergetic microorganismic convection
generated by Opercularia asymmetrica ciliates living in a colony as effective fluid transport on
the micro-scale”. In: Journal of biomechanics 42.14 (2009), pp. 2255–2262.
101



Appendix A
Notes on functions in spherical coordinates
A.1 Gegenbauer differential equation of degree ±
1
2
The Gegenbauer differential equation is
(1 − x
2
)y
′′ − (2α + 1)xy′ + n(n + 2α)y = 0. (A.1)
When α =
1
2
, we have the Gegenbauer differential equation of degree 1
2
, the Legendre differential equaiton,
(1 − x
2
)y
′′ − 2xy′ + n(n + 1)y = 0. (A.2)
The solution to Legendre differential equation is
y(x) = AnPn(x) + BnQn(x), (A.3)
where n is integer, An, Bn are constant, Pn(x)is Legendre function of the first kind, and Qn(µ)is Legendre
function of the second kind.
102



For Legendre polynomial satisfying the differential equation, we write the Legendre differential equation
as
d
dx
(1 − x
2
)
dPm
dx
+ m(m + 1)Pm = 0, (A.4)
where m is integer. Integral both sides of Legendre differential equation, we have
(1 − x
2
)
dPm
dx + m(m + 1) Z
Pmdx = 0. (A.5)
Set m = n − 1 and substitute into above equation, we have
(1 − x
2
)
dPn−1
dx + (n − 1)n
Z
Pn−1 = 0. (A.6)
Set function Fn(x) satisfying Fn(x) = −
R
Pn−1(x)dx, and substitute into above equation, we obtain
(1 − x
2
)
d
2Fn
d
2x
+ n(n − 1)Fn = 0. (A.7)
Notice that the above equation is Gegenbauer differential equation of degree −
1
2 when α = −
1
2
. The
solution is
y(x) = CnFn(x) + DnHn(x), (A.8)
where n is integer, Cn, Dn are constant, Fn(x) is solution function of the first kind, and Hn(x) is the
solution function of the second kind [51], with defining
F0(µ) = 1, F1(µ) = −µ,
H0(µ) = −µ, H1(µ) = −1.
(A.9)
103



Properties of the listed functions are listed below. The orthogonality property of Legendre polynomials is
Z 1
−1
Pn(µ)Pm(µ)dµ =
2
2n + 1
δnm. (A.10)
The differential relation and recurrence relation of Legendre polynomials are
(2n + 1)Pn(µ) = d
dµ(Pn+1(µ) − Pn−1(µ)),
(n + 1)Pn+1(µ) = (2n + 1)µPn(µ) − nPn−1(µ).
(A.11)
The orthogonality property of Fn(µ) is
Z +1
−1
Fn(µ)Fm(µ)
1 − µ2
dµ =
2
n(n − 1)(2n − 1)δnm. (A.12)
By using differential relations, we have the link to Legendre polynomials as
Fn(µ) = Pn−2(µ) − Pn(µ)
2n − 1
,
Hn(µ) = Qn−2(µ) − Qn(µ)
2n − 1
.
(A.13)
A.2 Envelope model with radial surface velocity
Consider the surface velocity of the swimming squirmer has both radial and tangential components. The
boundary conditions are
u(r = a) = X∞
n=1
AnPn(µ)er +
X∞
n=1
BnVn(µ)eθ,
u(r → ∞) = −Uez,
(A.14)
104



By substituting the far-field boundary condition into the general solution of fluid velocity (2.7), we have
C0 = −U. Substitute C0 into general solution, the constant coefficients Cn can be solved by solving the
following system of equations
−U +
C3
a
3
+
C1
a
= A1,
U +
C3
2a
3
−
C1
2a
= B1,
Cn+2
a
n+2 +
Cn
a
n
= An,
nCn+2
2a
n+2 +
(n − 2)Cn
2a
n
= Bn (n ≥ 2).
(A.15)
The coefficients Cn can be found as
C1 =
a
2
(3U + A1 − 2B1),
C3 =
a
3
2
(−U + A1 + 2B1),
Cn = (n
2
An − Bn)a
n
,
Cn+2 = (Bn −
n − 2
2
An)a
n+2
.
(A.16)
Thus, the fluid velocity can be obtained as
ur(r, µ) =
− U + (−U + A1 + 2B1)
a
3
2r
3
+ (3U + A1 − 2B1)
a
2r

P1(µ)
+
X∞
n=2

(Bn −
n − 2
2
An)
a
n+2
r
n+2 + (n
2
An − Bn)
a
n
r
n

Pn(µ),
uθ(r, µ) =
U + (−U + A1 + 2B1)
a
3
4r
3
− (3U + A1 − 2B1)
a
4r

V1(µ)
+
X∞
n=2

(Bn −
n − 2
2
An)
nan+2
2r
n+2 + (n
2
An − Bn)
(n − 2)a
n
2r
n

Vn(µ).
(A.17)
105



Substitute the above coefficients into the general hydrodynamic force expression with satisfying zero net
hydrodymaic force condition F =
R
σ · ndS = 0. We can obtain the swimming velocity as
U =
1
3
(2B1 − A1). (A.18)
A.3 Bessel functions and spherical Bessel functions
The Bessel differential equation is
xy′′ + xy′ + (x
2 − α
2
)y = 0. (A.19)
The solutions are Bessel function of the first kind Jα(x) and the second kind Yα(x). These two linearly
independent solutions to Bessel’s equation can be formulated as Hankel functions of the first and second
kind, defined as
H(1)
n
(x) = Jn(x) + iYn(x),
H(2)
n
(x) = Jn(x) − iYn(x).
(A.20)
These are also called Bessel function of the third kind.
When x is complex number, we have modified Bessel differential equation
xy′′ + xy′ − (x
2 + α
2
)y = 0. (A.21)
106



And the solutions are modified Bessel function of the first kind Kα(x) and the second kind Iα(x). Function
Kα can be expressed as Hankel functions
Kn(x) = π
2
i
α+1H(1)
α (ix), −π < arg x ≤
π
2
,
Kn(x) = π
2
(−i)
α+1H(2)
α (−ix), −
π
2
< arg x ≤ π,
(A.22)
In spherical coordinates, for example, from Helmholtz differential equation, we have spherical Bessel equation
xy′′ + 2xy′ + (x
2 − n(n + 1))y = 0. (A.23)
The solutions are spherical Bessel function of the first kind jn(x) and second kind yn(x). The relations to
ordinary Bessel functions are
jn(x) = r
π
2x
Jn+ 1
2
(x),
yn(x) = r
π
2x
Yn+ 1
2
(x).
(A.24)
Similarly, we have the spherical Hankel functions are
h
(1)
n
(x) = jn(x) + iyn(x),
h
(2)
n
(x) = jn(x) − iyn(x).
(A.25)
In heat equation, we have modified spherical Bessel differential equation
xy′′ + 2xy′ − (x
2 + n(n + 1))y = 0. (A.26)
107



The solutions are modified spherical Bessel function of the first kind kn(x) and second kind in(x). The
relation between kn(x) to spherical Hankel functions are
kn(x) = −
1
2
πinh
(1)
n
(ix), −π < arg x ≤
π
2
,
kn(x) = −
1
2
πi−nh
(2)
n
(−ix), −
π
2
< arg x ≤ π,
(A.27)
108



Appendix B
Numerical methods
B.1 Adjoint-based optimization
To search for optimal surface motions that maximize feeding in the sessile sphere model, we considered
an optimization method based on variational analysis and steepest ascent [90]. The problem consists of a
PDE-constrained optimization problem, where the goal is to find optimal Bn that maximize the Sherwood
number, subject to the concentration field c satisfying the advection-diffusion equation and surface velocity
satisfying the constant energy constraint,
max
Bn,c
Sh(Bn, c),
subject to L[c] = 0 and X
N
n=1
2B2
n
n(n + 1) = 1.
(B.1)
Here, the linear operator L = Peu · ∇ − ∇2
is that of the advection-diffusion equation along with the
corresponding boundary conditions.
We use a variational approach to derive an adjoint system of equation. Given a surface motion, we
consider small variations δBn in the coefficients associated with each mode. The corresponding small
109



variations in the velocity field and concentration field are given by δu and δc, and result in a variation in
Sh number 4.7
δSh = −
1
4π
Z
S
∇δc · nˆdS. (B.2)
The concentration variation δc must satisfy
Pe(u + δu) · ∇(c + δc) = ∇2
(c + δc), (c + δc)(r = 1) = 1, (c + δc)(r → ∞) = 0. (B.3)
Subtracting the PDE for c and keeping the leading order in δc, we arrive at
Pe(δu · ∇c + u · ∇δc) = ∇2
δc, δc(r = 1) = 0, δc(r → ∞) = 0. (B.4)
Multiplying the above equation with a test function g(r, µ) and integrating over the entire fluid domain,
we get
Pe Z
V
(gδu · ∇c + gu · ∇δc)dV =
Z
V
g∇2
δc dV. (B.5)
Using integration by parts and standard vector calculus identities, together with the appropriate boundary
conditions and continuity property of the fluid, we obtain
−Pe Z
V
cδu · ∇g dV +
Z
S
(g∇δc) · nˆdS =
Z
V
δc(Peu · ∇g + ∇2
g)dV, (B.6)
110



Following a standard argument, we get that the test function g must satisfy the following partial differential
equation and boundary conditions
Peu · ∇g + ∇2
g = 0, g(r = 1) = 1, g(r → ∞) = 0, (B.7)
and the consistency equation
Z
S
∇δc · nˆdS = Pe Z
V
cδu · ∇g dV. (B.8)
From (B.2) and (B.8), we get that the variation in Sh number is given by
δSh = −
1
4π
Z
S
∇δc · nˆdS = −
Pe
4π
Z
V
cδu · ∇g dV = −
Pe
2
Z ∞
1
Z 1
−1
c(δu · ∇g) r
2
drdµ. (B.9)
We expand the test function g(r, µ) in terms of Legendre polynomials in µ,
g(r, µ) = X∞
m=0
Gm(r)Pm(µ), (B.10)
and substitute back into (B.9) with the perturbation from surface velocity δu(r = 1) = (δBn)Vneθ. We
arrive at an expression for the gradient of nutrient uptake at each mode,
δSh
δBn
= −Pe X∞
m=0
X∞
k=0
1
2k + 1 Z ∞
1

CkfnrG
′
mEmnk −
CkGmfnθ
r
Fmnk
r
2
dr. (B.11)
We now consider a finite number of velocity modes and express the input surface velocity as βnVn(µ),
where βn =
q 2
n(n+1)Bn. The weighted coefficients βn associated with each velocity mode must satisfy
the constraint on the energy dissipation rate, that is, P
n
β
2
n = 1.
111



Starting from an initial vector β
0 = (β1, β2, . . . , βn, . . .) of weighted coefficients, our goal is to find the
optimal value of β that simultaneously maximize Sh and satisfy the constraint ∥β
(j)
∥ = 1 at each iteration
j in the optimization process. Thus, to get the value of β
(j)
at subsequent iterations j > 0, we project the
feeding gradient onto the constraint space ∥β
(j)
∥ = 1 using the linear projection (I −β
(j) ⊗β
(j)
), where
I is the identity matrix. That is, the steepest ascent direction of Sh with respect to weighted coefficients
βn at the j
th iteration is given by
∇dSh = ∇βSh − (β
(j)
· ∇βSh)β
(j)
. (B.12)
The optimization process consists of updating β
(j+1) following the gradient ascent direction, where α is
step size that can be adjusted in each iteration.
β
(j+1) =
β
(j) + α∇dSh
∥β
(j) + α∇dSh∥
. (B.13)
B.2 Legendre spectral method
B.2.1 Steady-state advection-diffusion equation
We solve (4.6) using the Legendre spectral method; see, e.g., [90]. To this end, we first expand the concentration field c(r, µ) in terms of the Legendre polynomials Pm(µ)
c(r, µ) = X∞
m=0
Cm(r)Pm(µ), (B.14)
where Cm(r) are unknown coefficients associated with Legendre basis functions Pm(µ). We substitute
(B.14) into (4.6) and project the governing equation onto the k
th Legendre polynomial Pk(µ), to arrive at
112



an infinite set of coupled boundary-valued ordinary differential equations for the unknown coefficients
Ck(r), k = 0, . . . , ∞
PeX∞
n=1
X∞
m=0
Bn

Emnkfnr
∂Cm
∂r − Fmnkfnθ
Cm
r

=
∂
2Ck
∂r2
+
2
r
∂Ck
∂r −
k(k + 1)
r
2
Ck,
Ck|
r=a=1 = δ0k, Ck|
r=∞ = 0.
(B.15)
Here, the coefficients Emnk and Fmnk are obtained by projection, using the orthogonality property of the
Legendre polynomials,
Emnk =
2k + 1
2
Z 1
−1
PnPmPkdµ, Fmnk =
(2k + 1)
2n(n + 1) Z 1
−1
(1 − µ
2
)P
′
nP
′
mPkdµ. (B.16)
The terms fnr and fnθ are the r components of the flow field
Sessile ciliated sphere:fnr =

1
r
n+2 −
1
r
n

, fnθ =

n
r
n+2 −
n − 2
r
n

,
Motile ciliated sphere: f1r =
2
3

−1 +
1
r
3

, f1θ =
2
3

2 +
1
r
3

,
fnr =

1
r
n+2 −
1
r
n

(n ≥ 2), fnθ =

n
r
n+2 −
n − 2
r
n

(n ≥ 2),
(B.17)
In the numerical calculation, we truncate the number of modes in the expansion (4.23) of the concentration field to account for a finite number M of modes. To reach far-field boundary condition, we used a nonuniform radial mesh such that the grid is denser near the sphere surface and more sparse in the far-field.
Specifically, we used the exponential function r = e
s(ζ)
, where ζ ∈ [0, 1] and s(ζ) = w1ζ
3 + w2ζ
2 + w3ζ
is a third order polynomial of ζ with constants w1, w2, w3 chosen for getting converged results. To express (B.15) in terms of the transformed variable ζ, we used the chain rule
dCk
dr =
dCk
dζ
dζ
dr ,
d
2Ck
dr2
=

dζ
dr 2
d
2Ck
dζ2
+
d
2
ζ
dr2
dCk
dζ . (B.18)
113



Considering N velocity modes and M concentration modes, the differential equations in (B.15) can be
rewritten as
PeX
N
n=1
X
M
m=0
Bn

Emnkfnr
dCm
dζ
dζ
dr − Fmnkfnθ
Cm
r

−(
dζ
dr )
2
d
2Ck
dζ2
−
d
2
ζ
dr2
dCk
dζ −
2
r
dζ
dr
dCk
dζ +
k(k + 1)
r
2
Ck = 0.
(B.19)
We discretized the spatial derivatives using the center difference scheme
dCm,j
dζ =
Cm,j+1 − Cm,j−1
2△ζ
,
d
2Cm,j
dζ2
=
Cm,j+1 − 2Cm,j + Cm,j−1
(△ζ)
2
.
(B.20)
We computed the concentration field and Sherwood number for various Péclet number and tested the
convergence of the results as a function of mesh size △ζ, computational domain R, and the number of
modes M in the concentration expansion at Pe = 100. To test the effect of mesh size on the convergence of
our simulations, we fixed the number of modes M and computational domain size R, and varied the mesh
size as δζ = [ 1
100 ,
1
200 ,
1
400 ,
1
800 ]. We computed the relative error of Sherwood number as a function of mesh
size δζ; the results converged as the mesh size got smaller. Consistent with the second order accuracy of
the discretization (B.20), we obtained a convergence rate close to 2. When testing the convergence as a
function of the number of modes M in the concentration expansion and the computational domain size R,
we found that a higher number of concentration modes and larger computational domain size are needed
for higher Péclet numbers. Also, because the concentration field becomes more front and back asymmetric
as Péclet increases, a denser mesh is required to capture the rapid change near the surface.
114



B.2.2 Numerical scheme
Use central difference, and choose M number of expansions in the spectrum and N number of modes in
the velocity field, and consider the number of grid points in r space is Nr. The governing equation is
PeX
N
n=1
X
M
m=0
Bn

Emnkfnr
Cm,j+1 − Cm,j−1
2△ζ
dζ
dr − Fmnkfnθ
Cm,j
r

− (
dζ
dr )
2 Ck,j+1 − 2Ck,j + Ck,j−1
(△ζ)
2
−
d
2
ζ
dr2
Ck,j+1 − Ck,j−1
2△ζ
−
2
r
dζ
dr
Ck,j+1 − Ck,j−1
△ζ
+
k(k + 1)
r
2
Ck,j = 0,
Ck,1 = δ1k, Ck,Nr = 0,
(B.21)
This equation can be solved by solving a linear system Ax = b, where A represents a matrix, and b and x
are vectors. The unknown vector is Ck, which has size M × Nr. The matrix A has M2 number of block
matrices, and each block matrix has size N2
r
. Write coefficients of every point in space along one row of
block matrices, the coefficients of block matrices in k row are
PeX
N
n=1
BnfnrEmnk
dζ
dr
1
2△ζ
−

1
(△ζ)
2
(
dζ
dr )
2 +
d
2
ζ
dr2
1
2△ζ
+
2
r
dζ
dr
1
2△ζ

δmk
Cm,j+1
−PeX
N
n=1
Bn
fnθ
r
Fmnk +

(
dζ
dr )
2
2
(△ζ)
2
+
k(k + 1)
r
2

δmk
Cm,j
−PeX
N
n=1
BnfnrEmnk
dζ
dr
1
2△ζ
−

1
(△ζ)
2
(
dζ
dr )
2 −
d
2
ζ
dr2
1
2△ζ
−
2
r
dζ
dr
1
2△ζ

δmk
Cm,j−1
(B.22)
For every k, we have same above expression for each block row. Thus, we have the matrix AM2×N2
r
. And
the b vector is
b = [−1, 0, . . . , 0]T (B.23)
115



The liner system Ax = b is




1
.
.
.
.
.
.
.
.
.
1


N2
r


0
.
.
.
.
.
.
.
.
.
0


N2
r
. . .


0
.
.
.
.
.
.
.
.
.
0


N2
r


1
.
.
.
.
.
.
.
.
.
1


N2
r
. . .
. . .


M2×Nr2




C0




C1


.
.
.


=




1
0
.
.
.
0




0
.
.
.
.
.
.


.
.
.
.
.
.


(B.24)
The size of the matrix in the linear system is (M × Nr)
2
; the size of the unknown vector and the RHS is
(M × Nr) × 1. Thus, all modes in concentration decomposition are solved simultaneously.
B.2.3 Unsteady advection-diffusion equation
We use a spectral method based to numerically solve 7.4 for the unsteady concentration field c
′ = c(t, r, µ)−
(1 + ϵz). Specifically, following [90, 79, 78], we expand the solution in terms of Legendre polynomials
c
′
(t, r, µ) = X∞
m=0
Cm(t, r)Pm(µ). (B.25)
116



We substitute this expansion into (7.4), and project on the k
th Legendre polynomial. Using properties of
Legendre polynomials, we get
∂Ck(t, r)
∂t =
1
r
2

∂
∂r

r
2 ∂Ck(t, r)
∂r

− n(n + 1)Ck(t, r)

,
Ck|t=0 = 0,
Ck|r=1 = −1δ0k − ϵδ1k,
dCk
dr ]



r→∞
= 0.
(B.26)
Only the first two terms C0(t, r) and C1(t, r) survive. We discrete the space r using a finite difference
method and time t using a Backward Euler method for C0(t, r) and C1(t, r). Finally, given c
′
in terms of
C0(t, r) and C1(t, r), we add the linear concentration background 1 + ϵz to recover the full concentration
field c(t, r, µ).
We use a similar spectral method to solve the advection-diffusion equations in (7.34) and (7.35).
Namely, we project ch and c˜ onto the Legendre polynomials Pm(µ) to obtain the expansions
ch(r, µ, t) = X∞
m=0
Cm(r, t)Pm(µ), (B.27)
c˜(r, µ, t) = X∞
m=0
C˜m(r, t)Pm(µ). (B.28)
Substitute (B.27) into (7.34) to get
∂Ck
∂t +
X∞
n=1
X∞
m=0
Bn

Emnkfnr
∂Cm
∂r − Fmnkfnθ
Cm
r

=
1
Pe
∂
2Ck
∂r2
+
2
r
∂Ck
∂r −
k(k + 1)
r
2
Ck

,
Ck|t=0 = 0, Ck|r=1 = −δ0k,
dCk
dr



r→∞
= 0.
(B.29)
117



Similarly, substitute (B.28) into (7.35) to get
∂C˜
k
∂t +
X∞
n=1
X∞
m=0
Bn
"
fnr
Emnk
∂C˜m
∂r + E1nk
− fnθ
Fmnk
C˜m
r
+ F1nk#
=
1
Pe
∂
2C˜
k
∂r2
+
2
r
∂C˜
k
∂r −
k(k + 1)
r
2
C˜
k
!
,
C˜
k|t=0 = 0, C˜
k|r=1 = −δ1k,
dC˜
k
dr



r→∞
= 0.
(B.30)
The terms fnr and fnθ related to the flow field are given by
fnr =

1
r
n+2 −
1
r
n

, fnθ =

n
r
n+2 −
n − 2
r
n

. (B.31)
The coefficients Emnk and Fmnk are given by
Emnk =
2k + 1
2
Z 1
−1
PnPmPkdµ, Fmnk =
(2k + 1)
2n(n + 1) Z 1
−1
(1 − µ
2
)P
′
nP
′
mPkdµ. (B.32)
Equations (B.29) and (B.30) are then discretized in space and time. Following [90, 79], we use a central
difference discretization in space and backward Euler in time. We compute the time evolution of the concentration field and associated nutrient uptake. The latter, using properties of the Legendre polynomials,
can be expressed as
J(t) = Jh(t) + ϵJ˜(t) =
∂C0
∂r + ϵ
∂C˜
0
∂r !
r=1
. (B.33)
Using a similar spatial discretization scheme and a Backward Euler method for the time-evolving term, the
unsteady concentration field can be solved numerically.
118 
Asset Metadata
Creator Liu, Jingyi (author) 
Core Title Flow physics and nutrients transport in microscopic living systems 
Contributor Electronically uploaded by the author (provenance) 
School Andrew and Erna Viterbi School of Engineering 
Degree Doctor of Philosophy 
Degree Program Mechanical Engineering 
Degree Conferral Date 2025-05 
Publication Date 04/28/2025 
Defense Date 03/13/2025 
Publisher University of Southern California. Libraries (digital), University of Southern California (Los Angeles, California, USA) (original) 
Tag advection-diffusion,ciliated microorganisms,OAI-PMH Harvest,Stokes flow 
Format theses (aat) 
Language English
Advisor Kanso, Eva (committee chair),  Luhar, Mitul (committee member), Bermejo-Moreno, Ivan (committee member), Haselwandter, Christoph (committee member) 
Creator Email liu027@usc.edu 
Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC11399KF04 
Unique identifier UC11399KF04 
Legacy Identifier etd-Liu-36919-47784 
Document Type Dissertation 
Format theses (aat) 
Rights Liu, Jingyi 
Internet Media Type application/pdf 
Type texts
Source 202504w05-usctheses (batch), University of Southern California Dissertations and Theses (collection), University of Southern California (contributing entity) 
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. 
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 uscdl@usc.edu
Abstract (if available)
Abstract Nutrient acquisition strategies are fundamental to the survival and ecological roles of microorganisms. In aquatic ecosystems, ciliated microbes have evolved specialized mechanisms to optimize nutrient uptake by actively modifying their surrounding fluid. These organisms either swim freely or attach to surfaces, generating feeding currents through the coordinated motion of cilia. The effectiveness of these strategies is shaped by hydrodynamic constraints, yet the extent to which flow physics governs feeding efficiency and ecological function remains an open question. To understand their foraging strategies, I explore three key questions: (1) the benefits and limitations of swimming versus attachment for nutrient acquisition, (2) the constraints imposed by flow physics and nutrient transport on ciliate morphology, and (3) the effects of environmental variations, such as nutrient gradients, on feeding efficiency. I solve the advection-diffusion equation governing nutrient transport in cilia-driven flows using numerical and analytical techniques, including spectral methods and adjoint-based optimization. The findings reveal that, despite their contrasting behaviors, swimming and attached ciliates can achieve hydrodynamically equivalent feeding rates, resolving a long-standing debate on the efficiency of these two strategies. Furthermore, flow physics drive the functional design of ciliate morphology, shaping their ecological adaptations. Finally, the results suggest that ciliates may employ chemical sensing to dynamically orient themselves for enhanced feeding in non-uniform nutrient environments. These findings open new directions for studying microswimmer locomotion in dynamic chemical landscapes and investigating adaptive motion in response to environmental cues. 
Tags
Stokes flow
advection-diffusion
ciliated microorganisms
Linked assets
University of Southern California Dissertations and Theses
doctype icon
University of Southern California Dissertations and Theses 
Action button