Modeling Brownian Microparticle Trajectories in Lab-on-a-Chip Devices with Time Varying Dielectrophoretic or Optical Forces

Abstract

Lab-on-a-chip (LOC) devices capable of manipulating micro/nano-sized samples have spurred advances in biotechnology and chemistry. Designing and analyzing new and more advanced LOCs require accurate modeling and simulation of sample/particle dynamics inside such devices. In this work, we present a generalized computational physics model to simulate particle/sample trajectories under the influence of dielectrophoretic or optical forces inside LOC devices. The model takes into account time varying applied forces, Brownian motion, fluid flow, collision mechanics, and hindered diffusion caused by hydrodynamic interactions. We develop a numerical solver incorporating the aforementioned physics and use it to simulate two example cases: first, an optical trapping experiment, and second, a dielectrophoretic cell sorter device. In both cases, the numerical results are found to be consistent with experimental observations, thus proving the generality of the model. The numerical solver can simulate time evolution of the positions and velocities of an arbitrarily large number of particles simultaneously. This allows us to characterize and optimize a wide range of LOCs. The developed numerical solver is made freely available through a GitHub repository so that researchers can use it to develop and simulate new designs.

1. Introduction

Recent decades have seen rapid advances in the field of biotechnology. Lab-on-a-chip (LOC) devices capable of manipulating micro- and nano-sized bio-samples have played a significant role in fueling this progress. The research on designing novel and complex LOC devices is ongoing. A computational-physics framework capable of simulating the trajectory of a micro- or nano-sample can be a very useful tool for designing and analyzing LOCs [1,2].

Although the main application of LOCs is in the field of biology, the design process of the devices requires modeling several physical phenomena. Trapping and manipulating small-sized samples are often a key feature of LOCs. On-chip manipulation of micro-sized samples (or microparticles) is usually accomplished by using viscous drag forces applied through fluid flow along with dielectrophoresis (DEP) [3,4,5] or optical [6,7,8]/ optoelectronic [9,10,11] techniques (DEP is more widely used). LOCs can have integrated micro-electrodes and microfluidic channels to carry out these functions [12]. Designing a device to achieve a specific functionality requires accurately predicting how a micro-object behaves under the influence of relevant forces. For example, microfluidic particle/cell sorter devices often utilize material-selective DEP forces to create separate trajectories for different types of particles (e.g., dead vs live cells) [13,14,15,16,17,18]. The design process hinges on engineering these trajectories. Similar arguments can be made for LOC devices used in other applications as well. A mathematical simulation model capable of predicting particle trajectories can help designers identify potential issues and optimize the device before proceeding to the fabrication step; consequently, saving time and cost.

Along with the effects of external forces, colloidal particles (we use the term particle throughout the paper to refer to any colloidal micro- or nano-specimen of roughly spherical shape) in LOCs exhibit Brownian motion which can affect their trajectories [19]. Brownian motion arises from random collisions between the particle and the molecules of the suspension fluid [20,21]. Although the phenomenon is usually not prominent for large particle sizes (tens of microns) [18], it can be notable when sizes are in sub-micron range. Therefore, it should be included in a generalized model. The Langevin equation, which takes into account the effects of external forces and Brownian motion, can be used for this purpose [20,22]. The viscous drag force experienced by the moving particle in the suspension fluid (which itself can be in motion) is integrated within the equation. This makes it ideal for modeling colloidal particle trajectories in a static or moving fluid.

Unlike colloids in a boundless isotropic fluid, particles in LOC devices are bounded within a finite-sized liquid chamber. As such, additional mechanics need to be taken into account for modeling particle behavior. Firstly, a particle close to the walls of the fluid chamber or a microfluidic channel would experience hindered diffusion instead of free diffusion [23,24,25]. This is due to the hydrodynamic interaction between the particle and the no-slip boundary layer [2] of the fluid near a solid surface. Another factor to consider is particle collision with the solid surfaces of the chip which alters its motion. In multi-particle systems, particle–particle collisions are also possible. The effect of collisions cannot be directly included within the Langevin equation and therefore, must be modeled separately. In this work, we present a numerical model that takes into account all the aforementioned physics to calculate particle trajectories inside LOCs. The numerical solver uses a finite-difference scheme to solve the Langevin partial-differential equation and obtain the time evolution of particle positions. We model hindered diffusion by using a particle-position dependent diffusion tensor in the Langevin equation. A collision detection algorithm is used at each time-step. If detected, an elastic collision model is used to modify the results of the Langevin mechanics. Our generalized model is capable of simulating heterogeneous multi-particle systems, where each particle can differ in size, mass, and material properties. In addition, a large number of particles can be simulated simultaneously. The numerical solver computes the particle positions in parallel and, therefore, does not require long computation time even for large number of particles. Thus, it is ideal for simulating and analyzing various different LOC environments.

The numerical solver is implemented in python and is made available publicly through a GitHub repository [26]. We chose python as it is open source and has many freely available scientific-computing libraries. The software and the mathematical libraries required to run the solver are freely available. Along with the numerical engine, we utilized the graphical libraries to create built-in visualization schemes for analyzing relevant data. This includes the ability to save the animation of the particle motion in the form of a video file and logging the position and velocity values of each particle.

To demonstrate the usefulness of the platform, we present simulation results of two LOC systems. The first example shows optical trapping of multiple microparticles. A dynamic nature is incorporated by activating and deactivating the trap at specific time instances. The corresponding simulation results show the trapping mechanics, as well as the dynamics of the particles when the trap is released. The second example covers a more complex LOC that uses DEP to sort and separate cells. Such LOCs have many applications [12,14,15]. A microfluidic device with flowing liquid is considered for this case. Material selective DEP force is used to direct viable (live) and non-viable (dead) yeast cells along different paths. These two distinctly different examples show the versatility of our solver for simulating various LOC configurations. Moreover, the solver is developed in modular form. Thus, switching from one device geometry to another requires modifying only a few modules and plugging in external force terms.

The rest of the paper is organized as follows: Section 2 introduces the Langevin equation and its discretization to the finite difference form. The particle position dependent diffusion tensor and collision mechanics are discussed in Section 3 and Section 4, respectively. Section 5 describes how all the physics were integrated in the numerical solver. The simulation results for the two different LOCs are discussed in Section 6. Finally, concluding remarks are made in Section 7.

2. Langevin Equation

The motion of a colloidal micro- or nano-particle under the influence of an external force-field can be modeled using the Langevin equation. For a colloidal particle in a low Reynolds number environment with shear-flow, the Langevin equation has the form [20,27,28,29]:

$$m\frac{\partial \mathbf{v}(\mathbf{r},t)}{\partial t}=\frac{k_{B}T}{\overleftrightarrow{\mathbf{D}}\left(\mathbf{r}\right)}\left[{\mathbf{v}}_f(\mathbf{r},t)-\mathbf{v}(\mathbf{r},t)+\sqrt{2}{\overleftrightarrow{\mathbf{D}}}_{\frac{1}{2}}\left(\mathbf{r}\right)\mathbf{W}\left(t\right)\right]+{\mathbf{F}}_{\mathrm{ext}}(\mathbf{r},t).$$

(1)

Here, t is the time variable, m, $\mathbf{r}$ and $\mathbf{v}(\mathbf{r},t)$ are the mass, position, and velocity of the particle, respectively, ${\mathbf{v}}_f(\mathbf{r},t)$ is the velocity of the fluid, $k_B$ is the Boltzmann constant, T is the temperature, ${\mathbf{F}}_{\mathrm{ext}}$ is the net external applied force acting on the particle, $\overleftrightarrow{\mathbf{D}}$ is the diffusion tensor, and $\mathbf{W}\left(t\right)$ is a vector white noise term. The tensor ${\overleftrightarrow{\mathbf{D}}}_{\frac{1}{2}}$ is defined as the element-wise square root of $\overleftrightarrow{\mathbf{D}}$. Each Cartesian component of $\mathbf{W}\left(t\right)$ is a random process with zero mean and unit variance. The $m\dot{\mathbf{v}}(\mathbf{r},t)$ term represents the inertia of the particle. This term can often be dropped for small particles where the mass is negligible [20]. The friction/drag force term, $[{\mathbf{v}}_f(\mathbf{r},t)-\mathbf{v}(\mathbf{r},t)]k_{B}T/\overleftrightarrow{\mathbf{D}}$, depends on the viscous properties of the fluid [17]. The white noise term $\sqrt{2}{\overleftrightarrow{\mathbf{D}}}_{\frac{1}{2}}\left(\mathbf{r}\right)\mathbf{W}\left(t\right)$ is the fluctuating force resulting from random collisions between the particle and fluid molecules. The ${\mathbf{F}}_{\mathrm{ext}}$ term encapsulates all external forces acting on the particle (e.g., optical trapping force, dielectrophoretic force, gravitational force, etc.). Whereas most research articles use the derivatives of the position vector, we express Equation (1) in terms of the velocity vector and its derivative. We do so for easier subsequent integration of collision physics into the model which depends on the velocity vectors. Thus, directly solving the velocity vectors simplifies the analysis.

Equation (1) is a stochastic differential equation. The presence of the discontinuous white noise term makes solving a stochastic differential equation more complex than solving an ordinary differential equation [30]. Our approach to numerically solve Equation (1) involves using the Euler–Maruyama method to convert it into a finite difference equation [30,31]. We use the following discretization scheme:

$$m\frac{{\mathbf{v}}_{i+1}-{\mathbf{v}}_i}{\Delta t}=\frac{k_{B}T}{\overleftrightarrow{\mathbf{D}}\left({\mathbf{r}}_i\right)}\left[{\mathbf{v}}_{f,i}-{\mathbf{v}}_{i+1}+\sqrt{\frac{2}{\Delta t}}{\overleftrightarrow{\mathbf{D}}}_{\frac{1}{2}}\left({\mathbf{r}}_i\right){\mathbf{w}}_i\right]+{\mathbf{F}}_{\mathrm{ext},i}\left({\mathbf{r}}_i\right).$$

(2)

Here, $i=0,1,\cdots N_t$ represents the time index and $\Delta t$ is the time step. The set of discrete time where the numerical solution is calculated is given by $t_i=i\Delta t$. The term $\mathbf{w}{}_i$ is a random vector whose Cartesian components are Gaussian random numbers with zero mean and unit variance. Solving for ${\mathbf{v}}_{i+1}$ from Equation (2) gives:

$${\mathbf{v}}_{i+1}=\frac{\overleftrightarrow{\Lambda }}{1+\overleftrightarrow{\Lambda }}{\mathbf{v}}_i+\frac{1}{1+\overleftrightarrow{\Lambda }}\left[{\mathbf{v}}_{f,i}+\sqrt{\frac{2}{\Delta t}}{\overleftrightarrow{\mathbf{D}}}_{\frac{1}{2}}\left({\mathbf{r}}_i\right){\mathbf{w}}_i+\frac{\overleftrightarrow{\mathbf{D}}\left({\mathbf{r}}_i\right)}{k_{B}T}{\mathbf{F}}_{\mathrm{ext},i}\left({\mathbf{r}}_i\right)\right].$$

(3)

Here, $\overleftrightarrow{\Lambda }\triangleq \frac{m\overleftrightarrow{\mathbf{D}}}{k_{B}T\Delta t}$ is a unit-less tensor quantity defined to make the equation compact. It should be noted that we use ${\mathbf{v}}_{i+1}$ instead of ${\mathbf{v}}_i$ in the right hand side of Equation (2). This ensures the stability of the numerical scheme for arbitrary $\Delta t$ values. This fact can be verified by setting the fluid velocity, the white noise term, and the external force to be zero and observing the simplified velocity update equation: ${\mathbf{v}}_{i+1}=\frac{\overleftrightarrow{\Lambda }}{1+\overleftrightarrow{\Lambda }}{\mathbf{v}}_i$. As the multiplicative factor is less than unity (every component of $\overleftrightarrow{\Lambda }$ is positive), the velocity will remain bounded over successive iterations and thus, stay stable.

After calculating the velocity vector at a given time step i, the position vector can be easily calculated using the Euler–Cromer method [32]:

$${\mathbf{r}}_{i+1}={\mathbf{r}}_i+{\mathbf{v}}_{i+1}\Delta t.$$

(4)

The position and velocity data of each particle is calculated using the same equations. The process is repeated for the next discrete time instance until the simulation end time is reached.

3. Diffusion Tensor and Hydrodynamic Interactions

The Langevin equation contains the diffusion tensor which is a representation of the fluid–particle interaction. The diffusion tensor in an isotropic homogeneous medium, ${\overleftrightarrow{\mathbf{D}}}_0$, depends on the material properties of the suspension medium, particle size, and the temperature. It is given by [23,24]:

$${\overleftrightarrow{\mathbf{D}}}_0=\frac{k_{B}T}{6\pi \eta r_o}\overleftrightarrow{\mathbf{I}},$$

(5)

where $\overleftrightarrow{\mathbf{I}}$ is a $3\times 3$ unit tensor, $r_o$ is the radius of the particle, and $\eta$ is the dynamic viscosity of the medium. Equation (5) is equivalent to using a scalar diffusion coefficient. For LOC devices, however, the fluid chamber or microfluidic channels are not large enough compared to the particle size for assuming isotropic conditions. In addition, particles are often located near the bottom surface during manipulation in such devices. So, the particle dynamics are different from the isotropic case. The hydrodynamic interactions between the particle and the static fluid layers near a surface hinder the diffusion process. Thus, the diffusion tensor becomes dependent on the relative position of the particle with respect to the solid surfaces of the device. The generalized diffusion tensor for such cases is given by [23,24]:

$$\overleftrightarrow{\mathbf{D}}\left(\mathbf{r}\right)=\frac{k_{B}T}{6\pi \eta r_o}\overleftrightarrow{\mathbf{H}}\left(\mathbf{r}\right).$$

(6)

$\overleftrightarrow{\mathbf{H}}\left(\mathbf{r}\right)$ is a $3\times 3$ unitless tensor that depends on the distance of the particle from the solid surfaces. For a geometry aligned with the coordinate system, only the diagonal elements of the tensor will be non-zero, i.e., $\overleftrightarrow{\mathbf{H}}\left(\mathbf{r}\right)=\left[H_{jk}\left(\mathbf{r}\right)\right]$, $H_{jk}\left(\mathbf{r}\right)=0$$\forall j\ne k$.

The majority of the hydrodynamic interactions in LOCs arise from the bottom/top surface along which particle manipulation takes place. Thus, we consider a case with a single solid surface located at the $z=0$ plane. The components of $\overleftrightarrow{\mathbf{H}}\left(\mathbf{r}\right)$ can be expressed as $H_{11}\left(\mathbf{r}\right)=H_{22}\left(\mathbf{r}\right)=H_{\Vert }\left(\mathbf{r}\right)$ and $H_{33}\left(\mathbf{r}\right)=H_{\perp }\left(\mathbf{r}\right)$ where [23,24,25]:

$$\begin{array}{c}{H}_{\Vert }\left(\mathbf{r}\right)=1-\frac{9r_o}{16z}+\frac{r_o^3}{8z^3}-\frac{45r_o^4}{256z^4}-\frac{r_o^5}{16z^5},\end{array}$$

(7)

$$\begin{array}{c}{H}_{\perp }\left(\mathbf{r}\right)=\frac{6z^2+2r_{o}z}{6z^2+9r_{o}z+2r_o^2}.\end{array}$$

(8)

Here, z is the z-position of the particle. In Figure 1, $H_{\Vert }\left(\mathbf{r}\right)$ and $H_{\perp }\left(\mathbf{r}\right)$ are plotted as functions of the normalized variable $z/r_o$. It can be observed that for particle position far from the $z=0$ surface (i.e., high $z/r_o$ values), the value of both terms approach unity and the diffusion tensor approaches the free space diffusion tensor. The closer the particle is to the surface, the larger the variation is from the free space case. Values less than unity suggests that the diffusion process is hindered near solid surfaces.

4. Collision Mechanics

In a multi-particle system, particle-particle collisions need to be taken into account when simulating particle trajectories. We consider a system with $N_p$ particles where particle p and particle q collide at time instance i. Their mass, position, and initial velocity are $m_p$, $m_q$, ${\mathbf{r}}_{p,i}$, ${\mathbf{r}}_{q,i}$, and ${\mathbf{v}}_{p,i}$, ${\mathbf{v}}_{q,i}$, respectively. The velocity vectors after the collision, ${\mathbf{v}}_{p,i}^{\mathrm{col},\mathrm{pp}}$ and ${\mathbf{v}}_{q,i}^{\mathrm{col},\mathrm{pp}}$, can be calculated using the theory of elastic collisions where the kinetic energy and momentum of the particles are conserved.

The collision mechanics in three-dimensional systems can be treated by considering the velocities along the normal and tangential directions separately. We note that the velocity along the normal direction of the colliding surfaces are altered after collision, whereas the tangential velocity components are unaffected [33]. First, we define the normal vector:

$$\widehat{\mathbf{n}}=\frac{{\mathbf{r}}_{p,i}-{\mathbf{r}}_{q,i}}{|{\mathbf{r}}_{p,i}-{\mathbf{r}}_{q,i}|}.$$

(9)

The normal and tangential component of the particle p velocity is then given by $\left({\mathbf{v}}_{p,i}·\widehat{\mathbf{n}}\right)\widehat{\mathbf{n}}$ and ${\mathbf{v}}_{p,i}-\left({\mathbf{v}}_{p,i}·\widehat{\mathbf{n}}\right)\widehat{\mathbf{n}}$, respectively. One-dimensional elastic collision formula can be applied for the normal directions. The velocity after the collision for particle p is given by:

$${\mathbf{v}}_{p,i}^{\mathrm{col},\mathrm{pp}}=\frac{\left({\mathbf{v}}_{p,i}·\widehat{\mathbf{n}}\right)\widehat{\mathbf{n}}(m_p-m_q)+2m_q\left({\mathbf{v}}_{q,i}·\widehat{\mathbf{n}}\right)\widehat{\mathbf{n}}}{m_p+m_q}+{\mathbf{v}}_{p,i}-\left({\mathbf{v}}_{p,i}·\widehat{\mathbf{n}}\right)\widehat{\mathbf{n}}.$$

(10)

It can be noted that the first term of Equation (10) is identical to the 1D elastic collision formula along the normal axis. Simplifying Equation (10) gives:

$${\mathbf{v}}_{p,i}^{\mathrm{col},\mathrm{pp}}={\mathbf{v}}_{p,i}-\frac{2m_p}{m_p+m_q}\left[\left({\mathbf{v}}_{p,i}-{\mathbf{v}}_{q,i}\right)·\widehat{\mathbf{n}}\right]\widehat{\mathbf{n}}$$

(11)

Swapping the subscripts p and q in Equation (11) gives the modified velocity equation for particle q.

In addition to particle–particle collisions, collisions with a wall/solid-surface may also occur in LOC systems (e.g., collision with bottom substrate, sidewalls of microfluidic channels, etc.). Assuming that the walls are much more massive than the particles, the velocity of a particle after collision can be calculated using a reflection operation. If particle p collides with a flat wall with surface normal ${\widehat{\mathbf{n}}}_{\mathrm{col},\mathrm{pw}}$, then the velocity is modified as:

$${\mathbf{v}}_{p,i}^{\mathrm{col},\mathrm{pw}}={\mathbf{v}}_{p,i}-2\left({\mathbf{v}}_{p,i}·{\widehat{\mathbf{n}}}_{\mathrm{col},\mathrm{pw}}\right){\widehat{\mathbf{n}}}_{\mathrm{col},\mathrm{pw}}$$

(12)

Using Equations (11) and (12), most common collision mechanics occurring in LOC environments can be modeled.

5. Physics Integration

We incorporate all the aforementioned physics into the particle dynamics solver. The solver starts with randomly populated position and velocity vectors for each particle. A vector white noise term is generated from a normal distribution. The particle position depended diffusion tensor is calculated from Equation (6). Then Equation (3) and (4) is used to find predicted particle velocity and position for the next time step (i.e., ${\mathbf{v}}_{i+1}$ and ${\mathbf{r}}_{i+1}$, respectively). As the Langevin equation does not take into account the presence of any obstacles in the particle path, the calculated position vector may represent an unphysical condition (i.e., two particles overlapping, particles going through walls, etc.). To correct for this, the Euclidean distances between each particle pair and each particle-wall pair are calculated to check for collisions. If a collision is detected, then position and velocity given by the Langevin equations needs to be adjusted. Using Equation (11) or Equation (12), the velocity vector is adjusted as:

$${\mathbf{v}}_{p,i+1}^{\mathrm{adj}}=\left\{\begin{array}{cc}{\mathbf{v}}_{p,i+1}^{\mathrm{col},\mathrm{pp}},\hfill & \mathrm{if}\mathrm{collision}\mathrm{with}\mathrm{particle}.\hfill \\ {\mathbf{v}}_{p,i+1}^{\mathrm{col},\mathrm{pw}},\hfill & \mathrm{if}\mathrm{collision}\mathrm{with}\mathrm{wall}.\hfill \end{array}\right.$$

(13)

The corresponding adjusted position vector is:

$${\mathbf{r}}_{p,i+1}^{\mathrm{adj}}={\mathbf{r}}_{p,i}+{\mathbf{v}}_{p,i+1}^{\mathrm{adj}}\Delta t.$$

(14)

The velocity and position vectors are replaced by the adjusted quantities (i.e., ${\mathbf{v}}_{p,i+1}\leftarrow {\mathbf{v}}_{p,i+1}^{\mathrm{adj}}$ and ${\mathbf{r}}_{p,i+1}\leftarrow {\mathbf{r}}_{p,i+1}^{\mathrm{adj}}$) before proceeding to the next time step. The Python code implements the entire process in parallel for all the particles (i.e., $p=0,1,\cdots N_p-1$) inside a single time loop. The implemented wall collision physics is general and works for any arbitrary polygon geometry. The flowchart shown in Figure 2 shows all the calculation steps.

6. Results

Using the developed algorithm, we present simulation results of two example LOC setups. The first example is a typical optical trapping setup with a Gaussian beam. The second example is a microfluidic LOC device that uses dielectrophoresis to sort or separate live and dead yeast cells.

6.1. Optical Trap

We consider an optical trapping experiment with a Gaussian beam. We implement a time dependent external force, modeling an optical beam that is turned ON and OFF at specific time instances. For simplicity, we assume a conservative gradient-force profile resulting from a Gaussian potential well, $u\left(\mathbf{r}\right)$, centered around the origin. The potential is defined as:

$$u\left(\mathbf{r}\right)=\frac{4}{3}\pi r_o^3{A}_d{e}^{-{\left(\frac{r}{w}\right)}^2}.$$

(15)

Here, $A_d=120k_{B}T$ $\mathsf{\mu }{\mathrm{m}}^{-3}$ is the depth of the volume density of the potential well (i.e., depth of the potential well normalized with respect to particle volume), $r_o$ is the radius of a particle, $w=10 \mathsf{\mu }\mathrm{m}$ is representative of the width of the well, and r is the radial distance from the trap center (assumed to be at the origin) to the particle position. The corresponding force profile, ${\mathbf{F}}_{\mathrm{opt}}\left(\mathbf{r}\right)$, is given by [7,34]:

$${\mathbf{F}}_{\mathrm{opt}}\left(\mathbf{r}\right)=-\nabla u\left(\mathbf{r}\right).$$

(16)

The force profile is shown in Figure 3. Due to the symmetry of the potential well, we only plot the radial component of the force, $F_r$. The x, y and z components of the force can easily be calculated from Equation (16) as well. Both the potential profile and the force profile are assumed to scale linearly with particle volume. This is a commonly used approximation for optical traps [35]. The conservative force-profile approximation is also commonly accepted for most cases. However, for some cases, especially involving near-field trapping, it is possible to have non-conservative force components [7]. Instead of Equation (16), more generalized approaches are necessary to model the force for those cases. For the current case, we define the ${\mathbf{F}}_{\mathrm{ext}}$ term in the particle dynamics model as:

$${\mathbf{F}}_{\mathrm{ext}}(\mathbf{r},t)=\left\{\begin{array}{cc}{\mathbf{F}}_{\mathrm{opt}}\left(\mathbf{r}\right),\hfill & \mathrm{if}1 \mathrm{s}<t<8 \mathrm{s}.\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(17)

This represents a case when the optical beam is ON during the time interval $1 \mathrm{s}<t<8 \mathrm{s}$. The interval boundaries are set arbitrarily. It should be noted that other possible external force terms (e.g., gravitational sedimentation, thermophoresis, etc.) can also be added here.

We consider $N_p=4$ polystyrene beads (density, $\rho =1055 \mathrm{kg}/{\mathrm{m}}^3$) suspended in water (dynamic viscosity, $\eta =8.9\times {10}^{-4} \mathrm{Pas}$, and temperature, $T=300 \mathrm{K}$) with radii $r_{o,0}=1.5 \mathsf{\mu }\mathrm{m}$, $r_{o,1}=2.5 \mathsf{\mu }\mathrm{m}$, and $r_{o,2}=r_{o,3}=2 \mathsf{\mu }\mathrm{m}$. The corresponding particle masses can be calculated from the volume and density values. A solid bottom surface is assumed to be located at the $z=0$ plane. With these conditions, the simulation was run from $t=0$ to $t=12 \mathrm{s}$ with a step size of $\Delta t=0.01 \mathrm{s}$. The simulations results are shown in Figure 4 and Figure 5. From Figure 4, we can see the particle positions at an arbitrary time instance ($t=6.3 \mathrm{s}$), when they are trapped on the optical spot. The associated multimedia file 1 is an animated video that shows the same results for all time instances. Figure 5 show the positions and velocities of the particles as a function of time. The position plots show that the particles start moving towards the origin at around $t=1 \mathrm{s}$, when the Gaussian optical spot (centered around the origin) is turned ON. As the particles move towards the optical spot, they start colliding with each other. This is represented by the spiky nature of the plots during the time interval from $t\approx 2 \mathrm{s}$ to $t\approx 8 \mathrm{s}$. After $t>8 \mathrm{s}$, the optical spot is turned OFF and the particles are no longer trapped. During this time the Brownian motion dominates as the particles move randomly. It should be noted that a particle does not always become trapped when the optical spot is active. Trapping occurs only when a particle is within the capturing range of the trap [36,37]. In other instances of the simulation, where the random initial position of a particle is significantly away from the optical spot, the particle does not become trapped.

6.2. Dielectrophoretic Cell Sorting or Separation

For our second example, we simulate a microfluidic device that uses DEP forces for cell sorting or separation which is a common application of LOCs [15,16,17]. Specifically, we focus on the separation of viable (live) and non-viable (dead) yeast cells using the device proposed by Doh et al. in their paper [15]. The device uses tapered electrodes embedded within a microfluidic chamber. The structure is shown in Figure 6a. The fluid walls are 50 $\mathsf{\mu }\mathrm{m}$ high. The electrodes having a height of $0.12$ $\mathsf{\mu }\mathrm{m}$ can be considered planar. The electrodes have a length of 350 $\mathsf{\mu }\mathrm{m}$ (along the x axis) and the gap between them is 16 $\mathsf{\mu }\mathrm{m}$. The top and bottom electrodes have a width (along y axis) of 20 $\mathsf{\mu }\mathrm{m}$, whereas the width of the tapered middle electrode spans $46–70 \mathsf{\mu }\mathrm{m}$. The top and bottom electrodes are electrically shorted externally. An AC voltage of 4 V at 5 MHz is applied between the center electrode and the other electrodes. A strong field is created between the electrodes generating dielectrophoretic forces. As a solution of cells flows from left to right, the electrodes exert different forces on live and dead cells, modifying their trajectories. Thus, specific cells can be diverted to specific microfluidic channels. We choose this device as an example as it incorporates some of essential particle manipulation mechanics (e.g., fluid flow, dielectrophoretic force, etc.) used in many LOC devices.

To use our numerical solver to simulate this device, first we have to model the forces acting on the cells. The time averaged DEP force, ${\mathbf{F}}_{\mathrm{DEP}}$, on the cells can be approximated as [3,15,17]:

$${\mathbf{F}}_{\mathrm{DEP}}=2\pi {ϵ}_0{ϵ}_m{r}_o^3\Re \left[{\tilde{C}}_M\right]\nabla {\left|\mathbf{E}\right|}^2.$$

(18)

Here, ${ϵ}_0$ is the permittivity of free space, ${ϵ}_m=80$ is the relative permittivity of the medium, $r_o=3 \mathsf{\mu }\mathrm{m}$ [38] is the cell radius, ${\tilde{C}}_M$ is the Clausius–Mossotti factor, $\Re \left[·\right]$ is the operator that outputs the real part of its argument, and $\mathbf{E}$ is the electric field. $\mathbf{E}$ depends on the geometry of the electrodes and the device. It can be noted that the sign of $\Re \left[{\tilde{C}}_M\right]$ determines whether ${\mathbf{F}}_{\mathrm{DEP}}$ is attractive (positive DEP) or repulsive (negative DEP). Since ${\tilde{C}}_M$ depends on the material properties of the cells, they differ from live cells to dead cells [39]. We use the multi-shell model to calculate ${\tilde{C}}_M$ for viable and non-viable yeast cells [3]. The procedure is discussed in detail in the Supplementary Information Document. We take the geometrical and material parameters of the cells and the device dimensions from the literature [15,40,41]. Plugging in these values, we find that at $f=5\mathrm{MHz}$, the viable and non-viable cells have relative permittivities of ${ϵ}_{\mathrm{viable}}=199.94$, and ${ϵ}_{\mathrm{nonviable}}=18.82$, respectively, and have conductivities of ${\sigma }_{\mathrm{viable}}=0.36 \mathrm{S/m}$ and ${\sigma }_{\mathrm{nonviable}}=0.013 \mathrm{S/m}$, respectively. These values correspond to $\Re \left[{\tilde{C}}_{M,\mathrm{viable}}\right]=0.945$ and $\Re \left[{\tilde{C}}_{M,\mathrm{nonviable}}\right]=-0.25$. Thus, the viable cells experience positive DEP, whereas the non-viable cells experience negative DEP.

The electric field term in Equation (18) can be calculated using any of the widely available numerical electromagnetic field solvers for any arbitrary device geometry. We use Comsol Multiphysics^® for this purpose. The calculated electric potential and field distributions are shown in Figure 6b,c. The corresponding force values on a viable yeast cell (calculated using Equation (18)) are shown in Figure 6d. It can be seen that the force pulls the viable cells near the gap regions between the electrodes where the field intensity is maximum. For a non-viable cell, the sign of the force would be reversed and scaled by a factor of $|\Re \left[{\tilde{C}}_{M,\mathrm{nonviable}}\right]|/|\Re \left[{\tilde{C}}_{M,\mathrm{viable}}\right]|\approx 0.26$.

Unlike the optical tweezers example, the cell sorter device works with a flowing fluid instead of a stationary pool of liquid. A left to right fluid flow with velocity of 150 $\mathsf{\mu }\mathrm{m/s}$ is considered. Due to the presence of the microfluidic channels, the fluid velocity does not maintain the purely left-to-right direction throughout the device. Considering the fact that the average fluid flow direction is tangential to the channel surface, we use the following simplified piece-wise model:

$${\mathbf{v}}_f(\mathbf{r},t)=\left\{\begin{array}{cc}150\widehat{\mathbf{x}}+30\widehat{\mathbf{y}},\hfill & \mathrm{if}x>170\mathsf{\mu }\mathrm{m},y>30\mathsf{\mu }\mathrm{m}.\hfill \\ 150\widehat{\mathbf{x}}-30\widehat{\mathbf{y}},\hfill & \mathrm{if}x>170\mathsf{\mu }\mathrm{m},y<-30\mathsf{\mu }\mathrm{m}.\hfill \\ 150\widehat{\mathbf{x}},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(19)

Here, ${\mathbf{v}}_f(\mathbf{r},t)$ has units of $\mathsf{\mu }\mathrm{m/s}$. The y component of the fluid velocity for a given x component is approximated from the slope of the channels so that the average flow is tangential to the fluid wall. A 3D moving average filter was applied to make the velocity profile smoother. It should be stated that this flow profile is an approximation. We also ignore any fluid flow that might be induced by the non-uniform electric field [42]. It is possible to plug in more accurate data from an external solver in our code. However, the piece-wise model is sufficiently accurate for the current demonstration problem. In the future, we plan to implement a numerical fluid flow solver and integrate it with our particle dynamics solver.

From Equations (18) and (19), we plug in the ${\mathbf{F}}_{\mathrm{ext}}={\mathbf{F}}_{\mathrm{DEP}}$ and the fluid flow velocity in Equations (3) and (4) to solve the particle dynamics in the time range $t=0$ to $t=6 \mathrm{s}$ with a time step of $\Delta t=5\mathrm{m}\mathrm{s}$. Like the optical trapping example, the particle–particle collision and particle–wall collision models are included in the simulation. We simultaneously simulate the trajectory of 17 viable yeast cells and 17 non-viable yeast cells with random initial positions located in the domain $x\in [-480,-200]\mathsf{\mu }\mathrm{m}$, $y\in [-280,280]\mathsf{\mu }\mathrm{m}$ and $z\in [6,80]\mathsf{\mu }\mathrm{m}$. The simulation results are shown in Figure 7. Four time instances are depicted here. The entire time evolution of the cell positions can be seen in the attached multimedia file 2. Although both live and dead cells experience the same fluid flow, they experience different DEP forces as they have different $\Re \left[{\tilde{C}}_{\mathrm{M}}\right]$ values. We note that as the cells are carried left-to-right by the fluid flow, the live cells (depicted in red) become deflected top/bottom (near the electrode gaps) whereas the dead cells (depicted in blue) are pushed towards the center electrode. This is consistent with the force profile shown in Figure 6d. As the cells travel to the right of the device near the output fluid channels, the cells are grouped as can be seen in $t=1.3 \mathrm{s}$ and $t=2.7 \mathrm{s}$ frames of Figure 7. The top and bottom output fluid channels end up having a much higher concentration of live cells and the center channel contains most of the dead cells.

The sorting mechanics can be analyzed further by plotting the trajectory of a few cells, as shown in Figure 8. During $0<t<0.8 \mathrm{s}$, the y position of the all the cells remain constant as they are not yet near the electrodes to experience the dielectrophoretic force. At this stage, driven by the fluid flow, their motion is limited to the x direction. During $0.8 \mathrm{s}<t<3.5 \mathrm{s}$, the cells are within the region of influence of the electrodes. The live cells move toward the high field region (the region between the electrodes) and dead cells move away from there. Since the $y=0$ line is at the middle of the two high field regions (note Figure 6c), the dead cells accumulate around this line. On the other hand, the live cells accumulate near the top or bottom electrodes. At $t>3.5 \mathrm{s}$, each cell enters one of the three output microfluidic channels. The fluid flows in the top and bottom channels change directions creating a shift in cell trajectory. This can be noted in Figure 8 as the live cells follow a trajectory parallel to the top or bottom channels. The dead cells follow a constant y path through the center channel. This is consistent with the experimental results reported by Doh et al. in their paper [15]. Thus, the simulation accurately demonstrates the cell sorting or separation mechanics.

7. Conclusions

A numerical model that can simulate behavior of particles in LOC devices was developed. Using the model, we successfully simulated the operations of an optical trap and a dielectrophoretic cell sorter. The numerical solver is generally applicable and can simulate most LOC devices once the external force terms are defined. It can be used for the development and optimization of new and existing LOC devices. The solver was coded using open-source Python libraries and can be run on freely accessible compilers. We also made our code for the solver publicly available [26]. Due to the general nature, the developed solver can be of interest to researchers working on biomicrofluidics and LOC devices.