Designing Synthetic, Pumping Cilia That Switch the Flow Direction in

Oct 11, 2008 - The George W. Woodruff School of Mechanical Engineering, Georgia ... Oxford, OX1 3NP, U.K., and Chemical Engineering Department, ...
7 downloads 0 Views 1MB Size
Download PDF
12102

Langmuir 2008, 24, 12102-12106

Designing Synthetic, Pumping Cilia That Switch the Flow Direction in Microchannels Alexander Alexeev,† J. M. Yeomans,‡ and Anna C. Balazs*,§ The George W. Woodruff School of Mechanical Engineering, Georgia Institute of Technology, Atlanta, Georgia 30332, The Rudolf Peierls Centre for Theoretical Physics, Oxford UniVersity, Oxford, OX1 3NP, U.K., and Chemical Engineering Department, UniVersity of Pittsburgh, Pittsburgh, PennsylVania 15261 ReceiVed June 17, 2008. ReVised Manuscript ReceiVed September 22, 2008 Using computational modeling, we simulate the 3D movement of actuated cilia in a fluid-filled microchannel. The cilia are modeled as deformable, elastic filaments, and the simulations capture the complex fluid-structure interactions among these filaments, the channel walls, and the surrounding solution. The cilia are tilted with respect to the surface and are actuated by a sinusoidal force that is applied at the free ends. We find that these cilia give rise to a unidirectional flow in the system and by simply altering the frequency of the applied force we can controllably switch the direction of the net flow. The findings indicate that beating, synthetic cilia could be harnessed to regulate the fluid streams in microfluidic devices.

In biological systems, microscopic cilia, which are hairlike protrusions that cover numerous cells, are highly effective at regulating the movement of fluids through confined geometries. The behavior of nodal cilia is particularly intriguing;1-5 these long filaments are slanted at an angle with respect to the underlying surface, and their rotating motion exhibits an asymmetry that is thought to be responsible for the unidirectional flow. Synthetic analogues of such nodal cilia could be very useful in microfluidics, allowing chemicals to be pumped in a specific direction within the microchannels. To test the potential utility of such synthetic cilia,6,7 we simulate the movement of elastic filaments actuated by a periodic external force, capturing the complex hydrodynamic interactions between the deforming filaments and the viscous liquid within a microchannel. Surprisingly, we find that by simply increasing the frequency of the driving force on the filament these cilia can actually alter the direction of the net fluid flow within the entire system. Thus, the synthetic cilia could be harnessed to switch the pumping direction controllably and repeatedly, providing both spatial and temporal control over the delivery of dissolved compounds to a region of a microfluidic device. There is rapidly growing interest in modeling the dynamics of the cilia-fluid interactions.2-5,8-20 For example, Kim and Netz8 * Corresponding author. E-mail: [email protected]. † Georgia Institute of Technology. ‡ Oxford University. § University of Pittsburgh. (1) Nonaka, S.; Yoshiba, S.; Watanabe, D.; Ikeuchi, S.; Goto, T.; Marshall, W. F.; Hamada, H. PLoS Biol. 2005, 3, 1467–1472. (2) Cartwright, J. H. E.; Piro, N.; Piro, O.; Tuval, I. J. R. Soc. Interface 2007, 4, 49–55. (3) Cartwright, J. H. E.; Piro, O.; Tuval, I. Proc. Natl. Acad. Sci. U.S.A. 2004, 101, 7234–7239. (4) Smith, D. J.; Blake, J. R.; Gaffney, E. A. J. R. Soc. Interface 2008, 5, 567–573. (5) Smith, D. J.; Gaffney, E. A.; Blake, J. R. Bull. Math. Biol. 2007, 69, 1477–1510. (6) Evans, B. A.; Shields, A. R.; Carroll, R. L.; Washburn, S.; Falvo, M. R.; Superfine, R. Nano Lett. 2007, 7, 1428–1434. (7) den Toonder, J.; Bos, F.; Broer, D.; Filippini, L.; Gillies, M.; de Goede, J.; Mol, T.; Reijme, M.; Talen, W.; Wilderbeek, H.; Khatavkar, V.; erson, P. Lab Chip 2008, 888, 533–541. (8) Kim, Y. W.; Netz, R. R. Phys. ReV. Lett. 2006, 96, 158101. (9) Khatavkar, V. V.; Anderson, P. D.; den Toonder, J. M. J.; Meijer, H. E. H. Phys. Fluids 2007, 19, 083605.

carried out 3D Brownian dynamics simulations to investigate the pumping efficiency of cilia that were modeled as chains of spherical beads. More recently, Khatavar et al.9 performed 2D simulations on an elastic filament to determine how artificial cilia could be exploited to promote mixing in microchannels. The 3D simulations carried out by Mitran10 were focused on capturing salient features of actual pulmonary cilia rather than determining how the material properties of synthetic cilia affect the collective behavior of the system. The computational studies performed by Gauger et al.20 could be very useful for guiding experimental efforts on magnetically actuated synthetic cilia6 because the researchers focused on modeling the response of a vertical string of beads to an asymmetrically rotating magnetic field. Because we model cilia as elastic filaments of finite width that are actuated by a simple applied force, our approach can provide more general guidelines for optimizing the performance of cilia fabricated from a variety of polymeric materials.21 To capture the complex fluid-structure interactions in our system, we use our hybrid LBM/LSM computational approach,22-24 which integrates the lattice Boltzmann model (LBM)25,26 for hydrodynamics and the lattice spring model (LSM)27-29 for the (10) Mitran, S. M. Comput. Struct. 2007, 85, 763–774. (11) Gauger, E.; Stark, H. Phys. ReV. E 2006, 74, 021907. (12) Guirao, B.; Joanny, J. F. Biophys. J. 2007, 921900-1917. (13) Gueron, S.; Levit-Gurevich, K.; Liron, N.; Blum, J. J. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 6001–6006. (14) Gueron, S.; Liron, N. Biophys. J. 1993, 65499-507. (15) Dillon, R. H.; Fauci, L. J.; Omoto, C.; Yang, X. Z. Reprod. Biomech. 2007, 1101, 494–505. (16) Dillon, R. H.; Fauci, L. J.; Yang, X. Z. Comp. Math. Appl. 2006, 52, 749–758. (17) Gebremichael, Y.; Ayton, G. S.; Voth, G. A. Biophys. J. 2006, 913640-3652. (18) Vilfan, A.; Julicher, F. Phys. ReV. Lett. 2006, 96, 058102. (19) Cosentino Lagomarsino, M.; Bassetti, B.; Jona, P. Eur. Phys. J. B 2002, 2681-88. (20) Gauger, E. M.; Downton, M. T.; Stark, H., arXiv:0805.3114v1. (21) Aksak, B.; Murphy, M. P.; Sitti, M. Langmuir 2007, 23, 3322–3332. (22) Alexeev, A.; Verberg, R.; Balazs, A. C. Macromolecules 2005, 38, 10244– 10260. (23) Alexeev, A.; Verberg, R.; Balazs, A. C. Phys. ReV. Lett. 2006, 96, 148103. (24) Alexeev, A.; Verberg, R.; Balazs, A. C. Langmuir 2007, 23, 983–987. (25) Succi, S. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond; Oxford University Press: Oxford, U.K., 2001; p xvi. (26) Ladd, A. J. C.; Verberg, R. J. Stat. Phys. 2001, 104, 1191–1251. (27) Ladd, A. J. C.; Kinney, J. H. Physica A 1997, 240, 349–360.

10.1021/la801907x CCC: $40.75  2008 American Chemical Society Published on Web 10/11/2008

Letters

Langmuir, Vol. 24, No. 21, 2008 12103

spring nodes that are situated at the solid-fluid interface impose their velocities on the surrounding fluids. In turn, these LS nodes experience forces due to the fluid pressure and viscous stresses at the interface; this force is applied as a load to the neighboring LS nodes. As in the case of nodal cilia, the cilium in our simulations is tilted at an angle R from the surface (Figures 1 and 2). The cilium has a square cross-section, with a side b equal to 4 LB units, and a length of L ) 10b.34 The cilium is firmly attached to the lower surface of our simulation box, which has dimensions of L × 0.5L in the horizontal directions and a height of H ) 1.5L. We impose no-slip and no-penetration conditions for the top and bottom walls. Because we impose periodic boundary conditions in the horizontal directions, we are effectively modeling an array of synchronously beating cilia within the fluid-filled channel. Figure 1. Schematic of a cilium attached to a rigid substrate with an angle R.

micromechanics of elastic solids. The LBM consists of two processes: the propagation of fluid “particles” to neighboring lattice sites and the subsequent collisions between particles when they reach a site. The system is characterized by a particle velocity distribution function, fi, describing the mass density of fluid particles propagating in the direction i with a velocity ci. The time evolution of the distribution function is governed by a discretized Boltzmann equation,25 and the hydrodynamic quantities are moments of the distribution function (i.e., the mass density F ) Σifi, the momentum density j ) Fu ) Σicifi, with u being the local fluid velocity, and the momentum flux Π ) Σicicifi). We set the fluid density equal to F ) 1 and the viscosity equal to µ ) νF ) 1/6. The compliant cilia are modeled via the LSM as a network of harmonic “springs”, which connect nearest- and next-nearestneighbor lattice nodes. To capture the dynamics of this solid material, we assign a mass Mi to each lattice node at position ri and integrate Newton’s equation of motion F(ri) ) Mi(d2ri/dt2), where F is the total force acting on the node. This total force consists of the force due to the interconnecting springs and the force exerted by the fluid on the solid at the solid-fluid boundary. We use the velocity Verlet algorithm to integrate the equation of motion. For the present study, we modified our method to capture the structure of the cilium, which is constructed from cubic elementary LSM units to form an elastic beam. For a simple cubic lattice, where springs of stiffness k connect the nearest- and next-nearestneighbor nodes, 30 the Young’s modulus of the solids is given by E ) 5k/2∆xLS, and the solid density is Fs ) M/∆x3LS, where ∆xLS ) 4/3 is the spacing between the LS nodes.29 We set the Young’s modulus E equal to 1 and the density Fs equal to 1 (in LB units). This simple lattice model is restricted to a Poisson’s ratio of 1/4, although more complicated many-body interactions can be included to vary the ratio.29,31-33 In our LBM/LSM simulations, the fluid and solid phases interact through appropriate boundary conditions.22,23 In particular, lattice (28) Ladd, A. J. C.; Kinney, J. H.; Breunig, T. M. Phys. ReV. E 1997, 55, 3271–3275. (29) Buxton, G. A.; Care, C. M.; Cleaver, D. J. Model. Simul. Mater. Sci. 2001, 9, 485–497. (30) We assign a value of k/2 to those springs that are located on the cilia surfaces and k/4 to the springs at the edges. We also set the LSM masses equal to m/2 and m/4 for nodes at the surfaces and edges, respectively. (31) Note that the Poisson’s ratio of polymers can vary from 0.3 to -0.7 and be tailored by blending.32 Thus, polymers blends with Poisson’s ratios of ∼0.25 can be synthesized.33 (32) Perez-Rigueiro, J.; Viney, C.; Llorca, J.; Elices, M. Polymer 2000, 41, 8433–8439. (33) Bureau, L.; Alliche, A.; Pilvin, P.; Pascal, S. Mater. Sci. Eng. A 2001, 308, 233–240.

A vertical, oscillatory driving force is applied to the cilium’s free end; this force follows a sinusoidal law with a frequency f and a force amplitude a and induces periodic oscillations of the cilium in its plane of inclination. Herein, we alter f and a to vary the oscillatory regimes of the beating cilia. Given that the cilium’s moment of inertia is I ) b4/12, the dimensionless force amplitude A ) 1/3aL2/EI represents the ratio between the amplitude of the driving force and the bending rigidity of a cilium. We consider relatively large amplitudes A g 1; consequently, for sufficiently slow oscillations where the viscous dumping is weak, the cilium touches the bottom wall when the force is directed downward and is stretched in the vertical direction when the force is directed upward (Figure 2). At these limiting positions, further increases in the force cause only slight changes in the cilium geometry. We find that the cilium geometry is strongly affected by the viscous forces from the surrounding fluid. These forces are proportional to the cilium velocity and thus depend on the oscillation frequency. The dimensionless sperm number Sp ) L(ς⊥ω/EI)0.25 characterizes the relative importance of the viscous force and the bending rigidity of the cilium EI; here, ς⊥ ) 4πµ, is the viscous drag coefficient and ω ) 2πf is the angular velocity of the cilium’s undulations. For larger Sp, the dynamics is dominated by the viscous forces, which can entirely suppress the cilium oscillations for sufficiently large Sp. For small Sp, however, the cilium follows the deformations due to the external force in a quasi-static manner (i.e., cilium geometry and position are completely defined by the magnitude of the force applied to the cilium at a given moment and does not depend on the fluid viscosity). Our aim is to identify the oscillatory regimes where the cilia can effectively pump the surrounding fluid. In a low-Reynoldsnumber environment, time-irreversible motion is required to induce a net fluid flow. Obviously, extremes in the sperm number are not favorable for inducing the pumping of fluid, with small Sp leading to periodic oscillations and large Sp completely suppressing the cilia’s motion. For intermediate Sp, however, the effects of the fluid viscosity and cilium bending rigidity are of comparable importance, and time-irreversible motion can emerge (Figure 2) that causes a net fluid flow in the channel. Our simulations indeed indicate that there is a net fluid flow in a range of Sp between 1.5 and 6 (Figure 3a).35 Here, we use a scale (34) Note that the cilia are constructed from a lattice with only four nodes across. To verify the accuracy of the model, we simulated the bending of a cilium due to a perpendicular force applied to its free end and found that an increase in lattice resolution effectively does not affect the deflection. Note also that for a 50% declination the difference between our LSM model and the linear theory prediction is about 10% and it increases to about 20% for a 100% declination (i.e., when declination is equal to the cilium length).

12104 Langmuir, Vol. 24, No. 21, 2008

Letters

Figure 2. Snapshots of cilium bending due to an oscillatory force applied to the cilium’s free end. The cilium is attached to the substrate with an angle of R ) 45°, the dimensionless force amplitude is A ) 10, and the sperm number is Sp ) 3. Snapshots a and c show two limiting positions of the cilium oscillation: stretched in the vertical direction and touching the substrate, respectively. Snapshots b and d show intermediate positions when the cilium moves in the downward and upward directions, respectively. The color shows the strain in the solid, and the arrows represent the direction and magnitude of the fluid velocity.

based on the fluid viscosity and the bending rigidity of the cilia, EIς⊥-1 L-3 ≡ ωLSp-4, to calculate the dimensionless magnitudes of the period-averaged flow rate, U, which are additionally divided by the dimensionless force amplitude A, so the magnitudes can be readily compared. Surprisingly, we found that the flow direction depends on the sperm number. Specifically, for lower Sp, the flow direction is positive (i.e., forward) and coincides with the direction of the cilium’s inclination, whereas for larger Sp the fluid flows in the opposite (backward) direction. We note that a similar effect was recently reported by Gauger et al.20 For A between 0.5 and 10, we found that both the forward and backward flows occur, although for the lower A the magnitude (35) A range of polymeric materials could be used to manufacture artificial cilia with these sperm numbers. As an example, let us consider a cilium that is 10 µm in length and 1 µm in diameter oscillating with a frequency of about 10 Hz in water, whose viscosity is µ ≈ 10-3 kg/m · s and density is F ≈ 103 kg/m3. For such cilia, a sperm number in the range from 1 to 10 can be obtained with polymers having a modulus of approximately 20 Pa-200 kPa, which is in the range of experimentally realistic values. These values yield a bending rigidity EI in the range of 10-24 to 10-20 N · m2 and a driving force amplitude of 0.1 pN < a < 1 nN. Alternatively, the relevant oscillation regimes can be excited by varying the driving frequency in the range of 10 Hz and 100 kHz.

of the backward flow is relatively weak (Figure 3a). Furthermore, we found that the fastest forward flow rate, U+, occurs for 2.5 e Sp e 3.5, whereas the fastest backward flow, U-, happens for 4.5 e Sp e 5. No net flow takes place at Sp ≈ 4. This distinctive behavior is solely governed by the interplay between the elastic forces of the deformed cilium and the viscous effects of the oscillating fluid (as indicated by the sperm number), whereas the magnitude of the applied force sets the ratio between the forward and backward flow rates. Figure 4a,b shows the oscillating patterns that cilia exhibit during one period and clearly reveals that the modes of motion for Sp ) 3 and 5 are different. Indeed, for larger Sp, the effect of fluid viscosity is more important and thus introduces a lag in the cilium motion. In fact, these shapes correspond to the first two modes of a flagella beating in a viscous fluid as described by Wiggins et al.36 These two oscillatory modes induce flows in opposite directions when a tilted cilium is attached to a substrate. (36) Wiggins, C. H.; Riveline, D.; Ott, A.; Goldstein, R. E. Biophys. J. 1998, 74, 1043–1060.

Letters

Figure 3. Net flow rate (a) and efficiency (b) as a function of the sperm number for a cilium tilted 45° with respect to the substrate.

To characterize the pumping performance of the cilia, we estimated the mechanical efficiency of pumping, , which is the ratio between the work required to pump a viscous fluid in the channel with the given flow rate and the work applied to drive the cilium. As seen in Figure 3b, the efficiency has two maxima that correspond to the fastest forward and backward flows (seen here most clearly for A ) 10). These maxima take place at Sp values that are somewhat lower that those leading to the maximum flow rates. Interestingly, not only the flow rate but also the pumping efficiency depends on the forcing amplitude. Typically, larger amplitudes result in better efficiency. The latter behavior is more pronounced for backward pumping, where a 10-fold increase in the amplitude (from 1 to 10) leads to 500× better efficiency. The increase is less dramatic for the forward flow, where the same change in the amplitude gives rise to approximately a 2.5-fold enhancement in the efficiency.

Langmuir, Vol. 24, No. 21, 2008 12105

We summarize our results for the maximum net flow rate and efficiency for the forward and backward flows in Figure 5a,b, respectively. Note that this data corresponds to slightly different Sp values, which lie in the ranges of 2.5 e Sp e 3.5 and 4.5 e Sp e 5 for the fastest forward and backward flows, respectively. The data is presented for three different cilium angles R and a range of forcing amplitudes A. For A < 3, there is relatively little difference between the angles for the values U+ and +, the respective maximum flow rate and efficiency of the forward flow, respectively. In this range of the amplitudes, U+ and + increase with increasing A and are roughly proportional to A1.5 for the flux, and to A for the efficiency. Nonetheless, the 45° angle provides up to 50% better efficiency than the other angles considered. For larger A, however, the three angles give significantly different magnitudes of the maximum forward flow rate and efficiency. We found that for all R, A ) 3 is close to a critical amplitude Acr for which the cilia can reach the limiting positions during an oscillation period for Sp ≈ 3, which corresponds to the fastest flow rate. Moreover, additional increases in the driving force do not change these limiting positions but reduce the part of the period that the cilia spend traveling between these positions. In other words, when oscillating with larger A, the cilia are “frozen” in the bounding positions for a part of the oscillation period. The backward pumping performance is also strongly coupled with the forcing amplitude. Both the maximum backward flow rate U- and efficiency - rapidly increase with increasing amplitude. For A > 6, the magnitude of the fastest backward flow becomes comparable (see Figure 5a for R ) 60°) or even exceeds (see Figure 5a for R ) 30 and 45°) that of the forward flow. In other words, for a sufficiently large but fixed amplitude, fluid flows of the same magnitude can be excited in either the forward or backward direction by changing the oscillation frequency by a factor of about 8 (i.e., changing the Sp between maximum forward and backward flow rates; recall that Sp ∼ ω0.25 and Sp are different for the fastest forward and backward flows by a factor of about 1.7). Furthermore, as seen in Figure 3a, the fluid velocity between these two extremes can be regulated solely by changing the driving frequency. As shown in Figure 5b, the pumping efficiency of the backward flow also increases rapidly with increasing amplitude, although

Figure 4. Cilia deformation during one oscillation period for A ) 10 and R ) 45°: (a) forward flow with Sp ) 3 and (b) backward flow with Sp ) 5. The blue lines show cilium centerlines when the force is directed upward, the green lines show the cilium when the force is downward, and the black dashed line is for the initial cilium position. Note that the oscillation of the cilium is asymmetric relative to its initial position as a result of the effect of the solid wall.

12106 Langmuir, Vol. 24, No. 21, 2008

Figure 5. Maximum net flow rate (a) and maximum pumping efficiency (b) are plotted as a function of the dimensionless force amplitude for different angles of cilium attachment to the substrate. The solid and empty symbols are for the forward and backward flows, respectively.

it is still relatively small compared to the forward flows with identical driving amplitudes. Indeed, backward flows take place for larger Sp (Figure 3) and therefore require a larger amount of work for the same flow rate when compared to the forward flows. We note that the forward flow is limited by the saturation effect when the cilia reach the limiting position. For the backward flow, however, the cilia did not reach these positions even for the largest amplitude considered in our study. Nonetheless, we found that the efficiency of the backward flow does decrease for A > 9. We relate the decrease in efficiency to the fact that for larger forces the curvature of the cilium becomes of the same order as the thickness b and this sufficiently reduces the flexibility

Letters

of the cilium. It appears that for the backward flow the saturation can be a result of a larger bending curvature and therefore can be regulated by choosing the appropriate thickness of the cilium. Another interesting effect observed in our simulations is that for A > 5 and Sp ≈ 4, which correspond to the transient regime between the forward and backward flows, the cilia exhibit a 3D out-of-plane motion. More specifically, the cilia produce circulatory motion. This 3D movement appears to arise from a buckling of the cilium due to the applied vertical force that develops even in the highly viscous environment. Indeed, the force applied to the cilium significantly exceeds the critical buckling force Fcr ) 0.25π2EI/L2, which gives a dimensionless amplitude of about 0.82. We did not observe, however, such a 3D motion for cilia oscillating with larger or smaller Sp. In summary, by adapting our LBM/LSM technique, we have formulated a new approach for simulating ciliated surfaces. This approach allows us to model the effect of the viscous forces on the moving, deformable filament and the concurrent effects of the filament on the fluid flow in the system. The cilium is tilted at an angle with respect to the surface, and thus, even though the actuating force has a simple sinusoidal form, the beating filament nonetheless leads to a complicated dynamic behavior within the system. In particular, the direction of the generated unidirectional flow can be switched by simply varying a single control parameter, the frequency of the applied driving force. By varying the tilt angle R, we could isolate the optimal geometry for maximizing the efficiency of the pumping process. The results provide guidelines for synthesizing microfluidic devices that do not require external pumps to move fluids through the microchannels. Consequently, the systems can be made more portable, expanding the applicability and utility of laboratory-on-a chip assemblies. Acknowledgment. We gratefully acknowledge helpful conversations with Halim Kusumaatmaja and Gareth Alexander. Funding from DOE (to A.C.B.) and ONR (to J.M.Y.) is also gratefully acknowledged. LA801907X