Isolating cells of interest from a heterogeneous population has been of critical importance in biological studies and clinical applications. In this study, a novel approach is proposed for utilizing an active ciliary system in microfluidic devices to separate particles based on their physical properties. In this approach, the bottom of the microchannel is covered with an equally spaced cilia array of various patterns which is actuated by an external stimuli. 3D simulations are carried out to study cilia-particle interaction and isolation dynamic in a microfluidic channel. It is observed that these elastic hair-like filaments can influence particle's trajectories differently depending on their biophysical properties. This modeling study utilizes immersed boundary method coupled with the lattice Boltzmann method. Soft particles and cilia are implemented through the spring connected network model and point-particle scheme, respectively. It is shown that cilia array with proper stimulation is able to continuously and non-destructively separate cells into subpopulations based on their size, shape, and stiffness. At the end, a design map for fabrication of a programmable microfluidic device capable of isolating various subpopulations of cells is developed. This biocompatible, label-free design can separate cells/soft microparticles with high throughput which can greatly complement existing separation technologies.

A large variety of biological organisms utilize the synchronized motion of beating cilia to capture food particles1 or prevent settlement of fouling agents onto their surfaces.2 In humans, cilia covering the upper generations of the lung airway are critical for expelling particulates out of the respiratory tract.3 Synthetic cilia are also being used to control fluid flow in microchannels for enhancing the performance of microfluidic devices. Specifically, artificial cilia is utilized in lab-on-chip devices for mixing and pumping purposes4 or to prevent the fouling of microfluidic devices.5 These ciliary systems are typically long, flexible grafted filaments which are anchored to the walls of microchannel and can be externally actuated by electromagnetic fields.4 

There are several methods to fabricate/actuate artificial cilia in micro-scale. For instance, double layer polyimide (PI) and chromium (Cr) microbeam can be electrostatically actuated. By applying a voltage difference between two layers, it can operate at 200 Hz and subsequently generate local mean fluid velocities up to 500 μm/s.6 Furthermore, magnetic-actuated cilia made from polydimethylsiloxane (PDMS) contain superparamagnetic nanoparticles (Fe3O4) which can create a net fluid flow of about 8 μm/s at frequencies up to 35 Hz.7 They can be fabricated by dissolving polycarbonate track etch (PCTE) filter membrane which is used as mold. There are also hydrogel-actuated8 and resonance-actuated9 artificial cilia which are fabricated with soft lithography approach. In the former case, the change in pH level results in a swelling or shrinking of the hydrogel and in the latter case, cilia oscillates in resonance with piezoelectric.

Theoretical and computational studies on cilia–microparticle interactions10–13 mostly focus on influence of active/passive harmonic motion of synthetic cilia on microparticles motion.11,14 It is shown that depending on the frequency of beating cilia, the microparticles can be either driven downward toward the substrate or upward into the fluid above the cilia layer.10 For instance, the hydrodynamic forces generated by active motion of ciliary systems favour capturing microparticles within the grafted filaments.15 It is also observed that even passive motions of synthetic cilia swayed by the ambient flow field can enhance the deposition16 or repel away the adhesive particles.17 Overall, particle trajectories converge to one of the released, propelled, or propelled states.18 It should be noted that these studies are based on cilia with sticky tips to regulate the motion of microscopic particles. Such adhesive tip provides desirable downward force for trapping particles.

Isolating cells/particles of interest from a heterogeneous population is critical in biological studies and clinical applications such as single cell sequencing, capturing circulating cancer cells, and rare cell isolation.19–21 As a result, a variety of methods has been developed to separate particles based on their size, shape, and stiffness Different physical mechanisms are employed to manipulate cells, such as magnetic fields,22,23 electric fields,24,25 optical forces,26,27 and acoustic fields.28,29 For instance, acoustic forces have been recently explored for alignment, separation, and enrichment of particles and cells.30 Specifically, standing surface acoustic waves are being used for continuous separation of particles based on their volume, density, and compressibility.29,31 Alternatively, labelling cells through specific binding of fluorescent antibodies are routinely used to separate cells in devices such as fluorescence activated cell sorters (FACS).32 This separation method is expensive and cell integrity during the cell-sorting process may be compromised. Other mechanisms such as hydrodynamic flows33 and valve-based switching34 can also be utilized to separate specific subpopulations of cells. For instance, Huang et al. proposed a continuous particle separation method through deterministic lateral displacement which use asymmetric bifurcation of laminar flow around obstacles.35 In a similar study, Takagi et al. use asymmetrically arranged multiple branches to continuously separate particles in a microchannel.36 It should be noted that these studies utilized size difference to perform isolation of target subpopulation through hydrodynamic filtration.37 Other biophysical characteristics were also being targeted. For instance, Seo et al. did membrane-free microfiltration study by asymmetric inertial migration.38 Furthermore, there are also studies which solely focus on stiffness based separation of cells in a microfluidic device.19 A very extensive review of various separation method based on size, shape, and stiffness can be found in a study by Pamme.39 Generally, low throughput, bulky instrument, and low biocompatibility can be identified as drawbacks of discussed separation methods. Therefore, a biocompatible label-free method that can separate cells with high throughput would greatly complement existing separation technologies.

Several studies have shown a reduction in cell stiffness with increasing metastatic efficiency in human cancer cell lines.40 Thus, mechanical stiffness can be used to identify abnormal cell populations in detecting cancer as well as identifying infectious disease.19 As discussed before, microfluidic methods were developed to classify and enrich cell populations utilizing mechanical stiffness.41,42 It is very challenging to isolate rare cells from a mixture when their mechanical stiffness slightly differs from regular cells. However, it is relatively easier to separate microparticles based on size and shape as proven by many hydrodynamic filtration methods.39 Inspired by cilia-microparticle interactions studies, we aim to harness the specific features of the active ciliary systems to propose a novel approach for sorting/separating cells based on size, shape, and stiffness. The hydrodynamic forces induced by beating cilia can regulate the trajectories of the suspended cells/particles based on their biophysical properties. To the best of our knowledge, the possibility of using active coordinated motion of cilia arrays to separate cells has not been extensively studied yet.

In our modelling study, the bottom wall of microfluidic channel is decorated by an equally spaced cilia assay. The spring connected network (SN) model is coupled with lattice Boltzmann through immersed boundary (IB) method to study hydrodynamic interactions of active ciliary array with soft particulates. In what follows, the model formulations are first discussed. Then through 3D simulations, we show how particles with different physical properties move at different heights above the ciliated wall. Because of that, soft microparticles may follow different lateral trajectories and therefore can be sorted along the width of the channel. Our approach is able to continuously and non-destructively separate cells into subpopulations with high throughput. We use our simulation results to find proper model parameters to efficiently separate cells with different size, shape, and mechanical stiffness. Ultimately, a design map for fabrication of a programmable microfluidic device capable of isolating various subpopulations of cells is provided.

The mesoscale Lattice-Boltzmann Method (LBM)43 is considered as an efficient solver for the Navier–Stokes equation due to its particulate nature and parallel computing capabilities.44–48 LBM is typically considered as a second order accurate method in space and time,49 and the local density distribution of fictitious particles of LBM is represented by gix,t. The dynamics of the LBM involves streaming and collision processes in phase space x,t.

(1)

where Ωi is the collision operator, ci denotes the lattice velocity, and Fi is the force term. Bhatnagar-Gross-Krook (BGK) collision scheme50 is used here

(2)

where τ is a relaxation constant. The population distribution at equilibrium, gieqx,t, can be calculated from the local macroscale fluid velocity

(3)

where wi is the weight coefficient. It should be noted that the existence of external bodies in the flow are introduced through a force term, Fi, expressed in terms of force density f and fluid velocity u

(4)

wheref is derived from solid models which is distributed on lattice nodes via Immersed Boundary method (IBM). D3Q19 lattice model are used in the current study.

FIG. 1.

(a) Spring connected network cell membrane model. (b) Kinematics for local stretching and bending modes of response.

FIG. 1.

(a) Spring connected network cell membrane model. (b) Kinematics for local stretching and bending modes of response.

Close modal

Spring connected network (SN) which is basically a particle based cell model is adopted due to its simplicity in mathematical description.51–53 The cell membrane consists of a network of vertices Xi,i1NV connected with nonlinear springs forming two-dimensional triangulated network as shown in Fig. 1. The potential energy of the membrane is defined as

(5)

The in-plane energy term characterizes the elastic energy stored in the membrane which includes attractive nonlinear wormlike chain (WLC) and repulsive power function (POW) potentials. It can be defined as

(6)

where x=l/lm, lm is the maximum spring extension, p is the persistence length, kBT is the energy per unit mass, and kp is the POW force coefficient. Bending potential corresponds to bending stiffness of erythrocyte which can be expressed as

(7)

where kb is the bending constant, θj is the instantaneous angle between two adjacent triangles having the common edge j, and θ0 is the spontaneous angle. Varea stands for local and global area stiffness of lipid bilayer which is defined as

(8)

where ka and kd are the global and local area constraint constants, respectively. Similarly, the last potential term conserves red blood cell (RBC) volume

(9)

where kv is volume constraint constant. The nodal forces corresponding to the each energy can be calculated as fi=Vxi/xi.

In this study, Immersed Boundary method (IBM) is used to transfer force density and velocity distribution data back and forth between LBM and solid models.54,55 In IBM, the fluid is solved on a Eulerian grid which does not conform to the solid mesh. Therefore, to couple fluid and immersed solid particle, the nodal forces at half time step are calculated using each energy terms as a parametric surface X(s,t) and then, the immersed boundary method spread this force density onto the fluid lattice to represent the effect of immersed solid body. To update the new position of immersed solid particles in time, similarly, the solid nodal velocity will be interpolated from the local fluid nodes using immersed boundary method. In other word, interpolated velocity field will be used to update the new position of immersed mesh.56 They can be written as

(10)

where Fs,t is the structural force at location s at time of t and ds is the discretized length of the immersed structure. The force spreads to the local fluid nodes and velocity is interpolated from the local fluid nodes through a delta function δr

(11)

where r=xXs,t, for instance, is the distance between the structural force/velocity at location s and lattice node.

Cilia diameter is negligible compared to its height. Therefore, Point-Particle LBM (PPM) scheme can be used to model the micromechanics of these elastic filaments.57 Via a network of harmonic springs, PPM characterizes an elastic body formed by a chain of regularly spaced beads where the first bead in the chain (i=0) is tethered to the wall (y=0). The elastic energy associated with the chain is11 

(12)

where θi is the angle between neighboring bonds and kb represents the cilium bending modulus. The instantaneous angle can be expressed as θi=cos1b̂i+1,i.b̂i1,i, where b̂i,j=rirj/rirj, with b̂0,1=j being the wall normal direction. Moreover, kb equals πErc4/4l where E is Young's modulus, l is cilium length, and rc is cilium radius. To use harmonic spring representation of elastic cilia, the nodal forces corresponding to the spring elastic energy should be calculated accordingly. There are also other forces applied on these cilia nodes. Specifically, the actuation force should be taken into account.

In biological cilia, the beating starts from the bottom to the tip. In biomimetic cilia, the actuation force can be originated from either a body force acting on whole cilium or local force at the tip of cilium where concentration of magnetic particle is highest. The difference is the small variation on cyclic motion of cilia which is not the focus of this study. Therefore, for simplicity, the cyclic actuating force is applied to the three points at the extreme tip of each cilium. Thus, a periodic force with functional form of Fextt=Fx0cosωti+Fz0sinωt is applied at the cilia tip, where ω is the frequency of actuation. The beating pattern is symmetric and the actuation force oscillates harmonically in XZ plane.

The third force on cilia nodes is hydrodynamic coupling force. An equal and opposite hydrodynamic coupling force is also interpolated onto the fluid lattice to conserve local and global momentum in the system. It should be noted that Lattice Boltzmann fluid model senses the presence of cilia through hydrodynamic coupling forces which are directly spread onto fluid lattice node using immersed boundary method. The hydrodynamic coupling between the fluid and the cilia is captured through a frictional force that is proportional to the slip velocity between the beads and fluid. Namely, the force on the i th cilium bead is given by FiH=6πηrcṙiu(ri) where ṙi is the velocity of ith particle and u(ri) is the interpolated fluid velocity at position ri of the particle i. Again, the interpolation of fluid lattice velocity is carried out by immersed boundary method.

Finally, these three forces (elastic, actuation, and hydrodynamic coupling forces) are input into Velocity Verlet algorithm to update the position and velocity of cilium beads to be used in the next time step.

To benchmark 3D soft particle model, we have followed optical tweezer experiment by Mills et al.58 and characterized the deformation of RBC undergoing a dynamic load of 0200pN as shown in Fig. 2(a). The bending stiffness, local area constraint, global area constraint, and volume constraint constants are chosen as 2.77×1018kgm2s2, 23.1×105kgs2, 4.72×106kgs2 and 249kgm1s2, respectively. The stretching test is performed when pair forces are applied at 2% of mesh nodes at two ends of a RBC. Our stretching test results for axial and transverse diameters of the deformed RBC agree well with experiments58 and other numerical studies.59,60

FIG. 2.

(a) The axial and transverse diameters of RBC deformed under stretching test. (b) Oscillation period over different shear rate for a RBC in shear flow, comparison of our model, and experiment.

FIG. 2.

(a) The axial and transverse diameters of RBC deformed under stretching test. (b) Oscillation period over different shear rate for a RBC in shear flow, comparison of our model, and experiment.

Close modal

To present test cases where the particle deformation is coupled with the fluid flow, the deformation of RBC under pure shear flow is investigated. In flow with low-shear rate, cells tumble while preserving their biconcave shape. As the fluid shear rate increases, the tumbling gradually reduces and RBCs deform into ellipsoidal shape. The oscillation period in low shear rates is plotted out in Fig. 2(b) where a good agreement with experimental study by Abkarian et al.61 is achieved.

To simulate cilia-flow interaction, we followed Bhattacharya et al.'s works11,12 where they used Point-Particle Method (PPM) to model micromechanics of these elastic filaments. They used this tool to study the motion of soft microscopic particles over active adhesive cilia. To ensure cilia can properly reproduce drag coefficient for flow normal to the filaments, at different flow rates (Q), the deflection of unactuated cilia is studied. Basically, cilium tip deflection is recorded as a function of imposed flow rate (Q) in a confined geometry. We validated our model with analytical solution proposed by Wexler et al.62 They developed a mathematical model showing a linear relationship between deflection of confined fibres and flow rate at low flow rates. The slope of flow rate-tip deflection profile can be calculated as

(13)

where H and D are is channel height and width, respectively; E is Young's modulus; d and I are cilia diameter and moment of inertia, respectively; and u1 is a function of c=h/H. In their experimental study, the longest rectangular shape cilia used had the dimension of 22μm×56μm×241μm, where the length to diameter aspect ratio is 6.

To benchmark our cilia model, we examined the flow of fluid past an elastic fibre in a channel with dimensions of Lx=80μm,Ly=10μm,andLz=40μm. It is observed that our cilia model agrees well with analytical solution at moderate flow rates as shown in Fig. 3. However, at low flow rates, our model predicts a nonlinear relation between cilia tip deflection and flow rate. Discrepancy is due to cilium aspect ratio of 20 used in our simulation as opposed to highest aspect ratio of 6 used in Ref. 62. Since PPM directly spread force density on lattice nodes through immersed boundary, it is only capable of simulating high aspect ratio cilia array. However, there is another method for modelling cilia, lattice spring model (LSM), which does not have this limitation. In LSM, no-slip conditions are explicitly imposed at the surface of the cilium. The reason we have chosen PPM over LSM is its considerably lower computation cost which is more desirable in our computationally intensive modelling study. For more detail information, readers can refer to the study of Bhattacharya et al.11,12 where PPM and LSM models have been carefully benchmarked and validated.

FIG. 3.

Cilia tip deflection, uy, versus flow rate, Q, for H=40μm, D=10μm, h=13.3μm, E=0.132MPa, d=0.66μm, and u113.3/40=6.1× 10−4.

FIG. 3.

Cilia tip deflection, uy, versus flow rate, Q, for H=40μm, D=10μm, h=13.3μm, E=0.132MPa, d=0.66μm, and u113.3/40=6.1× 10−4.

Close modal

In a straight channel, the synchronized motion of cilia actively changes the axial components of velocity near ciliated region. These beating hair-like filaments can switch the directionality of the net flow63 and even draw or repel away flow particles to/from ciliated region. In these flows, the sperm number quantifies the relative importance of the viscous force compared to the bending rigidity of oscillating cilia. The sperm number can be expressed as

(14)

where ζ=4πυρ is the viscous drag coefficient, ω=2πf is the angular velocity of the driving force, and EI is the bending rigidity of cilium. For relatively small Sp, the dynamic of cilium is governed by its elasticity and cilia array is unable to generate net flows at low Re.64 For high Sp numbers, no net fluid flow is generated due dominant effect of fluid viscosity. Once the cilium elasticity and fluid viscosity forces are of comparable magnitudes, coupling between elastic filaments and viscous fluid can create net flow.

The respective dimensions of computational box are Lx×Ly×Lz=80×50×50 lattice units. The kinematic viscosity in lattice unit is given by υLB=cs2τ12ΔtLB.65 The acceptable range for relaxation time is between 0.5 and 2. BGK LBM accuracy deteriorates as τ grows by introducing an error proportional to τ0.52. In many studies, depending on the Reynolds number, a relaxation time τ is typically recommended between 0.8 and 1.0.66,67 Due to small Reynolds number, 0.2, a relaxation time τ of 1.0 was used for all simulations for computational efficiency and larger time steps.68 Also, a simple case with τ=0.6 was simulated to see possible effect. However, no difference in particle's trajectories are observed between cases with τ=1 and τ=0.6. We also assume ΔtLB=1, ΔxLB=1, and cs=1/3. Therefore, the value of viscosity in lattice units can be derived as 0.16. Zou/He rule69 is applied at upper and lower walls to enforce non-slip wall boundary condition. Since the width of real channel is many times larger than its height, we assume that the effect of side walls on particle migration can be negligible. Therefore, to simplify the model, periodic boundary condition is used in Y directions. Moreover, a body force replaces necessary pressure gradient to create main flow. Thus, we further simplified model and treated inlet and outlet of the channel as periodic boundary condition. It should be noted that Mach number for all cases are kept below 0.05 in order to avoid deleterious compressibility effects. The length of cilium and channel height are l=13μm and L=40μm, respectively. Cilia spacing along channel length is 25 μm where rows are tilted 30° as shown in Fig. 4(a). The radius of each point particle has to be much less than 1LB unit for the model to be valid, thus rcilium/l=0.016.11 In the PPM simulations, the bond length between the point particles needs to be approximately 1LB unit, hence Ncilium=20. Bhattacharya et al.11,12 extensively benchmarked Particle Point Method (PPM) showing that with bond length of approximately 1 LB unit, best results with minor error can be achieved. The leakage is not a concern in PPM because cilium is not an enclosed boundary and the diameter of cilia is much smaller than lattice unit.

FIG. 4.

(a) 3D demonstration of particle-cilia interaction in microchannel. (b) Particles with distinct biophysical properties tend to flow downstream in different height of microchannel.

FIG. 4.

(a) 3D demonstration of particle-cilia interaction in microchannel. (b) Particles with distinct biophysical properties tend to flow downstream in different height of microchannel.

Close modal

Circulating tumor cells (CTC) are relatively stiffer and their deformability varies a lot.70 As cells become malignant, their stiffness decreases significantly.71 Various values of bending stiffness are proposed for blood cells in literature.52,72,73 To cover a target range of deformability, bending stiffness for soft and stiff particles are chosen as 2×1017 and 2×1019kgm2s2, respectively. Furthermore, to cover different sizes and shapes, 8μm RBC and 4μm, 8μm, and 12μm circular particles are chosen to study separation dynamic as shown in Fig. 4(a). 320, 1280, and 5120 triangular elements are used to mesh 4, 8, and 12 μm particle's membrane, respectively. With lattice spacing of 0.85μm, the ratios of solid to fluid spacing for 4, 8, and 12 μm particles are 0.64, 0.64, and 0.53, respectively. Although to completely avoid leakage, this ratio should be less than 0.5,74–77 others have shown that for ratios between 0.5 and 1, the leakage is small.78 To study the independence of the solution on this ratio, we carried out a simulation for a sphere settling in viscous fluid and found the difference between the theoretical and simulation-based terminal velocity prediction is within 2%, which indicates that the FSI code correctly reproduces the kinematics of the sphere in a viscous fluid.

The local area constraint, global area constraint, and volume constraint constants are chosen as 23.1×105kgs2, 4.72×106kgs2, and 249kgm1s2, respectively. To eliminate any residue membrane stress at equilibrium state, nodes are allowed to move on the membrane to minimize in-plane energy before starting the simulation. In this stress-free condition, the membrane is relaxed and equilibrium lengths of springs are recalculated. More detail information about spring network model can be found in Refs. 56 and 73. To start the simulation, cilia are actuated with 120Hz frequency without particles at first, and after reaching equilibrium, particles are added to remove the effect of initial conditions on particle transient motion. The dimensionless values for other relevant parameters used in the simulations are Re=0.2, Sp=6, Fx0l2/EI=64, and Fz0/Fx0=4. The Reynolds number's length-scale is based on cilia length scale. To give readers a sense of regimes involved, two other non-dimensional numbers is also provided. The Stokes numbers is fa2/ν=0.0012, where f is cilia beating frequency. Furthermore, the shear Stokes number is computed as ka2/ν=0.0076 where k is the typical shear rate in microchannel.

The biophysical properties of flow particles directly influence the dynamics of their interaction with ciliary system. Because cilia array is positioned in an angle respect to main flow, the position of particles along Z axis are stabilized at different height of channel as shown in Fig. 4(b). It is observed that bigger particles tend to be driven downward toward ciliated surface, while smaller particle are more effectively repelled away by harmonic motion of cilia. Similarly, particles with different stiffness and shapes also flow downstream at different heights of microchannel as shown in Fig. 4(b).

Similar to Herringbone structures,79 there will be a lateral drift in Y direction when unactuated cilia array are placed with an angle to the main flow. This lateral drift above ciliated region is due to cilia tilted configuration. Although in actuated cilia cases the flow will no longer be steady, the lateral drift in Y direction still exists. Particles above ciliated region get trapped in this lateral stream as they flow downstream. Utilizing this lateral drifting motion, we aim to separate/sort flow particles along the width of the channel. As shown in Fig. 5 and supplementary material Video 1, harmonic motion of tilted cilia array laterally displaces particles in Y direction. The magnitude of lateral drifting velocity is considerably smaller than the main flow. To provide a better insight into underlying physics of flow dynamics, the streamlines during one actuation cycle are demonstrated in supplementary material Video 2. It is observed in Fig. 5 that 12μm particles tend to laterally drift more compared to 8μm and 4μm particles. Also, the RBC is shown to displace even more than 12μm spherical particles by 10% because of its biconcave shape. The lateral drifting velocity also depends on the membrane bending stiffness as shown in Fig. 5. For instance, stiff 12μm particles drift 6% more compared to their softer counterpart with same biophysical properties. It is hypothesized that the difference in particle's drifting velocities are mainly due to their stabilized position along the Z axis.

FIG. 5.

Size, shape, and stiffness dependent trajectories of particles due to harmonic motion of tilted cilia array. As particles flow downstream in the X direction, they also displace laterally in the Y direction. After transition period, particle's position in Z direction does not change anymore. The yellow dashed and dotted black lines represent the case with soft and stiff 12 μm particle where cilia is not actuated. Separation efficiency is defined as the ratio of the separation distance between two trajectories to the length particles move downstream, Se=Δz/Δx.

FIG. 5.

Size, shape, and stiffness dependent trajectories of particles due to harmonic motion of tilted cilia array. As particles flow downstream in the X direction, they also displace laterally in the Y direction. After transition period, particle's position in Z direction does not change anymore. The yellow dashed and dotted black lines represent the case with soft and stiff 12 μm particle where cilia is not actuated. Separation efficiency is defined as the ratio of the separation distance between two trajectories to the length particles move downstream, Se=Δz/Δx.

Close modal

To check if the sorting process is not based on main flow and passive presence of cilia, a case with 12 μm particles in a channel with unactuated ciliary system is also simulated as shown in Fig. 5. Similar to Herringbone structures, unactuated tilted ciliary system can also laterally displace particles but cannot effectively separate particles. It is observed that separation efficiency of soft and stiff particles in unactuated system is 10 times smaller than actuated ciliary system. The separation efficiency can be defined as the ratio of the separation distance between two trajectories to the length particles move downstream, Se=Δz/Δx. Overall, using this setup, target particles with specific biophysical property can be isolated from a heterogeneous population because they follow unique lateral trajectories as shown in supplementary material Video 3. The dynamics of interaction between cilia and microparticles can be better realized by watching three multimedia files provided as electronic supplementary material.

It should be noted that the dynamics of cilia-particle interaction can be greatly influenced by some model parameters such as cilia array angle, the actuation frequency, and cilia length. The cilia array angle presumably is the most important parameter influencing the particle's lateral velocity. In this study, 20°, 30°, 45°, and 60° angles are studied to better understand the effect of this parameter on separation dynamics as shown in Fig. 6. Largest lateral displacement can be achieved when cilia array are placed in 45° angle. Furthermore, trajectories of particles with different sizes are more effectively separated when cilia array angle is 45°. However, 30° angle is shown to have higher separation efficiency in the case of particles with different bending stiffness.

FIG. 6.

The effect of ciliary angle on lateral displacements and separation efficiencies of soft particulates.

FIG. 6.

The effect of ciliary angle on lateral displacements and separation efficiencies of soft particulates.

Close modal

The actuation mode of external stimuli regulates beating frequency of ciliary system which also directly affect the lateral displacement of flow particles as shown in Fig. 7. It is observed that the higher the beating frequency, the bigger the lateral displacement of suspended particulates is. However, it is shown that 120 Hz leads to best separation between soft and stiff particles.

FIG. 7.

The effect of actuation frequency on lateral displacement of soft and stiff 12 μm particles.

FIG. 7.

The effect of actuation frequency on lateral displacement of soft and stiff 12 μm particles.

Close modal

The ratio of cilia length to channel height also is one other influential factor. It is demonstrated in Fig. 8 that the ciliary system with higher aspect ratio induce more lateral force on flow particulates. Similar to cilia array angle, the separation efficiency based on size, shape, and bending stiffness changes as cilia lengths increase or decrease in constant channel height. As observed in Fig. 8, the size based separation efficiency is 8.6% for cilia array with 0.33 aspect ratio while it is 6.5% for 0.4 aspect ratio. However, the ciliary system with 0.4 aspect ratio is found to be more capable of separating particles based on their bending stiffness.

FIG. 8.

The effect of ciliary aspect ratio on lateral displacement and separation efficiency of soft 8 μm and 12 μm particles. The zoom-in images demonstrate the oscillations in particle's trajectory because of synchronized motion of cilia.

FIG. 8.

The effect of ciliary aspect ratio on lateral displacement and separation efficiency of soft 8 μm and 12 μm particles. The zoom-in images demonstrate the oscillations in particle's trajectory because of synchronized motion of cilia.

Close modal

It is shown that soft microparticles follow unique trajectories based on their specific biophysical properties when flowing downstream in a microfluidic channel decorated with active synthetic ciliary system. In this section, we use our findings to propose a structural design for a ciliary system which can achieve high separation efficiency. As discussed before, the separation efficiency based on size, shape, and stiffness strongly depend on cilia array angle. This feature can be utilized to design an efficient ciliary pattern to effectively isolate target particles from a heterogeneous mixture.

For ciliary system with 30° angle, 0.33 length ratio, and 120Hz actuation frequency, it is observed that microparticles laterally displace 100μm as they move 10mm downstream (see Fig. 5). Subsequently, their trajectories will be separated by 2% to 15% which can ultimately lead to 2 to 15μm gap between their mass centres. We may be able to use longer channels to gain more separation gap, but it is necessary to use very wide channels. To address this issue, we can use ciliary arrays in blocks with different angles. Thus, by placing a 0.5mm block with 45° cilia angle after 1mm block with 30° cilia angle, particles can be re-directed in opposite direction to their original position and at the same time, we can exploit the fact that the separation efficiencies varies for patterns with different ciliary angles. This process can be repeated multiple times. The schematic demonstration of proposed design is shown in Fig. 9(a). Finally, by adding an expansion nozzle at downstream, we can further amplify separation gap between particle's trajectories.

FIG. 9.

(a) The schematics of a high-throughput, optimized particle sorter for efficient isolation of target particles with specific biophysical property. (b) Separation of soft and stiff 12μm in the proposed design for particle sorter.

FIG. 9.

(a) The schematics of a high-throughput, optimized particle sorter for efficient isolation of target particles with specific biophysical property. (b) Separation of soft and stiff 12μm in the proposed design for particle sorter.

Close modal

Even a slight difference in stiffness can be detected by active motion of tilted cilia array. Cancer cells are relatively stiffer compared to normal healthy cells.71 Assuming a separation efficiency as low as 0.05% between trajectories of cancer and healthy cells, they can be potentially distanced by 100μm at the end of expansion nozzle in a 20 cm microfluidic channel. As an example, we also carried out the simulation to demonstrate the separation process of stiff and soft 12μm particles as shown in in Fig. 9(b). It should be noted that the length of 30° and 45° blocks are 1mm and 0.62mm, respectively. It is observed that soft and stiff particles will be separated by 50μm as they move downstream of a 20mm channel. It should also be noted that these two particles move within 200μm width of ciliated region. Thus, utilizing active ciliary system with certain patterns, we can separate particles in just one step with high throughput as shown in Fig. 9(a). Assuming 15 mm/s average velocity of particles and 50 μm releasing interval, we estimate that this method can separate 300 particles per second.

Various methods have been suggested for fabrication of magnetically actuated cilia embedded in a closed micro channel.80–87 Among them, PCTE porous membrane,85 micro-holes drilled into acrylic block,84 silicon wafers prepared by photolithography method,83 or candy-based replica86 can be used as sacrificial or permanent mold to fabricate ciliary systems. In all these approaches, magnetic particles in solution are inserted into micro-structures first. Then, after pouring PDMS, degassing, and curing, the mold either get dissolved or removed. Therefore, fabrication of micro-channels based on this design seems practical. Additionally, the microfluidic device can be designed in a way that is suitable for different isolation scheme. As discussed before, the beating frequency also alter the lateral displacement as well as the size, shape, and stiffness based separation efficiency, refer to Fig. 7. Thus, one pattern of cilia array can be used to isolate various type of target particles just by tuning the actuation frequency. The scope of this work is focused on the conceptual design and numerical simulation of ciliary system based label-free isolation. The fabrication and experiment tests will be explored in our future work.

Simple and efficient particle separation methods are fundamentally important in biological and chemical analyses such as cancer cell detection, rare cell isolation, tissue engineering, and drug screening. In this study, a novel approach is proposed to isolate a specific subpopulation of microparticles using an active ciliary system. Our model utilizes immersed boundary (IB) method coupled with lattice Boltzmann. Soft particle and cilia are implemented through a spring connected network model and a point-particle scheme, respectively. It should be noted that unlike previous studies, we do not assume adhesive interaction between the particles and cilia. Therefore, antifouling characteristics of cilia array prevent any direct contact between cilia and sensitive particulates.

3D simulations are carried out to closely study cilia-particle interaction and isolation dynamic. It is observed that soft microparticles flow downstream in different trajectories based on their biophysical properties. It is hypothesized that depending on their closeness to beating cilia, they will experience different lateral forces. The goal of this study is to demonstrate that our approach is capable of continuously and non-destructively separating cells into subpopulations based on their size, shape and stiffness. At the end, proper structural parameters for an active ciliary system capable of isolating microparticles with high separation efficiency are provided. We have shown that our design is very sensitive to changes in stiffness. Overall, our programmable/high throughput microfluidic device is biocompatible which can greatly complement existing separation technologies.

See supplementary material for three multimedia files visualizing the dynamics of interaction between cilia and microparticles.

This work was partially supported by the National Science Foundation (Grant Nos. NSF 429IIP-1701136 and DMS-1516236), the National Institutes of Health (NIH) (Grant No. R01HL131750), and the Pennsylvania Infrastructure Technology Alliance (PITA) program.

1.
J.
Shimeta
and
M.
Koehl
, “
Mechanisms of particle selection by tentaculate suspension feeders during encounter, retention, and handling
,”
J. Exp. Mar. Biol. Ecol.
209
(
1
),
47
73
(
1997
).
2.
R.
Fox
,
Invertebrate Zoology
(
Wadsworth Publishing
,
2001
).
3.
Y.
Enuka
 et al., “
Epithelial sodium channels (ENaC) are uniformly distributed on motile cilia in the oviduct and the respiratory airways
,”
Histochem. Cell Biol.
137
(
3
),
339
353
(
2012
).
4.
J. M.
den Toonder
and
P. R.
Onck
, “
Microfluidic manipulation with artificial/bioinspired cilia
,”
Trends Biotechnol.
31
(
2
),
85
91
(
2013
).
5.
R.
Mukhopadhyay
, “
When microfluidic devices go bad
,”
Anal. Chem.
77
(
21
),
429A
432A
(
2005
).
6.
J.
den Toonder
 et al., “
Artificial cilia for active micro-fluidic mixing
,”
Lab Chip
8
(
4
),
533
541
(
2008
).
7.
A.
Shields
 et al., “
Biomimetic cilia arrays generate simultaneous pumping and mixing regimes
,”
Proc. Natl. Acad. Sci. U.S.A.
107
(
36
),
15670
15675
(
2010
).
8.
B.
Pokroy
 et al., “
Fabrication of bioinspired actuated nanostructures with arbitrary geometry and stiffness
,”
Adv. Mater.
21
(
4
),
463
469
(
2009
).
9.
K.
Oh
 et al., “
Characterization of mixing performance for bio-mimetic silicone cilia
,”
Microfluid. Nanofluid.
9
(
4–5
),
645
655
(
2010
).
10.
R.
Ghosh
 et al., “
Designing oscillating cilia that capture or release microscopic particles
,”
Langmuir
26
(
4
),
2963
2968
(
2010
).
11.
A.
Bhattacharya
 et al., “
Propulsion and trapping of microparticles by active cilia arrays
,”
Langmuir
28
(
6
),
3217
3226
(
2012
).
12.
A.
Bhattacharya
and
A. C.
Balazs
, “
Stiffness-modulated motion of soft microscopic particles over active adhesive cilia
,”
Soft Matter
9
(
15
),
3945
3955
(
2013
).
13.
M.
Ballard
 et al., “
Enhancing nanoparticle deposition using actuated synthetic cilia
,”
Microfluid. Nanofluid.
17
(
2
),
317
324
(
2014
).
14.
H.
Shum
 et al., “
Active ciliated surfaces expel model swimmers
,”
Langmuir
29
(
41
),
12770
12776
(
2013
).
15.
A.
Tripathi
,
A.
Bhattacharya
, and
A. C.
Balazs
, “
Size selectivity in artificial cilia–particle interactions: Mimicking the behavior of suspension feeders
,”
Langmuir
29
(
14
),
4616
4621
(
2013
).
16.
J.
Branscomb
and
A.
Alexeev
, “
Designing ciliated surfaces that regulate deposition of solid particles
,”
Soft Matter
6
(
17
),
4066
4069
(
2010
).
17.
A.
Tripathi
,
H.
Shum
, and
A. C.
Balazs
, “
Fluid-driven motion of passive cilia enables the layer to expel sticky particles
,”
Soft Matter
10
(
9
),
1416
1427
(
2014
).
18.
A. C.
Balazs
 et al., “
Designing bioinspired artificial cilia to regulate particle–surface interactions
,”
J. Phys. Chem. Lett.
5
(
10
),
1691
1700
(
2014
).
19.
G.
Wang
 et al., “
Stiffness dependent separation of cells in a microfluidic device
,”
PLoS One
8
(
10
),
e75901
(
2013
).
20.
M.
Eisenstein
, “
Cell sorting: Divide and conquer
,”
Nature
441
(
7097
),
1179
1185
(
2006
).
21.
S.
Sohrabi
 et al., “
Characterization of nanoparticle binding dynamics in microcirculation using an adhesion probability function
,”
Microvasc. Res.
108
,
41
47
(
2016
).
22.
N.
Pamme
and
C.
Wilhelm
, “
Continuous sorting of magnetic cells via on-chip free-flow magnetophoresis
,”
Lab Chip
6
(
8
),
974
980
(
2006
).
23.
N.
Xia
 et al., “
Combined microfluidic-micromagnetic separation of living cells in continuous flow
,”
Biomed. Microdevices
8
(
4
),
299
308
(
2006
).
24.
M. S.
Pommer
 et al., “
Dielectrophoretic separation of platelets from diluted whole blood in microfluidic channels
,”
Electrophoresis
29
(
6
),
1213
1218
(
2008
).
25.
I.
Doh
and
Y.-H.
Cho
, “
A continuous cell separation chip using hydrodynamic dielectrophoresis (DEP) process
,”
Sens. Actuators, A
121
(
1
),
59
65
(
2005
).
26.
M. M.
Wang
 et al., “
Microfluidic sorting of mammalian cells by optical force switching
,”
Nat. Biotechnol.
23
(
1
),
83
87
(
2005
).
27.
A. E.
Ekpenyong
 et al., “
Viscoelastic properties of differentiating blood cells are fate-and function-dependent
,”
PLoS One
7
(
9
),
e45237
(
2012
).
28.
F.
Petersson
 et al., “
Free flow acoustophoresis: Microfluidic-based mode of particle and cell separation
,”
Anal. Chem.
79
(
14
),
5117
5123
(
2007
).
29.
J.
Shi
 et al., “
Continuous particle separation in a microfluidic channel via standing surface acoustic waves (SSAW)
,”
Lab Chip
9
(
23
),
3354
3359
(
2009
).
30.
L.
Ren
 et al., “
A high-throughput acoustic cell sorter
,”
Lab Chip
15
(
19
),
3870
3879
(
2015
).
31.
A. A.
Nawaz
 et al., “
Acoustofluidic fluorescence activated cell sorter
,”
Anal. Chem.
87
(
24
),
12051
12058
(
2015
).
32.
K.
Takahashi
 et al., “
Non-destructive on-chip cell sorting system with real-time microscopic image processing
,”
J. Nanobiotechnol.
2
(
1
),
5
(
2004
).
33.
D.
Di Carlo
 et al., “
Continuous inertial focusing, ordering, and separation of particles in microchannels
,”
Proc. Natl. Acad. Sci. U.S.A.
104
(
48
),
18892
18897
(
2007
).
34.
A. Y.
Fu
 et al., “
An integrated microfabricated cell sorter
,”
Anal. Chem.
74
(
11
),
2451
2457
(
2002
).
35.
L. R.
Huang
 et al., “
Continuous particle separation through deterministic lateral displacement
,”
Science
304
(
5673
),
987
990
(
2004
).
36.
J.
Takagi
 et al., “
Continuous particle separation in a microchannel having asymmetrically arranged multiple branches
,”
Lab Chip
5
(
7
),
778
784
(
2005
).
37.
M.
Yamada
and
M.
Seki
, “
Hydrodynamic filtration for on-chip particle concentration and classification utilizing microfluidics
,”
Lab Chip
5
(
11
),
1233
1239
(
2005
).
38.
J.
Seo
,
M. H.
Lean
, and
A.
Kole
, “
Membrane-free microfiltration by asymmetric inertial migration
,”
Appl. Phys. Lett.
91
(
3
),
033901
(
2007
).
39.
N.
Pamme
, “
Continuous flow separations in microfluidic devices
,”
Lab Chip
7
(
12
),
1644
1659
(
2007
).
40.
W.
Xu
 et al., “
Cell stiffness is a biomarker of the metastatic potential of ovarian cancer cells
,”
PLoS One
7
(
10
),
e46609
(
2012
).
41.
B.
Lincoln
 et al., “
Deformability‐based flow cytometry
,”
Cytometry, Part A
59
(
2
),
203
209
(
2004
).
42.
S. C.
Hur
 et al., “
Deformability-based cell classification and enrichment using inertial microfluidics
,”
Lab Chip
11
(
5
),
912
920
(
2011
).
43.
S.
Chen
and
G. D.
Doolen
, “
Lattice Boltzmann method for fluid flows
,”
Annu. Rev. Fluid Mech.
30
(
1
),
329
364
(
1998
).
44.
D. A.
Reasor
,
J. R.
Clausen
, and
C. K.
Aidun
, “
Coupling the lattice-Boltzmann and spectrin-link methods for the direct numerical simulation of cellular blood flow
,”
Int. J. Numer. Methods Fluids
68
(
6
),
767
781
(
2012
).
45.
T.
Kruger
,
F.
Varnik
, and
D.
Raabe
, “
Efficient and accurate simulations of deformable particles immersed in a fluid using a combined immersed boundary lattice Boltzmann finite element method
,”
Comput. Math. Appl.
61
(
12
),
3485
3505
(
2011
).
46.
J. F.
Zhang
,
P. C.
Johnson
, and
A. S.
Popel
, “
Red blood cell aggregation and dissociation in shear flows simulated by lattice Boltzmann method
,”
J. Biomech.
41
(
1
),
47
55
(
2008
).
47.
J.
Tan
 et al., “
Characterization of nanoparticle dispersion in red blood cell suspension by the lattice Boltzmann-immersed boundary method
,”
Nanomaterials
6
(
2
),
30
(
2016
).
48.
J.
Tan
 et al., “
Numerical simulation of cell squeezing through a micropore by the immersed boundary method
,”
Proc. Inst. Mech. Eng., Part C: J. Mech. Eng. Sci.
232
,
502
514
(
2017
).
49.
J.
Latt
,
Hydrodynamic Limit of Lattice Boltzmann Equations
(
University of Geneva
,
2007
).
50.
Y.
Qian
,
D.
d'Humières
, and
P.
Lallemand
, “
Lattice BGK models for Navier-Stokes equation
,”
Europhys. Lett.
17
(
6
),
479
(
1992
).
51.
M.
Dao
,
J.
Li
, and
S.
Suresh
, “
Molecularly based analysis of deformation of spectrin network and human erythrocyte
,”
Mater. Sci. Eng., C
26
(
8
),
1232
1244
(
2006
).
52.
D. A.
Fedosov
 et al., “
Multiscale modeling of red blood cell mechanics and blood flow in malaria
,”
PLoS Comput. Biol.
7
(
12
),
e1002270
(
2011
).
53.
M.
Nakamura
,
S.
Bessho
, and
S.
Wada
, “
Spring-network-based model of a red blood cell for simulating mesoscopic blood flow
,”
Int. J. Numer. Methods Biomed. Eng.
29
(
1
),
114
128
(
2013
).
54.
J.
Tan
 et al., “
Coupled particulate and continuum model for nanoparticle targeted delivery
,”
Comput. Struct.
122
,
128
134
(
2013
).
55.
S.
Sohrabi
 et al., “
Numerical simulation of particle transport and deposition in the pulmonary vasculature
,”
J. Biomech. Eng.
136
(
12
),
121010
(
2014
).
56.
S.
Sohrabi
and
Y.
Liu
,
A Cellular Model of Shear‐Induced Hemolysis
(
Artificial Organs
,
2016
).
57.
G. A.
Buxton
,
C. M.
Care
, and
D. J.
Cleaver
, “
A lattice spring model of heterogeneous materials with plasticity
,”
Modell. Simul. Mater. Sci. Eng.
9
(
6
),
485
(
2001
).
58.
J.
Mills
 et al.,
Nonlinear Elastic and Viscoelastic Deformation of the Human Red Blood Cell with Optical Tweezers
(
MCB-TECH Science Press
,
2004
), Vol.
1
, pp.
169
180
.
59.
D.
Fedosov
,
B.
Caswell
, and
G. E.
Karniadakis
, “
Coarse-grained red blood cell model with accurate mechanical properties, rheology and dynamics
,” in
Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 2009 (EMBC 2009)
(IEEE,
2009
).
60.
I. V.
Pivkin
and
G. E.
Karniadakis
, “
Accurate coarse-grained modeling of red blood cells
,”
Phys. Rev. Lett.
101
(
11
),
118105
(
2008
).
61.
M.
Abkarian
,
M.
Faivre
, and
A.
Viallat
, “
Swinging of red blood cells under shear flow
,”
Phys. Rev. Lett.
98
(
18
),
188302
(
2007
).
62.
J. S.
Wexler
 et al., “
Bending of elastic fibres in viscous flows: The influence of confinement
,”
J. Fluid Mech.
720
,
517
544
(
2013
).
63.
A.
Alexeev
,
J.
Yeomans
, and
A. C.
Balazs
, “
Designing synthetic, pumping cilia that switch the flow direction in microchannels
,”
Langmuir
24
(
21
),
12102
12106
(
2008
).
64.
A.
Alexeev
 et al., “
Using actuated cilia to regulate motion of microscopic particles
,” in
ASME 2010 First Global Congress on NanoEngineering for Medicine and Biology
(
American Society of Mechanical Engineers
,
2010
).
65.
S.
Sohrabi
and
Y.
Liu
, “
Modeling thermal inkjet and cell printing process using modified pseudopotential and thermal lattice Boltzmann methods
,”
Phys. Rev. E
97
(
3
),
033105
(
2018
).
66.
T.
Krüger
,
F.
Varnik
, and
D.
Raabe
, “
Shear stress in lattice Boltzmann simulations
,”
Phys. Rev. E
79
(
4
),
046704
(
2009
).
67.
D. J.
Holdych
 et al., “
Truncation error analysis of lattice Boltzmann methods
,”
J. Comput. Phys.
193
(
2
),
595
619
(
2004
).
68.
A. J.
Ladd
, “
Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation
,”
J. Fluid Mech.
271
(
1
),
285
309
(
1994
).
69.
Q. S.
Zou
and
X. Y.
He
, “
On pressure and velocity boundary conditions for the lattice Boltzmann BGK model
,”
Phys. Fluids
9
(
6
),
1591
1598
(
1997
).
70.
M.
Lekka
 et al., “
Elasticity of normal and cancerous human bladder cells studied by scanning force microscopy
,”
Eur. Biophys. J.
28
(
4
),
312
316
(
1999
).
71.
S.
Byun
 et al., “
Characterizing deformability and surface friction of cancer cells
,”
Proc. Natl. Acad. Sci. U.S.A.
110
(
19
),
7580
7585
(
2013
).
72.
D. A.
Fedosov
,
J.
Fornleitner
, and
G.
Gompper
, “
Margination of white blood cells in microcapillary flow
,”
Phys. Rev. Lett.
108
(
2
),
028104
(
2012
).
73.
D. A.
Fedosov
,
B.
Caswell
, and
G. E.
Karniadakis
, “
A multiscale red blood cell model with accurate mechanics, rheology, and dynamics
,”
Biophys. J.
98
(
10
),
2215
2225
(
2010
).
74.
Y.
Wang
 et al., “
An immersed boundary-lattice Boltzmann flux solver and its applications to fluid–structure interaction problems
,”
J. Fluids Struct.
54
,
440
465
(
2015
).
75.
A.
De Rosis
and
E.
Lévêque
, “
Central-moment lattice Boltzmann schemes with fixed and moving immersed boundaries
,”
Comput. Math. Appl.
72
(
6
),
1616
1628
(
2016
).
76.
A. D.
Rosis
,
S.
Ubertini
, and
F.
Ubertini
, “
A comparison between the interpolated bounce-back scheme and the immersed boundary method to treat solid boundary conditions for laminar flows in the lattice Boltzmann framework
,”
J. Sci. Comput.
61
(
3
),
477
489
(
2014
).
77.
C. S.
Peskin
, “
The immersed boundary method
,”
Acta Numer.
11
,
479
517
(
2002
).
78.
Y.
Cheng
,
L.
Zhu
, and
C.
Zhang
, “
Numerical study of stability and accuracy of the immersed boundary method coupled to the lattice Boltzmann BGK model
,”
Commun. Comput. Phys.
16
(
1
),
136
168
(
2014
).
79.
C.-H.
Hsu
 et al., “
Microvortex for focusing, guiding and sorting of particles
,”
Lab Chip
8
(
12
),
2128
2134
(
2008
).
80.
A. R.
Lokanathan
 et al., “
Cilia-Mimetic hairy surfaces based on end-immobilized nanocellulose colloidal rods
,”
Biomacromolecules
14
(
8
),
2807
2813
(
2013
).
81.
Y.
Wang
 et al., “
Out of the cleanroom, self-assembled magnetic artificial cilia
,”
Lab Chip
13
(
17
),
3360
3366
(
2013
).
82.
J.
Hussong
 et al., “
Experimental investigation of the flow induced by artificial cilia
,”
Lab Chip
11
(
12
),
2017
2022
(
2011
).
83.
B.
Zhou
 et al., “
Design and fabrication of magnetically functionalized flexible micropillar arrays for rapid and controllable microfluidic mixing
,”
Lab Chip
15
(
9
),
2125
2132
(
2015
).
84.
C.-Y.
Chen
 et al., “
Magnetically actuated artificial cilia for optimum mixing performance in microfluidics
,”
Lab Chip
13
(
14
),
2834
2839
(
2013
).
85.
B.
Evans
 et al., “
Magnetically actuated nanorod arrays as biomimetic cilia
,”
Nano Lett.
7
(
5
),
1428
1434
(
2007
).
86.
C.
Moraes
 et al., “
Supersoft lithography: Candy-based fabrication of soft silicone microstructures
,”
Lab Chip
15
(
18
),
3760
3765
(
2015
).
87.
F.
Liu
 et al., “
An inverted micro-mixer based on a magnetically-actuated cilium made of Fe doped PDMS
,”
Smart Mater. Struct.
25
(
9
),
095049
(
2016
).

Supplementary Material