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.

## I. INTRODUCTION

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}

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 dynamics^{26–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/*μ*m^{2}^{23}). 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.

## II. MODEL AND METHOD

### A. Cilia dynamics

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 *N*_{f} beads, and arranged in a prism with a triangular base. In our implementation, *N*_{f} 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 *l*_{b} and harmonic constant *K*_{1}. 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 *K*_{2}. The harmonic constants are set to keep the bond length (and, therefore, the cilium length *L*_{c}) 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\alpha =12K\alpha (180\xb0\u2212\alpha i)2$, with *α* the angle defined by $AiAi+1Ai+2\u0302$, formed by three successive beads on a given filament.

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:

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 *i*_{0} is fixed at 5. During the backward stroke, *i*_{0} 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(|*l*_{i} − *l*_{b}(*i*)|) = *dl*_{max}.

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 *i*_{0} reaches the value of 12; afterward, we restore the equilibrium length *l*_{b}. 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 (*L*_{1}) with the length of the two remaining ones (*L*_{2} and *L*_{3}) by

The motion switches from the forward to the backward stroke when *dL* < *dL*_{min}, and vice versa when *dL* > *dL*_{max}. 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 *m*_{c} 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/cm^{3}.^{29} We then estimated the volume of the cilium by considering the experimental length *L*_{c} = 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 *m*_{c} = 5 × 10^{−18} kg for each bead. For the bonds, we chose a preferred length *l*_{b} = 0.3 *µ*m, to achieve a final cilium length of *L*_{c} = 6 *µ*m and width of 0.3 *µ*m. The values of the stiffness constants satisfy the following conditions *K*_{2} = *K*_{1}/10 and *K*_{α} = *K*_{1}/100.^{9} The value of *K*_{1} has been tuned to set the desired timescale for the beating cycle: *K*_{1} = 9500 × 10^{−2} (kcal/mol)/*μ*m^{2}. The propagation velocity of the pivot point *i*_{0} during the backward strokes [cf. Eq. (2)], $vi0$, is fixed at 3 ms^{−1}, meaning that the value of *i*_{0} 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 *i*_{0} before moving on to the next. Finally, the threshold values were empirically fixed at *dL*_{min} = −0.48 *µ*m and *dL*_{max} = 0.5 *µ*m to obtain the desired beating shape, and the stalling threshold is *dl*_{max} = 0.03*l*_{b}. 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/*μ*m^{2}, with the highest value approaching the experimentally reported density of 3 cilia/*μ*m^{2}.^{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 $\u223c8\mu m$; therefore, our systems correspond roughly to 1–4 MCC.

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 $\u223c600ms$ and $\u223c2000ms$, depending on the system.

In our setup, we used periodic boundary conditions for the *x*–*y* 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.

L (μm)
. | N_{cilia}
. | d_{c}/L_{c}
. | ρ_{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)
. | N_{cilia}
. | d_{c}/L_{c}
. | ρ_{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 |

### B. Lattice Boltzmann molecular dynamics

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

where $F\u20d7iC=\u2212\u2207Ui(r\u20d7i)$ corresponds to the conservative forces derived from the potential energy, and $F\u20d7iR$ is a random white noise accounting for thermal fluctuations. The term $F\u20d7iD=\u2212\gamma (v\u20d7i\u2212u\u20d7i)$ is a Stokes-like drag force, with $v\u20d7i$ the velocity of bead *i*, $u\u20d7i$ 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 *d*_{c}. 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.

## III. RESULTS

### A. Tuning the model

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.

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 *d*_{c}, 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 *d*_{c}, normalized by the cilium length *L*_{c} = 6 *µ*m. In the inset plot, we report the value of the beating period obtained while simulating a system at *d*_{c} = 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 *f*_{b} = 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}

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 $\u223c2.5$ times smaller than what is measured experimentally $\u223c0.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 *K*_{1}, 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 ⟨*v*_{tip}⟩ = 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

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}

### B. Metachronal wave characteristics

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\u20d7,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\u20d7,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\u20d7,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.

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

with $\delta B(r\u20d70,t)=B(r\u20d70,t)\u2212\u27e8B(r\u20d70,t)\u27e9r\u20d70,t$, where $r\u20d70$ is the position of a cilium and *t* the time. The average $\u27e8.\u27e9r\u20d70,t$ is taken over all cilia positions $r\u20d70$ and all frames *t*. Figure 8 shows an example for a cell of 20 × 20 cilia (*d*_{c}/*Lc* = 0.105). The correlation $C(r\u20d7,\tau )$ 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\u20d7,\tau )$ 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\u20d7,\tau )$ with an oscillatory function of the form

with *ω* the pulsation of the wave, and $k\u20d7$ 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 *T*_{wave} = 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}

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.

Δ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 |

### C. Effect of cilia on the fluid

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 *L*_{c}, see Fig. S6 in Ref. 9.

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.

### D. Model for the mucus layer

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 10^{3} 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 *L*_{x} = *L*_{y} = 12.5 *µ*m and *L*_{z} = 6.25 *µ*m, where we inserted about 360 peptide particles. These particles were free or interlinked via harmonic springs (*K*_{1}) to form short peptides of 4 or 8 beads.

An overview of the results is presented in Fig. 12(b). We obtain the smallest diffusion coefficient *D*_{min} ≈ 3.5 × 10^{−10} *μ*m^{2}/ns using the combination of 40 peptides of 8 beads, with each bead of mass 8*m*_{c}. 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 *L*_{x} = *L*_{y} = *L*_{z} = 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.

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).

### E. Pumping efficiency

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:

where 3*m*_{c} is the mass of one triangular cross section of the structure, *τ*_{pow} is the duration of the power stroke, and *v*_{b}, 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 *N*_{pow}, 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/*μ*m^{2}, averaging on the 100 cilia of the box and on 4–5 cycles, we obtain ⟨*P*_{in}⟩ ≈ 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:

where *m*_{pept} = 64*m*_{c} 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 ⟨*P*_{out}⟩ ≈ 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 *ɛ* = *P*_{out}/(*N*_{cilia}*P*_{in}) corresponds to the pumping efficiency of a system of *N*_{cilia}. To compute the average efficiency, we computed *P*_{out} 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 (*N*_{pow}) 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.

## IV. CONCLUSION

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 cilia^{48,49} and ultimately to design new therapeutic strategies for respiratory diseases.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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