1 Introduction
Cilia and flagella are contractile hair-like structures, put into motion by biochemical energy, and protruding on the free surfaces of eukaryotic or prokaryotic cells. While prokaryotes are single-celled organisms, eukaryote cells have membrane-bound organelles and are found in every mammal. Many living organisms, going from prokaryotic bacteria to mammals, use ciliary and/or flagellar propulsion as a swimming mechanism. Usually, flagella are external appendices used by microswimmers such as the alga Chlamydomonas reinhardtii for locomotion purposes, while cilia are generally internal appendices, shorter and more numerous, used for moving materials such as nutrients, dust or proteins into living organisms. Ciliary propulsion is a universal phenomenon, and many examples could be cited. For the particular case of the human body, cilia are responsible for the left–right asymmetry of the heart in early embryonic development, for the transport of nutrients in the brain and for the transport of mucus in the mucociliary clearance process, which is the background of the present work (see Satir & Christensen (Reference Satir and Christensen2007) for a review about the structure and function of mammalian cilia).
During the breathing process, a large number of foreign particles (bacteria, dust, pollutants or allergens) can penetrate the organism. The human body has then developed three mechanisms to protect itself from these particles: coughing, alveolar clearance and mucociliary clearance, which occurs on the epithelial surface of the respiratory system. In order to trap the particles, a layer of fluid called the airways surface liquid (ASL) covers the epithelial surface. Due to differences in the concentration of mucins inside the ASL, it is generally assumed to be the superposition of two different layers: the periciliary layer (PCL) and mucus. The PCL, of
$7~\unicode[STIX]{x03BC}\text{m}$
depth, is located between the epithelial surface from where the cilia protrude and the mucus phase just above it. Mainly composed of water and a few mucins of low molecular weight, this phase, often considered as being a Newtonian fluid similar to water, is a kind of lubricant allowing the mucus to slip on it and the cilia to beat without too much viscous resistance.
Mucus is a highly non-Newtonian fluid, with a strong viscoelastic behaviour. Additionally, the inner structure of the macromolecules composing it confers to the mucus a clear thixotropic behaviour, which is currently being studied and characterized to understand its inner rheological properties which exhibit a large variability (Lafforgue et al.
Reference Lafforgue, Poncet, Seyssiecq-Guarente and Favier2016). Mucus is indeed composed of 95 % water, but also contains macromolecules called the mucins (Lai et al.
Reference Lai, Wang, Wirtz and Hanes2009). It serves as a physical barrier against infectious agents and dust, but also to humidify the air flowing into the respiratory system and to catch particles. Its height varies between 5 and
$100~\unicode[STIX]{x03BC}\text{m}$
depending on many factors, including the position in the respiratory system (Widdicombe & Widdicombe Reference Widdicombe and Widdicombe1995), the pathology for a particular person and the quality of the air inhaled. Its viscosity can also vary by several orders of magnitude within the same day (Kirkham et al.
Reference Kirkham, Sheehan, Knight, Richardson and Thornton2002).
Cilia are the other protagonists of the mucociliary process. They are organized as tufts (around 200–300 cilia per tuft) at the epithelium surface. Formed by nine pairs of microtubules placed regularly in a circle and one pair of microtubules at the centre (which form the so-called axoneme), their purpose is to propel the surrounding fluid layers. Their diameter varies from 0.2 to
$0.3~\unicode[STIX]{x03BC}\text{m}$
and their length from 6 to
$7~\unicode[STIX]{x03BC}\text{m}$
(Sleigh, Blake & Liron Reference Sleigh, Blake and Liron1988).
The cilium motion can be decomposed into two steps: a stroke phase and a recovery phase. The stroke phase is characterized by almost straight cilia orthogonal to the flow in order to maximize the pushing effect, while the cilia are bending during the recovery phase in a more inclined plane to get closer to the epithelial surface in order to minimize the viscous resistance and therefore to reduce their impact on the flow. During the stroke phase, which takes around
$1/3$
of the total beating period, the tips of the cilia enter the mucus phase (Widdicombe & Widdicombe Reference Widdicombe and Widdicombe1995). Their beating frequency is estimated to vary between 10 and 20 Hz. It should be noted that the spatial asymmetry is essential for the cilia to generate propulsion in creeping flows, while the temporal asymmetry (recovery phase longer than stroke phase) is not necessary to induce a mucus motion (Khaderi et al.
Reference Khaderi, Baltussen, Anderson, den Toonder and Onck2010).
It has been experimentally observed that cilia usually do not beat randomly (Sleigh Reference Sleigh1962), but instead adapt their beating according to their neighbours, giving birth to the so-called metachronal waves (MCWs) observed at the tips of cilia. Metachronal waves occur when adjacent cilia beat with a constant phase lag
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
between one another. For
$0<\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}<\unicode[STIX]{x03C0}$
, the MCWs move in the opposite direction to the fluid propelled and are called antipleptic MCWs. On the contrary, for
$-\unicode[STIX]{x03C0}<\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}<0$
, the MCWs move in the same direction as the flow and are called symplectic MCWs. When the phase lag
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
is null, all cilia beat in a synchronized way. Finally, when the phase lag between neighbouring cilia is
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\pm \unicode[STIX]{x03C0}$
, a standing wave appears.
The universality of ciliary propulsion has intrigued scientists for decades, and several studies have been conducted in order to understand it. One actual objective is to be able to mimic this process in order to, for instance, create cilium-based actuators for mixing, and use them as flow regulators in microscopic biosensors or as micropumps for drug-delivery systems (Li, Tan & Zhang Reference Li, Tan and Zhang2009; Chen et al. Reference Chen, Chen, Lin and Hu2013). Moreover, diseases such as asthma and chronic obstructive pulmonary disease (COPD) are related to the mucociliary clearance process (Gardiner Reference Gardiner2005). The objectives are to understand the underlying mechanism that allows hundreds of cilia to act as a whole for the transport of mucus and how it affects the flow generated. The results could bring a deeper understanding to such pulmonary diseases.
Numerically, Gueron and co-authors (Gueron et al.
Reference Gueron, Levit-Gurevich, Liron and Blum1997; Gueron & Levit-Gurevich Reference Gueron and Levit-Gurevich1999) showed in the 1990s that two neighbouring cilia beating randomly will quickly synchronize within a few beating cycles due to hydrodynamic interactions and form antipleptic MCWs. This has been recently confirmed by Elgeti & Gompper (Reference Elgeti and Gompper2013), who observed the formation of MCWs in a 3D single-phase environment. A certain degree of freedom in the beating pattern is required, as shown in the theoretical work of Niedermayer, Eckhardt & Lenz (Reference Niedermayer, Eckhardt and Lenz2008). They modelled the beating pattern of cilia by circular trajectories. By allowing some flexibility in the radii, they managed to introduce the coupling leading to MCW formation. Among the different models used for the study of ciliary propulsion, the envelope model (Taylor Reference Taylor1951; Reynolds Reference Reynolds1965; Tuck Reference Tuck1968; Blake Reference Blake1971a
,Reference Blake
b
; Brennen & Winet Reference Brennen and Winet1977) assumes that cilia are so densely packed that it is possible to consider their tips as an oscillating surface. Nevertheless, such a configuration has only been observed in nature for symplectic metachrony. Moreover, this technique is limited to small-amplitude oscillations and imposes no-slip and impermeability conditions at the oscillating surface. In the sublayer (or stokeslets) model (Keller & Brennen Reference Keller and Brennen1968; Blake Reference Blake1972; Lighthill Reference Lighthill1976; Phan-Thien, Tran-Cong & Ramia Reference Phan-Thien, Tran-Cong and Ramia1987; Gueron et al.
Reference Gueron, Levit-Gurevich, Liron and Blum1997; Gueron & Levit-Gurevich Reference Gueron and Levit-Gurevich1999; Smith, Gaffney & Blake Reference Smith, Gaffney and Blake2007; Niedermayer et al.
Reference Niedermayer, Eckhardt and Lenz2008; Gauger, Downton & Stark Reference Gauger, Downton and Stark2009; Ding et al.
Reference Ding, Nawroth, McFall-Ngai and Kanso2014) the cilia are modelled by a distribution of stokeslets, which impose a force on the surrounding fluid. A proper mirror image of the stokeslets is required to impose the no-slip condition on the surface from where the cilia protrude. However, the presence of a wall is known to alter the nature of the far field of the stokeslets, and it is thought that this can have important consequences on the hydromechanics of the cilia near the wall (Blake & Chwang Reference Blake and Chwang1974). Moreover, stokeslets can only be used for fluids with constant viscosities. Results using this method tend to show that symplectic metachrony would be more efficient for mucus transport than synchronously beating cilia, and that antipleptic metachrony would induce a lower flow rate. The opposite result was nevertheless obtained by Gauger et al. (Reference Gauger, Downton and Stark2009), who found that antipleptic metachrony was more efficient than symplectic metachrony for a particular cilium spacing. For a distance between two cilia of
$1.5L$
,
$L$
being the length of the cilia, they obtained an increase in the pumping performance of 40 % relative to a single cilium. However, while the beating pattern used in (Gauger et al.
Reference Gauger, Downton and Stark2009) was realistic, they considered a slow stroke phase and a quick recovery phase, which is the opposite of what is observed in nature. Some authors (Gueron et al.
Reference Gueron, Levit-Gurevich, Liron and Blum1997; Gueron & Levit-Gurevich Reference Gueron and Levit-Gurevich1999; Kim & Netz Reference Kim and Netz2006) showed that the energy spent by a cilium would decrease in the presence of metachronal motion. Others tried to model the internal axoneme of a cilium. Among them, Mitran (Reference Mitran2007) used an overlapping fixed–moving grid formulation, coupling the finite volume and the finite element methods, to study the emergence of MCWs. His model was very detailed and had a lot of advantages (two-layer flow, viscoelastic mucus) but required many experimental parameters, and the impact of the MCWs on the flow has not been considered so far. Niedermayer et al. (Reference Niedermayer, Eckhardt and Lenz2008) observed MCW formation in 1D cilium arrays. Those waves were stable only if the wavelength was four times higher than the cilium spacing. For the interested reader, a clear review of these computational modellings of the internal axoneme can be found in Fauci & Dillon (Reference Fauci and Dillon2006). A different approach is to let the cilia adapt their motion in order to find the most energetically efficient ciliary beating pattern. In that context, Eloy & Lauga (Reference Eloy and Lauga2012) and Lauga & Eloy (Reference Lauga and Eloy2013) computed the shape and energy-optimal kinematics of cilia from an energetic point of view. They found that the optimal kinematics strongly depends on the cilium bending rigidity, and closely resembles the two-stroke ciliary beating pattern observed in natural cilia. Similarly, Osterman & Vilfan (Reference Osterman and Vilfan2011) computed the ciliary beating pattern with optimal pumping efficiency of isolated cilia and arrays. Recently, new methods have been introduced, such as the immersed-boundary (IB) method used by Dauptain, Favier & Bottaro (Reference Dauptain, Favier and Bottaro2008) to model the swimming of pleurobrachia. Lukens, Yang & Fauci (Reference Lukens, Yang and Fauci2010) used an IB method to study the mixing produced by a carpet of cilia in the context of the mucociliary clearance process. A coupled IB lattice Boltzmann method (LBM) was also used by Sedaghat et al. (Reference Sedaghat, Shahmardan, Norouzi, Jayathilake and Nazari2016) to study several parameters in a 2D configuration using an Oldroyd-B model for the mucus rheology. They found that the transport of mucus was maximized when considering the mucus as a Newtonian fluid. While many studies regarding the emergence of MCWs have been conducted, only a few have addressed 3D configurations of the mucociliary clearance process. Among them, Elgeti & Gompper (Reference Elgeti and Gompper2013) managed to observe symplectic and laeoplectic (perpendicular to the power-stroke direction) MCW formations. Ding et al. (Reference Ding, Nawroth, McFall-Ngai and Kanso2014) did not study the emergence of MCWs but performed a comparative study of antipleptic metachrony versus symplectic metachrony in terms of transport efficiency and mixing. Their results showed that both antipleptic and symplectic MCWs enhanced the fluid transport and the mixing, with the antipleptic waves being the most efficient ones. However, the two aforementioned studies only considered a single fluid layer.
In this work, by using the solver developed by Li et al. (Reference Li, Favier, D’Ortona and Poncet2016) and already validated in similar configurations, the focus is placed on the emergence of MCWs in a 3D two-phase flow configuration with a viscosity ratio of 20 and a large number of cilia. In particular, it is shown that a simple hydrodynamical feedback based on mechanical concepts can trigger the emergence of both symplectic and antipleptic MCWs, while usually only a single layer of fluid is considered and only antipleptic MCWs are seen to emerge. From a numerical point of view, the main advantage of the present method is its ease of implementation. Additionally, the local character of the collisions in the LBM allows an easy and straightforward parallelization of the code, making the simulation of a large number of cilia possible. The numerical method also possesses the following advantages: (i) viscosity ratios up to
$O(10^{2})$
can be achieved (Porter et al.
Reference Porter, Coon, Kang, Moulton and Carey2012) and (ii) the mucus–PCL interface emerges intrinsically from the model. To the best of the authors’ knowledge, this solver is the only one that combines all of these capabilities. The main contribution of this work is the thorough analysis of the advantages of the antipleptic and symplectic MCWs over synchronized beating by computing appropriate transport and efficiency ratios, mapping an inter-cilium spacing from 1.67 to 5 cilium lengths. Finally, and for the first time, both the PCL and the mucus layer have been taken into account in the present study. The inherent advantages of the MCWs for flow transport are studied by (i) considering the efficiency of the waves in transporting the flow, (ii) comparing the flow generated and the volume transported by the different kinds of metachrony, (iii) comparing the energetic cost of the MCWs and (iv) analysing the capacity of the waves to transport particles from an energetic point of view. The following analysis shows that antipleptic MCWs are always the most efficient in transporting mucus.
The remainder of this paper is organized as follows. In § 2, details about the numerical method are given, and the different quantities used are introduced. In § 3, the results are presented, starting from the emergence of MCWs by considering the fluid feedback onto the cilia; then, a parametric study is carried out to quantify the impact of both the MCWs and the cilium spacing on the transport of mucus. A summary of the results and the future perspectives concludes this paper in § 4.
2 Numerical method
From the numerical point of view, modelling of the complete problem under realistic conditions remains a huge challenge, and several assumptions need to be made to simplify the problem, while keeping the essential ingredients of the physical mechanisms in the model. In particular, it is assumed that there is no mass transfer at the wall, that both fluids are Newtonian and that the cilia are equally spaced filament-like structures. More assumptions are detailed in this section, together with the numerical set-up and geometry.
2.1 Geometrical modelling and beating pattern
The computational domain is a box with a regular mesh composed of
$N_{x}\times N_{y}\times N_{z}$
points. The cilia are equally spaced along the bottom (
$x$
–
$y$
) wall, such that their base points are located at
$z=0$
(see figure 1
a for a schematic view of the geometry). Their motion is imposed to be in the
$x$
-direction only. The spacing between two adjacent cilia is denoted
$a$
in the
$x$
- and
$y$
-directions. The wavelength
$\unicode[STIX]{x1D706}$
of the MCW is set such that
$\unicode[STIX]{x1D706}=N_{x}$
when the metachrony is forced. For the case of synchronously beating cilia,
$N_{x}$
is chosen to be
$8a$
. The length of the cilia
$L$
is set to 15 lattice units (lu), except when stated differently, and 200 snapshots in time per beating cycle are uniformly distributed to model their motion. The ratio
$h/H$
between the PCL thickness and the height of the domain is fixed to 0.27 for all simulations.
Figure 1. (a) Schematic view of the computational domain. The present case corresponds to an antipleptic MCW. The domain is filled with PCL (in blue) and mucus (in red). (b) Beating pattern of a cilium with the parametric equation used. Steps 1–6 correspond to the recovery phase and steps 7–9 to the stroke phase.
The equations of motion for the cilia are inspired from Chatelin (Reference Chatelin2013) and reproduce the beating pattern by resolving a 1D transport equation along a parametric curve. Let
$P(\unicode[STIX]{x1D701},t)$
be the position of the curve at time
$t$
and at a normalized distance
$\unicode[STIX]{x1D701}$
from the base point of a cilium. With appropriate boundary conditions, a realistic beating pattern can be obtained using the following transport equation:
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}P^{\prime }}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D708}(t)\frac{\unicode[STIX]{x2202}P^{\prime }}{\unicode[STIX]{x2202}\unicode[STIX]{x1D701}}=0,\quad \text{BC:}~\left\{\begin{array}{@{}l@{}}P(0,t)=(0,0,0),\\ P^{\prime }(0,t)=(2\cos (2\unicode[STIX]{x03C0}t/T),0,\cos (2\unicode[STIX]{x03C0}t/T)),\end{array}\right.\end{eqnarray}$$
where
$\unicode[STIX]{x1D708}(t)=[1+8\cos ^{2}(\unicode[STIX]{x03C0}(t+0.25T)/T)]/T$
is the viscosity of the surrounding fluid,
$T$
is the beating period and
$P^{\prime }=\unicode[STIX]{x2202}_{\unicode[STIX]{x1D701}}P$
. The resulting angular amplitude between the beginning and the end of a stroke phase is
$\unicode[STIX]{x1D703}=2\unicode[STIX]{x03C0}/3$
, which agrees well with experimental data (Sleigh et al.
Reference Sleigh, Blake and Liron1988). It should be noted that a 3D beating pattern would allow more realistic simulations to be achieved while being more computationally expensive. Hence, the choice has been made to use this 2D beating pattern. It captures the essential ingredients of the beating, and as cilia have a diameter smaller than a lattice unit, the difference in the induced flow between two cilia overlapping in 2D or slipping onto each other in 3D is very small. Figure 1(b) gives a view of the beating pattern obtained by resolving (2.1a,b
). It should be noted that with the present model, both phases take the same amount of time while moving in the same (
$x,y$
) plane. This choice has been made in order to only study spatially asymmetric motions, which are the only mechanisms effective at low Reynolds numbers. The temporal asymmetry is indeed a mechanism that can enhance the flow only when inertial forces are no longer negligible (Khaderi et al.
Reference Khaderi, Baltussen, Anderson, den Toonder and Onck2010). Nevertheless, when the feedback of the fluid is taken into account (as in § 3.1), a non-symmetrical motion will develop, with a stroke phase slower than the recovery phase. More details regarding this temporal asymmetry are given in § 2.3.
The PCL is set such that it fills the region going from the bottom of the domain (
$z=0$
) up to an altitude of
$h$
. In all of the simulations, the value of
$h=0.9L$
has been used in order to allow the tips of the cilia to emerge into the mucus layer during their stroke phase, as observed in real epithelium configurations.
Both the PCL and the mucus are considered to be Newtonian fluids. The kinematic viscosity of the mucus is
$\unicode[STIX]{x1D708}_{m}=10^{-3}~\text{m}^{2}~\text{s}^{-1}$
and the viscosity ratio
$r_{\unicode[STIX]{x1D708}}=\unicode[STIX]{x1D708}_{m}/\unicode[STIX]{x1D708}_{PCL}$
between the mucus and the PCL is set to 20. It has indeed been recently shown (Chatelin & Poncet Reference Chatelin and Poncet2016) that mucus transport is maximized for viscosity ratios ranging from 10 to 20 with a stiff transition between the two fluid layers. The beating period of the cilia is equal to
$N_{it}\times \text{d}t$
(with
$\text{d}t=1$
using the classical LBM normalization),
$N_{it}$
being the number of iterations for a cilium to perform a complete beating cycle. An oscillatory Reynolds number
$Re^{osc}$
, based on the velocity of the cilium tips,
$U_{cil}=2\unicode[STIX]{x1D703}L/T_{osc}$
, can now be defined:
$$\begin{eqnarray}Re^{osc}=\frac{U_{cil}L}{\unicode[STIX]{x1D708}_{mucus}}=\frac{{\displaystyle \frac{4\unicode[STIX]{x03C0}/3}{T_{osc}}}L^{2}}{\unicode[STIX]{x1D708}_{mucus}}=\frac{\unicode[STIX]{x1D714}L^{2}}{\unicode[STIX]{x1D708}_{mucus}},\end{eqnarray}$$
where
$\unicode[STIX]{x1D714}$
is the angular frequency of the cilium beating. With realistic physical quantities corresponding to the ciliated epithelium surface, the value of
$Re^{osc}$
is of the order of
$10^{-5}$
(using
$L\approx 10^{-5}$
m,
$\unicode[STIX]{x1D708}_{mucus}\approx 10^{-3}~\text{m}^{2}~\text{s}^{-1}$
and
$U_{cil}\approx 10^{-3}~\text{m}~\text{s}^{-1}$
).
To avoid running simulations at such a low Reynolds number, which would require a very high number of iterations using a lattice Boltzmann (LB) scheme, a higher Reynolds number was chosen:
$Re=20$
. Indeed, the achievement of simulations at
$Re<1$
using an LB scheme, while still describing the cilia well (
$L=15$
lu at least), requires a large amount of CPU resources. Thus, inertial effects are introduced in the model, but it has been carefully checked that the results remain the same for creeping flows (
$Re=O(10^{-2})$
) by comparing them with the results obtained with the LBM formulation designed for Stokes flows (Zou et al.
Reference Zou, Hou, Chen and Doolen1995; Guo & Shu Reference Guo and Shu2013). The differences are found to be less than 7 % of the transported velocity for all phase lags
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}\neq 0$
. The corresponding beating patterns of the cilia shown in the following are similar also in terms of vorticity generation, and the same qualitative coordination concerning the antipleptic/symplectic behaviour is observed. However, for the particular case of synchronously beating cilia (i.e.
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=0$
), the inertial effects will play a non-negligible role in the flow dynamics. Despite being weak, they will indeed cancel the reversal of the flow that should occur when the cilia are in the recovery phase. More details regarding the inertial effects will be given in § 3.2.7.
Finally, it is worth noticing that the lattice has been chosen such that the numerical diffuse diameter of the cilia due to the IB method corresponds to the diameter of real cilia (
${\approx}0.3~\unicode[STIX]{x03BC}\text{m}$
). Hence, a realistic drag is taken into account in the present study.
2.2 Algorithm
The numerical model is described in Li et al. (Reference Li, Favier, D’Ortona and Poncet2016), and validated on several configurations involving flexible and moving boundaries in multiphase flows, with a second-order accuracy. Briefly, the idea is to add a forcing term
$F_{i}^{\unicode[STIX]{x1D70E}}=\boldsymbol{F}_{\unicode[STIX]{x1D70E}}^{SC}+\boldsymbol{F}_{\unicode[STIX]{x1D70E}}^{IB}$
to the discrete LB equation for each
$\unicode[STIX]{x1D70E}$
fluid component. The term
$\boldsymbol{F}_{\unicode[STIX]{x1D70E}}^{SC}$
is an interparticle potential force that takes into account the fluid–fluid cohesion forces (Shan & Chen Reference Shan and Chen1994), and
$\boldsymbol{F}_{\unicode[STIX]{x1D70E}}^{IB}$
is an IB-related force to ensure the no-slip condition at the fluid–solid interface.
The fluid part is first solved on a Cartesian grid with the LBM using the Bhatnagar–Gross–Krook (BGK) operator and a D3Q19 scheme. The collision and streaming steps proper to the LBM are first performed. The model of Porter et al. (Reference Porter, Coon, Kang, Moulton and Carey2012) is used to model the two-phase flow as it allows minimization of the magnitude of spurious currents near the fluid–fluid interface. More importantly, it also allows the consideration of higher density or viscosity ratios. Then, values for the fluid velocity are interpolated at the Lagrangian points. This allows the computation of an IB force to be spread onto the neighbouring Eulerian fluid nodes in order to ensure the no-slip condition along the cilia. The macroscopic fluid velocity is then updated. The motion of each cilium is decomposed into a finite number of steps (snapshots) during a period. If necessary, an interpolation can be carried out in order to have the velocity values along the cilia in between two steps. It should be noted that the geometric shape of the beating is fixed in all simulations and is not impacted by the feedback law introduced in § 2.3, which only affects the duration of the recovery and stroke phases.
Since the model of Porter et al. (Reference Porter, Coon, Kang, Moulton and Carey2012) uses a Shan–Chen (SC) repulsive force (Shan & Chen Reference Shan and Chen1993, Reference Shan and Chen1994), surface tension effects emerge intrinsically at the PCL–mucus interface. Hence, a sharp interface between the mucus and the PCL can be maintained at any time. However, small diffusion effects might occur on the lattices bordering the interface. Additionally, the cilia that enter the PCL–mucus interface may temporarily induce a small mixing. This is, however, corrected by the SC force which ‘unmixes’ the two fluids.
Periodic boundary conditions are used in the
$x$
- and
$y$
-directions, while no-slip and free-slip boundary conditions are used at the bottom and top walls respectively. The size of the computational domain ranges from 50 lu to 400 lu depending on the configuration considered, except for the size along the
$z$
-direction which is always set to 50 lu.
Taking advantage of the local character of the LBM algorithm, the code is parallelized using MPI (message passing interface) libraries by splitting the full computational domain into nine subdomains of size
$(N_{x}/3,N_{y}/3,N_{z})$
. More details on the numerical model can be found in Li et al. (Reference Li, Favier, D’Ortona and Poncet2016).
2.3 Feedback of the fluids onto the cilia
The basic idea is to modulate the beating motion of the cilia as a function of the fluid motion. To do so, it is assumed that all cilia follow the same beating pattern; meanwhile, a feedback of the fluids, which consists of accelerating or slowing down the motion of the cilia, is introduced.
Each cilium is discretized with
$N_{s}=20$
Lagrangian points. Let
$s$
be the subscript corresponding to the
$s\text{th}$
Lagrangian point, starting from the base tip at
$s=1$
, and
$\boldsymbol{V}_{i}^{s}$
the velocity on the
$s\text{th}$
Lagrangian point of the
$i\text{th}$
cilium. For each cilium, we define the average velocity over all Lagrangian points
$\boldsymbol{V}_{i}$
, which is linked to the number of steps (snapshots) this cilium will skip during one iteration of the fluid solver. The fluid feedback onto the cilia thus consists of modifying the norm of the velocity vector
$\Vert \boldsymbol{V}_{i}\Vert$
, while its direction remains unchanged.
Figure 2. Schematic view of a cilium with the corresponding forces exerted on the fluids.
The feedback is computed in three steps. First, the IB forces corresponding to mucus and PCL are projected onto the corresponding velocity vectors for each Lagrangian point. Then, an estimate of the feedback is computed based on the torques of the forces for each Lagrangian point. Finally, the beating pattern of the cilia is adjusted at the beginning of the next time step. The different forces and geometrical variables used are illustrated in figure 2, and the three steps are explained below.
-
(i) The interpolated IB forces applied by the
$i\text{th}$
cilium onto the fluids – respectively
$\boldsymbol{F}_{m}^{i}$
for the force imposed on the mucus phase and
$\boldsymbol{F}_{PCL}^{i}$
for the force imposed on the PCL phase – are projected on
$\boldsymbol{V}_{i}^{s}$
: (2.3)where
$$\begin{eqnarray}\boldsymbol{F}_{\ast }^{i,proj}=\frac{\boldsymbol{F}_{\ast }^{i}\boldsymbol{\cdot }\boldsymbol{V}_{i}^{s}}{\Vert \boldsymbol{V}_{i}^{s}\Vert ^{2}}\boldsymbol{V}_{i}^{s},\end{eqnarray}$$
$\ast$
stands for ‘
$m$
’ or ‘
$PCL$
’ depending on the position of the Lagrangian node. In order to take into account the difference of viscosity between the two layers, the forces
$\boldsymbol{F}_{m}^{i,proj}$
and
$\boldsymbol{F}_{PCL}^{i,proj}$
are weighted by a term of the form
$\unicode[STIX]{x1D70F}_{\ast }/(\unicode[STIX]{x1D70F}_{m}+\unicode[STIX]{x1D70F}_{PCL})$
. The total projected force of the fluids onto the
$s\text{th}$
segments of the
$i\text{th}$
cilium
$\boldsymbol{F}_{fluids\rightarrow cilia}^{i}$
is written as (2.4)
$$\begin{eqnarray}\boldsymbol{F}_{fluids\rightarrow cilia}^{i}=-\frac{(\unicode[STIX]{x1D70F}_{m}\boldsymbol{F}_{m}^{i}+\unicode[STIX]{x1D70F}_{PCL}\boldsymbol{F}_{PCL}^{i})\boldsymbol{\cdot }\boldsymbol{V}_{i}^{s}}{(\unicode[STIX]{x1D70F}_{m}+\unicode[STIX]{x1D70F}_{PCL})\Vert \boldsymbol{V}_{i}^{s}\Vert ^{2}}\boldsymbol{V}_{i}^{s}.\end{eqnarray}$$
-
(ii) For each Lagrangian point
$s$
, the norm of the torque of
$\boldsymbol{F}_{fluids\rightarrow cilia}^{i}$
with respect to the base point
$O$
of the
$i\text{th}$
cilium is computed by (2.5)where
$$\begin{eqnarray}\Vert {\mathcal{M}}_{O}^{i}(\boldsymbol{F}_{fluids\rightarrow cilia}^{i})\Vert =\Vert \boldsymbol{F}_{fluids\rightarrow cilia}^{i}\Vert L_{p},\end{eqnarray}$$
$L_{p}=\Vert \mathbf{X}_{p}\otimes \boldsymbol{V}_{i}^{s}\Vert /\Vert \boldsymbol{V}_{i}^{s}\Vert$
is the lever arm. The total fluid feedback onto the
$i\text{th}$
cilium considered is then computed by summing the computed quantities over all Lagrangian points: (2.6)The velocity of each cilium is finally modified as follows:
$$\begin{eqnarray}T^{i}=\mathop{\sum }_{s=2}^{N_{segments}}\Vert \boldsymbol{F}_{fluids\rightarrow cilia}^{i}\Vert L_{p}.\end{eqnarray}$$
(2.7)where
$$\begin{eqnarray}\Vert \boldsymbol{V}_{i}\Vert =\Vert \boldsymbol{V}_{0}\Vert +\unicode[STIX]{x1D6FC}T^{i},\end{eqnarray}$$
$\boldsymbol{V}_{0}$
is the initial speed of the cilia and
$\unicode[STIX]{x1D6FC}$
is a coupling parameter that allows tuning of the feedback strength.
-
(iii) Then, at the beginning of the next iteration, the beating pattern of the
$i\text{th}$
cilium is adjusted: (2.8)where
$$\begin{eqnarray}N_{i}^{t}=\text{mod}~(N_{i}^{t-1}+\Vert \boldsymbol{V}_{i}\Vert ,N_{i}^{total}),\end{eqnarray}$$
$N_{i}^{total}$
is the total number of snapshots defining the beating pattern of the cilia,
$N_{i}^{t-1}$
is the previous snapshot and
$N_{i}^{t}$
is the new position of the
$i\text{th}$
cilium (in terms of number of snapshots) at the current iteration.
2.4 Injection of passive tracers
In the context of real epithelial systems in the human body, mucus acts as a barrier against particles and pollutants. To gain insight into how particles are dispersed and advected in the mucus and PCL, passive tracers are injected into the domain on a (
$y,z$
) plane, when the flow has reached an established regime (see figure 3 for a view of a domain filled with tracers). In figure 3, the tracers initially seeded into the PCL are coloured in green, and the ones initially seeded into the mucus are coloured in grey. Their displacements are then computed and averaged over several cilium beating cycles. Their equations of motion are solved by a second-order Runge–Kutta (RK2) scheme, using the interpolated fluid velocity at each time step and the same procedure as in the IB method.
By taking a closer look at figure 3, one can observe that the particles initially seeded into the PCL are greatly mixed while staying mostly in the PCL. On the contrary, the particles initially seeded into the mucus layer stay in the mucus layer and are not mixed. This is mainly due to the surface tension effects present at the mucus–PCL interface which prevent too strong a mixing at the interface. This shows how wetting particles that are deposited to the air–mucus interface might enter the mucus layer but never reach the PCL.
Figure 3. Domain filled with passive tracers. The tracers initially seeded into the PCL are coloured in green and the ones initially seeded into the mucus are coloured in grey. The present case corresponds to an antipleptic MCW with
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/4$
. An array of eight cilia with a cilium spacing
$a/L=1.67$
is considered on a computational domain of size (
$N_{x}=201$
,
$N_{y}=26$
,
$N_{z}=50$
). The mucus phase (in red) is located above the PCL phase (in blue), and the ratio of viscosity is set as
$r_{\unicode[STIX]{x1D708}}=10$
. The colour bar represents the dimensionless velocity magnitude of the fluids on the periodic boundary in the
$x$
-direction.
3 Results and Discussion
3.1 Emergence of MCWs
Using the feedback force introduced in § 2.3, randomly beating cilia are observed to synchronize with their immediate neighbours, giving birth to symplectic or antipleptic metachrony. The time for the synchronization to occur depends on the set of parameters used. Both local and global synchronizations can be observed.
The empirical parameter
$\unicode[STIX]{x1D6FC}$
plays a role in the emergence of the waves. When
$\unicode[STIX]{x1D6FC}$
is set to too low a value (
$|\unicode[STIX]{x1D6FC}|<0.5$
), the cilia will beat randomly, while if set to a too high value (
$|\unicode[STIX]{x1D6FC}|>7$
for the case
$a/L=1.67$
, for example), the cilia will fully synchronize with each other without any phase lag. Additionally, higher absolute values of
$\unicode[STIX]{x1D6FC}$
usually decrease the times needed for MCWs to emerge.
Figure 4 shows a symplectic MCW emerging from an initially random state (every cilium was initially beating at a random step of its beat cycle). In the present case, three rows of eight cilia on a computational domain of size (
$N_{x}=361$
,
$N_{y}=136$
,
$N_{z}=50$
) are considered. The PCL is set such that
$h=0.9L$
and the ratio of viscosity between the PCL and the mucus phase is 15. The feedback coefficient
$\unicode[STIX]{x1D6FC}$
used is
$\unicode[STIX]{x1D6FC}=-3.5$
. The spacing between two neighbouring cilia is
$a/L=3$
in the
$x$
- and
$y$
-directions, with
$L=15$
lu. Hence, cilia never collide and overlapping of kernels from neighbouring cilia never occurs. One can clearly see the symplectic MCW that emerged from the initially random state of the cilia, with a wavelength
$\unicode[STIX]{x1D706}=180$
lu and a phase lag
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}\approx -\unicode[STIX]{x03C0}/2$
, confirming that hydrodynamic interactions only suffice to account for the emergence of MCWs.
Figure 4. Symplectic MCW emerging from an initially random state of the cilia. Here, 24 cilia arranged in an
$8\times 3$
rectangle are considered on a computational domain of size (
$N_{x}=361$
,
$N_{y}=136$
,
$N_{z}=50$
) with a cilium spacing of
$a/L=3$
. The mucus phase is in red and the PCL phase is in blue. The colour bar indicates the phase of a particular cilium within one beating period, which is represented by a circle at its base. (a) Three-dimensional view of the system. (b) Two-dimensional view of the same system in an (
$x,y$
) plane to highlight the 3D modulation in the
$z$
-direction.
Figure 5 shows a similar configuration (
$h=0.60L$
,
$r_{\unicode[STIX]{x1D708}}=15$
,
$\unicode[STIX]{x1D6FC}=-3.5$
), where 1024 cilia arranged in a
$32\times 32$
square are considered on a computational domain of size (
$N_{x}=161$
,
$N_{y}=161$
,
$N_{z}=32$
). One can clearly see the formation of an antipleptic MCW with a wavelength of
$\unicode[STIX]{x1D706}=80$
lu and a phase lag of
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}\approx \unicode[STIX]{x03C0}/8$
. It should be noticed that, in the simulation presented in figure 5, the cilium spacing has been reduced (
$a/L=0.23$
and
$L=22$
lu), to take into account a larger number of cilia. Thus, neighbouring cilia may overlap, which can cause spurious numerical effects, due to the 2D beating pattern of the cilia, when two Lagrangian points with opposite velocities occupy the same mesh point. However, it has been verified in some simulations that these effects do not play any role in the emergence of the MCW by setting to zero the IB force contribution of these Lagrangian points, therefore cancelling the spurious effects that may occur at the overlapping points of neighbouring cilia. In other words, the fluid velocity at the Eulerian node surrounding these particular Lagrangian points was not modified by the IB method, proving that it is not the source of the cilium synchronization. As shown in figures 4 and 5, weakly 3D effects may be observed but do not play a key role in the physics of emerging waves. In all of the simulations presented, the synchronization along the
$x$
-direction is quite strong. In some simulations, particular cilia may lose synchronization over time, but quickly readjust their beating according to the others. It also appears that once an equilibrium is reached, the synchronization along the
$y$
-direction is stable too. A remarkable fact is that for particular configurations (like the one displayed in figure 5), the interface between the mucus and the PCL also forms a wave travelling in the same direction as the MCW. If the mucus viscosity is increased, the torques felt by the cilia become stronger. Hence, the cilia are either more accelerated (
$\unicode[STIX]{x1D6FC}>0$
) or decelerated (
$\unicode[STIX]{x1D6FC}<0$
), and the mucus will adapt its displacement accordingly. However, when the stationary regime is reached, the clearance velocity should not be strongly modified, as shown by Chatelin & Poncet (Reference Chatelin and Poncet2016). According to these authors, the mucus velocity is decreased by only 50 % when the viscosity ratio is increased by a factor of 100 000. Preliminary results (not shown) seem to indicate that, assuming a least-effort behaviour of the cilia (meaning that
$\unicode[STIX]{x1D6FC}<0$
), antipleptic waves are obtained for small cilium spacings (
$a/L\leqslant 1.5$
) while symplectic waves are seen to emerge for larger cilium spacings (
$a/L\geqslant 1.5$
). Since, in nature, antipleptic waves are often observed for densely packed cilia, this implies that natural cilia adopt a least-effort behaviour. It also suggests the existence of a critical value for the cilium density,
$\unicode[STIX]{x1D70C}_{c}=\unicode[STIX]{x1D70C}(\unicode[STIX]{x1D6FC}_{c})$
, where
$\unicode[STIX]{x1D6FC}_{c}$
is the corresponding critical value of the coupling parameter, from which the kind of waves emerging can be controlled. On the contrary, with the present model, assuming that the cilia beat faster when encountering a resistance (meaning that
$\unicode[STIX]{x1D6FC}>0$
), symplectic waves are seen to emerge for small cilium spacings (
$a/L\leqslant 1.5$
) while antipleptic waves emerge for larger cilium spacings (
$a/L\geqslant 1.5$
), which is not observed in nature. Movies of the MCWs are provided as online supplementary movies available at https://doi.org/10.1017/jfm.2017.352.
Figure 5. Antipleptic MCW emerging from an initially random state of the cilia. Here, 1024 cilia arranged in a
$32\times 32$
square are considered on a computational domain of size (
$N_{x}=161$
,
$N_{y}=161$
,
$N_{z}=32$
) with a cilium spacing of
$a/L=0.23$
. The mucus phase is in red and the PCL phase is in blue. The colour bar indicates the phase of a particular cilium within one beating period, which is represented by a circle at its base. (a) Three-dimensional view of the system. (b) Two-dimensional view of the same system in an (
$x,y$
) plane to highlight the 3D modulation in the
$z$
-direction.
As already explained in § 2.1, when the feedback is taken into account, a non-symmetrical motion develops, with a stroke phase slower than the recovery phase for both antipleptic and symplectic MCWs. For the antipleptic case,
$\Vert \boldsymbol{V}_{i}\Vert$
is smaller during the stroke phase, and since
$\unicode[STIX]{x1D6FC}$
is negative and the torques
$T^{i}$
are always positive, this means that the values of
$T^{i}$
computed during the stroke phase are larger than the values computed during the recovery phase. This results in a weaker velocity for the cilia which cover fewer snapshots during the stroke phase. It is the opposite for the symplectic case where
$\unicode[STIX]{x1D6FC}$
is positive:
$T^{i}$
takes larger values during the recovery phase as
$\Vert \boldsymbol{V}_{i}\Vert$
. In other words, the feedback depends on the clustering character of the motion: when cilia are clustered, the torque exerted by the fluids onto the cilia is weaker, and when they are far from each other (during the stroke phase for antipleptic motion and during the recovery phase for symplectic motion), the torques are stronger. This makes sense since the cilia encounter more viscous resistance when they are not clustered, as will be discussed later in § 3.2.4. The resulting motion of the cilia then differs from what is observed in nature, but it indicates that cilia experience stronger stresses during the stroke phase of antipleptic motion and during the recovery phase of symplectic motion. This supposes that the beating kinematics of real cilia is not dictated only by hydrodynamical interactions and suggests that other biological parameters or functions (such as sensing) may play a role. This is in agreement with Guo et al. (Reference Guo, Nawroth, Ding and Kanso2014), who compared the performance of pumping-specialized cilia and swimming-specialized cilia as a function of metachronal coordination, and found that the latter almost always outperform the pumping-specialized cilia. As will be further detailed in § 3.2.5, their results are also in accordance with the outcome of the present paper: antipleptic waves are the most efficient ones for the transport of fluids. Finally, it is worth noticing that the degree of asymmetry is much higher for antipleptic MCWs compared with symplectic MCWs.
3.2 Quantitative study of MCWs
This section presents the results obtained, once the flow is well established, for the three configurations studied: (i) synchronized case (all cilia beat together with no phase lag); (ii) symplectic MCW where two neighbouring cilia beat with a negative phase lag, i.e
$-\unicode[STIX]{x03C0}<\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}<0$
; and (iii) antipleptic MCW where two neighbouring cilia beat with a positive phase lag, i.e
$0<\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}<\unicode[STIX]{x03C0}$
. It should be noted that in all of the following results, contrary to § 3.1, the metachrony is imposed in order to study specific phase lags
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
. Hence, the size of the domains must be changed accordingly, since different phase lags
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
imply different wavelengths, and hence more or fewer cilia. However, the quantities presented here have been averaged over the domain, and there are no effects in changing the size of the box. In order to study only spatially asymmetric motions, the recovery and stroke phases now take the same amount of time in all of the following results. It should also be noted that the standard deviations are not displayed in any of the figures that will be presented as they are extremely small (less than 0.01 %) for all considered quantities. Hence, they do not give additional information. Schematic views of synchronized beating and metachronal motions are displayed in figure 6. Another configuration, corresponding to randomly beating cilia, is also represented in figure 6. One can observe that the synchronized motion of cilia creates vorticity only in the periciliary zone, whereas metachronal motions also induce vorticity into the mucus layer. It is worth noticing that for symplectic motion (case (b)), vortex trails emerge from the cilium tips performing their stroke phase. The same phenomenon is observed in the antipleptic case during the recovery motion of cilia. The presence of coherent vortices is clearly visible in both cases.
Figure 6. Different kinds of collective coordination for the beating cilia. Here, 32 cilia arranged in an
$8\times 4$
square are considered on a computational domain of size (
$N_{x}=241$
,
$N_{y}=121$
,
$N_{z}=50$
) in each case, with a cilium spacing of
$a/L=2$
. (a) Synchronized motion; (b) symplectic metachronal motion; (c) antipleptic metachronal motion; (d) no synchronization (random state of cilia). The figure shows contours of the magnitude of the dimensionless vorticity
$\Vert \unicode[STIX]{x1D74E}\Vert$
using a logarithmic scale. The black lines show the frontier between the PCL at the bottom and the mucus layer above.
3.2.1 Transport and mixing zone
Figure 7. Normalized average displacement in the
$x$
-direction as a function of
$z/L$
for (a)
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\pm \unicode[STIX]{x03C0}/9$
and
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\pm \unicode[STIX]{x03C0}/4$
with
$a/L=0.8$
, and (b)
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\pm \unicode[STIX]{x03C0}/4$
and
$a/L=1.67$
. As in Ding et al. (Reference Ding, Nawroth, McFall-Ngai and Kanso2014), the mixing zone extends from 0 to
$1L$
and the transport zone from
$1L$
to
$1.4L$
. In both cases (a) and (b), the mucus–PCL interface is located at
$z/L=0.9$
and is indicated by a dashed horizontal line.
The displacement field is calculated as
$\boldsymbol{d}(\boldsymbol{x})=\int _{0}^{\text{T}}\boldsymbol{u}(\boldsymbol{x}(t),t)\,\text{d}t$
, where
$\boldsymbol{x}$
is the position vector and
$\boldsymbol{u}$
is the fluid velocity. Its component in the
$x$
-direction as a function of the axial position is averaged over several beating cycles and displayed in figure 7. The transport and mixing zones defined in Ding et al. (Reference Ding, Nawroth, McFall-Ngai and Kanso2014) are clearly visible. This shows that the antipleptic MCWs perform better than the other kind of coordination in transporting the particles. It should be noted that, in figure 7(b), the two curves with
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\pm \unicode[STIX]{x03C0}/4$
have a similar shape, whereas in figure 7(a), symplectic MCWs induce a negative transport (i.e. a counter-flow) in the mixing zone for the same phase lag. For a given value of the phase lag, the transport depends on the density of cilia. One can also notice the importance of the phase lag by looking at figure 7(a): for a phase lag
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=-\unicode[STIX]{x03C0}/9$
, the symplectic MCWs induce no counter-flow in the mixing zone, while the opposite happens for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=-\unicode[STIX]{x03C0}/4$
. Finally, it is worth noticing that, as the cilium spacing increases, the influence of the kind of metachrony decreases. Eventually, for very large cilium spacings, one would expect that the three curves displayed in figure 7(b) would merge into a single one corresponding to the displacement generated by an isolated cilium. Additionally, even with a larger cilium spacing, the particles are better transported in case (b). As in Ding et al. (Reference Ding, Nawroth, McFall-Ngai and Kanso2014), the displacement reaches a plateau for an axial position of
$1.4L$
in both cases (a) and (b). The beating pattern used in Ding et al. (Reference Ding, Nawroth, McFall-Ngai and Kanso2014) is different, but similar trends are observed. Nevertheless, a major difference must be pointed out: for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=-\unicode[STIX]{x03C0}/4$
and
$a/L=1.67$
(figure 7
b), the present results show that the displacement induced by the symplectic wave has the same shape as the synchronized and antipleptic cases, i.e. it induces transport in the mucus phase even if it is smaller. On the contrary, Ding et al. (Reference Ding, Nawroth, McFall-Ngai and Kanso2014) obtained no transport at all for the symplectic wave in the mucus phase, but instead a peak of transport under the cilium tips in the PCL at
$z=0.7L$
due to the presence of a vortex-like structure both below and above the cilium tips for this particular value of the phase lag
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
. Since the main difference between the two studies is the beating pattern of the cilia (same phase lag and cilium spacing), this highlights the sensitivity of the system to the beating pattern.
The displacement in the
$z$
- and
$y$
-directions has also been analysed. The results show that the
$y$
-component of the displacement can be neglected and that the
$z$
-component of the displacement is almost null above the cilium tips, but not zero under them. As will be discussed later, the
$z$
-component of the displacement reaches its maximal value in regions where the shear rate is maximal too. This agrees with the existence of a mixing zone under the cilium tips.
By making an analogy with the strain-rate tensor classically used in solid mechanics, the gradient of the displacement field
$\unicode[STIX]{x1D735}d$
can be computed by considering each Eulerian node of the Cartesian grid as a passive tracer. Since the displacement in the transverse direction is very small, the
$y$
-component of the displacement field can be neglected. Henceforth, the gradient is computed over every (
$x,z$
) plane, and an average is then performed over the
$y$
-direction. Following the methodology described in Ding et al. (Reference Ding, Nawroth, McFall-Ngai and Kanso2014),
$\unicode[STIX]{x1D735}d$
is decomposed into an antisymmetric component
$\unicode[STIX]{x1D64D}=(\unicode[STIX]{x1D735}d-(\unicode[STIX]{x1D735}d)^{\text{T}})/2$
corresponding to rotation and a symmetric component
$\unicode[STIX]{x1D64E}=(\unicode[STIX]{x1D735}d+(\unicode[STIX]{x1D735}d)^{\text{T}})/2$
corresponding to shear deformation. The two eigenvalues of
$\unicode[STIX]{x1D64E}$
are of the form
$\pm \unicode[STIX]{x1D6FE}$
and indicate the rates of stretching (
$+\unicode[STIX]{x1D6FE}$
) and compression (
$-\unicode[STIX]{x1D6FE}$
). The unit eigenvector
$\boldsymbol{e}_{\unicode[STIX]{x1D6FE}}$
corresponding to the positive eigenvalue indicates the direction of stretching and the other eigenvalue the direction of compression. Because of the incompressibility condition, both eigenvectors are orthogonal, and so plotting only one of them is sufficient to have the complete set of information. In figure 8, the stretching rate and its direction are presented in an (
$x,z$
) plane for the three cases studied. The stretching rate is maximal near the upper part, and at small distance above the cilium tips in all three cases. For the synchronized motion (case (a)), there is almost no stretching away from the cilia during the recovery phase, whereas for cases (b) and (c), corresponding respectively to symplectic and antipleptic motion, a weak stretching can always be observed (see the right-hand side of cases (b) and (c) in figure 8). A complex shape of the stretching rate is observed at the cilium tips during the stroke phase in antipleptic motion (see cilium 3 for instance), during the stroke phase in symplectic motion (see cilia 6 and 7) and during the stroke phase in synchronized motion (results not shown). From their orientations, one can expect an enhancement in the mixing in this region. As in Ding et al. (Reference Ding, Nawroth, McFall-Ngai and Kanso2014), the stretching direction is a nonlinear function of space. Except for zones where the shear rate is maximal, the
$\unicode[STIX]{x1D6FE}\boldsymbol{e}_{\unicode[STIX]{x1D6FE}}$
field is then oriented at
$45^{\circ }$
compared with the
$x$
-direction, and almost uniform. This is reminiscent of a linear shear profile
$d_{x}=cz$
, where
$c$
is a constant. The nonuniform aspect of the stretching orientation indicates the presence of ‘folding’ in the displacement field
$\boldsymbol{d}$
, which also plays a role in the mixing, as explained by Kelley & Ouellette (Reference Kelley and Ouellette2011). Ding et al. (Reference Ding, Nawroth, McFall-Ngai and Kanso2014) obtained that the stretching rate was maximal for the antipleptic case for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/4$
, which is exactly what is obtained here, even if the beating patterns used are different. One can then expect that this value of the phase lag
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
is the most efficient in mixing fluids by enhancing the stretching near the upper part of the cilia.
Figure 8. Stretching rate and direction for (a) synchronized motion (
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=0$
), (b) symplectic MCW (
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=-\unicode[STIX]{x03C0}/4$
) and (c) antipleptic MCW (
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/4$
). In each case, the cilium spacing is
$a/L=1.67$
and the size of the computational domain is (
$N_{x}=201$
,
$N_{y}=26$
,
$N_{z}=50$
).
In order to assess the reliability and robustness of the solver, a comparison with experimental data is also performed. The average clearance velocity was computed on the plane (
$x$
,
$y$
,
$3.2L$
) for all simulations. The highest clearance velocity is reached for the case
$a/L=1.67$
with
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/4$
and equals
$33.47~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$
. This is in good agreement with the experimental results of Matsui et al. (Reference Matsui, Randell, Peretti, Davis and Boucher1998), who observed a clearance velocity of
$39.8\pm 4.2~\unicode[STIX]{x03BC}\text{m}~\text{s}^{-1}$
for mucus. Matsui et al. (Reference Matsui, Randell, Peretti, Davis and Boucher1998) also observed that the PCL flows in the same direction as the mucus, which is the case in the present study. Finally, the PCL–mucus interface remains approximately flat in most of the simulations presented, as is the case in the micrographs performed by Sanderson & Sleigh (Reference Sanderson and Sleigh1981) on rabbit tracheal epithelium.
3.2.2 Directional pushing efficiency
Inspired by the works of Kim & Netz (Reference Kim and Netz2006), Gauger et al. (Reference Gauger, Downton and Stark2009) and Khaderi, Den-Toonder & Onck (Reference Khaderi, Den-Toonder and Onck2011), a positive flux
$Q_{p}$
and a negative flux
$Q_{n}$
are defined as follows. The
$x$
-component of the velocity field is considered over the (
$N_{x},y,z$
) plane and, at each time step of a full beating cycle, the negative and positive velocity values are separated. The instantaneous negative and positive fluxes are then computed over a sufficiently large number of cycles. The difference
$(Q_{p}-Q_{n})$
gives the net instantaneous flux, and the directional pushing efficiency
$\unicode[STIX]{x1D716}_{PN}$
is defined as
$\unicode[STIX]{x1D716}_{PN}=(Q_{p}-Q_{n})/(Q_{p}+Q_{n})$
. It should be noted that here, in contrast to the previous works of the aforementioned authors, the domain is not restricted to the transport area but covers the whole flow region instead. Such a choice is justified by the fact that it has been experimentally observed, using confocal microscopy and fluorescent markers, that the PCL and mucus are in reality transported at approximately the same rate (Matsui et al.
Reference Matsui, Randell, Peretti, Davis and Boucher1998). As a matter of fact, PCL transport seems to depend on the presence of a mucus layer above it, and ciliary mixing is thought to be responsible for the diffusion of momentum from mucus to PCL. Taking into account the transport zone and the mixing zone in the computation of
$\unicode[STIX]{x1D716}_{PN}$
then allows one to collect information about how transport and mixing work together during a beating cycle.
The variation of
$\unicode[STIX]{x1D716}_{PN}$
over one beating cycle for five different phase lags
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
, all for an antipleptic MCW with a cilium spacing of
$a/L=2$
, is illustrated in figure 9. Let us recall that, although the quantities considered here are the results of the motion of all cilia, the calculations have been made through the beating cycle of a single cilium of reference. Therefore, the fact that the corresponding curves have a huge decrease around
$t=\unicode[STIX]{x03C0}$
is related to the choice of this particular cilium. Indeed, the plane used to compute the positive and negative fluxes is located close to a particular cilium. This cilium is in the recovery phase at
$t=\unicode[STIX]{x03C0}$
, inducing a small value of
$\unicode[STIX]{x1D716}_{PN}$
at this particular time. The choice of another cilium or another plane would horizontally shift the drop observed around
$t=\unicode[STIX]{x03C0}$
.
To interpret the evolution of
$\unicode[STIX]{x1D716}_{PN}$
in figure 9, it is necessary to examine the topology of the flow during this beating cycle, shown in figure 10. As cilium 1 begins its recovery phase (case (a) of figure 10), there is no reversal flow at the periodic boundary frontier. Hence,
$\unicode[STIX]{x1D716}_{PN}=1$
, as is visible in figure 9 at the point (a). Then, when cilium 1 carries on its recovery phase (case (b) of figure 10), no reversal flow can be observed since cilium 1 is still at the beginning of its recovery phase and cilium 8 is finishing its stroke phase. At
$t=\unicode[STIX]{x03C0}/2$
(case (c) of figure 10), both cilia 1 and 8 are in their early recovery phase and the directional pushing efficiency begins to decrease (see points (b) and (c) on figure 9). At
$t=3\unicode[STIX]{x03C0}/4$
(case (d) of figure 10), cilium 1 is reaching the maximal speed of its recovery phase. The negative flow generated is important, and the efficiency quickly drops: see point (d) on figure 9. Between
$t=3\unicode[STIX]{x03C0}/4$
and
$t=\unicode[STIX]{x03C0}$
, the directional pushing efficiency weakly increases. This is due to the positive velocity near the base of cilium 1 which is performing a ‘whip-like’ motion: while its free end is still going to the left, its base starts moving to the right. Case (e) of figure 10 corresponds to the end of the recovery phase of the first cilium, and at the same moment, cilium 8 is at its maximal negative velocity: henceforth the directional pushing efficiency is still decreasing. At
$t=5\unicode[STIX]{x03C0}/4$
(case (f) of figure 10), cilium 8 has finished its recovery phase and cilium 1 is in the middle of its stroke phase. It counteracts the reverse flow created earlier, and
$\unicode[STIX]{x1D716}_{PN}$
finally increases. At
$t=3\unicode[STIX]{x03C0}/2$
and
$t=7\unicode[STIX]{x03C0}/4$
(cases (g) and (h) of figure 10), both cilia 1 and 8 are in their stroke phases, and the flow generated is purely in the positive
$x$
-direction: see points (g) and (h) on figure 9.
Figure 9. Directional pushing efficiency
$\unicode[STIX]{x1D716}_{PN}$
over one beating cycle for three different phase lags
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
. The results are obtained for an antipleptic MCW and
$a/L=2$
.
Figure 10. Snapshots of the fluid velocity taken at eight different instants during the beating cycle of the cilium located at the left of the images. The cilia are beating in an antipleptic motion with a phase lag
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/4$
and a spacing of
$a/L=2$
. In (a),
$t=0$
, and the recovery phase of the first cilium begins; (b–e) recovery phase; (f–h) stroke phase. Once (h) is completed, a new cycle begins: (h)
$\rightarrow$
(a). The figure shows contours of the norm of the normalized fluid velocity
$\Vert \boldsymbol{V}_{f}\Vert$
such that
$\boldsymbol{V}_{f}=\boldsymbol{V}_{f}^{dim}/V_{f}^{ref}$
, where
$V^{ref}=\unicode[STIX]{x1D706}/T$
is the reference speed of the present system.
It is now interesting to compare the average directional pushing efficiency
$\langle \unicode[STIX]{x1D716}_{PN}\rangle$
for the different kinds of synchronization over a full beating cycle (see figure 11
a,b). Before going further into detail, let us recall the fact that the current efficiency
$\langle \unicode[STIX]{x1D716}_{PN}\rangle$
is a criterion qualifying the non-isotropy of the transport. Thus, it only gives an insight into the capacity of the cilia to transport the flow in a unidirectional way and does not give any quantitative information on the volume of fluid flow that is effectively displaced. Figure 11(a,b) shows an interesting phenomenon: systems with larger cilium spacing have a better capacity to transport the flow in the same direction. Another intriguing fact is that, for the smallest cilium spacing (
$a/L=1.67$
, figure 11), antipleptic MCWs seem to have a better ability to create a unidirectional flow compared with synchronized or symplectic motion, with a peak for the average efficiency around
$\unicode[STIX]{x03C0}/2$
. This agrees particularly well with the results of Gauger et al. (Reference Gauger, Downton and Stark2009), who reported that antipleptic MCWs are more efficient than the synchronized motion of cilia, which itself is more efficient that the symplectic case. This also partially agrees with the results of Ding et al. (Reference Ding, Nawroth, McFall-Ngai and Kanso2014), who obtained peaks of efficiency for phase lags around
$\pm \unicode[STIX]{x03C0}/2$
, with a stronger one for the antipleptic case. In the present study, the minimal efficiency for metachronal motion is reached for a symplectic MCW with a phase lag of approximately
$-\unicode[STIX]{x03C0}/2$
.
When the cilium spacing is increased to
$a/L=2$
(figure 11), a peak of efficiency appears for the antipleptic MCW for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=2\unicode[STIX]{x03C0}/3$
, and the lowest values are reached in the symplectic case for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}\approx -\unicode[STIX]{x03C0}/2$
and
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}\approx -\unicode[STIX]{x03C0}/3$
. These trends will be found again in the results shown in § 3.2.3.
Then, if the cilium spacing is increased above a limit value (
$a/L=3.33$
, figure 11), the directional pushing efficiency becomes equal to 1 for all cases (antipleptic, symplectic and synchronized). This agrees well with the results obtained by Khaderi et al. (Reference Khaderi, Den-Toonder and Onck2011), who concluded that ‘the amount of flow enhancement depends on the inter-cilia spacing [but] the efficiency is not significantly influenced’. Indeed, as the cilia are set away from each other, the influence of their neighbours becomes negligible. Nevertheless, the authors reported a very low efficiency for synchronized cilia, whereas the present results always show a positive flow in the mucus phase for synchronized motion. As mentioned in § 2.1, this is a direct consequence of the inertial effects (
$Re>1$
) in the present simulations. It is recalled that they only affect the synchronized case. It has been carefully checked that all other conclusions drawn for antipleptic and symplectic metachrony remain the same for
$Re<1$
and
$Re=20$
. For a detailed study of the inertial effects at
$Re=20$
, see § 3.2.7. The present results are different from the results of Gueron et al. (Reference Gueron, Levit-Gurevich, Liron and Blum1997), as it appears that even when the cilium spacing is higher than 2 cilium lengths, the influence of neighbouring cilia cannot be neglected (see the case
$a/L=2.53$
, for example).
It is important to remember that this efficiency does not take into account the actual net flow volume transported.
Figure 11. Mean directional pushing efficiency
$\langle \unicode[STIX]{x1D716}_{PN}\rangle$
over a beating cycle as a function of the phase lag
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
for different values of the cilium spacing
$a/L$
.
3.2.3 Comparison of hydrodynamic efficiency
After having investigated the capacity to transport the flow in the desired direction, it is now interesting to look at the actual flow volume displaced to see whether a better efficiency to transport the flow directionally corresponds to a better flow transport. To do so, the global volumetric flow rate
$Q_{v}$
over a unit volume of size (
$1\times 1\times N_{z}$
) is introduced,
$$\begin{eqnarray}Q_{v}=N_{z}\frac{U^{\ast }\,\text{d}x^{2}}{L^{2}},\end{eqnarray}$$
where
$\text{d}x=1$
using the classical LBM normalization, and
$U^{\ast }=U^{av}/U^{ref}$
, where
$U^{ref}=(\unicode[STIX]{x1D706}/N_{cil})/T$
is the reference velocity and
$U^{av}=(1/n_{i}n_{j}n_{k})\sum _{i,j,k}U_{ijk}$
is the average fluid velocity inside the whole domain.
In figure 12, the dimensionless volumetric flow rate
$Q_{v}$
is plotted over one beating cycle for arrays of cilia with different phase lags (
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=0$
,
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=-\unicode[STIX]{x03C0}/4$
and
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/4$
) and different cilium spacings. For cilium spacings equal to
$a/L=1.67$
and
$a/L=2$
, the antipleptic MCW is the most efficient for transporting fluids. However, it should be noted that as the cilium spacing is increased, the ability of the antipleptic wave to transport fluids exhibits a huge decrease compared with the symplectic wave. For a larger cilium spacing (
$a/L=3.33$
), the fluid transport is not impacted any more by metachrony. Finally, no flow reversal occurs for the synchronized cases; this is the consequence of working at
$Re=20$
.
Figure 12. Normalized volumetric flow rate
$Q_{v}$
as a function of time over a beating cycle for different phase lags
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
and different cilium spacings
$a/L$
. The recovery phase occurs for
$t\in [0,\unicode[STIX]{x03C0}]$
and the stroke phase for
$t\in [\unicode[STIX]{x03C0},2\unicode[STIX]{x03C0}]$
.
The total volume of fluid displaced during a beating cycle for the different phase lags is compared in figure 13(a,b). For a small cilium spacing (
$a/L=1.67$
), the efficiency of the antipleptic metachrony is obvious. This agrees well with the results of Khaderi et al. (Reference Khaderi, Den-Toonder and Onck2011), who observed a larger net flow produced by antipleptic metachrony for this value of cilium spacing. Symplectic waves appear to be less efficient than or at best equally efficient to antipleptic motion, except for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=-7\unicode[STIX]{x03C0}/8$
for
$a/L=1.67$
, where there is a peak in the total displaced volume of flow. On the contrary, negative peaks are found for
$a/L=1.67$
and
$a/L=2$
for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=-\unicode[STIX]{x03C0}/3$
. There are two neighbouring maxima at
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/4$
and
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/2$
for
$a/L=1.67$
and
$a/L=2$
respectively, indicating that specific phase lags are more able to generate a strong flow. At this point, it is worth remembering that in the present model, the recovery and stroke phases occur in the same plane, which is not the case in real configurations. According to Downton & Stark (Reference Downton and Stark2009), who studied both 2D and 3D beating patterns, it is expected that in the presence of 3D beating patterns, both the directional pushing efficiency
$\unicode[STIX]{x1D716}_{PN}$
and the total displaced volume of flow would increase. It is also expected that a small flow might occur in the
$y$
-direction.
Figure 13. Total dimensionless displaced flow volume generated by an array of cilia over a beating cycle for different phase lags and cilia spacings:
$+$
,
$a/L=1.67$
; ○,
$a/L=2$
; *,
$a/L=2.53$
; ▫,
$a/L=3.33$
.
3.2.4 Energetic cost of the MCWs
The previous sections have been dedicated to the study of the directional transport and the volume displacements of fluids. An energetic perspective is now introduced, to characterize the potential benefits of metachrony completely.
The average power spent by the cilia during a beating cycle is given by
$$\begin{eqnarray}P_{cil}=\frac{\displaystyle \mathop{\sum }_{s,i}\boldsymbol{V}_{i}^{s}\boldsymbol{\cdot }(\boldsymbol{F}_{m}^{i}+\boldsymbol{F}_{PCL}^{i})}{N_{cilia}},\end{eqnarray}$$
using the forces illustrated in figure 2. The power spent is averaged over several beating cycles. To have a dimensionless power
$P^{\ast }$
, the power
$P^{\infty }$
spent by an isolated cilium (
$a/L=10$
) is computed such that
$P^{\ast }=P_{cil}/P^{\infty }$
.
Before going any further, it is important to remember here that the only parameter that differs between each value of the phase lag
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
for a given cilium spacing
$a/L$
is the size of the domain
$N_{x}$
over the
$x$
-direction, and hence the number of cilia acting in one wavelength. Moreover, between different cilium spacings
$a/L$
, the space
$a$
between two cilia is modified in both the
$x$
- and
$y$
-directions, and the sizes
$N_{x}$
and
$N_{y}$
of the domain must be changed accordingly. As a consequence, if the cilium spacing is increased, the density of cilia in both the
$x$
- and
$y$
-directions is decreased. In all of the simulations presented in § 3.2, all of the other parameters (ratio of viscosity
$r_{\unicode[STIX]{x1D708}}$
, beating frequency
$f$
, length of the cilia
$L$
, etc.) are fixed.
Figure 14(a,b) shows the average dimensionless power as a function of
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
. For small values of the phase lag (
$|\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}|\leqslant \unicode[STIX]{x03C0}/2$
), the average power spent decreases to a minimal value for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}\approx \unicode[STIX]{x03C0}/4$
. For this particular value, the system spends less power than the synchronized case: the antipleptic MCW allows the system composed of all cilia to encounter a smaller resistance from the surrounding fluids.
Nevertheless, it is obvious, from a fluid mechanics point of view, that when cilia beat together in a synchronized way, the viscous resistance felt by each cilium is reduced.
In the case of metachrony, the behaviour is dual. Indeed, cilia in the stroke phase of the antipleptic motion encounter a stronger viscous resistance from the flow due to their added respective remoteness. Therefore, compared with the synchronized motion, they will transfer more energy to the flow, and the energy transferred will be entirely used to propel the fluids in the desired direction. The opposite is true during the recovery stroke of the antipleptic motion: the clustered behaviour of the cilia allows them to experience a lower viscous resistance from the flow, and therefore to limit the amount of power transferred during this phase. However, globally, a system with an antipleptic MCW and
$|\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}|<\unicode[STIX]{x03C0}/2$
encounters less viscous resistance than the synchronized motion of cilia, as can be seen from figure 14(a,b). This results in a better efficiency of the cilia with antipleptic metachronal motion in transferring their momentum to the flow, while in the meantime requiring less power. For the symplectic motion, a similar phenomenon occurss: cilia in their stroke phase are clustered and therefore unable to fully exert their pushing action, while cilia in their recovery phase are away from each other and generate a stronger reversal flow compared with the antipleptic motion. This results in a lower capacity of the symplectic MCW to transport mucus, while still allowing the system to spend less power than synchronously beating cilia for
$|\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}|<\unicode[STIX]{x03C0}/2$
. An energetic ratio will be introduced in § 3.2.5 to quantify, for this value of
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
, the capacity of cilia with antipleptic or symplectic metachrony to transfer their momentum to the flow.
Figure 14. Power spent by the system for different phase lags
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
and different cilium spacings:
$+$
,
$a/L=1.67$
; ○,
$a/L=2$
; ▫,
$a/L=3.33$
.
This is not true for large phase lags (
$|\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}|\geqslant \unicode[STIX]{x03C0}/2$
): a system with antipleptic or symplectic metachrony spends more power than the synchronous case. One can suppose that this is a direct consequence of the reduced number of cilia in a wavelength. Then, to explain the better efficiency of the antipleptic MCW over the symplectic one for large phase lags, an investigation of the flow topology is necessary. Figure 15 shows an antipleptic MCW (on the left) and a symplectic MCW (on the right) at the same time for a phase lag of
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\pm 2\unicode[STIX]{x03C0}/3$
. For this value of
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
, the antipleptic wave is more efficient in transporting the mucus (see figure 13
a,b), although the average amounts of power spent by both systems are relatively similar (see figure 14
a,b). One can easily see the main difference between the two cases. In the symplectic case, the cilium in stroke phase encounters the reversal flow generated by the other cilium in recovery phase. When these two secondary flows meet, vortices are generated and the global transport of fluid is less efficient, as a fraction of the energy transferred to the flow is used to cancel this reversal flow. On the contrary, in the antipleptic case, the cilium in stroke phase does not immediately feel the influence of the reversal flow created by the other cilium in recovery phase, which is behind it. Hence, this cilium is able to fully exert its pushing effect on the mucus phase. Moreover, for the antipleptic case, there is a suction effect due to the combined motion of the cilia in the stroke and recovery phases, maximizing the propulsion of mucus. One can then expect a weak blowing effect between the cilia in the stroke and recovery phases of symplectic motion.
Figure 15. Comparison of the flow generated by an antipleptic and a symplectic wave for the cilium spacing
$a/L=1.67$
. (a) Antipleptic MCW with
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=2\unicode[STIX]{x03C0}/3$
. (b) Symplectic MCW with
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=-2\unicode[STIX]{x03C0}/3$
. The plane is coloured with the magnitude of the dimensionless fluid velocity.
When the cilium spacing is large (see the case
$a/L=3.33$
of figure 14
a,b, for example), the power spent by the system is the same for all phase lags: the cilia are too far from each other to be impacted by the phase difference of their immediate neighbours. One can observe that for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=0$
(i.e. for synchronized beating), the average power spent by a cilium,
$P^{\ast }$
, is equal to 1 in figure 14(a) for the cilium spacings
$a/L=1.67$
and
$a/L=3.33$
, meaning that synchronous systems spend the same amount of power as an isolated cilium. This shows the beneficial cost of antipleptic metachrony for small phase lags.
3.2.5 Displacement ratio
Inspired by the work of Kim & Netz (Reference Kim and Netz2006), a displacement ratio, which can be seen as the transport efficiency of the waves, is introduced to quantify the capacity of a given system to transport particles, with respect to a given amount of power. In this context,
$\unicode[STIX]{x1D702}_{1}$
is defined by the mean displacement in the
$x$
-direction during one beating cycle, divided by the mean power that a cilium had to spend during this beating cycle. Since the main purpose of mucociliary clearance is to transport mucus, and since experimental data (Winters & Yeates Reference Winters and Yeates1997) report that the total thickness in the vertical direction of the mucus layer is in the range
$[1.4L;10L]$
, values for the displacement were taken on an arbitrary plane
$z/L=3.2$
near the extremity of the domain. To obtain a value for the displacement, the instantaneous average fluid velocity over the
$x$
-direction is computed, and the resulting value is then multiplied by the period of a full beating cycle, giving the mean displacement
$\langle d_{x}\rangle$
over one period on the (
$x,y,3.2L$
) plane.
By dividing this mean displacement by appropriate quantities, a dimensionless expression for the displacement ratio is obtained,
$$\begin{eqnarray}\unicode[STIX]{x1D702}_{1}=\frac{\langle d_{x}\rangle \displaystyle \frac{N_{cil}}{\unicode[STIX]{x1D706}}}{P^{\ast }}.\end{eqnarray}$$
In the synchronized case, i.e.
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=0$
,
$\unicode[STIX]{x1D706}$
is infinite and thus the size of the domain over the
$x$
-direction was used and divided by the number of cilia.
In figure 16(a,b), one can see, from an energetic point of view, the superiority of the antipleptic wave in transporting mucus for small cilium spacings, which confirms the previous findings. Moreover, for the smallest cilium spacing (
$a/L=1.67$
), a clear peak of efficiency can be seen for an antipleptic MCW with
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/4$
. If the cilium spacing is increased (
$a/L=2$
), similar results are found. Nevertheless, for
$a/L=2.53$
, a different behaviour occurs: the displacement ratio is found to be worst for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\pm \unicode[STIX]{x03C0}/2$
. Then, for cilia far away from each other, the displacement ratio
$\unicode[STIX]{x1D702}_{1}$
remains constant for all phase lags (the cilia do not feel the influence of the others any more).
Figure 16. Displacement ratio
$\unicode[STIX]{x1D702}_{1}$
as a function of the phase lag
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
for different cilium spacings
$a/L$
.
Similar results regarding the efficiency of the antipleptic MCW are found in Osterman & Vilfan (Reference Osterman and Vilfan2011), where the optimal beating pattern is investigated. In their study, they observed that antipleptic MCWs were often the most efficient. They also found, as is the case in the present study, that an increase in the density of cilia results in an increase of the efficiency up to a critical point where clustering becomes counterproductive.
3.2.6 Mixing
To investigate how the mixing can be enhanced by metachrony, the average stretching rate over the transport and mixing areas during a beating period is computed. Figure 17 shows the results obtained for all cilium spacings. Clearly, the antipleptic MCW is the most efficient for stretching fluids, and hence for mixing fluids. Clear peaks are visible for antipleptic waves with
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/4$
and
$a/L=1.67$
, with
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}\approx \unicode[STIX]{x03C0}/3$
and
$a/L=2$
, and with
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}\approx \unicode[STIX]{x03C0}/3$
and
$a/L=2.53$
. On the contrary, symplectic MCWs are almost always less efficient for mixing than antipleptic MCWs, except for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=-7\unicode[STIX]{x03C0}/8$
and
$a/L=1.67$
, where there is an enhancement in the mixing. This peak is also present for the opposite phase lag
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=7\unicode[STIX]{x03C0}/8$
, where it reaches approximately the same value. This is perfectly coherent with all previous results, and partially with the ones of Ding et al. (Reference Ding, Nawroth, McFall-Ngai and Kanso2014), who observed, for a unique layer of fluid with a uniform viscosity and a cilium spacing of
$a/L=1.67$
, the existence of two peaks in the shear rate: a weak one for symplectic waves for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=-\unicode[STIX]{x03C0}/2$
and a stronger one for antipleptic waves for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/2$
. The main difference here is that, in the present study, the symplectic MCWs with
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}\approx -\unicode[STIX]{x03C0}/2$
seem to have the worst capacity for mixing, while the antipleptic MCWs reach their full mixing capacity for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/4$
. The combined study of figures 14 and 17 gives a good insight into the mixing efficiency of the system.
Figure 17. (a) Average stretching rate in the transport and mixing areas as a function of the phase lag
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
for different cilium spacings:
$+$
,
$a/L=1.67$
; ○,
$a/L=2$
; *,
$a/L=2.53$
; ▫,
$a/L=3.33$
. (b) Three-dimensional view of the corresponding plot.
3.2.7 Quantification of inertial effects
In this section, the influence of the Reynolds number is considered. The present objective is to compare the qualitative behaviour of the MCWs between
$Re<1$
and
$Re>1$
. To reach Reynolds numbers of the order of
$10^{-2}$
, the size of the cilia has been reduced to half the size of real cilia. The geometry is then divided by a factor of 2 in the three spatial directions. The main consequence of this choice is that the computed quantities (fluxes, displaced volume of fluid, etc.) displayed here are approximately eight times smaller compared with the results previously shown (as in figure 13). However, the qualitative behaviour of the MCWs remains the same, as well as the drag exerted by the cilia. Thus, a comparison of the MCW behaviour between
$Re<1$
and
$Re>1$
is possible. In Figure 18, one can see the total displaced volume of fluid for
$Re=0.02$
,
$Re=1$
and
$Re=20$
, and a cilium spacing of
$a/L=1.67$
. A transition between
$Re<1$
and
$Re>1$
for the synchronous cases is expected, and can be seen in figure 18. As mentioned in § 3.2.3, for
$Re\leqslant 1$
, the synchronized cases are less efficient than the symplectic cases. Except for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=0$
, no other notable quantitative differences are present between the cases
$Re<1$
and
$Re=20$
. The general behaviour of both the antipleptic and symplectic MCWs remains similar for all the other phase lags
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
. In particular, a peak around
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/4$
is present for all of the Reynolds numbers tested (
$Re=0.047$
, 2, 5, 10, not shown). The similarity of the metachronal cases (i.e.
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}\neq 0$
) for this range of Reynolds numbers is certainly due to the fact that there are always cilia in the stroke phase. Hence, transport is observed during the whole beat cycle, and the inertial effects thus remain minor. This is no longer the case for fully synchronously beating cilia: despite being weak, inertial effects cancel the reversal of the flow.
Figure 18. Total displaced volume of fluid for
$Re=0.02$
,
$Re=1$
and
$Re=20$
as a function of the phase lag
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}$
. The results of these simulations are obtained with
$L=7$
lu and
$a/L=1.67$
.
4 Conclusions
By using an LB solver coupled to an IB method, and considering a purely hydrodynamical feedback from the fluid, symplectic and antipleptic MCWs can emerge in a 3D two-phase flow configuration with a viscosity ratio of 20. It is known that MCWs may emerge due to hydrodynamic interactions (Elgeti & Gompper Reference Elgeti and Gompper2013). However, for the first time, both types of metachrony are observed to emerge in a two-layer fluid with different viscosities, using a simple feedback law. This feedback depends on a coupling parameter
$\unicode[STIX]{x1D6FC}$
that can be used to tune the strength and direction of the waves. It is observed that cilia experience weaker torques when they are in a clustered configuration, i.e. during the recovery phase of antipleptic motion, as well as during the stroke phase of symplectic motion. The resulting beating pattern is a slow stroke phase and a fast recovery phase for both cases, suggesting that even if a simple hydrodynamic interaction is sufficient to let MCWs emerge, it is not sufficient to reproduce all features of the beating patterns observed in real ciliated surfaces, and other biological issues may play a role in the transport mechanism. The study of Hussong, Breugem & Westerweel (Reference Hussong, Breugem and Westerweel2011) has shown that metachronal motion can switch from symplectic to antipleptic with increasing inertial effects. Here, the influence of cilium density is highlighted: by assuming a least-effort behaviour for the cilia, the present model shows that the metachrony can switch from antipleptic to symplectic by lowering the cilium density. A more detailed study of the quality of synchronization along the
$x$
- and
$y$
-directions is the next step towards a better understanding of the emergence of MCWs.
A thorough comparative study of the antipleptic and symplectic MCWs has been performed, and the results show that antipleptic MCWs are the most efficient ones for transporting mucus. This is in accordance with most recent studies addressing this point (Khaderi et al.
Reference Khaderi, Den-Toonder and Onck2011; Osterman & Vilfan Reference Osterman and Vilfan2011; Ding et al.
Reference Ding, Nawroth, McFall-Ngai and Kanso2014), and is now confirmed in a two-layer environment. In the range of cilium spacing studied, and especially for small phase lags (
$|\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}|<\unicode[STIX]{x03C0}/2$
), the antipleptic MCWs have a better ability to (i) transport the flow in the same direction compared with symplectic MCWs with the same wavelength, (ii) generate a higher flow rate, (iii) advect particles at a given power input and (iv) generate a higher stretching and, hence, are more able to mix fluids. On the contrary, symplectic MCWs do not appear to have a great impact on the flow and are often less efficient than antipleptic MCWs. This is not in agreement with Elgeti & Gompper (Reference Elgeti and Gompper2013), who reported that symplectic waves are almost as efficient as antipleptic waves using an optimized-efficiency model. This difference may be due to several factors, including the two-layer flow character of the present work, the ratios of viscosity considered or even the cilium beat shape used. However, these results are in accordance with Khaderi et al. (Reference Khaderi, Den-Toonder and Onck2011), who explained this better transport efficiency as being due to vortices obstructing the flow for symplectic motion.
Among the quantities introduced earlier, some are more decisive for the transport of particles. The most important ones are the net volume of mucus displaced, the transport efficiency
$\unicode[STIX]{x1D702}_{1}$
and the mixing capacity of the system. It has been shown throughout this work that antipleptic MCWs with
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/4$
generate the largest flow while in the meantime showing the highest value of the transport efficiency (the power spent,
$P^{\ast }$
, for
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/4$
being minimal). A preliminary study of the mixing also strongly indicates that they might be more suitable for mixing fluids. The better transport efficiency of the antipleptic MCW is certainly due to the clusterized aspect of the cilia in recovery phase which minimizes their impact on the flow, while cilia in stroke phase are able to fully exert their pushing action. While being more able to advect particles than any other form of coordinated motion, antipleptic MCWs with
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}<\unicode[STIX]{x03C0}/2$
also require less power. It should be noted that the results presented throughout this paper concern a viscosity ratio
$r_{\unicode[STIX]{x1D708}}$
of 20. Chatelin & Poncet (Reference Chatelin and Poncet2016) have shown that mucus transport is maximized for viscosity ratios ranging from 10 to 20. With the present model, and within this interval of viscosity ratios, the same trends are always observed (results not shown), and the best transport capacities are always obtained for antipleptic MCWs with
$\unicode[STIX]{x0394}\unicode[STIX]{x1D6F7}=\unicode[STIX]{x03C0}/4$
. The transport properties of the waves also remain similar and do not vary significantly. However, it is worth noticing that the best transport capacity was obtained for
$r_{\unicode[STIX]{x1D708}}=20$
. Hence, the optimum is probably outside the range of viscosity ratios studied. This also strongly indicates that an increase in the mucus viscosity results in an increase of the transport properties. As in previous studies (Norton, Robinson & Weinstein Reference Norton, Robinson and Weinstein2011; Chatelin & Poncet Reference Chatelin and Poncet2016), an optimum viscosity for the transport of particles can certainly be found. A more detailed study of how different mucus viscosities affect the transport properties of the MCWs is the next step towards a better understanding of mucociliary clearance.
While in the present study only one prescribed beating pattern was investigated, the results can certainly be generalized to other beating patterns. An interesting path of investigation would be to study how this beat shape is influenced by different mucus viscosities, mucus rheologies or thicknesses of the PCL. It should be noted that a Reynolds number of 20 was used for computational reasons. While this does not affect the behaviour of either the antipleptic or symplectic MCW, weak inertial effects are seen to change the dynamics of the flow in the particular case of synchronously beating cilia. Finally, the mucus was considered as being a Newtonian fluid in the present work, while it is in reality highly non-Newtonian. It is expected that the clearance velocity of both the antipleptic and symplectic MCWs would be increased by considering the mucus as being a viscoelastic fluid. One could also expect the symplectic waves to become more efficient than the antipleptic MCWs for highly viscous mucus. Indeed, the clusterized aspect of the cilia during the stroke phase of symplectic motion could help them to overcome the viscous resistance of the mucus more easily (Knight-Jones Reference Knight-Jones1954). Nevertheless, it should be noted that studies (Norton et al. Reference Norton, Robinson and Weinstein2011; Chatelin & Poncet Reference Chatelin and Poncet2016) all tend to show, however, that the mucus cannot be too ‘solid’ or too ‘liquid’.
The numerical results obtained here constitute a useful basis of investigation to progress in the understanding of respiratory diseases linked to cilium beating disorders, such as COPD or severe asthma (Chanez Reference Chanez2005). They are also of interest for industrial purposes, such as the design of cilium-based actuators for mixing (Chen et al. Reference Chen, Chen, Lin and Hu2013) or flow regulators in microscopic biosensors.
Future works include the implementation of a realistic rheological model (Lafforgue et al. Reference Lafforgue, Poncet, Seyssiecq-Guarente and Favier2016) for the mucus in order to take into account its highly non-Newtonian behaviour, and the study of the mass transfer occurring at the epithelial surface. A long-term objective of this work is to build a numerical environment to predict the transport and mixing of drugs inside both the mucus layer and the PCL. With this tool, the effect of drugs could be tested virtually before proceeding to clinical trials. For instance, drugs acting on microscopic parameters such as the beating frequency, the viscosity of the mucus and the density of cilia could be tested and their effects on macroscopic quantities such as the mucus flow rate and mixing could be explored, in order to progress in the understanding and treatment of respiratory diseases.
Acknowledgements
The authors would like to thank Dr Z. Li and master intern J. Mercat for the development of the numerical code. S.C. and S.P. also acknowledge the Natural Sciences and Engineering Research Council of Canada for its financial support through a Discovery grant (RGPIN-2015-06512). This work was granted access to the HPC resources of Compute Canada and Aix-Marseille University (project Equip@Meso ANR-10-EQPX-29-01).
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2017.352.
