Propulsion and Trapping of Microparticles by Active Cilia Arrays

Article pubs.acs.org/Langmuir

Propulsion and Trapping of Microparticles by Active Cilia Arrays Amitabh Bhattacharya,† Gavin A. Buxton,‡ O. Berk Usta,§ and Anna C. Balazs*,† †

Department of Chemical and Petroleum Engineering, University of Pittsburgh, 1249 Benedum Hall, Pittsburgh, Pennsylvania 15261, United States ‡ Robert Morris University, 6001 University Boulevard, Moon Township, Pennsylvania 15108-1189, United States § The Center for Engineering in Medicine, 51 Blossom Street, Boston, Massachusetts 02114, United States ABSTRACT: We model the transport of a microscopic particle via a regular array of beating elastic cilia, whose tips experience an adhesive interaction with the particle’s surface. At optimal adhesion strength, the average particle velocity is maximized. Using simulations spanning a range of cilia stiffness and cilia−particle adhesion strength, we explore the parameter space over which the particle can be “released”, “propelled”, or “trapped” by the cilia. We use a lower-order model to predict parameters for which the cilia are able to “propel” the particle. This is the first study that shows how both stiffness and adhesion strength are crucial for manipulation of particles by active cilia arrays. These results can facilitate the design of synthetic cilia that integrate adhesive and hydrodynamic interactions to selectively repel or trap particulates. Surfaces that are effective at repelling particulates are valuable for antifouling applications, while surfaces that can trap and, thus, remove particulates from the solution are useful for efficient filtration systems.

1. INTRODUCTION In various living organisms, the coordinated motion of hairlike filaments, called “cilia”, is vital for performing such activities as movement, feeding, or pumping of fluids within the animal. For example, paramecia use cilia to undergo directed movement and marine suspension feeders use cilia to propel food into their bodies. In humans, the synchronized undulations of cilia in the respiratory tract help to pump viscous fluids away from the lungs, as well as propel dust and other particles out of the body.1 Inspired by the utility and efficiency of such biological cilia, researchers have begun to design synthetic analogues that could be harnessed to regulate flow in microfluidic devices.2−11 The synthetic cilia are fashioned from flexible filaments that are anchored to the floor of a microchannel and actuated by an applied field2−4,6,12−14 or forces.15These oscillating filaments not only provide effective pumping of fluids,2b,3 but can also promote the mixing of dissolved components 2b,3 and controllably switch the directionality of the net flow within microchannels.15 To date, however, the possibility of using active, synthetic cilia to direct the movement of microscopic particles, such as biological cells and microcapsules, within microchannels has not been extensively explored.16−18 Developing approaches for conveying cells or microcarriers to specified locations within microfluidic devices19,20 is vital for performing accurate microscale analysis or chemical synthesis. While the coordinated motion of the cilia is effective at propelling the surrounding fluid,1,3,15,21,22 it is also plausible that an adhesive interaction between the cilia and particulates is necessary for controlling the particle movement. In mammals, for example, beating cilia © 2012 American Chemical Society

located at the entrance of the oviduct cannot transport egg cells unless there is a critical level of adhesion between the cilia tips and cells.23,24 The combined effects of such cilial adhesion and generated flows might be responsible for transport in other biological environments.24 Notably, to the best of our knowledge, there have been no prior computational or theoretical studies aimed at probing the role of adhesion between motile cilia and particles in solution. Here, we use simulations and a lower-order analytical theory to model the interactions between microscopic particles and actuated cilia (see Figure 1) that encompass a “sticky” tip. As we show below, for a critical adhesion strength, active cilia can propel particles from one neighbor to the next, significantly increasing the particle’s velocity. If the adhesive interaction between a particle and this sticky tip is sufficiently high, the particle can become trapped within the cilia layer. These results can facilitate the design of synthetic cilia that integrate adhesive and hydrodynamic interactions to selectively repel or trap particulates. Surfaces that are effective at repelling particulates are valuable for antifouling applications, while surfaces that can trap and, thus, remove particulates from the solution are useful for efficient filtration systems. Furthermore, the findings can provide insight into physical factors that influence adhesive interactions between biological cilia and particulates. The paper is organized in the following manner. In section 2, we first describe the hybrid computational approach we developed to simulate the cilia−particle system, and then in section Received: December 8, 2011 Published: January 10, 2012 3217

dx.doi.org/10.1021/la204845v | Langmuir 2012, 28, 3217−3226

Langmuir

Article

The fluid velocity, u(x), is evolved using the lattice Boltzmann method (LBM),27 which can be viewed as an efficient solver for the Navier−Stokes equation. The hydrodynamic coupling between the fluid and the beads comprising the cilia is captured through a frictional force that acts on the beads.28,29 Namely, the force on the ith cilium bead is given by

FiH = −ζ(r i̇ − u(r i))

(2.2)

where ζ is the coupling constant, ṙ is the velocity of the ith particle, and u(ri) is the interpolated fluid velocity at the position, ri, of the particle i. The hydrodynamic coupling constant has the form ζ = 6πηa, where a = 0.15 is the input hydrodynamic radius of the particle in lattice units.27,28 To conserve local and therefore global momentum in the system, an equal and opposite force is interpolated onto the fluid grid. To actuate the cilia, we apply a cyclic force to the extreme tip of each cilium (i.e., at i = N). (We note that a force applied to the cilial tips has been considered as a mode of actuation in prior computational and theoretical studies.15,17,30,31) In particular, the force has the following functional form: i

Figure 1. Three-dimensional snapshot from the simulation. Black filaments denote cilia, and the red sphere is the particle.

3 we discuss the results we obtained via this approach. Also within section 3, we relate our results to typical experimental values for the ciliated surfaces. In section 4 we present the lower-order theoretical model we derived to describe the system and show that this model can reproduce a key feature of the phase map obtained via the simulations. Finally, we conclude by discussing the implications of these findings in section 5.

Fext(t ) = Fx0 cos(ωt )i + F y0 sin(ωt ) cos(ψ)j + F y0 sin(ωt ) sin(ψ)k

where ω is the frequency of actuation. The plane in which the actuation force oscillates harmonically is tilted by a small angle ψ with respect to the x−y plane (Figure 2). Hence, a cilium tip also moves in a slightly tilted plane; this ensures that a cilium tip does not collide with itself (i.e., with the lower half of the cilium) during the return or “recovery” stroke. The cilium tip is farther from the wall during the effective stroke than in the recovery stroke (Figure 3). This motion is

2. COMPUTATIONAL MODELING Each cilium has a length L and is formed from a chain of N beads of radius a (Figure 2). The coordinates of the N beads are specified as

Figure 2. Schematic of cilium and actuation force (blue arrow). The plane of actuation, tilted with respect to the x−y plane by an angle ψ, is indicated in light blue. Only red beads in the cilium stick to the particle.

Figure 3. (a) Superimposed 3D snapshots of an E = 3E0 cilium, taken at equal time intervals, during one cycle of forcing. The cilium is black when the tip moves in negative x direction (effective stroke) and red when the tip moves in positive x direction (recovery stroke). (b) Path of cilium tip over one cycle, projected onto the x−y plane, for E = E0 (red, dashed), E = 3E0 (blue, dot−dashed), and E = 5E0 (green, solid).

{ri; i = 1, ..., N}, and the first bead in the chain (i = 0) is tethered to the wall (y = 0). The elastic potential energy of the chain is25,26

U chain =

N

∑ (1/2)[ks(|ri − ri − 1| − l)2 + k b(θi − π)2 ] i=1

(2.3)

responsible for net fluid flow in the direction of the effective stroke (negative x direction). Given a fixed cilia geometry and low Reynolds number, the cilial motion depends solely on the Sperm number Sp = L[4πμω/EI]1/4, which is the ratio of the viscous to the elastic forces on a cilium.32 We first simulated actuated cilia without the particle for E ∈ {E0 → 5E0}, where Sp based on (E0,ω) equals 4.3. It is important to note that the imposed actuation force is taken to be proportional to cilia stiffness, i.e., Fx0 ∝ E and the ratio Fy0/Fx0 = 1/4 is fixed. The reason the tip forcing is proportional to cilia stiffness is as follows. To facilitate

(2.1)

In eq 2.1, θi is the angle between neighboring bonds and is given by θi = cos−1(bî +1,i·b̂i−1,i), where b̂i,j = [ri − rj]/|ri − rj|, with b̂0,1 = j being the wall normal direction. The constants ks, kb, and l = L/N represent the stretching modulus, bending modulus, and equilibrium distance between adjacent beads, respectively. Here, kb = EI/l, where E is Young’s modulus and I = πa4/4. Figure 1 shows the 15 × 12 cilia array we consider in our simulation; the intercilium separation is δx = δz = L/3. 3218


Langmuir

Article

when z = 1/2, so the maximum tension that can be sustained by the spring is equal to Fmax = Dλ/2. In other words, the bond will break at tensile forces greater than Fmax. The hydrodynamic interactions between the mobile particle and fluid are captured by introducing “bounce-back” boundary conditions at the particle−fluid surface within the LBM code.27 No-slip boundary conditions are applied to the fluid at the particle surface and the top and bottom walls. Additionally, periodic boundary conditions are applied in the wall-parallel directions (x and z). The box size was chosen with the following constraints in mind: the box size parallel to the wall should be much larger than 2(R + L) (R is the radius of the particle and L is the length of the cilia) so that a given cilium cannot interact with both the particle and the periodic images of the particle. The 15 × 12 lattice size was sufficiently large to meet that constraint and yet sufficiently small to allow the simulations to run in computationally reasonable time frames. (Given that the flow is primarily along x (due to the actuation force), we note that the hydrodynamic interactions along the z direction are weaker than along x and, hence, the box size could be a bit smaller in the z direction.) The respective dimensions of the computational box are Lx × Ly × Lz = 45 × 24 × 36 LB units, and the time period of actuation is T = 30 000 LB time units. The length of each cilium is L = 9 LB units. Here, 1 LB unit of length denotes the grid spacing and 1 LB unit of time denotes the time step used in the lattice Boltzmann scheme. The dimensionless values for the relevant parameters used in the simulations are N = 9, a/L = 0.016, l/L = 1/9, Re = ωL2/ν = 0.1, Sp(E0) = L[4πμω/E0I]1/4 = 4.3, Fx0L2/EI = 64, Lx/L = 5, Lz/L = 4, Ly/L = 2.7, R/L = 0.58, λL = 189, (rc−re)/L = 3.66 × 10−3, and re/L = 0.055.

comparisons between cilia of different stiffnesses, we attempted to maintain a similarity between the trajectories of the different cilial tips. If a cilium is stiff (i.e., low Sperm number), the elastic forces dominate, so the external force needs to be proportional to EI/L2 to ensure that the cilium trajectory remains the same with increasing stiffness. Here, we keep L and I fixed and make the external force proportional to cilium stiffness, E, in order to obtain similar tip trajectories for the different cilia. Figure 3b clearly shows that the trajectory of the cilium tip does not vary over a large range of stiffness when we make this choice. Into this system, we then introduce a spherical, neutrally buoyant particle that has a radius R and position rs. The three beads at the cilium tip (i ∈ {N − 2, N − 1, N}) are “sticky”; that is, whenever |rs − ri| − R < rc, a springlike interaction is introduced between the tip and nearest point on the particle surface, (say, rp). The latter interaction is modeled via an attractive Morse potential:

V a(|r p − r i|) = D[1 − exp{−λ(|r p − r i| − re)}]2

(2.4)

applied only when |rp − ri| > re, and the spring is untethered when |rp − ri| ≥ rc (see Figure 4). Here D is the adhesion strength and

3. RESULTS FROM THE COMPUTATIONAL MODEL 3.1. Validation of the Model. We validated the above point-particle LBM scheme by comparing the output from these simulations to findings from a previous study,15 in which we coupled the lattice spring model (LSM) to the lattice Boltzmann model (LBM) to simulate an elastic cilium that is immersed in a fluid and actuated by a sinusoidal force, which is applied to the tip of the cilium. The latter LBM/LSM approach involved resolving the cilium using a hexahedral mesh and explicitly imposing no-slip conditions at the surface of the cilium (i.e., at the bounding face elements of this hexahedral mesh). In particular, an interpolated “bounce-back” scheme was used to impose the no-slip condition. We refer to the latter simulation approach as “LSM”, while we refer to the current point-particle based approach as “PPM”. In the prior study,15 we specifically focused on a cilium that is tilted with respect to the bottom surface, and we found that in this configuration the actuated cilium gave rise to a nonzero bulk flow in the system. Moreover, by simply altering the Sperm number, Sp (by varying the frequency of the applied force), we could controllably switch the direction of the net flow. Figure 5 shows a plot of the bulk flow Ub as a function of Sp from both the LSM simulations and the PPM simulations. For both sets of simulations, the base of the cilium was titled at an angle of α = 45° with respect to the surface, along the positive x axis. Furthermore, a periodic force with amplitude Fy0 is applied in the vertical direction to the tip of the cilium, and the size of the periodic computational box is kept as Lx = L, Ly = 1.5L, and Lz = 0.5L. For the PPM simulations described in this section, L = 10 LB units. Based on the hydrodynamic radius of the point particles, the aspect ratio of the cilium (given by 2a/L) in the PPM simulations is approximately 0.05, while the aspect ratio used in the LSM simulations is 0.1. The latter mismatch cannot be avoided, since, in the PPM simulations, the bond length between the point particles has to be approximately 1 LB unit,

Figure 4. Schematic of bond formation of cilium tip with surface of sphere. (a) Example of bond formation (green line) between cilium tip and point of attachment on surface of spherical particle. (b) Details of a tacky bond. λ characterizes the range of the interaction. (The value of λ is chosen so that the length scale over which the Morse potential falls off is much smaller than the size of the particle, i.e., 1/λ ≪ R.) We choose rc = re + (ln 2)/λ, the separation at which the attractive force is maximized, i.e., d2Va/dr2|r=rc = 0. Each adhesive bead can form only one bond with the particle. The points of contact rp on the particle are evolved in time using quaternion equations.33,34 A repulsive exclusion force, arising from a similar Morse potential (but with a different value for D), acts on all cilium beads and particles when |rs − ri| − R < re. Similar exclusion forces are implemented for bead−bead and bead−wall interactions. The positions of the cilia and particles are evolved via a velocity Verlet algorithm. The bond strength for the cilia−particle interaction is determined by finding the maximum tension that can be sustained by the Morse potential spring. The Morse potential can be rewritten as V(r) = D[1 − z(r)]2, where z(r) = exp[−λ(r − re)]. The force arising from the Morse potential is then given as F(r) = −dV/dr = −2Dλz(1 − z). From the latter expression, one can see that F is minimized (i.e., dF/dz = 0) 3219


Langmuir

Article

ratio of the cilium should not affect our validation significantly, since slender body theory36 indicates that the hydrodynamics of a flagellum depends weakly (i.e., on the logarithm) on its aspect ratio. In our comparison of the PPM and LBM simulations (Figure 5), the net fluid velocity Ub = [1/(LxLyLz)]∫ ⟨u1(x)⟩ d3x is rescaled to yield a dimensionless number and A = Fy0L2/(3EI) = Fy0L2/(3kbl) = 10 is the dimensionless forcing amplitude. The PPM simulations are able to reproduce the main results in the prior study,15 namely, actuating the bent cilium leads to a net flow in the channel, and that, above a critical Sperm number, Sp = Sp0, the bulk velocity Ub changes sign. The values of U b from the point-particle-based simulations agree reasonably well with the values from the LSM simulations for 1 < Sp < 3 and 5 < Sp < 9 (see Figure 1). The agreement, however, breaks down closer to the critical Sperm number (i.e., for 3 < Sp < 5), and the values of Sp0 do not match for the two methods: Sp0 = 3.7 for LSM and Sp0 = 4.3 for PPM. This may

Figure 5. Plot of bulk velocity in the channel versus Sperm number due to actuated bent cilium, calculated using LBM coupled with (solid) LSM and (dashed) PPM schemes.

and the radius of each point particle has to be much less than 1 LB unit for the model to be valid.27,28 The mismatch in aspect

Figure 6. (a) Snapshots of LBM simulation at E = 3E0 and D = 1.33D0, taken over the effective stroke of a cycle. Each panel in (i)−(iv) shows the x−y projection of unattached cilia (black filaments), attached cilia (green filaments), particle (red sphere), and velocity field (blue arrows). The cross section of the velocity field is taken at the x−y plane cutting through the particle center, located at z = rzs . (b) Mean of streamwise fluid velocity plotted with respect to height, for E = 3E0. 3220


Langmuir

Article

be due to the fact that, in the LSM simulations,15 the cilium buckles close to Sp = Sp0. In the PPM simulations, the cilium does not have any cross-sectional moment of inertia, which is critical to capturing the effects of buckling, and therefore leads to the discrepancy between the two methods around Sp0. Nonetheless, the PPM approach also revealed the fundamental phenomenon of flow reversal by an actuated bent cilium, and yielded reasonable agreement with a highly resolved LSM simulation15 over a range of Sperm numbers. 3.2. Particle Trajectories and Velocities. We now focus on the 15 × 12 cilia array shown in Figure 1 and introduce a particle into the solution. Initially, the cilia is actuated, but the particle is fixed at a random position, within the interaction range of the cilia (i.e., rys (0) < L + R). After t > 10π/ω, the particle is freed, and the cilia−particle adhesion, with strength D, is activated. In the simulations, we varied only E and D within the following ranges: E ∈ {E0 → 6E0} and D ∈ {D0 → 10D0}, where D0 = Fx0/2λ ∝ E. Four independent simulations were carried out at each (E,D), for 60 actuation cycles. Figure 6a shows snapshots of the actuated cilia, the induced fluid motion, and the motile particle; the simulations are taken over the effective stroke of a cycle, with E = 3E0 and D = 1.33D0. Each panel shows the x−y projection of the cilia (green filaments) that are attached to the particle (red sphere), the unattached cilia (black filaments), and the velocity field (blue arrows). The cross section of the velocity field is taken at the x−y plane cutting through the particle center, located at z = rzs . The size of the blue arrows indicates the magnitude of the fluid velocity induced by the cilia. The average fluid velocity has a Couette flow profile UxFluid(y) ∝ (Ly − y), for y ≥ L, as can be seen in Figure 6b. The ensuing particle motion depends strongly on both adhesion strength and cilia stiffness. Our results show that any given particle trajectory converges to a “released”, “propelled”, or “trapped” state. Figure 7a displays a trajectory where the

effective stroke, and the particle is thus pushed/pulled forward; the particle detaches from the cilia during the recovery stroke. Effectively, over each actuation cycle, the particle is passed from one cilium to the next, in the x direction. Unlike the “released” states, a “propelled” state is characterized by a lower average particle height and significant variation in the particle height. Figure 7c shows a trajectory at an even higher adhesion level. At this level of adhesion, the cilia attach to and do not release a “trapped” particle. Consequently, the time-averaged velocity of the particle is equal to zero. The “trapped” state is characterized by negligible average velocity and large variation in particle height. Figure 8 shows the time-averaged particle velocity, U(E,D) = ⟨ṙxs⟩, plotted with respect to D/D0, for differing cilium stiffness

Figure 8. Average particle velocity over the trajectory, U = ⟨ṙxs⟩, versus adhesion strength, for different values of Young’s modulus, E. Error bars at each (E,D) show maximum and minimum values of U for four independent simulations.

E, at fixed ω. The averaging was performed over the last 25 actuation cycles of each simulation. U(E,D) was normalized with respect to the average particle velocity at zero adhesion strength U0(E) = U(E,D)|D=0; the relative velocity U*(E,D) = U(E,D)/U0(E) indicates the change in average velocity due to the effects of adhesion. For all curves, there is an optimal adhesion level where U* is maximized. Notably, at E = 6E0, the particle moves roughly 100% faster at this adhesion; significant increases also occur at other E values. The maximum in U* occurs due to the “propelled” states at these adhesion levels, where the cilium periodically attaches to the particle surface, transferring momentum in the process (see discussion below). During these “propelled” states, the adhesive cilium also pulls the particle closer to the lower wall, where the net streamwise fluid velocity is higher, leading to faster particle advection. As noted above, the average fluid velocity induced by the cilia has a Couette flow profile UxFluid(y) ∝ (Ly − y), for y ≥ L (Figure 6b). Therefore, during the “released” state, for which the average particle height is ⟨rys ⟩ ≈ L + R, the particle samples a lower fluid velocity than it would during the “propelled” state, for which ⟨rys⟩ ≈ L. It is also apparent from Figure 8 that increasing the adhesion strength above the optimal value eventually leads to the “trapped” state (U* ≈ 0). To gain further insight into the factors contributing to the increase in particle velocity seen in Figure 8, we decompose and average the forces acting on this propelled microscopic particle

Figure 7. Particle trajectories (solid black curves), starting at green circles, ending at red squares, projected on x−y plane, for E = 3E0. Filtered particle trajectories, with height h(t), shown by dashed red curves. Trajectory at (a) low adhesion, D = 0.66D0, shows “released” state, (b) medium adhesion, D = 1.33D0, shows “propelled” states, and (c) high adhesion, D = 5D0, shows particle “trapping”.

particle escapes the cilia and, eventually, attains the “released” state, at a height rys ≈ R + L with respect to the lower wall. This particle is transported almost completely by fluid advection. Figure 7b shows the particle trajectory for a higher adhesion level, corresponding to the “propelled” state. Here, cilia tips periodically attach themselves to the particle during the 3221


Langmuir

Article

the particle during the effective stroke gives rise to the increase in the velocity of the propelled particle at the optimal adhesion level. When the adhesion is increased beyond the optimal strength, however, the adhesive tips of the cilia are unable to detach from the particle during the return stroke. As a result, the net interaction force from the cilia during the return stroke begins to counteract the extra force from the cilia during the effective stroke. To illustrate this effect, we plot the forces acting on the particle during the effective and return strokes for a fixed stiffness and varying adhesion strength, averaged over several actuation cycles and random seeds (Figure 10). Specifically, we

from the nodes in the cilia during the effective stroke (i.e., the first half of the actuation cycle). The effective stroke is chosen for this analysis since, at optimal adhesion, the cilia do not exert large forces on the particle during the return stroke. We first denote the location of the ith node in the array as ri(t) at any given time t. (Unlike in section 2, i here denotes the global index of a cilia node in the simulation box.) The force from this node on the propelled particle, FMorse,i(t), is classified as a “pulling” force, FPull,i(t), or a “pushing” force FPush,i(t), depending on whether it points into (i.e., FMorse,i·(ri − rs) > 0), or away (i.e., FMorse,i·(ri − rs) ≤ 0) from the particle center, respectively. For any given time step, we then calculate the spatial averages, FxPull(δx,t) and FxPush(δx,t), of the x component of these forces over the y−z planes located at a distance δx with respect to the particle center. Specifically, FxPull/Push(δx,t) is the sum of the FxPull/Push,i(t) over cilia nodes i that lie in the y−z plane given by x = rxs + δx. In the propelled state, the particle undergoes approximately repetitive motion over several actuation cycles (Figure 7b), and therefore, we can average the integral of these forces over the effective stroke, as follows: Pull/Push

Fx̃

(δx) =

1 Nact

Nact − 1 ⎡

∑

j=0

⎢ ⎣

jT + T /2

∫jT

⎤ FxPull/Push(δx , t ) dt ⎥ ⎦

(3.1)

where Nact is the number of actuation cycles over which the average is being taken, and T = 2π/ω is the time period of the actuation force. F̃ xPull/Push(δx) gives the net pull/push forces acting on the particle, at an x-coordinate displaced by δx with respect to the particle center, during the effective stroke. Figure 9 Figure 10. Decomposition of forces acting on sphere, plotted with respect to adhesion strength. The force over the whole cycle is decomposed into interactions due to Morse potential FM and hydrodynamic forces FxH. The Morse potential force, FM, is further Eff and return decomposed into contributions from effective stroke FM Ret stroke FM . Error bars show maximum and minimum values of the force over four independent simulations.

calculate the total force on the particle from the cilia during the effective stroke (FxM,Eff), during the return stroke (FxM,Ret), and during the whole actuation cycle (FxM = FxM,Eff + FxM,Ret). We also calculate the net hydrodynamic force FxH acting on the particle over the whole actuation cycle. In Figure 10, we first note that FxM,Ret and FxM,Eff always have opposite signs. As expected, the total interaction force FxM is maximized at the optimal adhesion strength (D = 1.66D0). Above the optimal adhesion strength, FxM,Ret increases more rapidly, eventually overcoming FxM,Eff, leading to the particle slowing down and getting trapped. (Note that, at the highest adhesion strength (D = 5D0), where the particle has reached the trapped state, the net interaction force FxM is not zero, and in fact points in the positive x direction. This is not surprising, since, in the trapped state, FxM has to resist the net hydrodynamic force FxH acting on the particle due to the net fluid flow induced by the actuated cilia in the positive x direction.) 3.3. Phase Map for System. From the above simulations, we can generate the phase map shown in Figure 11. A measure that differentiates a trajectory in the “released” state from that in the “propelled” state is the variance h′ = [⟨(rys − h)2⟩]1/2 of the particle height rys (t), with respect to its moving average h(t) t+T/2 s ry (t′) dt′, where T = 2π/ω. The height h(t) gives an = ∫ t−T/2 estimate of the average height of particle trajectory, and will not depend on t if rys (t) becomes periodic, with time period T.

Figure 9. Spatial distribution of cilia−particle interaction forces acting on particle over effective stroke. The forces have been binned into different δx locations, where δx is the x displacement of the cilium bead with respect to the center of the particle. Error bars denote standard deviation of the forces, taken over several (∼25) actuation cycles.

shows the contribution from F̃xPull/Push(δx), as well as F̃xTotal = F̃ xPull + F̃ xPush, as a function of δx for the optimal adhesion level D = 1.66D0 at E = 3E0. The effective stroke points in the negative x direction, and therefore, negative δx corresponds to the “front” and positive δx corresponds to the “back” of the particle. As can be seen from the contributions to F̃xTotal, the force F̃xPull dominates at the front of the particle (δx/R < 0), while the force F̃ xPush dominates at the back of the particle. This combination of pulling (from the front) and pushing (from the back) on 3222


Langmuir

Article

of the cilia tips, per unit area.35 Our results indicate that, for E = E0 = 0.132 MPa, the I/II boundary will occur at σ = 293 Pa and γ = 6.1 × 10−5 J m−2 (which follows from Dc* = Dc/D0 = 2 at E = E0), while, for E = 6E0 = 0.794 MPa, the boundary will occur at σ = 967 Pa and γ = 2 × 10−4 J m−2 (which follows from Dc* = 1.1 at E = 6E0).

4. ANALYTICAL MODEL We now develop a lower-order model (with six degrees of freedom) to show that our prediction of a sharp transition at the I/II boundary is indeed robust. Using this lower-order model, we also show that this boundary can be captured by focusing solely on the particle−cilia interaction in the y direction. The key elements of this model are lumped-mass equations for cilia motion and the equation for particle motion in the y direction. The latter equations are then coupled with models for bond formation and breakage. The lumped-mass equation is derived for the position of the cilium tip, which is described by local coordinates r(t) = X(t)i + Y(t)i + Y(t)j = r(t)r̂ (Figure 12a), the origin r = 0 being the

Figure 11. Phase map of particle motion as a function of (D,E). Simulations were run at each point marked by a green symbol. Crosses mark phase I, circles mark phase II, and squares mark phase III. Solid magenta lines divide the phases. The yellow dashed line marks the I/II boundary predicted by the low-order model. The contour plot of the height variance of the particle from LBM simulations, h′(D,E), is shown in the background.

In Figure 11, the shaded grayscale contour map of h′ is plotted in the background. A trajectory in the “released” state is characterized by a small variance, i.e., h′/L ≪ εh (darkest shading), while a trajectory in “propelled” state has a large variance, h′/L > εh (lighter shading). We choose εh = 0.03 as the threshold for sensing a “propelled” state. A particle is “trapped” if its average velocity is low, i.e., U* < εu; here, we choose εu = 0.2. Using the above definitions, we separate the (E,D) plane into “released” (I), “propelled” (II), and “trapped” (III) states (Figure 11). The I/II boundary can be sensitive to the value of εh; therefore, in Figure 11, this boundary is also displayed via the contour plot of h′(E,D). Interestingly, we find that, for a given stiffness E, the increase in h′ with D is not gradual. Rather, for any E, the I/II boundary marks a critical adhesion level Dc(E), above which the particle trajectory makes a sharp transition in h′, going from being in a “released” state to a “propelled” state, regardless of initial particle depth in the array. It is worth recalling that the adhesion strength in the phase map is scaled by D0 = Fx0/2λ, which, as noted above, is proportional to E (since we took F0 ∝ E). The vertical nature of the phase boundaries between the different regions in the phase map indicates that this indeed was an appropriate scaling. 3.4. Comparison to Experimental Values. From recent experimental studies on synthetic cilia, we can specify the length of the cilia to be L = 20 μm and the beating frequency to be 1/T = 40 Hz.2a If we assume the surrounding fluid to be water, then we can take the kinematic viscosity of the fluid to be ν = 10−6 m2 s−1 and the density of the fluid to be ρ = 103 kg/m3. It follows from the dimensionless parameters used in the simulation (i.e., R/L = 0.58, a/L = 0.016, Sp(E0) = L[4πμω/ E0I]1/4 = 4.3, λL = 189, see section 2), that the radius of the particle is R = 12 μm, the cross-sectional radius of each cilium is a = 333 nm, the inverse of length scale over which Morse potential decays is λ = 9.6 × 106 m−1, the moment of inertia of the cilium cross section is I = πa4/4 = 9.7 × 10−27 m4, and the scale of Young’s modulus is E0 = 0.132 MPa. On the basis of these values, we can specify the I/II boundary in terms of surface energy γ = D/πa2, and interfacial tensile strength σ = Dλ/2πa2

Figure 12. (a) Schematic of cilia−particle interaction with cilia. Position vectors for interacting and not interacting cilia (inset) are represented by r ̅ (black vector) and r (green vector), respectively. (b) Different stages of particle−cilia interaction during an effective stroke. The width of the green shaded region gives the diameter (2d) of the contact area over which the cilia interact with the particle at that instant. The stages depicted in panels (i), (ii), (iii), and (iv) correspond to the snapshots Figure 6a, (i), (ii), (iii), and (iv), respectively.

point where the cilium is tethered to the wall. The tip experiences an external force Ftip(t) due to actuation forces and/or interaction with the particle. The viscous drag per unit length on the cilia (using approximations from resistive force theory36) is C∥ = 2παμ/ln(L/a) along the filament, and C⊥ = 2C∥ perpendicular to the cilia. For a rigid rod, the constant α equals 1 3223


Langmuir

Article

in a boundless fluid medium; we use α = 0.25, which we obtain from the simulations (see below). The lower value of α can be justified as an increase in the mobility of the cilium due to hydrodynamic interactions from closely spaced, neighboring cilia. Assuming overdamped motion, we equate viscous drag with forces and moments on the whole cilium to obtain the evolution equation for the tip:

r ̇ = r ṙ ̂ + rθθ̇ ̂ = μ−1[−∇r U E(r) + Ftip]

tips, resists motion due to viscous drag, leading to a nonzero contact area between the top of the cilia array and the particle (Figure 12b,ii). Once the free cilia tips start to recede from the particle (i.e., Ẏ ≤ ṙys ), the bonds stop forming (Figure 12b,iii); the cilia tips that are already attached to the particle stay attached, until the viscous drag on the particle overcomes the maximum tensile force (∼Dλ/2) that an average cilia−particle bond can sustain (Figure 12b,iv). Due to forces from the Morse potential, cilia forming bonds with the particle will be out of phase with the free cilia (Figure 12b,iv). The average position vector of the bound cilia tips is (i) described by local coordinates r(t) ̅ = (1/n)∑i∈C Δr (t), where Δr(i) is the local coordinate of the ith cilium tip and C is the set of n cilia bound to the particle. Given that d is the radius of the spherical cap over which the bonds are formed, the number of bonds is n = ρcπd2. Below, r(t) denotes the local position vector of the free cilia tips, while r(t) denotes the average local ̅ position vector of attached tips. Assuming R ≫ δx, and a uniform distribution of bonds over the contact area, the above summation becomes an integral:

(4.1)

Here, UE = [kr±(r − L)2/L + kθ(π/2 − θ)2]/2 is the elastic energy of this cilium, and μ−1 is its mobility matrix, defined as μrθ−1 = μθr−1 = 0, μrr−1 = 2/(C∥L), μθθ−1 = 3/(C⊥L). The constants kr+ and kr− are the respective stretching and compression moduli and kθ is the bending modulus, while L is the equilibrium length. Only the effective stroke of the cilium (i.e., 0 < t < π/ω) is calculated, since the particle transitions from the “released” state only during this time frame. The constants kr+ = 400EI/L2, kr− = 0.02EI/L2, kθ = 5EI/L, and α = 0.25 are calibrated with respect to a cilium simulated using the above point-particle LBM method, for the particular cilia grafting density ρc = 1/δx2 = 9/L2. As a validation of this lower-order model, we use eq 4.1 to determine the tip velocity and, as can be seen from Figure 13, find good agreement between the values generated from this approach and those obtained via our LBM simulations. Notably, the comparison is better for the vertical component of the motion than the lateral component. This is due to the fact that, despite the reduced value for α, the lower-order model underestimates the cooperative interactions between the cilia in the x direction. On the other hand, due to the presence of the walls, these cooperative interactions are damped in the y direction, leading to a better estimate for ṙy via the model. Since the transition from the “released” state to the “propelled” state occurs at a low threshold for the variance h′, we assume that the height of the particle is rys ≈ R + L, for which the particle forms bonds only with the cilia tips (see Figure 12a). These bonds cannot support a strong force parallel to the sphere’s surface (see Figure 4), since the length of a bond, rc, is close to the equilibrium distance for the Morse potential (i.e., rc − re ≪ re). Furthermore, we assume that the intercilium spacing is much smaller than the particle radius, i.e., δx/R ≪ 1. Hence, the Morse potential force per unit area (or the adhesive stress) in the lower-order model has primarily a vertical (y) component, and is approximately uniformly distributed throughout the cilia−particle contact area. During the effective stroke (see Figure 12b), cilia−particle bonds are formed when tips of free cilia approach the particle (i.e., Ẏ > ṙys ). The particle, while being pushed up by the cilia

2 d r ̅ (t ) = 2 Δr(ζ, t )ζ dζ d 0

∫

(4.2)

Here, Δr(ς,t) is the local coordinate of a cilium tip that is attached to the sphere at a radius ζ from the y axis of the sphere. Keeping in mind that both the contact radius d and Δr(ς,t) change with time, the time derivative of r(t) ̅ can be obtained from eq 4.2 as follows:

⎡ d ⎛ 2 ⎞⎤⎡ d r̅ = ⎢ ⎜ 2 ⎟⎥⎢ dt ⎣ d t ⎝ d ⎠⎦⎣

∫d

⎡ 2 ⎤⎡ d ⎛ + ⎢ 2 ⎥⎢ ⎜ ⎣ d ⎦⎣ d t ⎝

0

∫d

⎤ Δr(ς, t )ς dς⎥ ⎦ 0

⎞⎤ Δr(ς, t )ς dς⎟⎥ ⎠⎦

(4.3a)

⎡ 4 d d ⎤⎡ 0 ⎤ d r̅ = ⎢− 3 Δr(ς, t )ς dς⎥ ⎥⎢ ⎦ ⎣ d d t ⎦⎣ d dt ⎡ 2 ⎤⎡ d d + ⎢ 2 ⎥⎢ (Δr(d , t )d) ⎣ d ⎦⎣ d t

∫

+

∫d

0

⎤ ∂Δr(ς, t ) ς dς⎥ ⎦ ∂t

(4.3b)

To evaluate the second term in the right-hand side (RHS) of eq 4.3a, we use the chain rule to differentiate with respect to t, so that d/dt = (dd/dt)(∂/∂d) + ∂/∂t. Recalling that (d/dx)[∫ 0x f(t) dt] = f(x), we obtain the second term in eq 4.3b. Next, we assume Δr(ς, t)|ς=d = Δr(d,t) = r(t); a cilium at the edge of the

Figure 13. Comparison between LBM simulation (magenta) and lower-order model (blue) for cilium tip velocity in (a) x and (b) y directions for E = 5E0. 3224


Langmuir

Article

We use β = 4 to reflect the decrease in mobility due to the walls; this value is calibrated from the LBM simulations. Due to incompressibility and large particle size, we assume uy ≈ 0 for the fluid velocity around the particle. A separate exclusion force also acts on the particle to ensure that its height stays above y = ymin ≈ R + L/2. The variables r, r,̅ d, and rys are evolved over 0 ≤ t ≤ T using eqs 4.1, 4.4b, 4.5, 4.6, and 4.7, along with the following initial conditions at t = 0:

contact radius will have the same length and orientation as an unattached cilium. This is because, outside the contact area, only unattached cilia exist (see Figure 12b,iii) and we assume (for Ẏ > ṙs) that there is a continuous change in cilia height between free and attached cilia at the contact radius. The latter assumption is strictly true only when the free cilia tips are approaching the particle in the vertical direction (i.e., when Ẏ > ṙs). Now substituting eq 4.2 into the first term on the RHS of eq 4.3b, we obtain

⎡ 2 ⎤⎡ d d ⎡ 2 dd ⎤ d r̅ = ⎢− r ̅(t ) + ⎢ 2 ⎥⎢ (rd) ⎥ ⎣ d dt ⎦ ⎣ d ⎦⎣ d t dt 0 ∂Δr(ς , t ) ⎤ + ς dς⎥ ⎦ ∂t d

∫

r(0) = r ̅(0) = r 0 = L cos(θ0)i + L sin(θ0)j; d = 0;

∫

d(Y ̇ − ryṡ )

for all δ > 0 and Y ̇ − ryṡ ≥ 0

2δ

5. CONCLUSIONS In summary, we formulated both computational and theoretical approaches to model the interaction between active, adhesive cilia and microscopic particles in solution. In the process, we developed methods that provide a general framework for predicting the interaction of particles with a range of active biological and synthetic cilia. Using these approaches, we found that the particle velocity and trajectory depend strongly on the particle−cilia adhesion and cilia stiffness E (or, alternatively, the Sperm number). Notably, the simulations revealed that, at an optimal adhesion level, the cilia can significantly increase the velocity U at which a particle is propelled along the top of the layer. This increase in U is due to a concerted pushing and pulling exerted by the respective cilia in front and behind the particle. This increase in the particle velocity was most pronounced for the stiffest cilia considered here. Via the computational and theoretical models, we also determined a phase map as a function of E and D and thereby pinpointed the critical adhesion levels Dc(E), above which a particle goes from being in a “released” state to being in a “propelled” state, regardless of initial particle depth in the array. These results revealed that Dc/D0 decreases with increasing E, at a rate much slower than 1/E, implying (since D0 ∝ E) that Dc increases with increasing E. In other words, floppier cilia can capture particles at lower adhesion strength. This result is consistent with the understanding that particles get stuck more easily on softer passive adhesive surfaces.37 We also pinpointed the region in phase space where the particles become trapped within the active layer.

(4.5)

with ḋ = 0 otherwise. To derive the second term in the RHS of eq 4.4b, we approximate (following δ ≪ R) the local position vectors of all the attached cilia by the average local position vector r ̅ (Figure 12b,iii). Hence, the forces acting on these tips will be uniform, so using eq 4.1, we can state that, within the contact radius (i.e., 0 < ζ < d):

∂Δr(ς, t ) = μ−1[−∇r ̅ U E( r ̅) + Fext(t ) ∂t + F Morse(r ys , ry̅ )] 2 d ∂Δr(ς, t ) = 2 ς dς ∂t d 0

∫

(4.6)

where F = −∂V − R − ry̅ |)/∂ry̅ j is the Morse potential force acting per cilia−particle bond. The second equality in eq 4.6 follows from the first equality, which states that ∂Δr/∂t does not depend explicitly on ς. The particle undergoes overdamped motion in response to the equal and opposite Morse potential force from n bonds: Morse

ryṡ =

(rys ,ry̅ )

− nF yMorse(r ys , ry̅ ) 6πβR μ

a

(4.8)

with θ0 ≈ 0.2 rad. (At this configuration, the elastic and external forces balance each other, so Fx0L sin(θ0) = kθ(π/2 − θ0).) We seek a particle trajectory r̃ys (t) that is periodic, with the same time period T as the effective stroke, so that r̃ys (T) = r̃ys (0). This solution is obtained iteratively, by repeating the evolution over several cycles of the effective stroke, assigning rys (0) for a cycle to rys (T) calculated from the previous cycle; r̃ys (t) is reached when |rys (T) − rys (0)| converges to zero. Using the above formalism, we obtain, for a range of (E,D), the height variance h′ of r̃ys (t). As above, when h′/L < εh, the periodic trajectory is binned as a “released” state; otherwise it is identified as a “propelled” state. The yellow line in Figure 11 indicates that the lower-order model captures the I/II boundary, in good agreement with the predictions obtained from the simulations. Taken together, the different approaches corroborate our findings that the particle motion can be controlled by tailoring the stiffness of the cilia or the stickiness of its tips.

(4.4a)

2(r − r ̅) dd d r̅ 2 0 ∂Δr(ς, t ) = + 2 ς dς (4.4b) ∂t dt d dt d d In eq 4.4b, we need to determine ḋ and the second term in the RHS. To evolve d, we divide the effective stroke into Ẏ − ṙys > 0 and Ẏ − ṙys ≤ 0 (see Figure 12b). Bonds are formed when the top surface of the cilia array approaches the sphere (i.e., Ẏ − ṙys > 0), moving unattached tips within a distance of re + rc of the particle surface. We denote R′ = R + re + rc as the radius that is relevant for cilia−particle interaction. If the particle is submerged in the array (i.e., δ = Y − (rys − R′) > 0) and Ẏ − ṙys > 0, then we model the contact radius d in terms of the depth of the sphere in the array, i.e., d = (R′δ)1/2. When the top surface of the array recedes from the sphere (i.e., Ẏ − ṙys ≤ 0), we assume that the number of tips interacting with the particle stays constant. Using the definition of δ (the depth of the particle in the layer), along with the latter assumptions, we evolve d according to ḋ =

r ys(0) ≈ R + 0.9L

(|rys

(4.7) 3225


Langmuir

Article

(23) Norwood, J. T.; Hein, C. E.; Halbert, S. A.; Anderson, R. G. W. Proc. Natl. Acad. Sci. U.S.A. 1978, 75 (9), 4413−4416. (24) Talbot, P.; Riveles, K. Reprod. Biol. Endocrinol. 2005, 3, 52. (25) Kirkwood, J. G. J. Chem. Phys. 1939, 7, 506. (26) Keating, P. N. Phys. Rev. 1966, 145, 637. (27) Ladd, A. J. C. J. Fluid Mech. 1994, 271, 285−309. (28) Ahlrichs, P.; Dunweg, B. Int. J. Mod. Phys. C 1998, 9 (8), 1428− 1438. (29) Usta, O. B.; Ladd, A. J. C.; Butler, J. E. J. Chem. Phys. 2005, 122, 09490. (30) Vilfan, A.; Julicher, F. Phys. Rev. Lett. 2006, 96, 058102. (31) Khaderi, S. N.; Baltussen, M. G. H. M.; Anderson, P. D.; Ioan, D.; den Toonder, J. M. J.; Onck, P. R. Phys. Rev. E 2010, 82, 027302. (32) Lowe, C. P. Philos. Trans. R. Soc. London, B 2003, 358, 1543− 1550. (33) Omelyan, I. P. Phys. Rev. E 1998, 58, 1169. (34) Martys, N. S.; Mountain, R. D. Phys. Rev. E 1999, 59, 3733. (35) The adhesion scale D0 is itself a linear function of Young’s modulus E, so that D0(E) = Fxext/2λ = 32EI/L2λ (section 3.2). Using the values of E0, I, L, and λ in section 3.4, we obtain D0(E0) = 1.07 × 10−17 J and D0(6E0) = 6.42 × 10−17 J. The latter values of D0, along with the dimensional values λ = 9.6 × 106 m−1 and a = 333 nm, yield estimates of γ = D/πa2 and σ = Dλ/2πa2 for D = Dc(E0) = 2D0(E0) and D = Dc(6E0) = 1.1D0(6E0). (36) Brennen, C.; Winet, H. Annu. Rev. Fluid Mech. 1977, 9, 339− 398. (37) Jagota, A.; Bennison, S. J. Integr. Comp. Biol. 2002, 42, 1140− 1145.

These results provide guidelines for tailoring the functionality of the ciliated surface by tailoring the stiffness of the cilia and the “stickiness” of the cilia tips. Namely, for a range of E and D, the cilia are highly effective at preventing particles from penetrating the layer and, thus, can provide valuable antifouling properties to a range of surfaces. Conversely, there is a range of E and D where the beating cilia can be used to extract particles from solution and, thus, could act as an active filtration system.

■ ■ ■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected].

ACKNOWLEDGMENTS A.C.B. gratefully acknowledges the ARO for partial support of A.B. and the NSF for partial support of A.C.B.. REFERENCES

(1) Sleigh, M. Comp. Biochem. Physiol. 1989, 94A, 359. (2) (a) Evans, B. A.; Shields, A. R.; Carroll, R. L.; Washburn, S.; Falvo, M. R.; Superfine, R. Nano Lett. 2007, 7, 1428. (b) Shields, A. R.; Fiser, B. L.; Evans, B. A.; Falvo, M. R.; Washburn, S.; Superfine, R. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 15670−15675. (3) den Toonder, J.; Bos, F.; Broer, D.; Filippini, L.; Gillies, M.; de Goede, J.; Mol, T.; Reijme, M.; Talen, W.; Wilderbeek, H.; Khatavkar, V.; Anderson, P. Lab Chip 2008, 8, 533−541. (4) Fahrni, F.; Prins, M. W. J.; van IJzendoorn, L. J. Lab Chip 2009, 9, 3413−3421. (5) van Oosten, C. L.; Bastiaansen, C. W. M.; Broer, D. J. Nat. Mater. 2009, 8, 677−682. (6) Vilfan, M.; Potocnik, A.; Kavcic, B.; Osterman, N.; Poberaj, I.; Vilfan, A.; Babic, D. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 1844− 1847. (7) Liu, F.; Ramachandran, D.; Urban, M. W. Adv. Funct. Mater. 2010, 20, 3163−3167. (8) Sanchez, T.; Welch, D.; Nicastro, D.; Dogic, Z. Science 2011, 333, 456−459. (9) (a) Aizenberg, J.; Zarzar, L. D.; Kim, P. Adv. Mater. 2011, 23, 1442−1446. (b) Sidorenko, A.; Krupenkin, T.; Aizenberg, J. J. Mater. Chem. 2008, 18, 3841−3846. (10) Tabata, O.; Hirasawa, H.; Aoki, S.; Yoshida, R.; Kokufuta, E. Sens. Actuators, A 2002, 95, 234−238. (11) Dayal, P.; Kuksenok, O.; Bhattacharya, A.; Balazs, A. C. J. Mater. Chem. 2012, 22, 241−250. (12) Gauger, E. M.; Downton, M. T.; Stark, H. Eur. Phys. J. E 2009, 28, 231. (13) Pokroy, B.; Epstein, A.; Persson-Gulda, M. C. M.; Aizenberg, J. Adv. Mater. 2009, 21, 463. (14) Harris, K. D.; Cuypers, R.; Scheibe, P.; van Oosten, C. L.; Bastiaansen, C. W. M.; Lub, J.; Broer, D. J. J. Mater. Chem. 2005, 15, 5043−5048. (15) Alexeev, A.; Yeomans, J. M.; Balazs, A. C. Langmuir 2008, 24, 12102−1210. (16) Timonen, J. V. I.; Johans, C.; Kontturi, K.; Walther, A.; Ikkala, O.; Ras, R. H. A. Appl. Mater. Interfaces 2010, 2, 2226−2230. (17) Ghosh, R.; Buxton, G. A.; Usta, O. B.; Balazs, A. C.; Alexeev, A. Langmuir 2010, 26, 2963. (18) Oh, K.; Chung, J.; Devasia, S.; Riley, J. J. Lab Chip 2009, 9, 1561−1566. (19) Maeda, S.; Hara, Y.; Yoshida, R.; Hashimoto, S. Angew. Chem., Int. Ed. 2008, 47, 6690. (20) Lahann, J.; Mitragotri, S.; Tran, T. N.; Kaido, H.; Sundaram, J.; Choi, I. S.; Hoffer, S.; Somorjai, G. A.; Langer, R. Science 2003, 299, 371. (21) Kim, Y. W.; Netz, R. R. Phys. Rev. Lett. 2006, 96, 158101. (22) Wiggins, C. H.; Goldstein, R. E. Phys. Rev. Lett. 1998, 80, 3879. 3226


Propulsion and Trapping of Microparticles by Active Cilia Arrays

Recommend Documents