## 2 Numerical Method
The control points can be treated as force generators to the fluid, which are computed by using the following direct forcing method at Cartesian grid point \((i,j,k)\):  
\[\overrightarrow{\mathbf{P}}_{i,j,k}^{n+1}=\rho\Big{(}\frac{\overrightarrow{ \mathbf{U}}_{i,j,k}^{n+1}-\vec{\mathbf{u}}_{i,j,k}^{n}}{\Delta t}+(\vec{ \mathbf{u}}\cdot\
abla\vec{\mathbf{u}})_{i,j,k}^{n}\Big{)}+\
abla p_{i,j,k}^ {n}-\mu\
abla^{2}\vec{\mathbf{u}}_{i,j,k}^{n}, \tag{4}\]  
where \(\Delta t\) is the marching time step between the time levels \(n\) and \(n+1\), and \(\overrightarrow{\mathbf{U}}_{i,j,k}^{n+1}\) is the distributed cilium velocity which can be obtained as  
\[\overrightarrow{\mathbf{U}}_{i,j,k}^{n+1}=f(i,j,k,\overrightarrow{\mathbf{X }}_{l}^{n+1})\cdot\overrightarrow{\mathbf{U}}(\overrightarrow{\mathbf{X}}_{ l}^{n+1}). \tag{5}\]  
Since the cilium control points \(\overrightarrow{\mathbf{X}_{l}}\) do not coincide with the fluid grid points \(\vec{\mathbf{x}}_{i,j,k}\), the filter function \(f(i,j,k,\overrightarrow{\mathbf{X}}_{l})\) is employed to transfer the information between those two grids [40], which is defined by  
\[f(i,j,k,\overrightarrow{\mathbf{X}}_{l})=\frac{1-tanh\big{(}(r-1)/0.5\big{)}} {2},\;\;r=\frac{|\vec{\mathbf{x}}_{i,j,k}-\overrightarrow{\mathbf{X}}_{l}|} {h}, \tag{6}\]  
where \(h\) is the Eulerian grid spacing (see Figure 2a). It should be noted that \(\Delta S/h<1\) to avoid penetration of streamlines into the Lagrangian grid spacing. Three-dimensional kinematics of the cilia \(\mathbf{U}(\mathbf{X}_{l})\) are prescribed using the ciliary beat profiles shown in Figure 1b,c, which is introduced as the'standard beat' throughout this work. It is assumed that cilia beat with the frequency of \(CBF=1/T\). It can be clearly seen in Figure 1d that, through one beat cycle, the magnitude of the cilium velocity components in the \(x\)-direction \(U\) is significantly higher compared with the other two velocity components--\(V\) in the \(y\)-direction and \(W\) in the \(z\)-direction. This means that the driving force imparted by cilia is dominant in the \(x\)-direction. The packing among the cilia is also depicted in Figure 1f. The metachronal wave produced by actions of cilia travels approximately at an angle of 135\({}^{\circ}\) to the effective stroke, as reported in the experiment on ciliary activity of cultured rabbit trachea epithelium by [PERSON] and [PERSON] [3] (see Figure 1e). Hence, considering one full wave over the computational domain using the same propagation angle, the initial arrangements of the cilia are chosen such that they are in similar phase along the wave fronts, with a fixed phase difference of \(\Delta\Phi=\pi/6\) and the interciliary spacing being \(3\sqrt{2}\)\(\upmu\)m. The velocity of the cilia \(\mathbf{U}(\mathbf{X}_{l})\) is prescribed accordingly and the position of them at every time step \(\Delta t\) is calculated as \(\mathbf{X}^{n+1}=\mathbf{X}^{n}+\Delta t\times U(\mathbf{X}_{l})\). Once the force exerted to the fluid \(\mathbf{F}\) was computed explicitly, Equations (1) and (2) are then discretised and solved for the velocity field and pressure using the pressure-increment projection algorithm, analogous to that described in Reference [41].