Metachronal patterns by magnetically-programmable artificial cilia surfaces for low Reynolds number fluid transport and mixing

Abstract

Motile cilia can produce net fluid flows at low Reynolds number because of their asymmetric motion and metachrony of collective beating. Mimicking this with artificial cilia can find application in microfluidic devices for fluid transport and mixing. Here, we study the metachronal beating of nonidentical, magnetically-programmed artificial cilia whose individual non-reciprocal motion and collective metachronal beating pattern can be independently controlled. We use a finite element method that accounts for magnetic forces, cilia deformation and fluid flow in a fully coupled manner. Mimicking biological cilia, we study magnetic cilia subject to a full range of metachronal driving patterns, including antiplectic, symplectic, laeoplectic and diaplectic waves. We analyse the induced primary flow, secondary flow and mixing rate as a function of the phase lag between cilia and explore the underlying physical mechanism. Our results show that shielding effects between neighboring cilia lead to a primary flow that is larger for antiplectic than for symplectic metachronal waves. The secondary flow can be fully explained by the propagation direction of the metachronal wave. Finally, we show that the mixing rate can be strongly enhanced by laeoplectic and diaplectic metachrony resulting in large velocity gradients and vortex-like flow patterns.

---

Bio-inspiration in the design of dynamic surfaces to control fluid motion requires a deep understanding of the behaviour of living organisms. Biological cilia are hair-like slender structures on the micrometre scale. Due to their small size, a single cilium operates at low Reynolds numbers with viscous forces dominating over inertial forces.¹ During the effective stroke, the cilium moves fast and almost like a straight rod, while during the recovery stroke, it rolls close to the surface in a slow tangential motion.^2,3 In addition to this individual asymmetric motion, biological cilia are also engaged in collective metachronal waves, where neighbouring cilia beat out-of-phase rather than synchronously.^1,4–7 Cilia can be found abundantly in humans and microorganisms to sense and generate motion.^8–10 Motile cilia can produce a net fluid flow at very low Reynolds numbers, and play an essential role in many living systems, such as the self-propulsion of Paramecium,⁴ enhancement of feeding functions in coral reefs,¹¹ and the swimming and capturing of food by starfish larvae.¹² Also inside our body, cilia play an important role by pumping biofluids, such as mucus in the airways¹³ and cerebrospinal fluid in the brain ventricle.¹⁴

By mimicking this biological function, large arrays of collectively beating artificial cilia also provide possibilities for fluid manipulation in lab-on-chip devices,^10,15 mobile microrobots,^16,17 and bioengineering systems.^18–21 Much work has been done on cilia material design and fabrication, that can be actuation strategies, and cilia motion optimization.²² Cilia can be fabricated by various methods, such as self-assembly techniques,²³ casting processes,^24,25 and photolithography,²⁶ and actuated in many ways, for instance, by electric fields,²⁷ magnetic fields,^28–30 light^31,32 and pH.³³ The experimental realization of metachronal motion has recently become possible by using pneumatic³⁴ or magnetic cilia.^26,35 In these systems, the asymmetric motion of individual cilia and metachronal waves can be independently controlled.

Metachrony is widely observed in nature, from single-celled marine organisms to insects and crustaceans.³⁶ In cilia-driven flows, there are four metachronal wave patterns:³⁷ antiplectic metachrony (AM) and symplectic metachrony (SM), where the direction of the wave propagation and the effective stroke are in the opposite and in the same direction, respectively, and laeoplectic metachrony (LM) and diaplectic metachrony (DM), where the effective stroke are to the left and to the right, respectively, of the propagation direction of the metachronal wave. The metachronal waves that are generated by cilia have been studied by many research groups, e.g. ref. 3,38,39 using a wide range of theoretical and computational techniques. The envelope model, for instance, has been used in many studies, including Blake⁴⁰ and Michelin and Lauga.^41,42 Guo and co-workers used surface undulations to represent the envelope of the cilia tips to identify two optimal metachronal beating patterns for pumping and swimming.⁴³ A 3D lattice spring model (LSM) was used to model cilia motion with a bead-spring configuration,^15,44,45 also coupled to the lattice Boltzmann method for the fluid.⁴⁶ Chateau et al.⁴⁷ studied the formation of antipleptic and symplectic metachronal waves in 3D cilia arrays that were immersed in a two-fluid environment, while others^48–50 studied metachronal waves by employing slender body theory.

Many studies have shown that antiplectic metachronal waves are more efficient than symplectic metachronal waves in transporting fluid, and different explanations and views have been proposed.^{10,26,28,35,39,51–54} Some of these studies also include the metachronal effect on mixing.^10,26,53 In 2015 and 2020, Gorissen et al.⁵⁵ and Milana et al.⁵² used a pneumatically actuated array of artificial cilia to create metachronal motion, showing that the fluid pumping capability is improved compared to synchronous motion. In 2017, Marume et al.⁵⁶ experimental created a metachronal motion by designing an array of magnetic artificial cilia with different magnetizations by controlling the orientation of the magnetic particle distribution inside the cilia. In 2018, Hanasoge et al.⁵⁷ also reported metachronal motion of an array of magnetic artificial cilia in which the cilia length was varied. In 2020, both Gu et al.³⁵ and Dong et al.²⁶ demonstrated programmable metachronal waves of magnetic cilia and showed their capability of fluid pumping. In 2021, Zhang et al.⁵⁸ created metachronal motion of microscopic magnetic artificial cilia (μMAC) and demonstrated that the metachronal μMAC can create pronounced fluid pumping in a microfluidic device. However, most investigations of artificial cilia only considered fluid transport and mixing due to SM and AM while in biological systems also waves can develop in perpendiculer directions (LM and DM). The advent of experimental artificial cilia systems in which the cilia motion and the metachronal wave direction can be indepently controlled, opens the exciting opportunity to fully exploit the coupled effect of metachrony in two perpendicular directions (AM/SM versus LM/DM) on fluid transport and mixing. This is the aim of the current paper. To do so, we carry out numerical simulations based on programmable magnetic artificial cilia developed recently.²⁶

In this paper, we simulate the beating kinematics and induced fluid flow of nonidentical, pre-programmed magnetic cilia subjected to a uniformly rotating magnetic field and analyse the resulting mixing response and primary/secondary fluid flow. In Section 1, we summarize the coupled three-dimensional computational solid-fluid interaction model that is used to describe the coupled cilia deformation and fluid flow, subjected to magnetic actuation. In Section 2, the fluid transport and mixing due to one cilium is optimized, based on the initial cilium orientation and magnetization profile angle. In addition, we investigate the fluid transport and mixing due to an array of programmable magnetic cilia subject to a broad range of metachronal waves. The final conclusions are presented in Section 3.

1 Methods

To describe the cilia deformation induced by the rotating magnetic field and the accompanying fluid flow at low Reynolds numbers, a three-dimensional computational solid-fluid interaction model is used based on the finite element method.⁴⁸ In this model, the cilia are represented by shell elements accounting for large deflections using an updated-Lagrange framework. Since we focus on low Reynolds numbers, we modelled the fluid flow by using the Stokes equation, the solution of which can be written by Green's functions in a semi-infinite fluid.⁵⁹ The drag forces on the cilia surface are regarded as a distribution of point forces f(r). The fluid velocity u^f(r) at a point due to these point forces f(r) exerted on the fluid by the cilia surface S at a position r′ can be written as

u^f(r) = G(r − r′, h(r′))f(r′), (1)

with h(r′) being the distance of between the point force f(r) and the bottom surface, and G(r − r′, h(r′)) is the Green's function for a point force f(r) acting in a Stokes flow near the no-slip boundary.⁵⁹ The point forces are assumed to be distributed on the surface of the cilia in the form of a traction t(r′), that is varying linearly over each triangular surface element,Here “nelm” is the total number of cilia surface elements. Since equation 2 is also valid for every node on the cilium surface u^c(rⁱ), the equations can be assembled in matrix form giving u^c= Gt. The traction t exerted by the cilia on the fluid can be found by inverting this relation giving the cilia surface traction: t = G⁻¹u^c. These fluid tractions (drag forces) are then imposed as external forces to the solid mechanics model of the cilia. The fluid and solid mechanics models are then coupled to each other by imposing no-slip conditions on the cilia surface. Finally, magnetic body couples are applied as external loads to the solid mechanics model. These magnetic couples N are due to the non-uniform distribution of the cilia magnetization M(s) and the externally applied rotating magnetic field (B(t)) by using N = M × B. More details on the numerical model can be found in previous work.⁴⁸ The model has been validated in terms of the net fluid flow for different phase shifts and intercilia spacings by a careful comparison with experimental data.²⁶

The cilia pumping efficiency was quantified by the average fluid volume transported through specific planes. The average volumetric flow rate in the x-direction through the plane (y₁ < y < y₂ and z₁ < z < z₂) for x = x₀, is defined by

with T the period of the beating cycles. Similarly, the average volumetric flow rate in the y-direction, through the plane (x₁ < x < x₂ and z₁ < z < z₂) for y = y₀, is defined byTo characterize the cilia-driven flow mixing efficiency, the overall mixing number m was used⁵³ by analysing the shortest distance between different-colour tracer particles (red and blue, see Fig. 1d and 8)where x_i and x_j are the positions of tracer particles of different colour (red, blue) with the minimum taken over all distances squared between particle i of one colour and the j = 1…N particles of the other colour (N is the total number of same-colour particles and there are equal numbers of red and blue particles). In this paper, we study chaotic mixing generated by the cilia motion by analyzing m as a function of the number of cycles N. Using the method described in,⁶⁰ we approximate the exponential decay using

ln(m/m₀) = −βN_cycles, (6)

where m₀ is the initial value of the mixing number and where the fitted parameter β represents the mixing rate.

Fig. 1 Schematic representation of the one-cilium setup (a) Schematic showing the initial cilia orientation θ and magnetization profile M(s) that is subjected to a uniformly distributed rotating magnetic field B(t) with frequency 1 Hz and magnitude 40 mT. (b) Non-reciprocal ciliary motion consisting of an effective stroke and a recovery stroke; the grey area is the positive swept area (S_A) of the cilium tip in a full beating cycle T. (c) Snapshot of fluid transport passing through the selected x–z plane (y = 0.8L) and y–z plane (x = 2L). (d) Snapshot of tracer particle mixing generated by the cilium after one cycle. The particles are placed in the box between − 0.5L ≤ x ≤ 2.5L, −0.6L ≤ y ≤ 1.4L, 0 ≤ z ≤ 2.5L, with a spacing distance 0.1L. To calculate the mixing rate β, 50 periods simulations are performed.

The initial configuration of each cilium is the one shown in Fig. 1(a), with no-slip boundary conditions taken on the cilia surface and the substrate. Before the flow is evaluated, the simulations have run for three periods to make sure that the flow has reached a steady state. We also changed the initial state of the driving magnetic field to make sure that the results are not affected by the initial transient.

2 Results

2.1 Fluid transport and mixing by a single cilium

Subject to a rotating magnetic field, the motion of the cilium is determined by its magnetic and geometric properties. We analyze a cilium of length L = 1 mm, width W = 0.75 mm, thickness T = 100 μm with Young's modulus 144 kPa and Poisson ratio 0.5. In this paper, the cilia beat frequency f is set to 1 Hz, the viscosity μ and density ρ of glycerol are 0.876 Pa s and 1.260 g cm⁻³, so that the Reynolds number equals 0.001, comparable to the recent experimental study.²⁶ The cilium's initial orientation θ and pre-programmed magnetization profile M (s), depicted in Fig. 1a, determine the motion of a single cilium. The non-reciprocal ciliary motion consists of an effective and a recovery stroke (Fig. 1b), leading to a non-zero swept area S_A in one beating period.

The ciliary non-reciprocal motion can generate a non-zero net flow.⁴⁸ We calculate the average volumetric flow rate Q_x and Q_y in the x-direction and y-direction through two perpendicular planes, as shown in Fig. 1c. In addition, tracer particles were added in the box around the cilium to compute the mixing rate β (Fig. 1d). By changing the cilium's magnetic and geometric properties, we obtain different fluid transport and mixing rates. Specifically, we analyse the effect of different initial orientations of the cilium, θ, and different angles α in the magnetization profile M(s) = M₀(cos(αs),0,sin(αs)) (see Fig. 1a and Fig. S1, ESI†) when subjecting the cilium to an external rotating magnetic field B(t) with frequency 1 Hz and magnitude 40 mT. It should be noted that the results in Fig. 2 will also depend on the cilia stiffness and actuation frequency. However, since our simulations are based on the experimental system that we developed in our previous paper, only small variations in elastic properties are allowed to ensure robust devices without chemical degradation. We therefore did not vary the elastic properties in the optimization procedure shown in Fig. 2. In addition, to ensure low Reynolds number conditions in both the experiments as well as the simulations, we constrained the frequency to values that are in accordance with the low Re regime.

Fig. 2 Fluid transport and mixing due to a single cilium. (a) Swept area S_A, (b) volumetric flow rate Q_x, (c) volumetric flow rate Q_y and (d) fluid mixing rate β, as a function of the cilium's initial orientation θ and magnetization profile angle α.

For fluid transport, we explore the volumetric flow rate through the selected planes y = 0.8L and plane x = 2L in the x-direction (Q_x) and the y-direction (Q_y), respectively, see Fig. 1c. From Fig. 2a and b we can deduce that there can be a positive swept area A (Fig. 1b), which leads to a positive flow in the x-direction, and a negative swept area B (Fig. S2, ESI†) which leads to a negative flow in the x-direction. In previous work,^26,48,61 we found that the fluid transport in the x-direction closely correlates with the swept area S_A, which is confirmed in Fig. 2a and b. The swept areas in Fig. 2a shows the existence of an optimal ciliary motion that maximizes fluid transport Q_x around area C (Fig. 2b). From Fig. 2c we see that the fluid transport is much smaller than Q_x and maximizes around area D. It should be noted that the non-zero fluid flow in the y-direction is due to the fact that we report the flow through a side-plane (y = 0.8L in Fig. 1c). The flow through the symmetry plane (y = 0.5L in Fig. 1c) is exactly equal to zero. Furthermore, the flow in the y-direction is suppressed if the cilium position is close to the surface with θ < 0.2π and when the swept area is zero (θ > 0.2π). Fig. 2d shows the fluid mixing rate β from which we can see that the optimal mixing happens in area E, which is very close to the maximum positive swept area A. At the same time, there is a second mixing area F, caused by the maximum negative swept area B, showing that the mixing effect is related to the amplitude of the swept area. Clearly, the reciprocal motion of the cilium can not only generate net flow but also mix the fluid (and possibly particles) around it.

For the study of the cilia arrays in the next section, we choose α = 1.85π and θ = 0.125π for each cilium as the optimal parameters for the swept area S_A which are close to the values used in the experimental set-up.²⁶ Note that in that case we have a positive swept area and the effective stroke is along the positive x-direction, see Fig. 1b.

2.2 Fluid transport induced by a 6 × 6 cilia array

Biological cilia typically form arrays and usually exhibit a collective wavelike metachronal motion. In this section, we study the dynamics of a 6 × 6 cilia array (Fig. 3). The period of metachronal waves is the period of the magnetic driven field, which is 1s. The wavelength can be expressed as . The cilia are geometrically identical and magnetically different due to a variation in phase shift in their magnetization distribution (M_ij(s) = M₀(cos(αs + Δϕ_ij),0,sin(αs + Δϕ_ij))) where Δϕ_ij = Δϕ_x·i + Δϕ_y·j for the i-th cilium in the primary direction along the x-axis and the j-th cilium in the secondary direction along the y-axis. Here, the primary direction is the beating direction of each cilium, while the secondary direction is orthogonal to the primary direction. When Δϕ_x = 0 and Δϕ_y = 0, all the cilia beat in phase, i.e. we have synchronous beating (see movie ESI,† 1). When Δϕ_x ≠ 0, Δϕ_y = 0, antiplectic metachrony (AM) (see movie ESI,† 2) and symplectic metachrony (SM) (see movie ESI,† 3) occur, when the wave propagation and the effective stroke are in the opposite and same direction, respectively.³⁷ On the other hand, with Δϕ_x = 0, Δϕ_y ≠ 0, laeoplectic metachrony (LM) and diaplectic metachrony (DM) occur, when the effective stroke is to the left and right, respectively, of the propagation direction of the metachronic wave (see movies ESI,† 4, 5 and 6 for laeoplectic metachrony for different values of Δϕ_y). The laeoplectic metachrony (LM) and diaplectic metachrony (DM) are symmetric in the plane y = 0. If both Δϕ_x ≠ 0 and Δϕ_y ≠ 0, antiplectic-laeoplectic metachrony occurs with Δϕ_x < 0 and Δϕ_y > 0 (see movie ESI,† 7 and 9) and symplectic-laeoplectic metachrony occurs when Δϕ_x > 0 and Δϕ_y > 0 (see movie ESI,† 8).

Fig. 3 The simulation box of the 6 × 6 cilia array arrangement. A snapshot is shown for an antiplectic-laeoplactic metachronal wave with Δϕ_x = −π/6 and Δϕ_y = π/6 at a fixed spacing L_x = 1.75L and L_y = 0.75W (see also movie ESI,† 7). The dimensions of the box are −6L ≤ x ≤ 6L, −3L ≤ y ≤3L, 0 ≤ z ≤ 5L. The cilium (i, j) = (1,1) is marked by the red triangle.

As readouts we compute the net fluid flow in the x-direction (Q_x) by subtracting the flow through plane P3 from the flow through P1. Similarly, we compute the net fluid flow in the y-direction (Q_y) by subtracting the flow through plane P4 from the flow through plane P2. Furthermore, we note that for the ciliary array studied, the system is fully symmetric with respect to laeoplectic and diaplectic metachrony: Q_x in the primary direction will remain unchanged upon a transition from LM to DM, while for Q_y in the secondary direction the flow will change sign but its magnitude will remain equal. In this work, we therefore only compute Q_x and Q_y for LM and use the inherent symmetries to generate the results for DM.

2.2.1 Fluid transport in the primary direction. For the various metachronal motions of the cilia array, we varied the phase shift Δϕ_x in the primary direction from −π to π, and the phase shift Δϕ_y in the secondary direction from 0 to π, while keeping a fixed spacing L_x = 1.75L (along the cilia length) and L_y = 0.75W (along the cilia width), see Fig. 3. Different beating patterns are shown in Fig. S3 (ESI†) as a function of phase shifts Δϕ_x and Δϕ_y.

For the effect of Δϕ_x, antiplectic metachronal waves produce a larger directional average volumetric flow rate Q_x when Δϕ_x ∈ [−π/2, −π/6] than the symplectic waves (Fig. 4a and c). Our simulation results are in close agreement with the experiments that studied the effect of Δϕ_x with Δϕ_y = 0.²⁶ We know from previous studies^26,28,52,54 that, despite the different cilia kinematics, the average flow for antiplectic waves is larger than from symplectic waves. The reason for this phenomenon is the shielding effect, meaning that the flow that is generated by each cilium is obstructed by the cilia surrounding it. Both the positive flow (Q_xp) caused by the effective stroke and the negative flow (Q_xn) caused by the recovery stroke decrease compared to synchronous waves²⁸ due to shielding. When the shielding effect is more prominent for the negative flow than for the positive flow, the net flow will be larger and vice versa. Fig. 5a shows that for antiplectic waves, the negative flow is more obstructed than for the symplectic waves, resulting in larger net flow for antiplectic metachrony. It is interesting to note that for a kinematically different motion²⁸ also AM was found to be more efficient, but with the shielding effect leading to an enhanced positive flow instead of a reduced negative flow as we found here.

Fig. 4 The average volumetric flow rate Q_x and Q_y as a function of Δϕ_x and Δϕ_y. (a) Contours of Q_x as a function of Δϕ_x and Δϕ_y. (b) Contours of Q_y as a function of Δϕ_x and Δϕ_y. (c) Q_x as a function of Δϕ_x for 13 values of Δϕ_y ranging from −π to π. (d) Q_x as a function of Δϕ_y for 13 values of Δϕ_x ranging from −π to π. (e) Q_y as a function of Δϕ_x for 13 values of Δϕ_y ranging from −π to π. (f) Q_y as a function of Δϕ_y for 13 values of ranging from −π to π. The average Q_x and Q_y are shown as solid lines separately.

Fig. 5 The volumetric flow rate in the primary direction as a function of (a) Δϕ_x for Δϕ_y = 0 and (b) Δϕ_y for Δϕ_x = 0. The red lines show the positive flow generated during the effective stroke (Q_xp). The blue lines show the absolute value of the negative flow rate generated during the recovery stroke (Q_xn). The black lines show the average volumetric flow rate Q_x = Q_xp − |Q_xn|.

In Fig. 4a and d, we see that diaplectic metachrony and laeoplectic metachrony can increase the flow rate Q_x further for Δϕ_y ∈ [−π/2, −π/6] ∪ [π/6, π/2]. This phenomenon can also be explained by the shielding effect, as shown in Fig. 5b, where the negative flow is more obstructed than the positive flow, especially in the region Δϕ_y ∈ [−π/2, −π/6] ∪ [π/6, π/2] where all six cilia in a lateral row have different phases. This fact results in stronger flow for Δϕ_y ∈ [−π/2, −π/6] ∪ [π/6, π/2]. By combining these two factors, we found a maximum Q_x for antiplectic-laeoplectic metachronal waves with and Δϕ_y ∈ [π/6, π/2] and, similarly, for antiplectic-diaplectic metachronal waves with Δϕ_x ∈ [−π/2, −π/6] and Δϕ_y ∈ [−π/2, −π/6].

2.2.2 Fluid transport in the secondary direction. For the secondary volumetric flow rate Q_y, as shown in Fig. 4b, it is interesting to observe that similar flow rates can be generated as in the primary direction. The maximum value of Q_y can be found for antiplectic-laeoplectic metachronal waves with Δϕ_x = −5π/6 and Δϕ_y = π/3 (see also movie ESI,† 9) and, similarly, antiplectic-diaplectic metachronal waves with Δϕ_x = −5π/6 and Δϕ_y = −π/3. It should be noted, however, that the effect of Δϕ_x on Q_y is almost symmetric, see Fig. 4b and e. Moreover, the smaller the absolute value of Δϕ_x, the weaker the effect of Δϕ_x on the volumetric flow rate Q_y. Finally, as shown in Fig. 4f, Q_y peaks at small Δϕ_y values, after which it decreases with an increase of |Δϕ_y|.

To further investigate the effect of metachrony on the secondary flow Q_y, we analyze the direction of the metachronal wave. The metachronal wave direction is described by the angle of the wave with the positive x-direction, denoted as γ ∈ [0,2π] with tanγ = (−Δϕ_y)/Δϕ_x·c. Here, c = L_x/L_y is a geometrical factor that incorporates the relative spacing of the cilia, see Fig. 3. The fluid volume Q_y for all cases studied is plotted as a function of γ in Fig. 6a (movies of a selection of these cases can be found in the ESI,†i.e. ESI,† 10–14). We observe that Q_y is zero when γ = 0 or π, i.e. when Δϕ_y = 0 and the wave only travels in the horizontal (primary) direction. For the cases with γ = π/2 and γ = 3π/2, the metachronal wave travels along the y-axis with Δϕ_x = 0. As can be observed in Fig. 6b, the flow pattern is almost circular, generating only a very small flow travelling in the y-direction. The small flow in these cases occurs because of a peristalsis-like mechanism, pushing the fluid in the direction of the metachronal wave. We see from Fig. 6b that the secondary flow direction coincides with the metachronal wave direction γ and peaks at wave directions that are at an angle of around π/4 with the principal axes.

Fig. 6 (a) The volumetric flow rate Q_y as a function of the metachronal wave propagation angle γ (c = 2.27). The lines in (a) are for systems with different Δϕ_x. (b) Fluid velocity amplitude in the x–y plane (z = L) at phases Δϕ_x = 0 and Δϕ_y = π/6 (above, red star in (a)) and Δϕ_x = π/3 and Δϕ_y = π/6 (bottom, blue star in (a)). The white arrows represent the streamlines. The red arrows represent the wave propagation direction γ.

2.2.3 Spacing effects on fluid transport. The cilia spacing is another critical parameter that influences the efficiency of cilia-driven fluid transport. Fig. 7 shows the effect of L_x and L_y on the flow. The inter-cilium spacing L_x varies from 1.5L to 4L and L_y varies from 0.5L to 3L for the optimal values Δϕ_x = −2π/3 and Δϕ_y = π/6 and keeping all other parameters the same. The box size (see Fig. 3) scales with the cilia spacing, keeping the cilia arrays' minimal distance at 0.8L and 1.6L to the box walls, respectively.

Fig. 7 The volumetric flow rate Q_x (a) and Q_y (b) as a function of cilia spacing L_x and L_y at a fixed phase shift and Δϕ_y = π/6. The red point indicates the reference spacing used in Fig. 3–6.

Fig. 7a shows that Q_x decreases with an increase in spacing L_x due to a reduction of cilia density in the flow direction. However, here the lateral spacings between cilia L_y also has an effect on primary flow Q_x, with a pronounced minimum for spacings around L_y = 1.5W. To investigate this further, we plot the negative and positive flow as a function of L_y for a fixed value of L_x (=4L) in Fig. S4 (ESI†). It can be observed that both Q_xp and Q_xn increase with spacing L_y due to a reduction in shielding. However, since the dependences of Q_xp and Q_xn on L_y differ, a maximum total flow is reached around L_y = 1.5W. Fig. 7b shows a similar effect of cilia density on flow, with a decreasing flow Q_y with increasing spacing L_y. The effect of L_x on Q_y is small.

2.3 Fluid mixing induced by a 6 × 6 cilia array

In addition to fluid transport, fluid mixing is another function that can be achieved by metachronal cilia beating. Similar to the single cilium mixing analysis, we use the forward evolution of tracer particles⁶² with an initial distribution in patterns of two colours, as shown in Fig. 8a. The propagating tracer particles displace in time, see Fig. 8b for the particle distribution after one cycle. To quantitatively explore the effect of the metachronal phase lag Δϕ on fluid mixing, the fluid mixing rate β due to different ciliary motions is studied for various cases with Δϕ_x ranging from −π to π and Δϕ_y from 0 to π (see Fig. 9). As discussed above, the system is fully symmetric with respect to laeoplectic and diaplectic metachrony, so we only compute β for laeoplectic metachrony and use the inherent symmetry to generate the results for diaplectic metachrony.

Fig. 8 3D tracer particle distribution for Δϕ_x = π/6 and Δϕ_y = π/6 at a fixed spacing L_x = 1.75L and L_y = 0.75W at initial time t = 0 (a) and after one beating cycle t = T (b).

Fig. 9 The fluid mixing rate β distribution after 3 cycles in the cilia arrays for (a) different phase shifts Δϕ_x and Δϕ_y at a fixed spacing L_x = 1.75 and L_y = 0.75W and (b) different cilia spacing L_x and L_y at fixed phase shift Δϕ_x = − 2π/3 and Δϕ_y = π/6. (c) The fluid mixing rate β as a function of Δϕ_x for 10 values of Δϕ_y ranging from 0 to π. (d) The fluid mixing rates β as a function of Δϕ_y for 13 values of Δϕ_x ranging from −π to π.

We observe in Fig. 9a and c that both antiplectic metachronal beating as well as symplectic beating can enhance the fluid mixing rate compared to synchronous motion (Δϕ_x = 0 and Δϕ_y = 0), in accordance with previous studies.^53,63 However, where the previous studies only focused on metachronal waves in the primary direction (Δϕ_y = 0), here, we also show that the mixing rate β is strongly enhanced by laeoplectic (and diaplectic) metachrony, increasing it by more than twice when Δϕ_y∈ [π/6, π/2] ∪ [−π/6, −π/2] compared to Δϕ_y = 0 (see Fig. 9d). From a comparison between Fig. 9c and d, it follows that the fluctuations due to a variation of Δϕ_x are rather shallow and discrete, but that the principal effect on the mixing rate is due to variations in Δϕ_y, creating trends that are much more pronounced than the local variations in the plateau region shown in Fig. 9a. The velocity field in the plane z = L is shown in Fig. 10. It is almost uniform for Δϕ_y = 0 (Fig. 10a), but for the velocity field is strongly inhomogeneous with pronounced velocity gradients (Fig. 10b). Because of the existence of both Δϕ_x and Δϕ_y, all the 6 cilia in a row and column have different phases. The cilia in the effective stroke generate a strong forward flow while those in the recovery stroke generate a backward flow. This leads to clockwise and anticlockwise vortex-like flow patterns formed between the neighbouring cilia which results in a largely increased mixing process. This phenomenon underlines the importance of combining LM (and DM) in the secondary direction with AM and SM in the primary direction to enhance mixing.

Fig. 10 Velocity field in the plane z = L at three snapshots in time t = 0, T/3 and 2T/3 from top to bottom generated by (a) Δϕ_x = π/2, Δϕ_y = 0 and (b) Δϕ_x = π/2, Δϕ_y = π/3. Background color represents flow speed. The cilia configurations are shown in black colour. The white arrows represent the streamlines.

To study the effect of cilia spacing, Fig. 9b shows the effect of L_x and L_y on the mixing rate for Δϕ_x = −2π/3 and Δϕ_y = π/6. The initial positions of the 20 × 10 × 10 particles in the x, y, and z directions are scaled with the box dimensions. We observe that the fluid mixing rate β increases when the cilia spacing L_x and L_y decreases. Clearly, by placing cilia closer together, the local flow profiles become more non-uniform, enhancing the propensity to form vortex-like flow patterns, which enhances mixing.

3 Conclusions

We computationally studied fluid transport and mixing by magnetically-programmable artificial cilia surfaces at low Reynolds number, using the finite element method. The beating kinematics of magnetically-actuated cilia with pre-programmed magnetization profiles was studied to explore the induced flow and mixing due to one single cilium. We showed that the induced primary volumetric flow rate, in the beating direction, Q_x strongly correlates with the swept area, which can be tuned through the magnetization profile to be positive as well as negative. Flow in the secondary direction Q_y (perpendicular to the primary direction) was found to be two orders of magnitude smaller than in the primary direction. The single cilium mixing rate maximizes around the regions where the swept area is maximially positive or negative.

To mimic the collective metachronal waves of biological cilia, a 6 × 6 cilia array was investigation. Flow in the primary direction was found to be larger for antiplectic metachrony compared to symplectic metachrony due to the enhanced shielding of negative flow in case of antiplectic metachrony. The effect of laeoplectic metachrony (and diaplectic metachrony) on primary flow was found to be symmetric, resulting in a small increase of flow compared to synchronous beating in the secondary direction. Also here, this increase was found to be due to enhanced shielding of negative flow. The effect of metachrony on secondary flow is of a completely different origin and can be ascribed to the metachronal wave direction γ leading to a strongly-enhanced flow for wave directions that are at an angle of approximately π/4 with the principal (primary and secondary) axes. Finally, we found that laeoplectic metachrony and diaplectic metachrony strongly enhance the mixing rate compared to synchronous, antiplectic and symplectic metachronal waves which is caused by a very inhomogeneous flow field with large velocity gradients inducing local vortex-like flow patterns. The results of this study can be used as a guideline for the design of programmable ciliated surfaces with optimal fluid transport and mixing rates.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We would like to thank Xiaoguang Dong and Metin Sitti (Max Planck Institute for Intelligent Systems, Stuttgart, Germany) for insightful discussions. Rongjing Zhang is financially supported by the Chinese Scholarship Council. The authors would like to thank the Centre for High-Performance Computing at the University of Groningen for providing HPC resources and support.

Notes and references

1. E. M. Purcell, Am. J. Phys., 1977, 45, 3–11.
2. J. Gray and J. S. Gardiner, Proc. R. Soc. London, Ser. B, 1922, 93, 104–121.
3. S. Gueron and K. Levit-Gurevich, Proc. Natl. Acad. Sci. U. S. A., 1999, 96, 12240–12245.
4. S. L. Tamm, J. Cell Biol., 1972, 55, 250.
5. M. Fliegauf, T. Benzing and H. Omran, Nat. Rev. Mol. Cell Biol., 2007, 8, 880–893.
6. W. Gilpin, M. S. Bull and M. Prakash, Nat. Rev. Phys., 2020, 2, 74–88.
7. J. Elgeti, R. G. Winkler and G. Gompper, Rep. Prog. Phys., 2015, 78, 056601.
8. P. Satir and S. T. Christensen, Annu. Rev. Physiol., 2007, 69, 377–400.
9. G. J. Pazour and G. B. Witman, Curr. Opin. Cell Biol., 2003, 15, 105–110.
10. A. Shields, B. Fiser, B. Evans, M. Falvo, S. Washburn and R. Superfine, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 15670–15675.
11. O. H. Shapiro, V. I. Fernandez, M. Garren, J. S. Guasto, F. P. Debaillon-Vesque, E. Kramarsky-Winter, A. Vardi and R. Stocker, Proc. Natl. Acad. Sci. U. S. A., 2014, 111, 13391–13396.
12. W. Gilpin, V. N. Prakash and M. Prakash, Nat. Phys., 2017, 13, 380–386.
13. M. R. Knowles and R. C. Boucher, et al. , J. Clin. Invest., 2002, 109, 571–577.
14. R. Faubel, C. Westendorf, E. Bodenschatz and G. Eichele, Science, 2016, 353, 176–178.
15. A. Alexeev, J. Yeomans and A. C. Balazs, Langmuir, 2008, 24, 12102–12106.
16. M. Sitti, Nat. Rev. Mater., 2018, 3, 74–75.
17. M. Cianchetti, C. Laschi, A. Menciassi and P. Dario, Nat. Rev. Mater., 2018, 3, 143–153.
18. R. M. Arco, J. R. Vélez-Cordero, E. Lauga and R. Zenit, Bioinspiration Biomimetics, 2014, 9, 036007.
19. V. Cacucciolo, J. Shintake, Y. Kuwajima, S. Maeda, D. Floreano and H. Shea, Nature, 2019, 572, 516–519.
20. H. A. Stone, A. D. Stroock and A. Ajdari, Annu. Rev. Fluid Mech., 2004, 36, 381–411.
21. C.-Y. Lee, C.-L. Chang, Y.-N. Wang and L.-M. Fu, Int. J. Mol. Sci., 2011, 12, 3263–3287.
22. X. Zhang, J. Guo, X. Fu, D. Zhang and Y. Zhao, Adv. Intelligent Systems, 2020, 2000225.
23. M. Vilfan, A. Potočnik, N. Osterman, I. Poberaj, A. Vilfan and D. Babič, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 1844–1847.
24. S. Zhang, Z. Cui, Y. Wang and J. M. den Toonder, Lab Chip, 2020, 20, 3569–3581.
25. S. Zhang, R. Zhang, Y. Wang, P. R. Onck and J. M. den Toonder, ACS Nano, 2020, 14, 10313–10323.
26. X. Dong, G. Z. Lum, W. Hu, R. Zhang, Z. Ren, P. R. Onck and M. Sitti, Sci. Adv., 2020, 6, eabc9323.
27. W. Jiang, L. Wang, H. Liu, H. Ma, H. Tian, B. Chen, Y. Shi, L. Yin and Y. Ding, RSC Adv., 2014, 4, 42002–42008.
28. S. N. Khaderi, J. M. J. den Toonder and P. R. Onck, J. Fluid Mech., 2011, 688, 44–65.
29. E. M. Gauger, M. T. Downton and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2009, 28, 231–242.
30. N. Banka, Y. L. Ng and S. Devasia, J. Med. Dev., 2017, 11, 031003.
31. H. Zeng, P. Wasylczyk, D. S. Wiersma and A. Priimagi, Adv. Mater., 2018, 30, 1703554.
32. A. H. Gelebart, M. Mc Bride, A. P. Schenning, C. N. Bowman and D. J. Broer, Adv. Funct. Mater., 2016, 26, 5322–5327.
33. L. D. Zarzar, P. Kim and J. Aizenberg, Adv. Mater., 2011, 23, 1442–1446.
34. E. Milana, B. Gorissen, S. Peerlinck, M. De Volder and D. Reynaerts, Adv. Funct. Mater., 2019, 29, 1900462.
35. H. Gu, Q. Boehler, H. Cui, E. Secchi, G. Savorana, C. De Marco, S. Gervasoni, Q. Peyron, T.-Y. Huang and S. Pane, et al. , Nat. Commun., 2020, 11, 1–10.
36. Y. Collard, G. Grosjean and N. Vandewalle, Commun. Phys., 2020, 3, 1–10.
37. S. Childress, Mechanics of Swimming and Flying, Cambridge University Press, 1981.
38. J. Blake, J. Fluid Mech., 1972, 55, 1–23.
39. J. Elgeti and G. Gompper, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 4470–4475.
40. J. R. Blake, J. Fluid Mech., 1971, 46, 199–208.
41. S. Michelin and E. Lauga, Phys. Fluids, 2010, 22, 111901.
42. S. Michelin and E. Lauga, Phys. Fluids, 2011, 23, 101901.
43. H. Guo, J. Nawroth, Y. Ding and E. Kanso, Phys. Fluids, 2014, 26, 091901.
44. A. C. Balazs, A. Bhattacharya, A. Tripathi and H. Shum, J. Phys. Chem. Lett., 2014, 5, 1691–1700.
45. A. Tripathi, A. Bhattacharya and A. C. Balazs, Langmuir, 2013, 29, 4616–4621.
46. Z. Guo, C. Zheng and B. Shi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 65, 046308.
47. S. Chateau, J. Favier, U. D’Ortona and S. Poncet, J. Fluid Mech., 2017, 824, 931–961.
48. S. Khaderi and P. Onck, J. Fluid Mech., 2012, 708, 303.
49. S. Gueron, K. Levit-Gurevich, N. Liron and J. J. Blum, Proc. Natl. Acad. Sci. U. S. A., 1997, 94, 6001–6006.
50. D. Smith, E. Gaffney and J. Blake, Bull. Math. Biol., 2007, 69, 1477–1510.
51. N. Osterman and A. Vilfan, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 15727–15732.
52. E. Milana, R. Zhang, M. R. Vetrano, S. Peerlinck, M. De Volder, P. R. Onck, D. Reynaerts and B. Gorissen, Sci. Adv., 2020, 6, eabd2508.
53. Y. Ding, J. C. Nawroth, M. J. McFall-Ngai and E. Kanso, J. Fluid Mech., 2014, 743, 124–140.
54. R. Zhang, J. den Toonder and P. R. Onck, Phys. Fluids, 2021, 33, 092009.
55. B. Gorissen, M. De Volder and D. Reynaerts, Lab Chip, 2015, 15, 4348–4355.
56. R. Marume, F. Tsumori, K. Kudo, T. Osada and K. Shinagawa, Jpn. J. Appl. Phys., 2017, 56, 06GN15.
57. S. Hanasoge, P. J. Hesketh and A. Alexeev, Soft Matter, 2018, 14, 3689–3693.
58. S. Zhang, Z. Cui, W. Ye and J. den Toonder, ACS Appl. Mater. Interfaces, 2021, 13, 20845–20857.
59. J. R. Blake, Math. Proc. Cambridge Philos. Soc., 1971, 70, 303–310.
60. J. B. Weiss and A. Provenzale, Transport and mixing in geophysical flows, Springer, 2007, vol. 744.
61. S. Khaderi, M. Baltussen, P. Anderson, D. Ioan, J. Den Toonder and P. Onck, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 046304.
62. Z. Stone and H. A. Stone, Phys. Fluids, 2005, 17, 063103.
63. S. Chateau, U. d'Ortona, S. Poncet and J. Favier, Front. Physiol., 2018, 9, 161.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d1sm01680f

---