1. Introduction
Eukaryotic cilia are found in large numbers across various environments, including microorganisms, coral reefs, mucociliary layers and fallopian tubes (Smith, Gaffney & Blake Reference Smith, Gaffney and Blake2008; Shapiro et al. Reference Shapiro, Fernandez, Garren, Guasto, Debaillon-Vesque, Kramarsky-Winter, Vardi and Stocker2014; Pedley, Brumley & Goldstein Reference Pedley, Brumley and Goldstein2016; Mathijssen et al. Reference Mathijssen, Doostmohammadi, Yeomans and Shendruk2016), where they perform a suite of functions such as swimming, feeding and transportation of dissolved and particulate matter. A great many models have been employed to capture both the transient and steady-state behaviour of cilia in different geometric configurations, which vary from detailed descriptions of beating filaments (Gueron et al. Reference Gueron, Levit-Gurevich, Liron and Blum1997; Elgeti & Gompper Reference Elgeti and Gompper2013; Goldstein et al. Reference Goldstein, Lauga, Pesci and Proctor2016) and minimal models for the associated flows (Vilfan & Jülicher Reference Vilfan and Jülicher2006; Niedermayer, Eckhardt & Lenz Reference Niedermayer, Eckhardt and Lenz2008; Brumley et al. Reference Brumley, Polin, Pedley and Goldstein2012, Reference Brumley, Wan, Polin and Goldstein2014), through to the celebrated envelope, or ‘squirmer’, description of collective ciliary motion (Lighthill Reference Lighthill1952; Blake Reference Blake1971; Pedley Reference Pedley2016). The latter abstracts a ciliated surface as a no-slip boundary undergoing small-amplitude deformations, so it fails to resolve the flow field of individual cilia. On the other hand, by representing the ciliary tip as a sphere traversing a circular orbit above a no-slip wall, minimal models can capture the transient flow field generated by the cilium, but this calculation becomes relatively time-consuming for a large ciliary array. A computationally efficient means of calculating both near- and far-field ciliary flows was provided recently by Selvan et al. (Reference Selvan, Duck, Pihler-Puzović and Brumley2023) who demonstrated that the action of a point torque accurately resolves the total steady-state flow field of a single cilium.
Ciliary flows routinely transport both particulate and dissolved matter, as well as microswimmers in a range of model systems, including corals (Shapiro et al. Reference Shapiro, Fernandez, Garren, Guasto, Debaillon-Vesque, Kramarsky-Winter, Vardi and Stocker2014; Pacherres et al. Reference Pacherres, Ahmerkamp, Schmidt-Grieb, Holtappels and Richter2020; Murthy, Picioreanu & Kühl Reference Murthy, Picioreanu and Kühl2023), volvocine green algae (Short et al. Reference Short, Solari, Ganguly, Powers, Kessler and Goldstein2006), starfish larvae (Gilpin, Prakash & Prakash Reference Gilpin, Prakash and Prakash2017) and airways (Grotberg Reference Grotberg1994; Ramirez-San Juan et al. Reference Ramirez-San Juan, Mathijssen, He, Jan, Marshall and Prakash2020). In conjunction with ciliary flows, externally imposed fluid flows also help to transport dissolved and particulate matter, thereby influencing photosynthesis (Mass et al. Reference Mass, Genin, Shavit, Grinstein and Tchernov2010), respiration (Patterson, Sebens & Olson Reference Patterson, Sebens and Olson1991; Smith Reference Smith2006) and reproduction (Kölle et al. Reference Kölle, Dubielzig, Reese, Wehrend, König and Kummer2009; Poon et al. Reference Poon, Westwood, Laeverenz-Schlogelhofer, Brodrick, Craggs, Keaveny, Jékely and Wan2023) in different species. Given the scale of collective transport by ciliary arrays, it is common to consider coarse-grained flow for the purposes of mass transport, for example using squirmer models, rather than resolving the flow around individual cilia (Magar, Goto & Pedley Reference Magar, Goto and Pedley2003; Shapiro et al. Reference Shapiro, Fernandez, Garren, Guasto, Debaillon-Vesque, Kramarsky-Winter, Vardi and Stocker2014; Pacherres et al. Reference Pacherres, Ahmerkamp, Koren, Richter and Holtappels2022; Ahmerkamp et al. Reference Ahmerkamp, Jalaluddin, Cui, Brumley, Pacherres, Berg, Stocker, Kuypers, Koren and Behrendt2022). However, the latter is necessary for understanding the nutrient transport within the ciliary layer, and therefore essential for overall flux calculations.
In the present study, we utilise the point torque model by Selvan et al. (Reference Selvan, Duck, Pihler-Puzović and Brumley2023) to swiftly, but accurately, calculate the fluid flow around an individual cilium and use this to quantify mass transport around it. This is done by comparing two different models. The first model, based on Langevin dynamics, assumes that matter is an ensemble of Brownian particles subject to both deterministic and fluctuating forces, and has been employed widely to understand, for example, the transport of bacteria, eggs in an oviduct and particles in porous media (Verdugo et al. Reference Verdugo, Lee, Halbert, Blandau and Tam1980; Tartakovsky, Tartakovsky & Meakin Reference Tartakovsky, Tartakovsky and Meakin2008; Rusconi, Guasto & Stocker Reference Rusconi, Guasto and Stocker2014; Farago & Pontrelli Reference Farago and Pontrelli2020). We also develop a continuum model for examining the nutrient flux between a no-slip substrate and a surrounding fluid, which models transport using an advection–diffusion equation. This approach has been used for investigating mass transport by ciliary flows around corals (Shapiro et al. Reference Shapiro, Fernandez, Garren, Guasto, Debaillon-Vesque, Kramarsky-Winter, Vardi and Stocker2014), the capture rates of diffusive molecules on a cilium (Hickey, Vilfan & Golestanian Reference Hickey, Vilfan and Golestanian2021) and convection in the suspension of non-motile bacteria (Dunstan et al. Reference Dunstan, Lee, Hwang, Park and Goldstein2018). We validate the Langevin model of mass transport in the absence of a background flow using the continuum model, and proceed to investigate the mass transport due to a single cilium in the presence of a linear shear flow using the Langevin model.
The framework of this paper is as follows: in § 2 we derive the integral equation describing the trajectory of a Brownian particle in the presence of both ciliary and background flows, using the Langevin approach; in § 3, we examine the particle distribution and nutrient transport in the absence of a background flow; in § 4, we extend this to study the effect of a background linear shear flow on the particle trajectories and nutrient transport; finally, this is followed by a discussion in § 5.
(a) Schematic diagram of the Langevin model for the mass transport from a shallow disk-shaped source of concentration
$C_0$
in the presence of a ciliary flow (i.e. the flow generated by a point torque
$\boldsymbol{\Omega }$
) and linear shear flow above the rigid wall. (b) The corresponding sketch of the continuum model showing the regularised torque (i.e.
$\boldsymbol {\Omega }=|\boldsymbol {\Omega }|\hat {\boldsymbol {e}}_y$
) of radius
$\kappa$
(
$\approx 0.001 d$
) positioned above the rigid wall inside the box (
$l_b\times w_b \times h_b$
). A circular source of concentration
$C_0$
(situated at
$z=0$
) and a planar sink at
$z=h_b$
are applied in this case.
2. Langevin model
Consider Brownian particles of radius
$a_b$
, emitted from a disk-shaped source of constant concentration
$C_0$
, radius
$r_s$
and height
$h_s$
, which approximates a ciliary patch, in the presence of a single cilium and a background linear shear flow
$\boldsymbol {u}_b$
(figure 1
a). The finite height,
$0\lt h_s \ll 2r_s$
, is chosen so that a constant concentration boundary condition can be maintained numerically, but is sufficiently small so as to represent a surface concentration condition. Following the work of Selvan et al. (Reference Selvan, Duck, Pihler-Puzović and Brumley2023), we capture the time-averaged flow
$\boldsymbol {u}$
generated by a single cilium using a stationary point torque
$\boldsymbol {\Omega }$
positioned at
$\boldsymbol {x}_r=(0,0,d)$
, so that
\begin{align} &u_i=\frac {1}{8\pi \mu }\left[\frac {\epsilon _{ijk} r_k}{r^3}-\frac {\epsilon _{ijk}R_k}{R^3}+2\,\rm d\it \epsilon _{kj3} \left(\frac {\delta _{ik}}{R^3}-\frac {3R_iR_k}{R^5}\right)+\frac {6\epsilon _{kj3}R_iR_kR_3}{R^5}\right] \Omega _j, \end{align}
where
$i,j,k \in \{1,2,3\}$
,
$r=|\boldsymbol {r}|$
,
$\boldsymbol {r}=(r_1,r_2,r_3)=(x,y,z-d)$
,
$\boldsymbol {R}=(R_1,R_2,R_3)=(x,y,z+d)$
,
$R=|\boldsymbol {R}|$
,
$\delta_{ij}$
and
$\epsilon_{ijk}$
are the Kronecker delta and the Levi-Civita symbol, respectively, and
$\mu$
is the dynamic viscosity of the fluid. Then the dynamics of a Brownian particle at a point
$\boldsymbol {x}=(x,y,z)$
emitted from the source in the presence of
$\boldsymbol {u}$
and
$\boldsymbol {u}_b$
are given by
\begin{equation} \frac {{\text{d}}\boldsymbol {x}}{{\text{d}}t}=\boldsymbol {u}+\boldsymbol {u}_b+\sqrt {\frac {2k_BT}{\boldsymbol {f}}} \cdot \boldsymbol {\eta }(t). \end{equation}
Here,
$\boldsymbol {u}_b=\dot {\gamma }z\hat {\boldsymbol {e}}_x$
with
$\dot {\gamma }$
being the shear rate of the background flow and
$\hat {\boldsymbol {e}}_x$
denoting the unit normal vector along the
$x$
-axis. The white noise
$\boldsymbol {\eta }(t)$
satisfies the conditions
$\langle \boldsymbol {\eta }(t) \rangle =0$
and
$\langle \boldsymbol {\eta }(t) \boldsymbol {\eta }(t^{\prime})\rangle =\delta (t-t^{\prime}) \mathbf {I}$
, where
$\mathbf {I}$
is the identity matrix (Ramirez-San Juan et al. Reference Ramirez-San Juan, Mathijssen, He, Jan, Marshall and Prakash2020). The parameters
$k_B$
and
$T$
are the Boltzmann constant and absolute temperature, respectively. Given that the particles could experience additional drag due to the presence of the wall, the friction tensor
$\boldsymbol {f}$
is defined as
$\boldsymbol {f}=6\pi \mu a_b [\boldsymbol {\mathbf {I}}+{(9a_b)}/{(16 z(t))} (\boldsymbol {\mathbf {I}}+\hat {\boldsymbol {e}}_z \hat {\boldsymbol {e}}_z) ]$
, where
$\hat {\boldsymbol {e}}_z$
and
$z(t)$
denote the unit normal vector along the
$z$
-axis, and the height of a Brownian particle above the rigid wall at a given time
$t$
, respectively. The particle’s radius,
$a_b$
, and diffusivity,
$D_0$
, are related through the Stokes–Einstein equation, which can be approximated as
$D_0= k_B T / (6\pi \mu a_b)$
(Miller Reference Miller1924). For a particle that can move by Brownian motion, and through ciliary and linear shear flows, we impose the reflecting boundary condition on the plane rigid wall (
$z=0$
) and maintain open boundaries in all other directions.
Let us assume that the trajectory of a single particle traverses the following positions
$\{\boldsymbol {x}_0=\boldsymbol {x}(t_0),\,\boldsymbol {x}_1=\boldsymbol {x}(t_1),\,\boldsymbol {x}_2=\boldsymbol {x}(t_2),\dots \}$
in the time steps
$\{t_0,\,t_1=t_0+\delta t,\,t_2=t_0+2\delta t,\dots \}$
with
$\boldsymbol {x}(t_0)=\boldsymbol {x}_0$
being the initial position. The particle trajectory can then be obtained on each subinterval
$[t_{n-1},\,t_n]$
for
$n\in \mathbb {N}$
using the following equation derived from (2.2):
\begin{align} \boldsymbol {x}_n=\int _{t_{n-1}}^{t_n} \bigg [\boldsymbol {u}(\boldsymbol {x}_r,\boldsymbol {x}(t)) + \boldsymbol {u}_b(\boldsymbol {x}(t)) \bigg ] \,\rm d\it t+ \left(\frac {2k_BT}{\boldsymbol {f}}\right)^{1/2} \cdot \mathcal {N}\big(0,{\delta t}\big), \end{align}
where
$\mathcal {N}(0, {\delta t})$
is a three-dimensional vector whose components are normally distributed with mean
$0$
and variance
$\delta t$
. The numerical integration was performed using the MATLAB solver ode15s (with the relative tolerance of
$10^{-12}$
); on each subinterval of length
$\delta t$
we superimpose the particle movement due to the white noise expressed in (2.3) with its flow-induced movement which is solved for by the numerical integration of (2.2) without the noise term. We use the same procedure to simulate a system of
$N_s$
particles in the source of concentration
$C_0=N_s/V$
, where
$V=\pi r_s^2h_s$
is the disk’s volume. If any particles depart the source during a given time step, they are replaced at random positions within the disk to maintain the constant concentration boundary condition. Particles that return to the disk due to either Brownian motion or ciliary advection are also accounted for in the source; however, the ratio of particles entering to exiting the source is very small (i.e.
$\ll 1$
).
3. Absence of background flow (
$\dot {\gamma }=0$
)
We start by investigating the effect of ciliary flows on particle transport in the absence of a linear shear flow (
$\dot {\gamma }=0$
), specifically examining the dependence on both the particles’ diffusivity and the torque magnitude. We study the transport of large dissolved molecules in the vicinity of the cilium through to micron-sized particulate matter, which corresponds to
$D_0\approx 0.1 \,-\,2400\,\unicode{x03BC}\text {m}^2\,\text {s}^{-1}$
(Shapiro et al. Reference Shapiro, Fernandez, Garren, Guasto, Debaillon-Vesque, Kramarsky-Winter, Vardi and Stocker2014; Ramirez-San Juan et al. Reference Ramirez-San Juan, Mathijssen, He, Jan, Marshall and Prakash2020). Our reference case is that of
$|\boldsymbol{\Omega }|=\Omega _s=0.01 \ \mbox {fNm}$
for
$d=8\ \unicode{x03BC}\text{m}$
corresponding to the eukaryotic cilium (length
$l_c=12.75\,\unicode{x03BC}\text{m}$
; frequency
$f_c=16.9\,\rm s^{-1}$
) of P. damicornis obtained using Selvan et al. (Reference Selvan, Duck, Pihler-Puzović and Brumley2023), where the total steady-state flow field of a cilium was fitted with the point torque positioned above the rigid wall. Furthermore, the following parameters remain unchanged throughout our analysis:
$\mu =10^{-3} \ \mbox {Pa s}$
,
$T=298 \ \mbox {K}$
,
$r_s/d=1/\sqrt {\pi }$
,
$h_s/d=0.25$
and
$k_B=1.38\times 10^{-23} \ \mbox {J K}^{-1}$
. For the assumed parameters, the contribution to the drag
$\boldsymbol {f}$
due to the wall has a negligible effect on the reported results.
3.1. Validation with a continuum model
In order to provide macroscopic validation of the microscopic Langevin model described in § 2, we introduce a continuum model which solves an advection–diffusion equation in the presence of ciliary flows. The steady concentration of particles,
$C^*$
, is obtained from
\begin{align} D_0 \nabla ^2 C^*-(\boldsymbol {u} \cdot \nabla )C^*=0, \end{align}
in a box of dimension (
$l_b\times w_b \times h_b$
), with the circular concentration source of radius
$r_s$
located at the centre of the bottom face (figure 1
b), and an absorbing boundary condition (sink) located at
$z=h_b$
.
In contrast to the Langevin model, which requires only one boundary condition on the plane rigid wall and open boundaries in all other directions, here we require the following boundary conditions to be imposed along the box surfaces in order to solve (3.1):
\begin{align} &\frac {\partial C^*}{\partial x}=0\,\,\mbox {at}\,\,x=-\frac {l_b}{2},\frac {l_b}{2}, \qquad \frac {\partial C^*}{\partial y}=0\,\,\mbox {at}\,\,y=-\frac {w_b}{2},\frac {w_b}{2}, \end{align}
\begin{align} &\frac {\partial C^*}{\partial z}=0\,\,\mbox {at}\,\,z=0\,\,\mbox {(except at source)}, \end{align}
\begin{align} &C^*=kC_0 \,\,\mbox {at source,}\,\,\,\mbox {and}\,\,\,C^*=0\,\,\mbox {at}\,\,z=h_b\,\mbox {(sink)}. \end{align}
Equations (3.2a
) specify zero diffusive flux, but permit advective flux across the side boundaries. In general, the position of the side boundaries is chosen to be sufficiently large so that the diffusive flux is negligible. We use
$l_b=10d$
,
$w_b=10d$
and
$h_b=12d$
throughout our study for particles with diffusivity
$D_0=0.65\, \unicode{x03BC}\text{m}^2\,\rm s^{-1}$
, and have determined that increasing the box dimensions further has no bearing on the reported result. The factor
$k$
in (3.2c
) takes into account the different form of the source boundary condition, which is flat (i.e. imposed at
$z=0$
) for the continuum model, but has a finite volume in the Langevin model. The fluid in the immediate vicinity of a surface concentration condition would exhibit a slightly lower mean concentration, and thus representing the boundary condition through a finite volume requires a proportionally lower seeding density (with factor
$k$
).
Notably, the advective flow generated by the point torque in (2.1) is singular at the point of application. To overcome numerical issues while solving equation (3.1), we desingularise the velocity field (2.1) using the method of regularisation employed in the work of Cortez & Varela (Reference Cortez and Varela2015). A detailed description of this process is given elsewhere (Cortez Reference Cortez2001; Cortez & Varela Reference Cortez and Varela2015; Park, Kim & Lim Reference Park, Kim and Lim2019). Using the same notation as before, the modified velocity field
$\boldsymbol {u}^{\kappa }$
due to the regularised torque is
\begin{align} u_i^{\kappa } & =\frac {1}{2} Q(r) \varepsilon _{ijk} \Omega _j r_k-\frac {1}{2} Q(R) \varepsilon _{ijk} \Omega _j R_k+d\left\{D_1^{\phi }(R)\varepsilon _{ij3}\Omega _j+D_2^{\phi }(R) \varepsilon _{\alpha j3}\Omega _j R_{\alpha } R_i\right\}\nonumber \\ & \quad -\{\varepsilon _{\alpha j3}\Omega _j R_{\alpha } \delta _{i3}+\varepsilon _{ij3} R_3 \Omega _j\} H_3(R)-\varepsilon _{\alpha j3}\Omega _j R_{\alpha } R_3 R_i D_2^{\phi }(R)-{\rm d}\varepsilon _{ijk} \Omega _j R_k H_3(R)\nonumber \\ & \quad -{\rm d} \varepsilon _{ijk} \Omega _j R_3 R_k H_4(R)-\varepsilon _{ijk} (\Omega _j-\Omega _3 \delta _{j3}) R_k H_3(R)+d^2 \varepsilon _{ijk} \Omega _j R_k, \end{align}
where
$\alpha \in \{1,2\}$
; the functions
$Q$
,
$D_1^{\phi }$
,
$D_2^{\phi }$
,
$H_3$
and
$H_4$
are defined as follows (Park et al. Reference Park, Kim and Lim2019):
\begin{eqnarray} &&Q(r)=\frac {5\kappa ^2+2r^2}{8\pi (r^2+\kappa ^2)^{5/2}},\,\, D_1^{\phi }(R)=\frac {R^2-2\kappa ^2}{4\pi (R^2+\kappa ^2)^{5/2}},\,\,D_2^{\phi }(R)=\frac {-3}{4\pi (R^2+\kappa ^2)^{5/2}},\,\,\nonumber \\ &&H_3(R)=\frac {-3\kappa ^2}{8\pi (R^2+\kappa ^2)^{5/2}},\,\, H_4(R)=\frac {15\kappa ^2}{8\pi (R^2+\kappa ^2)^{7/2}}, \end{eqnarray}
and
$\kappa$
is the regularisation parameter (i.e. the radius of the ‘blob’) which is assumed to be very small compared with the reference length scale (here
$\kappa /d=0.001$
).
The regularised velocity expression in (3.3) is used within the advection–diffusion equation (3.1), which is solved with a customised finite element method via MATLAB’s solvepde solver (element, tetrahedron; maximum element size,
$3.5d$
; minimum size,
$0.17d$
; geometric order, quadratic) for the boundary conditions given in (3.2a
)–(3.2c
). Since the velocity decays like
$1/r^2$
from the regularised torque, the advective flux along the side boundaries of the domain is extremely small. However, in the absence of explicit boundary conditions for the advective term, the solvepde solver does not restrict the advective flux at these boundaries, instead handling them as open boundaries for advection (Padilla, Secretan & Leclerc Reference Padilla, Secretan and Leclerc1997).
An example of a (scaled) concentration profile computed using this model and averaged over the thickness
$|y| \leq r_s$
as
\begin{eqnarray} C^*(x,z)=\frac {1}{2r_s}\int _{-r_s}^{r_s} C^*(x,y,z) \,\rm d\it y, \end{eqnarray}
is shown in figure 2(b); here,
$D_0=0.65 \, \unicode{x03BC} \text {m}^2\,\rm s^-{^1}$
and
$C_0=0.2\,\unicode{x03BC} \text {m}^{-3}$
.
Scaled average concentration of particles with diffusivity
$D_0=0.65 \, \unicode{x03BC} \text {m}^2\,\text {s}^{-1}$
in the presence of ciliary flows of magnitude
$|\boldsymbol {\Omega }|=0.01\,\, \mbox {fNm}$
calculated using (a) the Langevin model simulated for
$n\delta t=350\, \text {s}$
to attain the steady-state distribution (
$C/C_0$
) and (b) the continuum model (
$C^*/C_0$
) with factor
$k=1.3$
(see appendix A) for the source concentration of
$C_0=0.2\,\unicode{x03BC}\text{m}$
$^{-3}$
. Lines with arrows indicate the streamlines of the flow field.
In figure 2(a), we show the corresponding average concentration of particles evaluated using the Langevin model from § 2 (in the presence of the ciliary flow only). This concentration was calculated by dividing the domain of thickness
$|y| \leq r_s$
into smaller boxes of equivalent thickness, i.e. of volume
${\rm d}V = 2r_s \times {\rm d}x\times {\rm d}z$
where
${\rm d}x={\rm d}z=0.1d\ll d$
so that the variations in average concentration profile close to the point torque were resolved. The reported results did not change significantly when taking other values of
${\rm d}x={\rm d}z$
provided they were much smaller than
$d$
. Then in each box, we found the concentration in the mth time interval
$[m\delta t,\,(m+1) \delta t]$
as
$\displaystyle C^{(m)}(x,z)=N_s^{(m)}/{\rm d}V$
, where
$m\in \mathbb {Z}^+$
and
$N_s^{(m)}$
corresponds to the number of particles in the box during that time interval. Hence, the average concentration profile in each box over time was found as
\begin{eqnarray} C(x,z)=\frac {1}{n}\sum _{m=1}^n C^{(m)}(x,z), \end{eqnarray}
where
$n\delta t$
is the total simulation time using the Langevin model. The results in figure 2(a) can be compared directly with results in figure 2(b). The predictions for the average concentrations found using the two models agree to within 10% when the factor
$k=1.3$
is used (see appendix A for detailed description). The continuum model presented above serves to validate the results of the Langevin model for particles with a diffusivity
$D_0=0.65\, \unicode{x03BC} \text {m}^2\,\text {s}^{-1}$
. For other diffusivity values, the lateral size of the box should be optimised to ensure negligible diffusive flux far from the source. We have displayed one such result for higher diffusivity (with
$D_0=1400 \, \unicode{x03BC} \text {m}^2\,\text {s}^{-1}$
) with an optimised box size in appendix A. We emphasise here that although the continuum model would be broadly useful for various ciliary arrangements and flow regimes, the remainder of this paper utilises exclusively the Langevin model to study ciliary mass transport.
3.2. Particle trajectories and distribution
Next, we investigate how the particles’ trajectories, and their steady-state distributions, are affected by the ciliary flows and particle diffusivity. The displacement of particles with diffusivity
$D_0=0.65 \ \unicode{x03BC} \text {m}^2\,\rm s^-{^1}$
exhibits no bias in the
$x$
–
$y$
plane when the ciliary flow is absent (i.e.
$|\boldsymbol{\Omega }| = 0$
, no advection, see figure 3
a). However, when the cilium generates a flow, particles are advected along the direction of the flow (figure 3
b). Examining the corresponding trajectories for
$10$
example particles (figure 3 aii, bii, cii) reveals that the strong near-field vortical flow circulates particles before they escape by diffusion. This is comparable to the particle transport in the presence of time-dependent ciliary motion obtained with the colloidal rotor model (Vilfan & Jülicher Reference Vilfan and Jülicher2006; Selvan et al. Reference Selvan, Duck, Pihler-Puzović and Brumley2023), where the oscillating particles can be trapped in closed loops (appendix D). Further, we have also examined the current Langevin formulation (2.3) in the presence of ciliary flow at the limit of zero diffusivity (appendix C), where the noise along the particle trajectories is negligible. However, increasing the diffusivity, e.g. by a factor of 10 in figure 3(c), diminishes the downstream bias of the particle distribution due to the relatively larger role of diffusion. The motion of Brownian particles is evidently controlled by the interplay between the particle diffusivity and strength of the point torque, encapsulated in a Péclet number. This is discussed in more detail in the following subsection.
Particles from the steady-state distribution and corresponding two-dimensional projection of the concentration map (a i,b i,c i), along with example trajectories of
$10$
particles (a ii,b ii,c ii) emitted from the source. Results are shown for (a)
$|\boldsymbol{\Omega }|=0$
with
$D_0=0.65 \, \unicode{x03BC} \text {m}^2\,\rm s^-{^1}$
(no-advection), (b)
$|\boldsymbol{\Omega }|=\Omega _s$
with
$D_0=0.65 \, \unicode{x03BC} \text {m}^2\,\rm s^-{^1}$
(advection) and (c)
$|\boldsymbol{\Omega }|=\Omega _s$
with
$D_0=6.5 \, \unicode{x03BC} \text {m}^2\,\rm s^-{^1}$
(advection, larger diffusivity), for source concentration of
$C_0=0.2\, \unicode{x03BC} \text {m}^{-3}$
.
3.3. Nutrient transport
Here, we further investigate the transport of particles of various diffusivities by evaluating their integrated flux, or current, between the no-slip surface and the bulk fluid. This quantity is the net emission rate of particles from the source under the influence of both advection and diffusion and is used to measure the effect of ciliary flows on mass transport (Saito Reference Saito1968; Berg Reference Berg1993). The current is defined as
\begin{align} I=\frac {N_e}{T_e}, \end{align}
where
$N_e$
is the number of particles emitted from the source in a given time
$T_e$
. In the absence of any flows, the diffusion current emanating from a disk-shaped source with surface concentration,
$C=C_0$
, and far-field concentration,
$C=0$
, is directly proportional to the diffusivity of the particles and the radius of the disk (Berg Reference Berg1993):
\begin{align} I=4r_sC_0 D_0. \end{align}
We plot the current
$I$
obtained in Langevin simulations as a function of
$D_0$
in figure 4(a). In the absence of advection (i.e. for
$|\boldsymbol {\Omega }|=0$
), it is in close agreement with (3.8) and prediction from the continuum model, i.e.
\begin{eqnarray} I=D_0\int _{0}^{2\pi }\int _{0}^{r_s} \frac {\partial C^*}{\partial r}\bigg |_{z=0} r \rm d\it r \rm d\theta . \end{eqnarray}
Meanwhile, the interplay between the ciliary flows and Brownian motion results in a deviation from this linear profile in the presence of advection (i.e.
$|\boldsymbol {\Omega }|\neq 0$
).
(a,b) Dimensional current
$I$
(defined for
$T_e=5\,\mbox {min}\{\tau _a,\tau _d\}$
when
$C_0=0.1 \, \unicode{x03BC} \text {m}^{-3}$
) as a function of particle diffusivity
$D_0$
showing the effects of (a) advection due to the point torque and (b) changing the torque magnitude. (c) Scaled current (
$I/C_0U_c{L}^2$
) as a function of scaled diffusivity (
$1/Pe_r$
) for different torque magnitudes.
For small
$D_0$
, the current is enhanced by the presence of a ciliary flow (
$|\boldsymbol {\Omega }| \gt 0$
), which serves to transport particles away from the source. However, at larger values of
$D_0$
, the relative enhancement due to advection diminishes. A comparison with the time-dependent ciliary motion (appendix D) revealed that the steady rotlet model tends to slightly underestimate particle current due to the reliance on a time-averaged flow approximation. Figure 4(b) shows
$I$
as a function of
$D_0$
for various magnitudes of the point torque. These results reveal that the current strongly depends on
$|\boldsymbol {\Omega }|$
. The extent of the deviation from the linear profile of (3.8) is determined by the relative sizes of diffusive,
$\tau _d$
, and advective,
$\tau _a$
, time scales, discussed in more detail in appendix B. For
$\tau _d \leq \tau _a$
, diffusion is the dominant transport mechanism, and the current scales as
$I \propto D_0$
.
The family of curves in figure 4(b) can be readily rescaled using the rotlet Péclet number,
\begin{eqnarray} Pe_r=\frac {\tau _a}{\tau _d}=\frac {U_c{L}}{D_0}, \end{eqnarray}
where
$U_c$
and
$L$
are the characteristic advection velocity and length scales, respectively (see appendix B). Plotting these data after rescaling in figure 4(c) shows that the current remains unchanged as a function of the scaled diffusivity for different point torque magnitudes. Furthermore, the transition from the nonlinear to the linear profile occurs around
$Pe_r=1$
, where the ciliary advection balances diffusion. To investigate the relative influence of advection compared with diffusion in the observed mass transport, we utilised the Sherwood number (Friedlander Reference Friedlander1957; Acrivos & Taylor Reference Acrivos and Taylor1962; Magar et al. Reference Magar, Goto and Pedley2003),
\begin{eqnarray} Sh=\frac {I}{C_0D_0L}, \end{eqnarray}
which measures the ratio of the particle current,
$I$
, to the diffusion rate,
$C_0 D_0 L$
. This is plotted as a function of
$Pe_r$
in figure 5(a), with the corresponding log–log plot shown in figure 5(b). In the limit of high rotlet Péclet number, we observe a scaling of
$Sh\sim Pe_r^{3/4}$
, which suggests that the relative effect of diffusion diminishes with increasing
$Pe_r$
. The rate of this decrease is higher compared with other flows reported in the literature, e.g. around a no-slip cylinder for which the exponent of
$Sh$
versus
$Pe_r$
is
$1/3$
(Friedlander Reference Friedlander1957). Taken together, our results imply that nutrient transport by an individual cilium can be characterised using the rotlet Péclet number.
(a) Sherwood number (
$Sh$
) as a function of rotlet Péclet number (
$Pe_r$
) for different torque magnitudes. (b) Data from (a) shown on a log–log scale, which demonstrates that
$Sh\sim Pe_r^{3/4}$
for
$Pe_r\gg 1$
.
4. Presence of background flow (
$\dot {\gamma }\neq 0$
)
4.1. Particle trajectories
In this section, we investigate mass transport in the presence of both a ciliary flow and an externally imposed linear shear flow. We start by focusing on the effect of ciliary orientation on the trajectories of particles with low diffusivity (i.e.
$D_0=0.65 \ \unicode{x03BC} \text {m}^2\,\text {s}^{-1}$
). Figure 6 shows particle trajectories in the presence of ciliary flows with two different orientations,
$\boldsymbol {\Omega }= \pm |\boldsymbol {\Omega }|\hat {\boldsymbol {e}}_y$
. There is a significant qualitative difference in the trajectories of particles, depending on the orientation of the model cilium. When the external flow coincides with the far-field flow direction of ciliary beating (
$\boldsymbol {\Omega }=|\boldsymbol {\Omega }|\hat {\boldsymbol {e}}_y$
, figure 6
a), the reverse flow generated by the cilium near the source traps the emitted particles in closed loops, before they ‘escape’ due to diffusion. Conversely, for the cilium with opposite orientation to the external flow direction (
$\boldsymbol {\Omega }=-|\boldsymbol {\Omega }|\hat {\boldsymbol {e}}_y$
, figure 6
b), source particles are immediately swept downstream by the external flow due to the constructive nature of both ciliary and external flows close to the source. In appendix D we show that when time-dependent ciliary motion is accounted for, particles oscillate in the near-field, but are still trapped within closed loops characterised by the same length scales as in figure 6, which results in the same average distribution of particles.
Representative trajectories of 10 particles with diffusivity
$D_0= 0.65 \ \unicode{x03BC} \text {m}^2\,\rm s^-{^1}$
(a i,b i) and the corresponding flow field due to both ciliary and external flows (a ii,b ii). Results are shown for (a)
$\boldsymbol {\Omega }=\Omega _s\hat {\boldsymbol {e}}_y$
and (b)
$\boldsymbol {\Omega }=-\Omega _s\hat {\boldsymbol {e}}_y$
and the external shear flow with
$\dot {\gamma }=0.5\,\text {s}^{-1}$
. The black and red arrows in (a i,b i) plots indicate the direction of linear shear and ciliary flows, respectively. The black lines with arrows in (a ii,b ii) are the streamlines of the flow field.
4.2. Nutrient transport
We now explore the nutrient transport due to ciliary flows in the presence of a background flow by computing the current as defined in (3.7). In figure 7(a–c), we plot
$I$
as a function of diffusivity
$D_0$
for different strengths of advection due to ciliary (
$\boldsymbol {\Omega }=\Omega _s \hat {\textbf {e}}_y$
) and external flows. As shown in figure 7(a), when the external flow is non-zero (
$\dot {\gamma }\neq 0$
), the current still varies linearly with
$D_0$
for sufficiently large diffusivity (diffusive region), but the relationship between the two is otherwise nonlinear (advective region). Compared with the analogous case of no external flow (i.e.
$\dot {\gamma }= 0$
), the emission of nutrients increases considerably in the latter region. In addition to this observation, we also observed that the current measured in time-dependent ciliary motion for the external flow with low shear rate tends to slightly overestimate the current calculated using the steady rotlet model (appendix D).
Particle current
$I$
as a function of
$D_0$
showing the effect of (a) advection due to both ciliary (
$\boldsymbol {\Omega }=\Omega _s \hat {\textbf {\textit{e}}}_y$
) and background flow, and (b,c) shear rate for different torque magnitudes and orientations. (d) Scaled current (
$I/C_0U_c{L}^2$
) as a function of scaled diffusivity (
$1/Pe$
) for different torque magnitudes and orientations and for different shear rates. (e) The same data as in (d) but focusing on the limit of
$1/Pe\to 0$
to show the effect of ciliary orientation in the presence of background flow. The remaining parameters are the same as in § 3.1.
In figures 7(b) and 7(c), we study the effect of shear rate on the current
$I$
as a function of
$D_0$
for different torque magnitudes and orientations. When the shear rate is increased, the current also increases for both torque magnitudes considered,
$|\boldsymbol{\Omega }|=\Omega _s$
and
$|\boldsymbol{\Omega }|=5\Omega _s$
. By contrast, the change in the ciliary orientation compared with the external flow appears not to affect the nutrient transport for the range of diffusivities considered here, though a closer examination of the region of small diffusivity suggests a more interesting picture (see below). To probe this further, we scale the family of curves in figures 7(b) and 7(c) as in § 3.3, but using the total (or ‘effective’) Péclet number
\begin{eqnarray} Pe=\frac {U_c{L}}{D_0}=Pe_r+Pe_{\dot {\gamma }}, \end{eqnarray}
where
$Pe_r=U_rL/D_0$
and
$Pe_{\dot {\gamma }}=U_eL/D_0$
are the rotlet and shear Péclet numbers, respectively;
$U_r$
and
$U_e$
are the average velocity experienced by the particles in the presence of ciliary and external flows, respectively (see B2). When rescaled, all data align well irrespective of the shear rate, ciliary orientation and torque magnitudes, as shown in figure 7(d). This effective Péclet number is not the only control parameter that governs the current across all flow conditions, because the structure of, and therefore transport by, the external shear and rotlet flows are fundamentally different. However, we have examined the particle current as the relative contributions of
$Pe_r$
and
$Pe_{\dot {\gamma }}$
are varied for fixed total
$Pe$
, and found that it does not depend strongly on the composition of
$Pe$
. In varying
$1/Pe$
from
$0.1$
to
$1.5$
, the scaled particle current varied by a factor of
${\sim}10$
. However, for each value of
$1/Pe$
, we found that varying the relative proportion of
$Pe_r$
and
$Pe_{\dot {\gamma }}$
resulted in just
${\sim}1\,\%$
difference in the scaled particle current. This suggests that
$Pe$
adequately characterises the total effect of flow relative to diffusion.
We now return to the effect of orientation on the less diffusive particles by examining the scaled current when
$1/Pe \to 0$
(see figure 7
e). When the far-field flow from the point torque is aligned with the external flow (solid line), the particle emission rate is diminished due to competition between the ciliary and external flows near the no-slip boundary. Conversely, for the point torque oriented so that its far-field flow opposes the external shear (red dash–dot line), the particle emission rate increases linearly with the shear rate due to the superposition of ciliary and external flows around the source. Increasing
$1/Pe$
weakens the importance of the orientation.
5. Discussion
We have explored the mass transport of ciliary flows which are modelled using a rotlet. We employ the Langevin model of mass transport to calculate the movement of particulate matter, which is emitted from a disk-shaped source positioned below the point torque. We validate the Langevin model with a continuum description, by solving the advection–diffusion equation in the absence of an external flow. In this case, we observe that ciliary beating results in a pronounced enhancement in the nutrient current for sufficiently small particle diffusivity (i.e. high Péclet number).
In the presence of an externally imposed linear shear flow, the ciliary orientation plays a crucial role in governing particle transport, whenever the ciliary flow is relatively strong compared with the external flow and Brownian noise. When the external flow coincides with the direction of ciliary pumping, the particles tend to become trapped in a closed loop before being released due to diffusion. In such cases, the particle emission rate decreases due to the trade-off between the strengths of ciliary and external flows. For ciliary and external flows that oppose one another, the superposed flow leads to downstream particle release, which, in turn, increases the particle emission rate.
Importantly, these findings reveal that average transport of dissolved or particle matter is driven by the flows within the typical ciliary length scale – and therefore associated with the recovery stroke of a cilium – rather than the far-field characteristics of the ciliary beating. This result hinges on being able to resolve the near-field ciliary flows. Our model is therefore expected to find utility in applications involving ciliary carpets, in which heterogeneity of ciliary alignment can further enhance the three-dimensional mixing effects. Patches of cilia with the same orientation could be modelled by increasing torque strengths, which requires further research into the relationship between torque magnitudes, the ciliary beating strength and the number of cilia within patches. Such work is on its way.
Lastly, we note that although this work has focused on the transport of nutrients away from a no-slip substrate – for example modelling transport of excess oxygen away from a coral surface – this framework could also be adapted to consider nutrient uptake by the surface by reversing role of the source and sink (Masoud & Stone Reference Masoud and Stone2019) together with the flow velocities, and noting the symmetry of the configuration.
Acknowledgements.
The work was enabled by the joint PhD programme between the University of Melbourne and the University of Manchester. S.A.S. is grateful to the University of Melbourne for providing funding through its Xing Lei Scholarship and the Melbourne Research Scholarship. S.A.S. performed simulations using The University of Melbourne’s High-Performance Computer Spartan (https://doi.org/10.4225/49/58ead90dceaaa). D.R.B. was supported by a Symbiosis in Aquatic Systems grant from the Gordon and Betty Moore Foundation (Grant 9351).
Declaration of interests.
The authors report no conflict of interest.
Appendix A. Factor ‘
$k$
’
The source geometries used in the Langevin and continuum models are slightly different. To compensate for this, we introduce a factor in the source boundary condition of the continuum model, see (3.2c
). For the reference case shown in figure 2, the largest concentration of particles was in the region of interest
$x \in [-3d/2,6d]$
and
$z \in [0,2d]$
. In this region, we compare the average absolute difference in the average concentration profiles obtained using the two models
\begin{eqnarray} \langle \Delta C \rangle =\frac {1}{p} \sum _p |C^*-C|, \end{eqnarray}
where
$p$
is the number of grid points in the region of interest. When
$k=1.3$
,
$\langle \Delta C \rangle$
attains the minimum value as shown in figure 8(a), resulting in less than
$10\,\%$
difference in the average particle concentration obtained using the two models that were quoted in § 3.1. The favourable comparison is also demonstrated in figure 8(b), which plots the corresponding heat map of the absolute difference,
$\Delta C=|C^*-C|$
, when
$k=1.3$
.
(a) Average absolute difference
$\langle \Delta C \rangle$
as a function of
$k$
for the particles discussed in § 3.1. The red cross indicates the minimum
$\langle \Delta C \rangle$
attained. (b) The corresponding absolute difference in the average particle concentration obtained using the two models when
$k=1.3$
. (c) The near-field average concentration of particles (scaled by
$C_0$
) computed using the Langevin model for parameters from § 3.1. Blue lines with arrows denote the flow streamlines, and the annular region
$\mathcal {R}$
discussed in appendix B is marked with black solid lines.
Scaled average concentration of particles with diffusivity
$D_0=1400 \, \unicode{x03BC} \text {m}^2\,\text {s}^{-1}$
in the presence of ciliary flows of magnitude
$|\boldsymbol {\Omega }|=5\Omega _s$
calculated for the source concentration of
$C_0=0.2\,\unicode{x03BC}\text{m}$
$^{-3}$
using (a) the Langevin model that was simulated for
$n\delta t=70$
s to attain the steady-state particle distribution (
$C/C_0$
) and (b) the continuum model (
$C^*/C_0$
) using
$k=1.3$
. Lines with arrows indicate the streamlines of the flow field. (c) The corresponding absolute difference in the average particle concentration
$\Delta C$
computed using the models.
Finally, we also compare the two models for larger values of particle diffusivity,
$D_0=1400\,\unicode{x03BC} \text {m}^2\,\text {s}^{-1}$
, under the influence of flow generated by the torque of magnitude
$|\boldsymbol {\Omega }|=5\Omega _s$
, still using
$k=1.3$
. The box dimensions along lateral directions are fixed as
$l_b=30d$
and
$w_b=30d$
to ensure that the diffusive flux becomes negligible, while the height is maintained as
$h_b=12d$
. The corresponding concentration maps are plotted in figures 9(a) and 9(b), respectively, and result in less than
$10\,\%$
disagreement as shown in figure 9(c) where we plot the absolute difference in the average particle concentrations. This illustrates that the level of agreement between the models is similar for all values of particle diffusivity and torque used in § 3. While
$k=1.3$
has proven effective for the highly diffusive particles (
$D_0=1400\, \unicode{x03BC} \text {m}^2\,\text {s}^{-1}$
) under large torque magnitude (
$|\boldsymbol {\Omega }|=5\Omega _s$
), the optimal value of
$k$
differs slightly as shown in figure 8(a). The optimal value of
$k$
depends on the dimension of the source in the Langevin model, particle diffusivity and torque magnitudes. However, we emphasise that there is little difference between the optimal values of
$k$
shown in figure 8(a), and note that the continuum model is only used here to compare with the Langevin model. The latter approach is used for all subsequent calculations and conclusions.
Appendix B. Characteristic scales
Here, we define the characteristic scales of the mass transport due to a single cilium. Consider the near-field average concentration of particles with diffusivity
$D_0=0.65\, \unicode{x03BC} \text {m}^2\,\text {s}^{-1}$
in the presence of ciliary flows shown in figure 8(c). These data were obtained from the Langevin model as discussed in § 3.1. It suggests that on average, the particles are advected along the vortex within the annular region,
$\mathcal {R}=\{(x,0,z)|r=\sqrt {x^2+(z-d)^2},2d/5\lt r\lt d\}$
. Hence, we choose the characteristic length and velocity scales of advection as the vortex length and the average velocity experienced by the particles due to the cilium in the annular region
$\mathcal {R}$
,
\begin{eqnarray} {L}=2d \quad \text {and} \quad U_c=U_r\sim \frac {0.4352|\boldsymbol {\Omega }|}{\unicode{x03BC} {L}^2}, \end{eqnarray}
respectively. In the presence of external linear shear flow, the characteristic velocity scale is modified by the combined effect of ciliary motion and external shear,
\begin{eqnarray} U_c=U_r+U_e\sim \frac {0.4352|\boldsymbol {\Omega }|}{\unicode{x03BC} {L}^2}+\dot {\gamma }{L}. \end{eqnarray}
From this, the advection time scale is given as
$\tau _a={L}/U_c$
. On the other hand, the time taken to traverse the same vortex length
${L}=2d$
in the absence of advection is given by the diffusive time scale,
$\tau _d={L}^2/D_0$
. In figure 8(c), the flow due to the point torque (
$d=8\, \unicode{x03BC} \text {m}$
and
$|\boldsymbol {\Omega }|=\Omega _s$
) transports the particles with the characteristic velocity
$U_c \sim 16\, \unicode{x03BC} \text {m}\,\text {s}^{-1}$
in the vortex of length
$L=16\, \unicode{x03BC} \text {m}$
. The results of the Langevin model were obtained using
$\displaystyle \delta t = \mbox {min}\{\tau _a,\tau _d\}/100$
.
Appendix C. Transition to non-diffusive tracers (
$D_0 \to 0$
)
In this section, we examine the effect of modifying the Brownian particles’ diffusivity from a finite value (
$D_0\gt 0$
), towards zero. Figure 10 displays example trajectories of particles, calculated using the Langevin formulation in the presence of ciliary flow, for several values of decreasing particle diffusivity. Unsurprisingly, for extremely small diffusivity – e.g.
$D=0.0025\, \unicode{x03BC} \text {m}^2\,\text {s}^{-1}$
as shown in figure 10(c) – the trajectories executed by the particles deviate only very slightly from the streamlines associated with the ciliary flow (Selvan et al. Reference Selvan, Duck, Pihler-Puzović and Brumley2023). This recirculation of particles is carefully considered when calculating the net emission rate from the source at the surface.
Example trajectories of
$10$
particles emitted from the source are shown for (a)
$D_0=0.25 \, \unicode{x03BC} \text {m}^2\,\text {s}^{-1}$
, (b)
$D_0=0.025\, \unicode{x03BC} \text {m}^2\,\text {s}^{-1}$
and (c)
$D_0=0.0025 \, \unicode{x03BC} \text {m}^2\,\text {s}^{-1}$
with
$|\boldsymbol {\Omega }|=\Omega _s$
for source concentration of
$C_0=0.2 \, \unicode{x03BC} \text {m}^{-3}$
.
Appendix D. Predicting particle transport using the colloidal rotor model
Here we use the colloidal rotor (or ‘minimal’) model, which captures the transient flow field of the single cilium, to predict the particle mass transport. Following Vilfan & Jülicher (Reference Vilfan and Jülicher2006); Selvan et al. (Reference Selvan, Duck, Pihler-Puzović and Brumley2023), we model the ciliary tip as a sphere of radius
$a$
translating with constant angular velocity
$\omega$
in a circular trajectory of radius
$r_0$
, whose centre is positioned at the distance
$h$
above the rigid wall (see figure 11
a). Hence, the evolution of Brownian particle position
$\boldsymbol {x}=(x,y,z)$
due to time-dependent ciliary
$\boldsymbol {u}^s$
and background linear flows
$\boldsymbol {u}_b$
is given by
\begin{eqnarray} \frac {\,\rm d\boldsymbol {x}}{\,\rm dt}=\boldsymbol {u}^s+\boldsymbol {u}_b+\sqrt {\frac {2k_BT}{\boldsymbol {f}}} \cdot \boldsymbol {\eta }(t), \end{eqnarray}
where the fluid velocity generated by the transient sphere located at
$\boldsymbol {X}(t)=(x_s(t),y_s(t),z_s(t))=(r_0\cos \omega t,0,{h-r_0}\sin \omega t)$
is
\begin{eqnarray} u_i^s=\frac {1}{8\pi \mu }S_{ij}(\boldsymbol {X}(t),\boldsymbol {x}) F_j(t), \end{eqnarray}
and
\begin{align} S_{ij}(\boldsymbol {X}(t),\boldsymbol {x}) & =\frac {\delta _{ij}}{r}+\frac {r_ir_j}{r^3}-\left(\frac {\delta _{ij}}{R}+\frac {R_iR_j}{R^3}\right)\nonumber\\& \quad +2z_s(t) \rho _{jk3}\left[\frac {\partial }{\partial R_k}\left(\frac {z_s(t) R_i}{R^3}-\left[\frac {\delta _{i3}}{R}+\frac {R_iR_3}{R^3}\right]\right)\right], \end{align}
with
$\boldsymbol {r}=(r_1,r_2,r_3)=(x-x_s(t),y-y_s(t),z-z_s(t))$
,
$\boldsymbol {R}=(R_1,R_2,R_3)=(x-x_s(t),$
$y-y_s(t),z+z_s(t))$
,
$\rho _{jk3}=\delta _{j\alpha }\delta _{\alpha k}-\delta _{j3}\delta _{3 k}$
and the force exerted on fluid by the sphere is
$F_j(t)=6\pi \mu a (\delta _{jk}+\frac {9a}{16 z_s(t)} (\delta _{jk}+\delta _{j3} ) )X^{\prime}_k(t)$
. All other notation is the same as in § 2.
(a) Diagram of the Langevin model of mass transport when the ciliary flow is modelled using a colloidal rotor. (a–c) Typical trajectories of 10 particles with diffusivity
$D_0= 0.65 \ \unicode{x03BC} \text {m}^2\,\text {s}^{-1}$
in the presence of ciliary and external flows. Results are shown for (b) clockwise ciliary orientation and
${\dot {\gamma }}=0$
, (c) clockwise ciliary orientation and
${\dot {\gamma }}=0.5$
s
$^{-1}$
, and (d) anticlockwise ciliary orientation and
${\dot {\gamma }}=0.5$
s
$^{-1}$
. As before, the black and red arrows indicate the direction of linear shear and ciliary flows, respectively.
Particle current
$I$
(defined for
$T_e=5\,\mbox {min}\{20/f_c,\tau _d\}$
when
$C_0=0.1 \, \unicode{x03BC} \text {m}^{-3}$
) as a function of
$D_0$
for clockwise ciliary orientation comparing the time-averaged rotlet and time-dependent colloidal rotor models for different values of shear rate (a)
$\dot {\gamma }=0$
, (b)
$\dot {\gamma }=0.5 \ \text {s}^{-1}$
, (c)
$\dot {\gamma }=1\ \text {s}^{-1}$
and (d)
$\dot {\gamma }=3 \ \text {s}^{-1}$
.
In our simulation, we used the following parameters for the colloidal rotor model:
$a=3 \,\unicode{x03BC}$
m,
$r_0=3\,\unicode{x03BC}$
m,
$h=6\,\unicode{x03BC}\text{m}$
and
$\omega =2\pi f_c$
; they have been optimised for the eukaryotic cilium of P. damicornis (Shapiro et al. Reference Shapiro, Fernandez, Garren, Guasto, Debaillon-Vesque, Kramarsky-Winter, Vardi and Stocker2014; Selvan et al. Reference Selvan, Duck, Pihler-Puzović and Brumley2023). We investigate the effect of time-dependent ciliary and external linear shear flows on the trajectories of particles with diffusivity
$D_0=0.65\, \unicode{x03BC} \text {m}^2\,\text {s}^{-1}$
by plotting 10 representative particle trajectories in figure 11(b–d). Particles oscillate along individual trajectories, but their average behaviour is akin to that predicted using the steady rotlet model. When the external flow is zero (
$\dot {{\gamma }}=0$
), the particles emitted from the source follow oscillatory trajectories forming closed loops due to the time-dependent ciliary motion, only escaping due to diffusion (figure 11
b). When the direction of far-field ciliary flow and the external flow of strength
$\dot {\gamma }=0.5$
s
$^{-1}$
coincide (figure 11
c), the recovery stroke generated by the cilium traps particles in closed loops before they are released due to both diffusion and external flow. By contrast, when the shear flow of strength
$\dot {\gamma }=0.5$
s
$^{-1}$
has the opposite direction to the far-field ciliary flow (figure 11
d), particles are advected downstream due to flow superposition.
While there is a qualitative resemblance between the time-averaged rotlet (figure 6) and the time-dependent colloidal rotor (figures 11
c and 11
d) models, we proceed to quantify the particle currents across a range of input parameters and external flow strengths. We measure the particle current,
$I$
, as a function of diffusivity,
$D_0$
, for various shear rates as shown in figure 12. In the absence of external shear (figure 12
a), the particles with higher diffusivity are more likely to escape during one period of a sphere as compared with particles with lower diffusivity. This slightly overestimates the particle current compared with the rotlet model due to the time-averaged approximation of ciliary flow. When external shear is present (figure 12
b–d), a larger emission of particles during one period of the sphere is observed, particularly at lower shear rates. As the shear rate increases, particles are predominantly advected by the external flow, resulting in closer agreement in the particle current profiles between the time-averaged and time-dependent flow models. Taken together, these results demonstrate that the steady (rotlet) model captures the key features of nutrient transport by unsteady flows across a range of different flow conditions and particle diffusivities.
Open access
