## 2 Numerical Method  
In this work, the time-dependent Newtonian PCL flow is solved by a projection method and the force exerted by cilia on the fluid is computed via the IB method. Figure 1a shows a three-dimensional bounded domain \(\Omega=[-24,24]\times[-24,24]\times[0,8]\)\(\upmu\)m\({}^{3}\) that contains immersed massless boundaries (cilia). The ciliated domain is considered as \(\Omega^{\prime}=[-16.5,16.5]\times[-16.5,16.5]\)\(\upmu\)m\({}^{2}\). For this fluid-structure interaction problem, the incompressible Navier-Stokes equations are considered which can be written as  
\[\rho(\partial\mathbf{u}/\partial t+\mathbf{u}\cdot\
abla\mathbf{u})=-\
abla p +\mu\
abla^{2}\mathbf{u}+\mathbf{F}, \tag{1}\]  
\[\
abla\cdot\mathbf{u}=0, \tag{2}\]  
where \(\mathbf{u}=(u,v,w)\) is the fluid velocity, \(p\) is the fluid pressure, and \(\rho\) and \(\mu\) are constant density and dynamic viscosity of the fluid, respectively. Figure 1: Model schematics. (**a**) Ciliary array of motile cilia beating within the periciliary liquid layer (PCL). (**b**) Side view and (**c**) top view of beating kinematics of an individual cilium reported by [PERSON] [34], which shows the changes in shape of each cilium \(\mathbf{X}=(X_{1},X_{2},X_{3})\) at intervals of \(T/12\) through one beat cycle (\(0<t<T\)). Cilia move in a planner motion during their rapid effective stroke (phase No. \(\sim\) 1–7 ), but the recovering cilia (phase No. \(\sim\)8–12 ) are not performing a planner beat. When viewed from above, the ciliary tips follow a anticlockwise oval path. (**d**) The cilium velocity components. (**e**) Interpretation of ciliary activity on cultured rabbit tracheal epithelium from the model of [PERSON] and [PERSON] [3]. Arrows \(e\) and \(m\) indicate the directions of effective stroke and metachronal wave, respectively. (**f**) Adopted ciliary pack for all the simulation cases unless stated otherwise. Since one portion of an infinite problem is simulated in the present work, periodic boundary conditions are applied to the side walls normal to the \(x\) and \(y\) axes. [PERSON] et al. [35] observed that the fluid velocity at the epithelium boundaries remain static due to the slow movements of the basal bodies (a protein structure found at the base of a cilium) at the cell boundaries of cultured mouse tracheal multiciliated cells. Hence, the bottom boundary at the epithelium is treated as no-slip wall in our simulations. Given that the boundary condition at the PCL-mucus interface is influenced by the surface tension force between the ML and PCL, as well as the complex fluid rheology of the mucus, it is reasonable to assume the top boundary as the free slip wall in the present mucus-free model. The mentioned boundary conditions can be presented in the mathematical form as:  
\[ \begin{cases}\mathbf{u}=0&\text{epithelium}\\ \frac{\partial u}{\partial z}=0,\frac{\partial v}{\partial z}=0,w=0&\text{top boundary}\\ \mathbf{u}|_{x=-24\mu m}=\mathbf{u}|_{x=+24\mu m},\mathbf{u}|_{y=-24\mu m}= \mathbf{u}|_{y=+24\mu m}&\text{periodic boundaries in the $x$- and $y$-directions}.\end{cases}  \tag{3}\]  
Note that in this work the geometry is described using terms such as top and bottom, but within the respiratory tracts, the orientation may be in any direction. Gravity is not considered in the present model. To study the flow-structure interaction in biological systems, the immersed boundary (IB) methods have been applied successfully. Employing this method [36; 37; 38; 39] (see Figure 2b), each cilium can be represented by a set of \(N\) control points \(X_{l},l=1,\dots,N\) as functions of Lagrangian coordinate \(S\) (\(0<S<L\)) and time \(t\) (\(0<t<T\)) embedding in the Eulerian fluid domain. Thus, the cilium length can be given by \(L=\Delta S.N\) where \(\Delta S\) is the Lagrangian grid spacing along the IB, as depicted in Figure 2a. Here, \(T\) is the period of cilium beat cycle.