Mucociliary clearance is the first defense mechanism of the respiratory tract against inhaled particles. This mechanism is based on the collective beating motion of cilia at the surface of epithelial cells. Impaired clearance, either caused by malfunctioning or absent cilia, or mucus defects, is a symptom of many respiratory diseases. Here, by exploiting the lattice Boltzmann particle dynamics technique, we develop a model to simulate the dynamics of multiciliated cells in a two-layer fluid. First, we tuned our model to reproduce the characteristic length- and time-scales of the cilia beating. We then check for the emergence of the metachronal wave as a consequence of hydrodynamic mediated correlations between beating cilia. Finally, we tune the viscosity of the top fluid layer to simulate the mucus flow upon cilia beating, and evaluate the pushing efficiency of a carpet of cilia. With this work, we build a realistic framework that can be used to explore several important physiological aspects of mucociliary clearance.

The process of mucus clearance is the first defense mechanism of the respiratory tract against inhaled particles, including viruses and bacteria. This mechanical shield relies on the collective motion of the cilia of ciliated cells lining the respiratory tract. Cilia are mostly immersed in a water-like environment, the periciliary fluid layer (PCL), whereas the overlying layer is made of a viscoelastic fluid, called mucus, whose role is to trap the inhaled particles (see Fig. 1).1,2

FIG. 1.

Schematic representation of the epithelium of the respiratory tract, where multiciliated cells are covered by cilia, whose motion allows the clearance of the mucus layer.

FIG. 1.

Schematic representation of the epithelium of the respiratory tract, where multiciliated cells are covered by cilia, whose motion allows the clearance of the mucus layer.

Close modal

Cilia in the respiratory tract are motile, hair-like filaments exhibiting a periodic asymmetric beating motion divided in two phases. During the power stroke (PS), the almost fully extended cilia move forward, penetrate slightly into the mucus layer, and push it forward. During the recovery stroke (RS), they unfold backward, without making contact with the mucus layer above. This cyclic pattern induces a net forward displacement of the mucus, and allows the clearance of the airway.1 Measurements obtained from high-speed video of differential interference contrast (DIC) microscopy show the frequency of this cycle to be 10–20 Hz.3–6 Moreover, scanning electron microscopy images showed that the motion of a cilium has a clockwise, out-of-plane component during the recovery stroke.3,7

The collective motion of respiratory cilia is highly coordinated, as observed on images and videos of DIC microscopy.3,7 The phase of the beat behaves like a progressive wave: on lines parallel to the wave propagation, successive cilia are phase-shifted, while cilia aligned perpendicularly are in phase and beat synchronously. This peculiar coordination is called the metachronal wave (MCW). It emerges from fluid mediated interactions between neighboring cilia and increases the efficiency of mucus transport.8,9 The wavelength for this wave has been estimated using different experimental approaches, with values ranging from 6 µm obtained from video transmission microscopy of human sphenoid sinus mucosa,10 up to 20 µm found by DIC imaging of human bronchial mucosa,7 and 50 µm using video reflection microscopy on several mammals.11 Finally, the velocity of tracheal mucus transport was measured by scintigraphy and found to be 5–10 mm/min for healthy subjects.12,13

Cilia malfunction leads to the impairment of mucociliary clearance—a hallmark of several respiratory diseases. Primary ciliary dyskinesia (PCD, also known as immotile cilia syndrome) is a congenital disease causing defects in the cilium ultrastructure, which, in turn, may lead to an abnormal motility or complete immotility of cilia, an absence of metachronal coordination, or even loss of cilia.1 Bacterial and viral infections may also cause cilia dismotility or immotility, and the presence of deciliated areas.14,15 Patients with cystic fibrosis, asthma, or chronic obstructive pulmonary disease present with reduced ciliary beating frequencies and suffer from dehydration of the airway surface fluid and hypersecretion of mucus, which results in an increased viscoelasticity of the mucus layer.16,17 Generally, the dysfunction of the clearance mechanism is associated with a higher risk of respiratory infection.

The detailed study of the complex mucociliary clearance mechanism is still a challenge from the experimental point of view, and powerful support might come from computational modeling. A wide range of models have been implemented and deployed to investigate, with a varying degree of resolution, the different aspects of the problem, from the cilia beating motion, to the coordination mechanism in cilia carpets and the associated fluid displacement.18–20 

Different mechanical models describing cilia movement are used to investigate the effect on the surrounding fluid using different techniques: slender body theory,21 lattice Boltzmann,8,22 and multiparticle collision dynamics.9 Also, continuous representation of the cilia carpets, and their spatial distribution, can be used as input to solve the surrounding fluid dynamics.23–25 

In this work, we use the computational technique of lattice Boltzmann molecular dynamics (LBMD) to investigate the collective behavior of cilia motion and the associated mucus displacement. The study is based on the original mechanical model introduced by Elgeti and Gompper in Ref. 9. In this model, cilia are discrete and independent entities, whose motion is controlled to reproduce the two phases of the beating dynamics. While originally the model was coupled to multi-particle collision dynamics26–28 to account for hydrodynamic effects, here, we implemented the coupling with the lattice Boltzmann technique. Our major effort was to tune the model in order to reproduce the physiological, characteristic time- and length-scales of cilia motility. Moreover, we explored the collective behavior for cilia carpets with a surface density close to the physiological one (ρ = 3.05 cilia/μm223). Finally, we accounted for the different viscosities of the mucus and PCL layers, and estimated the clearance efficiency. The final result is a realistic 3D framework, ready for the investigation at pertinent scales of the different aspects of the mucociliary clearance mechanism, from particle infections to cilia immotility.

In order to simulate the cilia dynamics, we implemented and adapted the phenomenological model introduced by Elgeti and Gompper in Ref. 9. In this model, a cilium is a crane-like structure, made of a network of springs, as represented in Fig. 2. The structure can be viewed as three interconnected vertical filaments—each formed by Nf beads, and arranged in a prism with a triangular base. In our implementation, Nf is set to 20. The bonds between the beads are described by harmonic potentials, with the constant K taking a different value depending on the bond type. Vertical bonds joining successive beads in a given filament and horizontal bonds linking beads belonging to the same cross section have the same equilibrium length lb and harmonic constant K1. The cross bonds making up the diagonals of the square cells on the vertical faces of the prism have a preferred length lc=lb2 and harmonic constant K2. The harmonic constants are set to keep the bond length (and, therefore, the cilium length Lc) mostly constant. The inter-filament bonds help avoid the twisting of the structure around the vertical axis. Finally, the bending rigidity of the cilium is modeled using an angular constraint Hα=12Kα(180°αi)2, with α the angle defined by AiAi+1Ai+2̂, formed by three successive beads on a given filament.

FIG. 2.

Representation of the cilium structure. Panel (a) represents the cilium’s triangular cross section and the elementary inter-filament lateral network, along with the characteristic bond lengths. Panel (b) represents the cilium during the power stroke. The structure is mostly straight above the pivot point, shown by the converging arrows, where it is bent forward. Panel (c) represents the cilium during the recovery stroke. The structure is bent backward at the base (see diverging arrows), and the point of maximal curvature shown by the converging arrows propagates to the tip, thus causing the cilium to unfold progressively. In panels (b) and (c), red bonds correspond to the active filament that controls the in-plane motion of the cilium. The green bonds, which extend during the backward stroke, ensure the out-of-plane motion. Arrows indicate the contraction or extension of the bonds during the different phases of the motion.

FIG. 2.

Representation of the cilium structure. Panel (a) represents the cilium’s triangular cross section and the elementary inter-filament lateral network, along with the characteristic bond lengths. Panel (b) represents the cilium during the power stroke. The structure is mostly straight above the pivot point, shown by the converging arrows, where it is bent forward. Panel (c) represents the cilium during the recovery stroke. The structure is bent backward at the base (see diverging arrows), and the point of maximal curvature shown by the converging arrows propagates to the tip, thus causing the cilium to unfold progressively. In panels (b) and (c), red bonds correspond to the active filament that controls the in-plane motion of the cilium. The green bonds, which extend during the backward stroke, ensure the out-of-plane motion. Arrows indicate the contraction or extension of the bonds during the different phases of the motion.

Close modal

To generate the active motion of the structure, the bonds of one of the three vertical filaments, which we refer to as the active filament (see red filament in Fig. 2), are extended or contracted in a time-dependent manner so as to obtain the cyclic motion typical of a cilium. This mechanism, devised originally in Ref. 9, aims to mimic the molecular action of the dynein motor proteins, responsible for the bending of the cilia. The resulting motion is characterized by two phases: the forward (power stroke) and the backward (recovery stroke) dynamics. The equilibrium length of each bond in the active filament is modified according to Eq. (1) [Eq. (2)] during the forward (backward) stroke:

lb(i)=lb11.2(ii0)2+1,
(1)
lb(i)=lb1+A(Nf1i)/Nfβfori<i01,lb11.2(ii0)2+1forii01,
(2)

where A = 0.3, β = 2.5 were originally set to reproduce the motion of the rabbit’s tracheal cilia.9 The index i refers to the position of the bond along the active filament (i = 1 being the first bond counting from the base). In the forward stroke, see Eq. (1), the pivot point i0 is fixed at 5. During the backward stroke, i0 starts at 7 and moves up the active filament at a constant velocity vi0, thus allowing the point of maximum curvature to propagate up to the tip of the cilium. Vertical bonds of variable preferred lengths are also subjected to a stalling mechanism. It is inspired by the stall force of dynein motors corresponding to the maximum external load that they can move against before stopping. This means that max(|lilb(i)|) = dlmax.

Additionally, it has been observed experimentally that while the forward stroke of a cilium is planar, an out-of-plane component is present during the backward stroke.3 This is accounted for by extending the bottom bonds (1 to 7) of a second filament (green filament in Fig. 2) from the start of the recovery stroke until i0 reaches the value of 12; afterward, we restore the equilibrium length lb. This allows the cilium to slowly go back in its plane before the next forward stroke starts. In order to prevent unwanted twisting of the body of the cilium, possibly exacerbated by the interactions with the surrounding fluid, the three bottom beads of each filament were kept frozen.

Finally, the switch between the forward and backward strokes is controlled by a shape-dependent condition.9 We compare the length of the active filament (L1) with the length of the two remaining ones (L2 and L3) by

dL=2L1L2L3.
(3)

The motion switches from the forward to the backward stroke when dL < dLmin, and vice versa when dL > dLmax. The time evolution of dL and the switching mechanism are shown in Fig. S1 in the supplementary material.

The described model was tuned in order to reproduce the main characteristic length- and time-scales of a physiological beating motion. First, we assign a mass mc to each bead of the structure. We assume that the real cilium has a mass density that matches that of its protein constituents, 1.35 g/cm3.29 We then estimated the volume of the cilium by considering the experimental length Lc = 6 µm and cross section (with radius r = 0.15 µm). We, therefore, obtained the volume of each bead in the model. We have adjusted the bead mass to ensure numerical stability and correct time scales, and obtained mc = 5 × 10−18 kg for each bead. For the bonds, we chose a preferred length lb = 0.3 µm, to achieve a final cilium length of Lc = 6 µm and width of 0.3 µm. The values of the stiffness constants satisfy the following conditions K2 = K1/10 and Kα = K1/100.9 The value of K1 has been tuned to set the desired timescale for the beating cycle: K1 = 9500 × 10−2 (kcal/mol)/μm2. The propagation velocity of the pivot point i0 during the backward strokes [cf. Eq. (2)], vi0, is fixed at 3 ms−1, meaning that the value of i0 is kept constant for 3 ms before increasing by 1. This switching kinetics influences the duration of the beating cycle of the cilium as well as the duration ratio between forward and backward strokes. It also ensures that the cilium has enough time to unfold for each value of i0 before moving on to the next. Finally, the threshold values were empirically fixed at dLmin = −0.48 µm and dLmax = 0.5 µm to obtain the desired beating shape, and the stalling threshold is dlmax = 0.03lb. The specific aspects of the beating motion will be discussed later on in the results section.

Using this model, we built square carpets of independent cilia. First, we simulated systems with different surface densities, by considering patches of 10 × 10, 14 × 14, and 20 × 20 cilia, see Fig. 3, top panel. For these systems, the surface densities were 0.63–2.5 cilia/μm2, with the highest value approaching the experimentally reported density of 3 cilia/μm2.23 Additionally, for the lowest density, we extended the size of the patch to check for finite size effects (see Fig. 3, bottom panels). Experimentally, a multi-ciliated cell (MCC) has characteristic size of 8μm; therefore, our systems correspond roughly to 1–4 MCC.

FIG. 3.

Schematic representation of the different systems simulated in this work. In the top layer, we represent the grid where the carpet cilia are anchored for different surface densities. In the bottom layer, for a given density, we represent systems of different sizes.

FIG. 3.

Schematic representation of the different systems simulated in this work. In the top layer, we represent the grid where the carpet cilia are anchored for different surface densities. In the bottom layer, for a given density, we represent systems of different sizes.

Close modal
FIG. 4.

Schematic view of the embedded particle/fluid systems. Particles live in the continuous space, while the LB fluid lives in the lattice space. The coupling is performed via an interpolation/extrapolation procedure of the fluid and particle velocities.

FIG. 4.

Schematic view of the embedded particle/fluid systems. Particles live in the continuous space, while the LB fluid lives in the lattice space. The coupling is performed via an interpolation/extrapolation procedure of the fluid and particle velocities.

Close modal

Each simulation is started with perfectly synchronized and aligned cilia, identically oriented, to have their power stroke parallel to the y direction. The chosen time-step for the molecular dynamics was Δt = 20 ns (a value that ensured numerical stability), and the simulations were extended to sample 10–30 beating cycles (depending on the density of the patch) and reach a steady state. The duration of simulations varied between 600ms and 2000ms, depending on the system.

In our setup, we used periodic boundary conditions for the xy directions. The simulation box in the z direction was set to 12.6 µm in all cases. An overview of the simulated systems is given in Table I.

TABLE I.

List of the systems simulated in this work. The box size L refers to the x and y directions. The cilia length is Lc = 6 µm.

L (μm)Nciliadc/Lcρcilia (μm−2)
12.6 100 0.21 0.63 
12.6 196 0.15 1.23 
12.6 400 0.11 2.52 
16.8 169 0.22 0.60 
21.7 289 0.21 0.61 
25.2 400 0.21 0.63 
L (μm)Nciliadc/Lcρcilia (μm−2)
12.6 100 0.21 0.63 
12.6 196 0.15 1.23 
12.6 400 0.11 2.52 
16.8 169 0.22 0.60 
21.7 289 0.21 0.61 
25.2 400 0.21 0.63 

The cilia dynamics was coupled to the dynamics of a surrounding fluid simulated by the lattice Boltzmann (LB) technique.30,31 The coupling between LB and MD introduces long-range hydrodynamic interactions into the MD simulation of the cilia, thus allowing for fluid mediated correlations (Fig. 4). The motion of the cilia is described by a Langevin-like equation of the form

miai=FiC+FiD+FiR,
(4)

where FiC=Ui(ri) corresponds to the conservative forces derived from the potential energy, and FiR is a random white noise accounting for thermal fluctuations. The term FiD=γ(viui) is a Stokes-like drag force, with vi the velocity of bead i, ui the velocity of the fluid at the position of bead i, and γ a friction coefficient that acts as a coupling parameter between the fluid dynamics and Molecular Dynamics (MD). We have worked with γ = 0.1 ns−1, chosen for numerical stability. For the water–mucus composite system, the fluid kinematic viscosity was set to the value for bulk water at ambient condition. In our simulation, the LB implementation uses the Bhatnagar–Gross–Krook (BGK) collision operator.32 The fluid and particles dynamics are evolved synchronously. For the sake of comparison, in our simulations, we used different values for the grid spacing Δx, in the range 0.4–0.7 µm, a size comparable to the cilium cross section, and the inter-cilia distance dc. The effect of the particle motion on the fluid is accounted for in the LB equation so that the coupling between fluid and cilia is a two-way exchange of momentum.33 

The simulations were performed using the code Muphy,34 already exploited for the simulations of a great variety of biological systems.33,35–39 The cilia dynamics was implemented as an independent Python snippet, optimized by Numba and Mpi for Python, and coupled to the main engine controlling the fluid and particle dynamics.

In this section, we discuss the tuning of some simulation parameters, as well as the effect of system surface density and cell size, on the beating dynamics.

We first tested the implementation by simulating an individual cilium embedded in a single phase water-like fluid and varying the lattice spacing Δx. The characteristic forward and backward motions are only mildly affected by the choice of Δx, which controls the coupling with the fluid, see Fig. 5. For instance, we note that when using the lowest resolutions of the grid (Δx = 0.6 and 0.7 µm), the direction of the power stroke slightly deviates from the y direction. The measured beating period slightly increases with Δx and is in the range 43–48 ms. The time evolution of the tip projection for different Δx is also shown in Fig. S2 in the supplementary material.

FIG. 5.

Time evolution of the projection of the tip of the cilium during a beating cycle and for different values of the grid spacing Δx. Panels (a) and (b) show the projection in the yz and xy planes, respectively. Panel (c) shows a characteristic 3D beating. The cilium is colored blue during the power stroke and cyan during the recovery stroke.

FIG. 5.

Time evolution of the projection of the tip of the cilium during a beating cycle and for different values of the grid spacing Δx. Panels (a) and (b) show the projection in the yz and xy planes, respectively. Panel (c) shows a characteristic 3D beating. The cilium is colored blue during the power stroke and cyan during the recovery stroke.

Close modal

We next consider a carpet of cilia and test the effect of the lattice spacing for different surface densities ρs. To each ρs value corresponds an inter-cilia spacing dc, see Table I. In Fig. 6, we plot the value of the beating period τb, obtained as an average over all the cilia in the carpet, as a function of the system inter-cilia spacing dc, normalized by the cilium length Lc = 6 µm. In the inset plot, we report the value of the beating period obtained while simulating a system at dc = 0.63 µm for different lattice spacings Δx. First, we note that the average beating period in the cilia carpet is basically independent of the surface density and increases only weakly with the grid spacing size. Averaging all simulations, we obtain a mean beating period of τb = 44.6 ± 2.1 ms and a mean frequency of fb = 22.6 ± 1.0 Hz for the motion of individual cilia. This value corresponds to the upper boundary of the usual range of frequencies given for respiratory cilia in humans, going from 10 to 20 Hz.1,4–6

FIG. 6.

Beating period as a function of the inter-cilia distance dc, normalized by the cilium length Lc. In the inset plot, for the system with dc = 0.63 µm, we report the beating period as a function of the lattice spacing Δx.

FIG. 6.

Beating period as a function of the inter-cilia distance dc, normalized by the cilium length Lc. In the inset plot, for the system with dc = 0.63 µm, we report the beating period as a function of the lattice spacing Δx.

Close modal

The recovery stroke takes about τrec = 38.9 ± 0.9 ms to complete before switching to the power stroke, which lasts τpow = 5.7 ± 2.0 ms. Thus, the recovery stroke takes up 87% of the full beating cycle, while the power stroke takes the remaining 13%, and their ratio is τpow/τrac = 0.15. When compared to data reported in the literature, this ratio is 2.5 times smaller than what is measured experimentally 0.4.5,7 This is mostly due to the unfolding mechanism used for the recovery stroke: It cannot be sped up since the active springs need to relax at each step to obtain the desired backward stroke shape. As a result, as discussed in Sec. II A, we decided to speed up the power stroke phase, mainly via tuning K1, to achieve the total beating period that matches the experimental value.

We estimated the pushing velocity of the cilium by computing the average velocity of its tip during the power stroke. We find ⟨vtip⟩ = 0.62 ± 0.02 mm s−1, with a maximum of 1.75–2 mm s−1, reached during the second part of the power stroke, when the cilium has passed its fully extended position. This velocity compares well with the value 0.76 mm s−1 reported for a beat frequency of 20 Hz.3 

The resulting Reynolds number is

Re=ρLcvtipη=4×103,
(5)

where ρ and η are the density and the dynamic viscosity of the fluid, respectively. This is the expected order of magnitude for cilia as reported in the literature.40 

During a beating cycle, a cilium reaches its maximum height in the middle of the power stroke, when it is almost perpendicular to the epithelial surface and fully extended. During the recovery stroke, it folds forward at the tip, so as to move back to the starting position of the power stroke without making contact with the mucus layer above. This asymmetry is an essential feature of the beating motion, which allows cilia to induce a net forward motion in the mucus layer. In our simulation [see Fig. 5(a)], the difference in height between the end of the recovery stroke and the highest point is between 0.5 and 1 µm, which makes it possible to achieve a penetration of 0.5 µm into the mucus layer without dragging it backward during the recovery stroke.3 

It has been observed experimentally that ensembles of cilia arranged in carpets and pushing on fluid exhibit a form of coordination of their beating cycles called the metachronal wave (MCW).3,41 In these conditions, cilia on lines perpendicular to the direction of the wave beat synchronously, while cilia on lines parallel to the direction of the wave exhibit a phase shift. It has been shown that hydrodynamic interactions are the key ingredient, allowing cilia to form the MCW.41 

In the following, we show that our simulated cells exhibit such a behavior. To quantitatively measure the loss of synchronicity and characterize the emerging MCW, we computed the observable B(r,t) corresponding to the projection of the base-to-tip vector (linking the bottom and the top beads of the active filament) in the direction of the power stroke.9 In Fig. 7(a), we report B(r,t) computed for a line of cilia as a function of the position in the y direction and as a function of time. We observe the progressive loss of synchronicity between the cilia as time progresses. In Fig. 7(b), we illustrate the propagation of the wave by showing two snapshots of a given cell at an interval of 10 ms. In the left part of the panel, we show two maps of the variable B(r,t) corresponding to the two snapshots. In these maps, we can visualize how wavefronts (i.e., lines of constant B) propagate in space. From this representation, it is possible to graphically estimate the wavelength; see the white segment in the same plot.

FIG. 7.

Simulation of a multi-ciliated cell containing 20 × 20 cilia in a box of dimensions Lx = Ly = Lz = 12.6 µm, with density ρs = 2.52 µm−2. Panel (a): Time evolution of B for a line of cilia. Panel (b): 2D map of the projection B(r,t) at two separate times in the same wave period (left). Snapshots of the system associated with the represented projection maps, colored according to the height of the beads (white for the lower beads and blue for the higher ones).

FIG. 7.

Simulation of a multi-ciliated cell containing 20 × 20 cilia in a box of dimensions Lx = Ly = Lz = 12.6 µm, with density ρs = 2.52 µm−2. Panel (a): Time evolution of B for a line of cilia. Panel (b): 2D map of the projection B(r,t) at two separate times in the same wave period (left). Snapshots of the system associated with the represented projection maps, colored according to the height of the beads (white for the lower beads and blue for the higher ones).

Close modal

We measure the normalized correlations in space and time of the variable B:

C(r,τ)=δB(r0,t)δB(r0+r,t+τ)r0,tδB(r0,t)2r0,t,
(6)

with δB(r0,t)=B(r0,t)B(r0,t)r0,t, where r0 is the position of a cilium and t the time. The average .r0,t is taken over all cilia positions r0 and all frames t. Figure 8 shows an example for a cell of 20 × 20 cilia (dc/Lc = 0.105). The correlation C(r,τ) has the structure of a progressive wave, with alternating highly correlated and highly anti-correlated regions. The structure is caused by the progressive phase shift between successive cilia in the wave propagation direction. The regions of constant C(r,τ) correspond to regions where cilia are in phase. They are thus called lines of synchrony and extend along lines perpendicular to the direction of propagation of the wave. We can extract characteristic parameters of the wave by fitting C(r,τ) with an oscillatory function of the form

C(r,τ)=Acos(ωτ+kr+ϕ),
(7)

with ω the pulsation of the wave, and k the wave vector such that the wavelength is λ = 2π/k. The details of the fitting procedure are given in the supplementary material, see also Figs. S3 and S4. Note that due to the small size of our system, we are not able to extract any relevant correlation length to characterize the overall decay. The period Twave = 2π/ω extracted from the fit is reported in Fig. 9(a) for different systems. It is mostly unaffected by the system density and is about 44.4 ms, which corresponds to a frequency around f = 22.5 Hz. These values match the period and frequency of the beating motion of individual cilia, as expected. From Fig. 9(b), we note that the wavelength seems equally independent of the density of cilia. We find an average wavelength of 10.9 µm, with a range of values going from 8.7 to 13.4 µm, where the highest value corresponds to a plateau reached when increasing the size of the box by a factor of 1.3 to 2 (see inset). This finding agrees with the experimental estimates for human bronchial cilia,7,λ ≃ 20.8 µm and f ≃ 15.6 Hz. Our values also agree with earlier calculations on a carpet of similar density and using the same model.9 However, note that in a recent study that compares characteristics of tracheal cilia in different mammals, a much higher value is reported, λ ≃ 50 µm.11 

FIG. 8.

Snapshot of the spatio–temporal correlations of the projection B in space, taken for τ = 50 ms, for a cell of 20 × 20 cilia (with ρs = 2.52 µm−2).

FIG. 8.

Snapshot of the spatio–temporal correlations of the projection B in space, taken for τ = 50 ms, for a cell of 20 × 20 cilia (with ρs = 2.52 µm−2).

Close modal
FIG. 9.

Period [Panel (a)] and wavelength of the metachronal wave [Panel (b)] as a function of the systems’ inter-cilia distance dc/Lc, and for systems of different sizes (insets).

FIG. 9.

Period [Panel (a)] and wavelength of the metachronal wave [Panel (b)] as a function of the systems’ inter-cilia distance dc/Lc, and for systems of different sizes (insets).

Close modal

The data we discussed were obtained from simulations with lattice spacing Δx = 0.7 µm. Data obtained for other lattice spacings are reported in Table II. We have also verified that, for a simulation that does not include the fluid dynamics, there is no collective coordination, and the MCW simply does not take place, see Fig. S5.

TABLE II.

Wave parameters for different values of the grid spacing Δx, obtained from the simulations of the system with box size L = 12.5 µm and ρs = 2.52 cilia/μm2.

Δx (μm)λ (μm)τb (ms)f (Hz)ω (s−1)
0.4 8.7 45.2 22.1 147.7 
0.5 8.8 43.3 23.1 145.2 
0.6 11.9 43.4 23.1 144.8 
0.7 9.9 44.0 22.7 142.6 
Δx (μm)λ (μm)τb (ms)f (Hz)ω (s−1)
0.4 8.7 45.2 22.1 147.7 
0.5 8.8 43.3 23.1 145.2 
0.6 11.9 43.4 23.1 144.8 
0.7 9.9 44.0 22.7 142.6 

We now investigate the pumping effect of the cilia. Once the steady state is reached, we recorded the fluid velocity fields spanning 1 to 2 beating cycles, and computed the fluid velocity profiles by averaging in time and in the xy-plane, see Fig. 10. The obtained asymmetric profile is reminiscent of a Poiseuille-like flow that would be generated by replacing the cilia action by a constant force density applied on a layer of height Lc, see Fig. S6 in Ref. 9.

FIG. 10.

Average fluid velocity profiles extracted from simulations of multi-ciliated cells of different surface densities ρs. The black dashed line represents a Poiseuille-like, parabolic velocity profile, resulting from a body force exerted in the layer [0, Lc], represented by the black arrows.

FIG. 10.

Average fluid velocity profiles extracted from simulations of multi-ciliated cells of different surface densities ρs. The black dashed line represents a Poiseuille-like, parabolic velocity profile, resulting from a body force exerted in the layer [0, Lc], represented by the black arrows.

Close modal

It is worth noting that when the density of the carpet is increased, the fluid velocity profile changes differently depending on the distance from the cilia, see the different curves in Fig. 10. In the upper layer (where mucus would be localized in physiological conditions), the velocity is higher for higher densities, showing that the higher the density, the stronger the pushing action is. On the contrary, in the bottom part (0–6 µm, corresponding to the PCL), the higher packing and crowding induces a screening effect, which reduces the fluid velocity.42 

In Fig. 11, we represent the action of the cilia beating on the fluid by showing the evolution of the velocity streamlines in the mucus region. We observe the heterogeneous evolution of the power stroke phase in the system in a time interval of 6 ms. The fluid is progressively pushed forward by the collective motion of the cilia.

FIG. 11.

Successive snapshots of the fluid streamlines in the region, corresponding to the mucus layer in physiological conditions, for a simulation of a multi-ciliated cell of size 14 × 14 in a box of size L = 12.5 µm, corresponding to a cilia density of ρs = 1.23 cilia/μm2.

FIG. 11.

Successive snapshots of the fluid streamlines in the region, corresponding to the mucus layer in physiological conditions, for a simulation of a multi-ciliated cell of size 14 × 14 in a box of size L = 12.5 µm, corresponding to a cilia density of ρs = 1.23 cilia/μm2.

Close modal

In this section, we introduce an effective model to describe the mucus layer in our simulations. The mucus is a viscoelastic gel, whose dynamic viscosity is several orders of magnitude higher than that of water. It ranges from 10−2 to 103 Pa s, whereas that of water, and, therefore, that of the PCL, is about 10−3 Pa s.16 

Our goal is to create an environment of higher viscosity in the top part of the simulation box in order to mimic the effect of mucus on the cilia dynamics and monitor the transport efficiency. From the lattice Boltzmann perspective, we have chosen a single phase fluid. We, therefore, alter the local viscosity by introducing short peptides (by analogy, with the presence of mucin glycoproteins in physiological mucus) that would alter the fluid response to cilia motion, see Fig. 12(a). An alternative approach would have been to define two separate fluid phases with different viscosities and particle coupling.43,44 Since our goal is to introduce an effect in the fluid, and not an explicit molecular description of the mucus, we did not introduce any inter-particle interactions. Therefore, all inter-particle interactions are mediated by the coupling with the fluid. To adjust the effect of the added peptides on the environment viscosity, we have first tuned some parameters (mass, peptide length, and density). This is done by considering how these parameters influence the peptide diffusivity. We built systems of dimension Lx = Ly = 12.5 µm and Lz = 6.25 µm, where we inserted about 360 peptide particles. These particles were free or interlinked via harmonic springs (K1) to form short peptides of 4 or 8 beads.

FIG. 12.

Panel (a): Snapshot of a simulation of a multi-ciliated cell with an explicit model of the mucus layer. Panel (b): Diffusion coefficient of the peptides in the mucus layer as a function of the mass and length of the peptide.

FIG. 12.

Panel (a): Snapshot of a simulation of a multi-ciliated cell with an explicit model of the mucus layer. Panel (b): Diffusion coefficient of the peptides in the mucus layer as a function of the mass and length of the peptide.

Close modal

An overview of the results is presented in Fig. 12(b). We obtain the smallest diffusion coefficient Dmin ≈ 3.5 × 10−10 μm2/ns using the combination of 40 peptides of 8 beads, with each bead of mass 8mc. This diffusivity is one order of magnitude smaller than the diffusion coefficient obtained for the center of mass of an isolated peptide in solution. Therefore, assuming a simple Stokes–Einstein relationship D ∝ 1/ν, we estimate that the presence of the peptides increases the viscosity of the environment by one order of magnitude, therefore approaching physiological conditions.

We performed a simulation of a ciliated cell with the explicit presence of the mucus. We considered three systems at different surface densities in a box of dimensions Lx = Ly = Lz = 12.5 µm. We measured the displacement of the mucus particles under the action of the beating cilia. In the top layer of Fig. 13, we report the velocity of the mucus peptides (calculated for the center of mass of each) as a function of time, distinguishing the contribution from peptides lying at different positions along the z axis. First, we note a cyclic behavior, with alternating high- and low-displacement times. The peptides closer to the cilia carpet (blue colors) are the most affected by the beating cycle and respond the most to the power stroke of the cilia. Interestingly, the clear cyclic pattern is blurred when the surface density increases.

FIG. 13.

Top: Velocity of the center of mass of each mucus peptide as a function of time and colored according to the position above the cilia, from blue for those close to the cilia, to red for those closer to the top of the box. Bottom: Average velocity of the mucus peptides (blue) and fraction of cilia in the power stroke (orange) as a function of time.

FIG. 13.

Top: Velocity of the center of mass of each mucus peptide as a function of time and colored according to the position above the cilia, from blue for those close to the cilia, to red for those closer to the top of the box. Bottom: Average velocity of the mucus peptides (blue) and fraction of cilia in the power stroke (orange) as a function of time.

Close modal

The correlation between cilia beating and mucus particle displacement is explicitly shown in the bottom part of Fig. 13, where we overlay the average velocity of the peptides and the fraction of cilia in the power stroke as a function of time. When the fraction of cilia in the power stroke increases, the peptides’ displacement reaches the maximum.

The non-constant fraction of cilia in the power stroke reflects the fact that the wave is not propagating uniformly in the carpet. This heterogeneous dynamics causes an oscillation of the fraction of cilia in the power stroke, with a characteristic period being half of that of the MCW propagation. The dynamics of cilia can be observed in Fig. 14, where snapshots of the carpet are displayed over the duration of two cycles (i.e., two peaks from Fig. 13, bottom left panel).

FIG. 14.

Successive snapshots of the wave over a duration of two peaks (see Fig. 13) for a ciliated cell with 10 × 10 cilia in a box of side L = 12.5 µm (ρc = 0.63 µm−2), colored according to the position of the beads along the z axis (white for the lower beads, blue for the higher ones).

FIG. 14.

Successive snapshots of the wave over a duration of two peaks (see Fig. 13) for a ciliated cell with 10 × 10 cilia in a box of side L = 12.5 µm (ρc = 0.63 µm−2), colored according to the position of the beads along the z axis (white for the lower beads, blue for the higher ones).

Close modal

To compute the power transferred by one cilium to the fluid, we compute its kinetic energy per unit time during the power stroke as follows:

Pin=3mcτpowb=1Nfvb2,
(8)

where 3mc is the mass of one triangular cross section of the structure, τpow is the duration of the power stroke, and vb, the average velocity of bead b along the active filament, during the power stroke. The power generated by a whole carpet of cilia is obtained by multiplying with the number of cilia that are simultaneously in the power stroke Npow, when it reaches a local maximum (i.e., for the peaks on the bottom left panel of Fig. 13). For a density of 0.63 cilia/μm2, averaging on the 100 cilia of the box and on 4–5 cycles, we obtain ⟨Pin⟩ ≈ 1.9 × 10−20 ± 7.7 × 10−20 W for a single cilium, with values being in a wide range 3.0 × 10−22–1.36 × 10−18 W.

To compute the power output, which effectively moves mucus peptides forward, we compute the average kinetic energy of the center of mass of a mucus peptide per unit time during a given acceleration phase (i.e., a given peak on the bottom left panel of Fig. 13) and multiply it by the number of peptides:

Pout=NpeptmpeptτaccelΔv2,
(9)

where mpept = 64mc is the mass of one peptide, Δv is the velocity gain during the acceleration phase, denoted by the positive slope of the peaks on the bottom layer of Fig. 13, and τaccel is the duration of this acceleration phase.

Then, we find ⟨Pout⟩ ≈ 3.7 × 10−20 ± 2.0 × 10−20 W, by averaging over all peaks displayed on the bottom left panel of Fig. 13.

The ratio ɛ = Pout/(NciliaPin) corresponds to the pumping efficiency of a system of Ncilia. To compute the average efficiency, we computed Pout for each peak and the corresponding efficiency. By considering all cilia in the system, we get an efficiency of about 2%. This order of magnitude is close to the upper boundary of the usual range of efficiency, between 0.1% and 1%, that is found for models of cilia driving the fluid.42,45,46 However, if we consider that only a fraction of the cilia is in the power strokes (Npow) during the peptide flow peak, and averaging all the peaks, we obtain a larger efficiency, ɛ ≈ 8%. As the order of magnitude of the fluid velocity is similar for all densities of cilia, we expect efficiency to remain in this order of magnitude for higher cilia densities.

We have implemented a working model of a respiratory cilium and simulated carpets of cilia using the LBMD simulation technique to account for hydrodynamic interactions. Special emphasis has been placed on tuning the model to reproduce the physiological characteristics of airway cilia and multiciliated cells. In particular, we have managed to reproduce the period and frequency of the beating cycle of a cilium, as well as obtain the emergence of the metachronal wave with a physiological wavelength from the collective coordination of a carpet of cilia. The simulated systems correspond, for size and density, to 1–4 multiciliated cells. Finally, we modeled the mucus layer, in order to account for its higher viscosity with respect to the PCL layer and explicitly quantify the clearance efficiency.

In conclusion, we have built a realistic and scalable framework for the investigation at the mesoscale, various aspects of the mucociliary clearance. As a next step, we will account for different levels of disorders in cilium orientation or in the spatial distribution of multiciliated cells to understand their effects on fluid transport and efficiency.23,47 Additionally, we will implement patho-physiological mechanisms, such as the removal or immotility of cilia,15 in order to quantify their impact on the clearance process, and how they impact the penetration of infectious particles in the PCL. From the methodological point of view, the implementation of a two-phase fluid to distinguish the response of the PCL and mucus layers could be explored.43,44

Finally, this framework could be used to help in the design of arrays of artificial motile cilia48,49 and ultimately to design new therapeutic strategies for respiratory diseases.

See supplementary material for cilium motion and analysis of the correlation function.

We acknowledge the financial support by the “Initiative d’excellence” program from the French State (Grant Nos. “DYNAMO,” ANR-11-LABX-0011-01, and “CACSICE,” ANR-11-EQPX-0008). Part of this work was performed using the HPC resources of GENCI (TGCC, IDRIS) (Grant No. x20226818) and LBT. We thank Geoffrey Letessier for technical support at the LBT.

The authors have no conflicts to disclose.

Emeline Laborie: Data curation (equal); Formal analysis (equal); Investigation (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Simone Melchionna: Methodology (equal); Software (supporting); Writing – review & editing (supporting). Fabio Sterpone: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Supervision (equal); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
E.
Houtmeyers
,
R.
Gosselink
,
G.
Gayan-Ramirez
, and
M.
Decramer
,
Eur. Respir. J.
13
,
1177
(
1999
).
2.
X. M.
Bustamante-Marin
and
L. E.
Ostrowski
,
Cold Spring Harbor Perspect. Biol.
9
,
a028241
(
2017
).
3.
M. J.
Sanderson
and
M. A.
Sleigh
,
J. Cell Sci.
47
,
331
(
1981
).
4.
M. A.
Chilvers
and
C.
O’Callaghan
,
Thorax
55
,
314
(
2000
).
5.
M.
Rautiainen
,
S.
Matsune
,
S.
Shima
,
K.
Sakamoto
,
Y.
Hanamure
, and
M.
Ohyama
,
Acta Oto-Laryngol.
112
,
845
(
1992
).
6.
R.
Hard
,
S. R.
Besch
,
D.
Tristram
,
J.
Han
, and
W.
Hicks
,
Laryngoscope
109
,
103
(
1999
).
7.
M. R.
Marino
and
E.
Aiello
,
Cell Motility
2
,
35
(
1982
).
8.
S.
Chateau
,
J.
Favier
,
U.
D’Ortona
, and
S.
Poncet
,
J. Fluid Mech.
824
,
931
(
2017
).
9.
J.
Elgeti
and
G.
Gompper
,
Proc. Natl. Acad. Sci. U. S. A.
110
,
4470
(
2013
).
10.
W.
Yi
,
K.
Park
,
C.
Rhee
, and
C.
Lee
, in
Annual International Conference of the IEEE Engineering in Medicine and Biology - Proceedings
(
IEEE
,
2000
), Vol. 3, p.
1783
null ; Conference date: 23-07-2000 Through 28-07-2000.
11.
A.
Burn
,
M.
Schneiter
,
M.
Ryser
,
P.
Gehr
,
J.
Rička
, and
M.
Frenz
,
Eur. Biophys. J.
51
,
51
(
2022
).
12.
L.
Morgan
,
M.
Pearson
,
R.
de Iongh
,
D.
Mackey
,
H.
van der Wall
,
M.
Peters
, and
J.
Rutland
,
Eur. Respir. J.
23
,
518
(
2004
).
13.
P.
Satir
and
M. A.
Sleigh
,
Annu. Rev. Physiol.
52
,
137
(
1990
).
14.
B. A.
Afzelius
,
J. Pathol.
204
,
470
(
2004
).
15.
R.
Robinot
,
M.
Hubert
,
G. Dias
de Melo
,
F.
Lazarini
,
T.
Bruel
,
N.
Smith
,
S.
Levallois
,
F.
Larrous
,
J.
Fernandes
,
S.
Gellenoncourt
,
S.
Rigaud
,
O.
Gorgette
,
C.
Thouvenot
,
C.
Trébeau
,
A.
Mallet
,
G.
Duménil
,
S.
Gobaa
,
R.
Etournay
,
P.-M.
Lledo
,
M.
Lecuit
,
H.
Bourhy
,
D.
Duffy
,
V.
Michel
,
O.
Schwartz
, and
L. A.
Chakrabarti
,
Nat. Commun.
12
,
4354
(
2021
).
16.
S. K.
Lai
,
Y.-Y.
Wang
,
D.
Wirtz
, and
J.
Hanes
,
Adv. Drug Delivery Rev.
61
,
86
(
2009
).
17.
Adivitiya
,
M. S.
Kaushik
,
S.
Chakraborty
,
S.
Veleri
, and
S.
Kateriya
,
Biology
10
,
95
(
2021
).
18.
L.
Xu
and
Y.
Jiang
,
Cells
8
,
736
(
2019
).
19.
D. J.
Smith
,
E. A.
Gaffney
, and
J. R.
Blake
,
Respir. Physiol. Neurobiol.
163
,
178
(
2008
).
20.
S. M.
Vanaki
,
D.
Holmes
,
S. C.
Saha
,
J.
Chen
,
R. J.
Brown
, and
P. G.
Jayathilake
,
J. Biomech.
99
,
109578
(
2020
).
21.
B.
Chakrabarti
and
D.
Saintillan
,
Phys. Rev. Fluids
4
,
043102
(
2019
).
22.
J.
Hall
and
N.
Clarke
,
J. Fluid Mech.
891
,
A20
(
2020
).
23.
G. R.
Ramirez-San Juan
,
A. J. T. M.
Mathijssen
,
M.
He
,
L.
Jan
,
W.
Marshall
, and
M.
Prakash
,
Nat. Phys.
16
,
958
(
2020
).
24.
J.
Hussong
,
W.-P.
Breugem
, and
J.
Westerweel
,
J. Fluid Mech.
684
,
137
(
2011
).
25.
Y.
Ding
and
E.
Kanso
,
Phys. Fluids
27
,
121902
(
2015
).
26.
A.
Malevanets
and
R.
Kapral
,
J. Chem. Phys.
112
,
7260
(
2000
).
27.
R.
Kapral
, in
Advances in Chemical Physics
, edited by
S. A.
Rice
(
John Wiley & Sons, Inc.
,
Hoboken, NJ, USA
,
2008
), pp.
89
146
.
28.
G.
Gompper
,
T.
Ihle
,
D. M.
Kroll
, and
R. G.
Winkler
, in
Advanced Computer Simulation Approaches for Soft Matter Sciences III
, edited by
C.
Holm
and
K.
Kremer
(
Springer
,
Berlin, Heidelberg
,
2009
), pp.
1
87
.
29.
H.
Fischer
,
I.
Polikarpov
, and
A. F.
Craievich
,
Protein Sci.
13
,
2825
(
2004
).
30.
P.
Ahlrichs
and
B.
Dünweg
,
J. Chem. Phys.
111
,
8225
(
1999
).
31.
S.
Succi
,
The Lattice Boltzmann Equation: For Fluid Dynamics and beyond
(
Oxford University Press
,
2001
).
32.
R.
Benzi
,
S.
Succi
, and
M.
Vergassola
,
Phys. Rep.
222
,
145
(
1992
).
33.
F.
Sterpone
,
P.
Derreumaux
, and
S.
Melchionna
,
J. Chem. Theory Comput.
11
,
1843
(
2015
).
34.
M.
Bernaschi
,
S.
Melchionna
,
S.
Succi
,
M.
Fyta
,
E.
Kaxiras
, and
J. K.
Sircar
,
Comput. Phys. Commun.
180
,
1495
(
2009
).
35.
M.
Chiricotto
,
S.
Melchionna
,
P.
Derreumaux
, and
F.
Sterpone
,
J. Chem. Phys.
145
,
035102
(
2016
).
36.
M.
Chiricotto
,
S.
Melchionna
,
P.
Derreumaux
, and
F.
Sterpone
,
J. Phys. Chem. Lett.
10
,
1594
(
2019
).
37.
F.
Sterpone
,
P.
Derreumaux
, and
S.
Melchionna
,
J. Phys. Chem. B
122
,
1573
(
2018
).
38.
O.
Languin-Cattoën
,
E.
Laborie
,
D. O.
Yurkova
,
S.
Melchionna
,
P.
Derreumaux
,
A. V.
Belyaev
, and
F.
Sterpone
,
Polymers
13
,
3912
(
2021
).
39.
S.
Timr
and
F.
Sterpone
,
J. Phys. Chem. Lett.
12
,
1741
(
2021
).
40.
G. R.
Fulford
and
J. R.
Blake
,
J. Theor. Biol.
121
,
381
(
1986
).
41.
E. R.
Brooks
and
J. B.
Wallingford
,
Curr. Biol.
24
,
R973
(
2014
).
42.
N.
Osterman
and
A.
Vilfan
,
Proc. Natl. Acad. Sci. U. S. A.
108
,
15727
(
2011
).
43.
A.
Tiribocchi
,
F.
Bonaccorso
,
M.
Lauricella
,
S.
Melchionna
,
A.
Montessori
, and
S.
Succi
,
Soft Matter
15
,
2848
(
2019
).
44.
A.
Tiribocchi
,
M.
Lauricella
,
A.
Montessori
,
S.
Melchionna
, and
S.
Succi
,
Int. J. Mod. Phys. C
30
,
1941004
(
2019
).
45.
C.
Eloy
and
E.
Lauga
,
Phys. Rev. Lett.
109
,
038101
(
2012
).
46.
H.
Guo
and
E.
Kanso
,
Phys. Rev. E
93
,
033119
(
2016
).
47.
A. V.
Kanale
,
F.
Ling
,
H.
Guo
,
S.
Fürthauer
, and
E.
Kanso
,
Proc. Natl. Acad. Sci. U. S. A.
119
,
e2214413119
(
2022
).
48.
W.
Wang
,
Q.
Liu
,
I.
Tanasijevic
,
M. F.
Reynolds
,
A. J.
Cortese
,
M. Z.
Miskin
,
M. C.
Cao
,
D. A.
Muller
,
A. C.
Molnar
,
E.
Lauga
,
P. L.
McEuen
, and
I.
Cohen
,
Nature
605
,
681
(
2022
).
49.
H.
Gu
,
Q.
Boehler
,
H.
Cui
,
E.
Secchi
,
G.
Savorana
,
C.
De Marco
,
S.
Gervasoni
,
Q.
Peyron
,
T.-Y.
Huang
,
S.
Pane
,
A. M.
Hirt
,
D.
Ahmed
, and
B. J.
Nelson
,
Nat. Commun.
11
,
2637
(
2020
).

Supplementary Material