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.
INTRODUCTION
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 6 Furthermore, magnetic-actuated cilia made from polydimethylsiloxane (PDMS) contain superparamagnetic nanoparticles () which can create a net fluid flow of about 8 at frequencies up to 35 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 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.
COMPUTATIONAL MODEL
Lattice Boltzmann fluid model
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 . The dynamics of the LBM involves streaming and collision processes in phase space .
where is the collision operator, denotes the lattice velocity, and is the force term. Bhatnagar-Gross-Krook (BGK) collision scheme50 is used here
where is a relaxation constant. The population distribution at equilibrium, , can be calculated from the local macroscale fluid velocity
where is the weight coefficient. It should be noted that the existence of external bodies in the flow are introduced through a force term, , expressed in terms of force density and fluid velocity
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.
(a) Spring connected network cell membrane model. (b) Kinematics for local stretching and bending modes of response.
(a) Spring connected network cell membrane model. (b) Kinematics for local stretching and bending modes of response.
Soft particle model
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 connected with nonlinear springs forming two-dimensional triangulated network as shown in Fig. 1. The potential energy of the membrane is defined as
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
where , is the maximum spring extension, is the persistence length, is the energy per unit mass, and is the POW force coefficient. Bending potential corresponds to bending stiffness of erythrocyte which can be expressed as
where is the bending constant, is the instantaneous angle between two adjacent triangles having the common edge j, and is the spontaneous angle. stands for local and global area stiffness of lipid bilayer which is defined as
where and are the global and local area constraint constants, respectively. Similarly, the last potential term conserves red blood cell (RBC) volume
where is volume constraint constant. The nodal forces corresponding to the each energy can be calculated as .
Immersed boundary coupling scheme
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 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
where is the structural force at location at time of and 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
where , for instance, is the distance between the structural force/velocity at location and lattice node.
Cilia array
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 is tethered to the wall . The elastic energy associated with the chain is11
where is the angle between neighboring bonds and represents the cilium bending modulus. The instantaneous angle can be expressed as , where , with being the wall normal direction. Moreover, equals where is Young's modulus, is cilium length, and 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 is applied at the cilia tip, where is the frequency of actuation. The beating pattern is symmetric and the actuation force oscillates harmonically in 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 th cilium bead is given by where is the velocity of th particle and is the interpolated fluid velocity at position of the particle . 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.
Model validation
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 as shown in Fig. 2(a). The bending stiffness, local area constraint, global area constraint, and volume constraint constants are chosen as , , and , 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
(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.
(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.
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 (), the deflection of unactuated cilia is studied. Basically, cilium tip deflection is recorded as a function of imposed flow rate () 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
where and are is channel height and width, respectively; is Young's modulus; and are cilia diameter and moment of inertia, respectively; is a function of . In their experimental study, the longest rectangular shape cilia used had the dimension of , where the length to diameter aspect ratio is .
To benchmark our cilia model, we examined the flow of fluid past an elastic fibre in a channel with dimensions of . 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 used in our simulation as opposed to highest aspect ratio of 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.
Cilia tip deflection, , versus flow rate, , for , , , , , and 10−4.
CILIA-PARTICLE INTERACTION
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
where is the viscous drag coefficient, is the angular velocity of the driving force, and is the bending rigidity of cilium. For relatively small , the dynamic of cilium is governed by its elasticity and cilia array is unable to generate net flows at low 64 For high 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 lattice units. The kinematic viscosity in lattice unit is given by 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 . 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 was simulated to see possible effect. However, no difference in particle's trajectories are observed between cases with and . We also assume , , and . 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 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 and , respectively. Cilia spacing along channel length is 25 where rows are tilted 30° as shown in Fig. 4(a). The radius of each point particle has to be much less than unit for the model to be valid, thus .11 In the PPM simulations, the bond length between the point particles needs to be approximately unit, hence . 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.
(a) 3D demonstration of particle-cilia interaction in microchannel. (b) Particles with distinct biophysical properties tend to flow downstream in different height of microchannel.
(a) 3D demonstration of particle-cilia interaction in microchannel. (b) Particles with distinct biophysical properties tend to flow downstream in different height of microchannel.
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 and , respectively. Furthermore, to cover different sizes and shapes, RBC and , , and 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 particle's membrane, respectively. With lattice spacing of , the ratios of solid to fluid spacing for 4, 8, and 12 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 , , and , 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 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 , , , and . 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 where is cilia beating frequency. Furthermore, the shear Stokes number is computed as where 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 particles tend to laterally drift more compared to and particles. Also, the RBC is shown to displace even more than 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 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 axis.
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, .
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, .
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, . 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.
PARAMETER STUDY
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.
The effect of ciliary angle on lateral displacements and separation efficiencies of soft particulates.
The effect of ciliary angle on lateral displacements and separation efficiencies of soft particulates.
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.
The effect of actuation frequency on lateral displacement of soft and stiff 12 μm particles.
The effect of actuation frequency on lateral displacement of soft and stiff 12 μm particles.
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.
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.
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.
SEPARATION DEVICE DESIGN
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 actuation frequency, it is observed that microparticles laterally displace as they move downstream (see Fig. 5). Subsequently, their trajectories will be separated by 2% to 15% which can ultimately lead to to 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 block with 45° cilia angle after 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.
(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 in the proposed design for particle sorter.
(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 in the proposed design for particle sorter.
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 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 particles as shown in in Fig. 9(b). It should be noted that the length of 30° and 45° blocks are and , respectively. It is observed that soft and stiff particles will be separated by as they move downstream of a channel. It should also be noted that these two particles move within 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.
CONCLUSION
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.
SUPPLEMENTARY MATERIAL
See supplementary material for three multimedia files visualizing the dynamics of interaction between cilia and microparticles.
ACKNOWLEDGMENTS
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.