Computer simulation plays a growing role in the design of microfluidic chips. However, the particle tracers in some existing commercial computational fluid dynamics software are not well suited for accurately simulating the trajectories of particles such as cells, microbeads, and droplets in microfluidic systems. To address this issue, we present a microfluidics-optimized particle simulation algorithm (MOPSA) that simulates the trajectories of cells, droplets, and other particles in microfluidic chips with more lifelike results than particle tracers in existing commercial software. When calculating the velocity of a particle, MOPSA treats the particle as a two-dimensional rigid circular object instead of a single point. MOPSA also checks for unrealistic interactions between particles and channel walls and applies an empirical correcting function to eliminate these errors. To validate the performance of MOPSA, we used it to simulate a variety of important features of microfluidic devices like channel intersections and deterministic lateral displacement (DLD) particle sorter chips. MOPSA successfully predicted that different particle sizes will have different trajectories in six published DLD experiments from three research groups; these DLD chips were used to sort a variety of different cells, particles, and droplets. While some of these particles are not actually rigid or spherical, MOPSA's approximation of these particles as rigid spheres nonetheless resulted in lifelike simulations of the behaviors of these particles (at least for the particle sizes and types shown here). In contrast, existing commercial software failed to replicate these experiments. Finally, to demonstrate that MOPSA can be extended to simulate other properties of particles, we added support for simulating particle density to MOPSA and then used MOPSA to simulate the operation of a microfluidic chip capable of sorting cells by their density. By enabling researchers to accurately simulate the behavior of some types of particles in microfluidic chips before fabricating the chips, MOPSA should accelerate the development of new microfluidic devices for important applications.

## I. INTRODUCTION

The field of lab-on-a-chip and micro total analysis devices has been growing rapidly for nearly 40 years,^{1} and microfluidic chips have found important applications in biological and chemical analyses. Many of these chips contain particles such as cells, microbeads, or droplets. For example, microfluidic cell sorters are emerging as powerful tools for point-of-care testing and biological research.^{2,3}

As the variety of microfluidic devices increases, computer-based simulations of microfluidics have become more important.^{4–6} Microfluidics researchers are increasingly using finite element analysis (FEA) software such as COMSOL Multiphysics (COMSOL Inc., Burlington, MA), Autodesk CFD (Autodesk Inc., Venice, CA), and Fluent (ANSYS Inc., Canonsburg, PA) to simulate and optimize their microfluidic chip designs. These software tools excel at calculating the fluid velocity field in a chip. However, these tools are not well suited for simulating the behavior of particles like cells, beads, and droplets flowing in microfluidic chips. For example, Fig. 1 shows the results from using COMSOL Multiphysics to simulate the behavior of a deterministic lateral displacement (DLD) cell sorter. DLD sorters use arrays of micron-scale posts to sort cells based on their size. As cells flow through the chip and interact with the posts, large and small cells follow different paths and are therefore separated by the chip.^{7,8} However, in the simulation results from COMSOL Multiphysics in Fig. 1, the larger 10 *μ*m cells (red) and smaller 3 *μ*m cells (green) follow exactly the same trajectory. If the simulation were accurate, the cells would have followed different trajectories. Additionally, by zooming in to the simulation results in Fig. 1(a), we see that the 10 *μ*m cell actually overlaps with the wall of a triangular post (an impossible situation). These observations confirm that the particle tracing algorithm in COMSOL Multiphysics treats the particle as a mathematical point with no area or volume, an assumption that can be acceptable in macro-scale physics but not in microfluidics. This assumption leads to additional problems: if the simulated single-point particle in COMSOL Multiphysics is located on a mesh node with a fluid velocity magnitude of zero, the particle will remain stuck at that location forever [Fig. 1(b)]. This sticking behavior remains regardless of which boundary condition is chosen for the wall in COMSOL Multiphysics (*freeze, stick, bounce, disappear*, or *reflect diffusely*).

The lack of commercial software tools that accurately simulate particles in microfluidics has profound consequences. For example, researchers developing new cell sorters like the DLD cell sorter shown in Fig. 1 typically have to fabricate and test these chips to determine whether they sort the desired cells. While some empirical models for predicting the sorting behavior of DLD chips do exist for chips with cylindrical posts,^{9} these models do not easily extend to posts with arbitrary shapes like the triangles shown in Fig. 1. And with each chip design iteration requiring fabrication and testing in the lab, it can take months for researchers to develop a functional device. This slows the progress of research and keeps valuable research tools and lifesaving medical diagnostics out of the hands of the people who need them.

To address this problem, we developed a microfluidics-optimized particle simulation algorithm (MOPSA) that is capable of simulating lifelike trajectories for some particles in situations where existing software tools fail. MOPSA treats a particle as a 2D rigid circular object instead of a single point when calculating the particle's velocity. MOPSA also checks whether the particle overlaps a solid object (a post or channel wall) at each simulated time step. If overlap is detected, MOPSA applies an empirical correcting function, *wallEffect*, to shift the particle to a new position where it will not overlap walls. By providing researchers with accurate simulations of some types of particles in microfluidic devices, MOPSA should accelerate the development of new microfluidic devices for important applications in biological research, medical diagnostics, and beyond.

Assuming that a particle is a rigid circle means that MOPSA cannot directly simulate the deformability or non-spherical shapes of some particles. In spite of this limitation, by assuming that the particles are rigid circles, MOPSA nonetheless predicts particle trajectories that are consistent with the experimentally observed behavior of published DLD chips sorting some types of deformable and non-spherical cells and droplets. However, undoubtedly there are particles that cannot be approximated as rigid circles and simulated successfully using MOPSA; for these particles other simulation approaches may be necessary and could even be incorporated into future versions of MOPSA (as discussed in the Conclusions).

## II. SIMULATION

Simulating a microfluidic chip using MOPSA requires two steps:

Calculate the fluid velocity field in the chip using conventional finite element analysis software

Calculate the trajectory followed by a particle through the chip using MOPSA

In this proof-of-concept demonstration, we used COMSOL Multiphysics to calculate the fluid velocity field and MATLAB to implement MOPSA. We used this approach to simulate three different microfluidic chip models:

The

*cross channel model*is a typical microfluidic cross-shaped intersection with 200*μ*m channel widths. Junctions like this are extremely common in microfluidic chips,^{10–13}so assessing the performance of MOPSA in these intersections is very important.The

*triangular DLD model*was chosen because of its complexity. Arrays of hundreds of triangular posts offer many opportunities for interactions between particles and the posts and challenge the robustness of MOPSA.The

*cylindrical DLD model*enables comparisons between MOPSA and published experimental results. The vast majority of DLD chips in the literature use arrays of cylindrical posts, and we use MOPSA to simulate the paths followed by particles in six published DLD experiments by three different research groups.^{7,14,15}

While we chose these three models for this proof-of-concept demonstration, we note that, in principle, MOPSA can simulate virtually any microfluidic chip design containing flowing particles, as long as the chip can be modeled in two dimensions (MOPSA is currently a two-dimensional simulation, although it could easily be extended to three dimensions as discussed in the Conclusions) and as long as the particles can be approximated as rigid circles or spheres. Finally, to compare MOPSA to existing commercial software tools for particle simulation in microfluidics, we also simulated these chip models entirely in COMSOL Multiphysics (without using MOPSA).

### A. Calculating fluid velocity field in COMSOL

For each microfluidic chip model simulated here, the fluid velocity fields were calculated using finite element analysis software (COMSOL Multiphysics). For accuracy, we used the *Laminar Flow* physics module and a customized free triangle mesh with a maximum mesh size that was equal to or less than the diameter of the smallest particle that would be simulated using that mesh. In the *cross channel model*, the maximum mesh size was 8 *μ*m (for simulating a particle diameter of 50 *μ*m). Inlets were assigned an inlet boundary condition of 5 mm/s normal inflow velocity. In the *triangular DLD models*, the maximum mesh size was 1 *μ*m (for simulating particle diameters of 1, 3, 10, and 12 *μ*m). In the *cylindrical DLD models*, the maximum mesh size was a quarter of the difference between two simulated particles' diameters (e.g., for simulating a DLD chip separating 8 and 9 *μ*m beads, the maximum mesh size was 0.25 *μ*m). Inlets were assigned an inlet boundary condition of 1 mm/s normal inflow velocity. In all models, outlets were assigned an outlet boundary condition of 0 Pa pressure, the remaining boundaries were walls (no-slip boundary condition), and the material filling the channels was water under incompressible flow. A stationary solver was used for calculating the fluid velocity field.

Require: Velocity profile in x and y direction of model (x.csv, y.csv), diameter of particle (D), initial position of particle (_{p}position), model geometry file (geometry.dxf), resolution of generating streamlines (_{i}resolution), vx, _{pre}vy, _{pre}β. |

Ensure: Simulated particle trajectory (trajectory) |

1: $polylines\u2190geometry.dxf$ ▷ Load microfluidic chip geometry from DXF file. |

2: $x\u2190x.csv$ ▷ Load pre-simulated fluid velocity profiles. |

3: $y\u2190y.csv$ |

4: $node\u2190x.csv$ ▷ Load mesh information from x.csv. |

5: node ← sortrows(_{sorted}node) ▷ Sort nodes based on x-direction. |

6: $geometryboundary[xmin\u2009xmax\u2009ymin\u2009ymax]\u2190[$min$(node(:,1))\u2003$ max$(node(:,1))\u2003$ min$(node(:,2)\u2003$ max$(node(:,2))]$ ▷ Define region of interest. |

7: $positionc\u2190positioni$ ▷ Initialize current position of particle. |

8: $trajectory\u2190positioni$ ▷ Initialize trajectory of particle. |

9: while $positionc\u2208geometryboundary$ do ▷ Move to next time step if particle is still in region of interest. |

10: $nodescovered\u2190$ findCoveredNodes$(Dp,positionc,nodesorted,node)$ ▷ Identify mesh nodes covered by particle. |

11: $vx\u2190$ getMean$(nodescovered,x)$ ▷ Set particle velocity to average fluid velocity of covered mesh nodes. |

12: $vy\u2190$ getMean$(nodescovered,y)$ |

13: $vy\u2190vy\xb7\beta $ ▷ Scale lateral (y) component of particle velocity by β if necessary. |

14: if v = 0 and _{x}v = 0 _{y}then ▷ If particle has zero velocity (is stuck)… |

15: $vx\u2190vxpre$ ▷…use the particle's velocity from the previous time step instead. |

16: $vy\u2190vypre$ |

17: end if |

18: $positionn\u2190positionc+resolution\xb7[vx\u2009vy]$ ▷ Calculate position of particle at next time step. |

19: $positionn\u2190$ wallEffect$(Dp/2,positionn,polylines)$ ▷ Apply wallEffect (see Algorithm 2). |

20: $trajectory\u2190[trajectory;positionn]$ |

21: $positionc\u2190positionn$ ▷ Record new particle position. |

22: $vxpre\u2190vx$ ▷ Record new “previous” particle velocity for use in next iteration. |

23: $vxpre\u2190vy$ |

24: end while |

Require: Velocity profile in x and y direction of model (x.csv, y.csv), diameter of particle (D), initial position of particle (_{p}position), model geometry file (geometry.dxf), resolution of generating streamlines (_{i}resolution), vx, _{pre}vy, _{pre}β. |

Ensure: Simulated particle trajectory (trajectory) |

1: $polylines\u2190geometry.dxf$ ▷ Load microfluidic chip geometry from DXF file. |

2: $x\u2190x.csv$ ▷ Load pre-simulated fluid velocity profiles. |

3: $y\u2190y.csv$ |

4: $node\u2190x.csv$ ▷ Load mesh information from x.csv. |

5: node ← sortrows(_{sorted}node) ▷ Sort nodes based on x-direction. |

6: $geometryboundary[xmin\u2009xmax\u2009ymin\u2009ymax]\u2190[$min$(node(:,1))\u2003$ max$(node(:,1))\u2003$ min$(node(:,2)\u2003$ max$(node(:,2))]$ ▷ Define region of interest. |

7: $positionc\u2190positioni$ ▷ Initialize current position of particle. |

8: $trajectory\u2190positioni$ ▷ Initialize trajectory of particle. |

9: while $positionc\u2208geometryboundary$ do ▷ Move to next time step if particle is still in region of interest. |

10: $nodescovered\u2190$ findCoveredNodes$(Dp,positionc,nodesorted,node)$ ▷ Identify mesh nodes covered by particle. |

11: $vx\u2190$ getMean$(nodescovered,x)$ ▷ Set particle velocity to average fluid velocity of covered mesh nodes. |

12: $vy\u2190$ getMean$(nodescovered,y)$ |

13: $vy\u2190vy\xb7\beta $ ▷ Scale lateral (y) component of particle velocity by β if necessary. |

14: if v = 0 and _{x}v = 0 _{y}then ▷ If particle has zero velocity (is stuck)… |

15: $vx\u2190vxpre$ ▷…use the particle's velocity from the previous time step instead. |

16: $vy\u2190vypre$ |

17: end if |

18: $positionn\u2190positionc+resolution\xb7[vx\u2009vy]$ ▷ Calculate position of particle at next time step. |

19: $positionn\u2190$ wallEffect$(Dp/2,positionn,polylines)$ ▷ Apply wallEffect (see Algorithm 2). |

20: $trajectory\u2190[trajectory;positionn]$ |

21: $positionc\u2190positionn$ ▷ Record new particle position. |

22: $vxpre\u2190vx$ ▷ Record new “previous” particle velocity for use in next iteration. |

23: $vxpre\u2190vy$ |

24: end while |

### B. Simulating particle trajectory with MOPSA

MOPSA is summarized in pseudocode in Algorithm 1. We implemented MOPSA in MATLAB to leverage MATLAB's strengths in matrix calculation. Using the built-in MATLAB function *INPOLYGON* significantly simplifies the MOPSA code. In principle, the pseudocode shown here could be implemented in any programming language.

MOPSA first imports the geometry of the microfluidic chip (the design of the chip in DXF format^{16}), the fluid velocity field calculated by COMSOL Multiphysics, and the description of the mesh used by COMSOL (lines 1–4 of Algorithm 1). Instead of treating the particle as a single mathematical point, MOPSA uses the average velocity of all of the mesh nodes covered by the particle to determine the particle's velocity profile (lines 11–12). *v _{x}* represents the magnitude of the velocity vector in the

*x*direction, and

*v*represents the velocity in the

_{y}*y*direction. As shown in Algorithm 1 line 13, MOPSA uses an empirically obtained parameter

*β*to “weight” the components of the particle velocity vector in direction orthogonal to fluid flow (the

*y*direction). For the

*cross channel model*,

*β*= 1 (meaning that the velocities calculated by COMSOL are used as-is). For the

*triangular*and

*cylindrical DLD*models,

*β*= 1.45. Details on how this value was determined are provided in Results and Discussion.

MOPSA next considers the possibility of particles with zero velocity. In theory, as long as there is fluid flow in a microfluidic channel, mesh nodes with zero fluid velocity should exist only at channel walls. However, in practice, a finite mesh size makes it possible for software like COMSOL Multiphysics to predict that nodes *near* channel walls may also have zero fluid velocity. A particle located at a mesh node with zero fluid velocity can become permanently trapped. To detect this situation, MOPSA utilizes a conditional statement (line 14 in Algorithm 1) to determine if the velocity of the particle in the time step will be zero. If so, MOPSA uses the particle's non-zero velocity from the *previous* time step for the current time step calculation (lines 15 and 16). As long as the amount of time between steps [*resolution*, defined in Eq. (1) below] is sufficiently small, this substitution does not seem to affect the accuracy of the simulation. In this manner, MOPSA ensures that the particle will not get permanently stuck. Alternatively, for applications in which particle sticking is desired (for example, using antibody-coated posts to capture cells of interest in a microfluidic chip^{17}), a threshold can be added to allow particle sticking below a user-specified particle velocity.

where *position _{c}* (

*μ*m) is the current position of the particle,

*position*(

_{n}*μ*m) is the position of the particle in the next time step,

*resolution*(s) is the amount of time between each calculation step, and

*v*and

_{x}*v*(

_{y}*μ*m/s) are the average fluid velocities of the mesh nodes covered by the particle, as defined earlier. The variable

*resolution*is a parameter that can be customized by the user. A smaller

*resolution*makes the simulation more accurate but also more time consuming. If the flow through a chip is relatively fast, then a small value for

*resolution*is needed to accurately predict the path followed by a particle in the fluid. If the flow is relatively slow, then a larger

*resolution*can be used without sacrificing simulation accuracy. To find an optimal value for

*resolution*, a MOPSA simulation can be repeated with successively smaller values for

*resolution*. Eventually the predicted particle trajectories will converge to a single trajectory that will not change with additional decreases in

*resolution*; this is the optimal value for

*resolution*for a given application.

Eq. (1) assumes that:

The density of the particle is the same as its surrounding fluid.

The net force applied to the particle is always zero.

Particles do not interfere with each other or channel walls.

Particles have smooth surfaces, so friction between fluid and particle surfaces is negligible.

The material of the particle is homogeneous.

If these assumptions are unsuitable for a particular application, MOPSA can be modified by replacing Eq. (1) with a more complicated drag law such as Stokes' Law^{19} or the work of Haider and Levenspiel^{20} (we demonstrate a modification of MOPSA later in this work).

Require: microfluidic chip wall geometry information (polylines), radius of particle (R), current position of particle (_{p}position), particle mapping parameter (_{c}k) |

Ensure: shifted position of particle (position) _{s} |

1: function wallEffect(R, _{p}position, _{c}polylines) |

2: $xp\u2190Rp\xb7\u2009cos\u2009(0:(2\pi /k):2\pi )+positionc(1,1)$ ▷ Calculate locations of the k points representing the particle boundary. |

3: $yp\u2190Rp\xb7\u2009sin\u2009(0:(2\pi /k):2\pi )+positionc(1,2)$ |

4: for $i=0\u2192totalpolylines$ do |

5: $polylinec\u2190polylines[i]$ ▷ Obtain the next polyline representing a channel feature (post, wall, etc.) |

6: if $polylinec(1,:)=polylinec(end,:)$ then ▷ Confirm polyline is closed. |

7: $xpolyline\u2190polylinec(:,1)$ |

8: $ypolyline\u2190polylinec(:,2)$ |

9: $in\u2190$ inpolygon$(xp,yp,xpolyline,ypolyline)$ ▷ Determine how many boundary points of the particle are inside the polyline. |

10: $countin\u21900$ |

11: $pointsin\u2190[]$ |

12: for $j=0\u2192totalin$ do |

13: if in(j) = 1 then |

14: $countin++$ |

15: $pointsin\u2190[pointsin;xp\u2009yp]$ |

16: end if |

17: end for |

18: if count > _{in}k or count = 0 _{in}then ▷ If the particle does not overlap with the wall… |

19: $positions\u2190positionc$ ▷…keep the current position of the particle. |

20: else ▷ However, if the particle does overlap with the wall: |

21: $pointmiddle\u2190[$mean$(pointsin(1,:),pointsin(end,:))]$ |

22: $direction\u2190$norm$(positionc\u2212pointmiddle)$ ▷…calculate the shift direction as a normal vector… |

23: $positions\u2190positionc+direction\xb7(countin/k)\xb7Rp\xb7w$ ▷…and apply the shift to the particle location. |

24: BREAK () |

25: end if |

26: end if |

27: end for |

28: return position ▷ Return the new particle location to the main MOPSA algorithm. _{s} |

29: end function |

Require: microfluidic chip wall geometry information (polylines), radius of particle (R), current position of particle (_{p}position), particle mapping parameter (_{c}k) |

Ensure: shifted position of particle (position) _{s} |

1: function wallEffect(R, _{p}position, _{c}polylines) |

2: $xp\u2190Rp\xb7\u2009cos\u2009(0:(2\pi /k):2\pi )+positionc(1,1)$ ▷ Calculate locations of the k points representing the particle boundary. |

3: $yp\u2190Rp\xb7\u2009sin\u2009(0:(2\pi /k):2\pi )+positionc(1,2)$ |

4: for $i=0\u2192totalpolylines$ do |

5: $polylinec\u2190polylines[i]$ ▷ Obtain the next polyline representing a channel feature (post, wall, etc.) |

6: if $polylinec(1,:)=polylinec(end,:)$ then ▷ Confirm polyline is closed. |

7: $xpolyline\u2190polylinec(:,1)$ |

8: $ypolyline\u2190polylinec(:,2)$ |

9: $in\u2190$ inpolygon$(xp,yp,xpolyline,ypolyline)$ ▷ Determine how many boundary points of the particle are inside the polyline. |

10: $countin\u21900$ |

11: $pointsin\u2190[]$ |

12: for $j=0\u2192totalin$ do |

13: if in(j) = 1 then |

14: $countin++$ |

15: $pointsin\u2190[pointsin;xp\u2009yp]$ |

16: end if |

17: end for |

18: if count > _{in}k or count = 0 _{in}then ▷ If the particle does not overlap with the wall… |

19: $positions\u2190positionc$ ▷…keep the current position of the particle. |

20: else ▷ However, if the particle does overlap with the wall: |

21: $pointmiddle\u2190[$mean$(pointsin(1,:),pointsin(end,:))]$ |

22: $direction\u2190$norm$(positionc\u2212pointmiddle)$ ▷…calculate the shift direction as a normal vector… |

23: $positions\u2190positionc+direction\xb7(countin/k)\xb7Rp\xb7w$ ▷…and apply the shift to the particle location. |

24: BREAK () |

25: end if |

26: end if |

27: end for |

28: return position ▷ Return the new particle location to the main MOPSA algorithm. _{s} |

29: end function |

After the particle's position at the next time step has been calculated, MOPSA's *wallEffect* function (Algorithm 2) determines if the particle will collide or overlap with solid walls. As shown in Fig. 2, *wallEffect* uses a large number of points (*k* = 200) to map the boundary of the particle located at position *a*. If any of these *k* points are located inside a solid wall (e.g., a triangular post, a channel wall, or any other device feature), then the algorithm recognizes that the particle is overlapping the wall and compensates by shifting the particle's position away from the wall to a new location $a\u2032$. The direction of this shift is determined by first identifying the two points of contact between the particle perimeter and the wall (*b* and *c*), then defining a line segment *bc* between these two points, and finally computing the midpoint *d* of the line segment. The direction of the particle shift is perpendicular to line *bc* (roughly pointing in the opposite direction of the fluid flow) and is calculated using Eq. (2)

where *position _{s}* is the shifted position of the particle,

*position*is the initial (wall-overlapping) position of the particle,

_{c}*direction*is a normal vector describing the shifting direction,

*k*is the total number of points used to map the perimeter of the particle,

*count*is the number of those perimeter mapping points located inside the solid wall,

_{in}*R*is the radius of the particle, and

_{p}*w*is a weight factor that can be adjusted for different models (

*w*= 1 was used in this work). Once this algorithm has corrected the overlapping position of the particle, MOPSA moves on to the next simulation time step.

The parameters *resolution* in Eq. (1) and *w* in Eq. (2) provide a mechanism for fine-tuning the *wallEffect* algorithm. For example, if the algorithm predicts an unrealistically large shift distance, reducing *resolution* will decrease the amount of overlap between the particle and wall (and thus decrease the amount of correction applied by the *wallEffect*. If the algorithm predicts an unrealistically small shift distance, increasing *w* will increase the shift distance. Finally, if the geometry of a wall has concave features, *wallEffect* may yield unexpected results because the particle perimeter could intersect with the wall in more than two points. In this case, changing *w* could compensate for unexpected results from concave wall geometries. In any case, *w* = 1 was used for all of the simulations presented here.

### C. Comparison with existing commercial particle tracers

To compare the performance of MOPSA to an existing commercial particle tracer, we repeated all of our particle simulations using the built-in *Particle Tracing for Fluid Flow* physics module in COMSOL Multiphysics. The simulations used a time-dependent solver, and a “drag force” boundary was added into the physics to use Stokes' Law for particle trajectory calculation. For the *triangular DLD model*, the inlet was assigned an inlet boundary condition of 100 particles per release with a uniform distribution. The particles were assigned diameters of 1, 3, 10, or 12 *μ*m (although, as described above, these diameters do not seem to affect the simulation results in COMSOL Multiphysics). The outlet was assigned an outlet boundary condition of *freeze wall*. The six *cylindrical DLD model* chips we simulated used the same conditions as the *triangular DLD model* but with 10 particles per release and the particle diameters shown in Table I. We also tested a range of boundary conditions for the post walls in the DLD models, including *bounce, diffuse scattering*, and *general reflection*. Finally, for the *cross channel model*, the inlet was assigned an inlet boundary condition of 10 particles per release (50 *μ*m particle diameter) with uniform distribution and the outlet was assigned an outlet boundary condition of *freeze wall*.

Figure . | λ (μm)
. | Δλ (μm)
. | G (μm)
. | Sorted particle types and sizes . | Figure in reference . |
---|---|---|---|---|---|

7(a) | 80 | 6 | 25 | 8 μm and 9 μm beads | Fig. 2(a) in Li et al.^{14} |

7(b) | 90 | 8 | 25 | 9 μm and 10 μm beads | Fig. 2(b) in Li et al.^{14} |

7(c) | 80 | 2 | 20 | 2 μm platelets and 6 μm red blood cells | Fig. 4(a) in Li et al.^{14} |

7(d) | 80 | 5 | 20 | 6 μm red blood cells and 10 μm white blood cells | Fig. 4(b) in Li et al.^{14} |

7(e) | 8 | 0.8 | 1.6 | 0.4 μm and 1.0 μm beads | Fig. 2(a) in Huang et al.^{7} |

7(f) | 120 | 6 | 60 | 11 μm and 30 μm droplets | Fig. 2(c) in Joensson et al.^{15} |

Figure . | λ (μm)
. | Δλ (μm)
. | G (μm)
. | Sorted particle types and sizes . | Figure in reference . |
---|---|---|---|---|---|

7(a) | 80 | 6 | 25 | 8 μm and 9 μm beads | Fig. 2(a) in Li et al.^{14} |

7(b) | 90 | 8 | 25 | 9 μm and 10 μm beads | Fig. 2(b) in Li et al.^{14} |

7(c) | 80 | 2 | 20 | 2 μm platelets and 6 μm red blood cells | Fig. 4(a) in Li et al.^{14} |

7(d) | 80 | 5 | 20 | 6 μm red blood cells and 10 μm white blood cells | Fig. 4(b) in Li et al.^{14} |

7(e) | 8 | 0.8 | 1.6 | 0.4 μm and 1.0 μm beads | Fig. 2(a) in Huang et al.^{7} |

7(f) | 120 | 6 | 60 | 11 μm and 30 μm droplets | Fig. 2(c) in Joensson et al.^{15} |

## III. RESULTS AND DISCUSSION

### A. Simulating the cross channel model

The fluid velocity field of the *cross channel model* was solved using COMSOL Multiphysics and is shown in Fig. 3(a). This model simulates a cross-channel with fluid entering from the left and exiting out the top, right, and bottom. For rigid spherical particles originating at many locations across the left-hand entrance channel, COMSOL's particle tracer and MOPSA predict very similar particle trajectories. However, for particles that start close to a channel wall, MOPSA's simulation is more realistic than COMSOL's. Figure 3(b) shows trajectories predicted by COMSOL (red) and MOPSA (green) for a 50 *μ*m diameter particle flowing adjacent to the channel wall. In the trajectory predicted by COMSOL, the particle overlaps with the channel wall as it goes around the corner. In contrast, MOPSA detects these channel walls and keeps the particle from overlapping them. The realistic wall interactions simulated by MOPSA result in a particle trajectory that is 18 *μ*m (36% of the particle's size) farther to the right than the trajectory calculated by COMSOL. For simple microfluidic chips with a single intersection, the particle tracing errors introduced by software like COMSOL Multiphysics may be acceptable and may not affect the overall accuracy of the simulation. However, for microfludic chips with hundreds of intersections in parallel (like the hydrodynamic filter of Yamada *et al.*^{10}), small errors in each intersection could combine to form a large net error associated with the path followed by particles flowing across the entire chip. In these cases, MOPSA should provide a much more accurate particle trajectory than existing software tools.

### B. Simulating deterministic lateral displacement (DLD)

We also tested MOPSA on the *triangular deterministic lateral displacement (DLD) model* cell sorter to determine how the algorithm performs in a chip with more complex geometry. Figure 4 shows the predicted particle trajectories for a 10 *μ*m diameter rigid spherical particle flowing through an array of triangular posts, calculated using COMSOL Multiphysics (red) and MOPSA (green). In both simulations, the particle starts at the same location. Since COMSOL treats the particle as a mathematical point, COMSOL's trajectory allows the particle to overlap with a triangular post, which is physically impossible. Additionally, later in the COMSOL simulation the particle encountered a zero-velocity mesh node and permanently stopped. In contrast, the trajectory calculated by MOPSA is much more consistent with the known mechanism of DLD. Repeated interactions with triangular posts laterally displace the particle, and no overlapping or sticking with posts is observed (even after repeating the simulation shown in Fig. 4 hundreds of times with different starting locations for the particle). The design of this chip in DXF format, fluid velocity field, and particle trajectory data are all available for download in online supplementary material.

The goal of DLD is to separate different cells and other particles based on their size: smaller particles follow the fluid streamlines straight through the device, while larger particles are “bumped” laterally by the posts and follow a diagonal path through the device. To confirm that MOPSA can simulate the particle sorting capabilities of a DLD chip, in Fig. 5 we used it to simulate the trajectories of 1, 10, and 12 *μ*m diameter rigid spherical particles flowing through an array of triangular posts. Although particles of all three sizes start in the same location, MOPSA predicts that the particles will exit the chip in different locations: the small 1 *μ*m particle follows the fluid streamlines and moves straight through the chip, and the large 12 *μ*m particle is “bumped” by each post and moves diagonally through the chip. The mid-sized 10 *μ*m particle initially follows a diagonal path but then switches to a straight-through path, leading to a trajectory that lies between the other two particles' trajectories. This suggests that the critical diameter for this DLD chip design (the size of the smallest particle that is laterally displaced in the chip) is close to 10 *μ*m. In contrast, particles of different sizes follow exactly the same trajectory in COMSOL's particle tracer (Fig. 1). These results confirm that MOPSA can successfully simulate size-based DLD separations that are difficult or impossible to simulate using existing commercial software, at least for rigid spherical particles. The design of this chip in DXF format, fluid velocity field, and particle trajectory data are all available for download in online supplementary material.

### C. Simulating published DLD experiments

We next used MOPSA to predict the paths followed by different sized particles in several published DLD chips. We identified six experiments in three published papers^{7,14,15} in which the chip designs and experimental conditions are described in sufficient detail to allow us to simulate them using MOPSA. These chips use arrays of cylindrical posts with dimensions defined in Fig. 6 and summarized in Table I.

We first simulated the DLD chip described by Li *et al.*^{14} who used this chip to sort 8 *μ*m and 9 *μ*m diameter rigid spherical beads. We found that with a value of *β* = 1.45, the MOPSA predicts that the 8 *μ*m and 9 *μ*m beads will follow different paths through the DLD chip [Fig. 7(a)]; this is consistent with the observed bead sorting behavior of this DLD chip.^{14}

We then left the value of *β* = 1.45 unchanged as we used MOPSA to predict the paths followed by particles through five additional published DLD chips from three research groups [Figs. 7(b)–7(f)]. In their original publications,^{7,14,15} these chips were used to sort a wide variety of particle sizes (from 0.4 *μ*m to 30 *μ*m diameters) and types (beads, droplets, red blood cells, white blood cells, and platelets). Some of these particles are not rigid or spherical—droplets are not rigid, and red blood cells, white blood cells, and platelets are neither rigid nor spherical—so one may rightly expect that MOPSA (which assumes rigid spherical particles) may not be able to accurately simulate the trajectories of these particles. However, MOPSA nonetheless successfully predicted that the different-sized particles follow different paths through the chip. This prediction is consistent with the particle-sorting behavior reported by the creators of these chips. In contrast, the built-in particle tracer in COMSOL Multiphysics predicted no differences in the paths followed by the different particles (see supplementary material Figs. 2–7). Undoubtedly there are some particles whose deformability and non-spherical shapes will render MOPSA unable to accurately simulate the behavior of these particles, but at least for the particle types and sizes we studied, MOPSA's assumption of rigid and spherical particles did not seem to adversely affect the quality of the simulation results. The chip design DXF files, fluid velocity profiles, and predicted particle trajectories from both MOPSA and COMSOL for each of the six simulations in Fig. 7 are available in supplementary material.

In some of the DLD chip simulations in Fig. 7, MOPSA predicted that the smaller particles will follow a somewhat jagged path when the particles flow vertically between two rows of posts. For example, in Fig. 7(c) the 2 *μ*m diameter platelets seem to make two sharp turns as they pass from above the third row of posts to below the row. This behavior was unexpected; to determine its source, we examined the fluid velocity fields calculated by COMSOL for these DLD simulations. Close inspection of the gaps between columns of posts (supplementary material Fig. 8) reveals that COMSOL predicted asymmetric fluid velocity fields in these gaps in the cases where the MOPSA-predicted particle paths were most jagged. For example, the fluid velocity field calculated by COMOSL for the chip in Fig. 7(c) [shown in supplementary material Fig. 8(c)] has a diagonal region of low flow between every two pillars; when MOPSA uses this fluid velocity field to predict the path followed by a particle between the pillars, the asymmetry in the field results in a jagged predicted particle path. These particle trajectories do not seem to affect the overall simulation results (that different-sized particles follow different paths through the DLD chips).

Why was it necessary to increase *β* to successfully simulate DLD chips in MOPSA? One explanation could be our assumption that the presence of the particle will not affect the fluid velocity profile in the chip. In actuality, the presence of a particle *could* increase the hydrodynamic resistance of the DLD chip.^{21} This could reduce the fluid velocity in the direction of flow (the *x* direction) and make the fluid velocity in the direction perpendicular to flow (the *y* direction) more significant. Increasing *β* from 1 to 1.45 adds additional weight to the *y* component of the fluid velocity and enables MOPSA to accurately simulate particles in DLD devices, at least for the particle sizes and types considered here. While this value of *β* may not be suitable for all DLD devices or particle sizes and types, it is noteworthy that all six previously published DLD chips we replicated in this work were successfully simulated using the same value for *β* (1.45).

### D. Extending MOPSA to simulate particles with different densities

In its current form, MOPSA is capable of simulating rigid spherical particles in a variety of different microfluidic applications. However, in some cases it may be necessary to add additional physics modules to MOPSA to simulate certain particle properties. For example, microfluidic chips have been demonstrated that sort cells and other particles by their densities.^{22,23} Since the current MOPSA assumes that all particles have the same density (and the particles' density equals the density of the fluid around the particles), MOPSA cannot currently be used to simulate the behavior of these density sorter chips.

To demonstrate that additional particle properties can be readily added to MOPSA when needed, we added support for particle density to MOPSA. Our approach uses the work of Haider and Levenspiel^{20} to predict the effects of buoyancy on particles. We then used our modified MOPSA to predict the paths followed by rigid spherical particles of four different densities (0.8, 1.0, 1.2, and 1.5 g/cm^{3}) in a microfluidic chip filled with water (1.0 g/cm^{3}) and oriented on its edge relative to Earth's gravity. In Fig. 8, all four particles start at the same location in a lower inlet. As the particles flow along the horizontal channel, our modified MOPSA correctly predicts that the particles with density greater than the fluid density will sink downward and exit via the lower outlet, and the particles with density greater than the fluid density will float upward and exit via the upper outlet. Since different cell types often have different densities,^{24} our modified MOPSA may be used to help design chips for sorting cells by their type, an important capability in biological research and medical diagnostics.

## IV. CONCLUSIONS

In this work, we presented MOPSA, an algorithm for optimized particle simulation in microfluidic chips. By treating particles as two-dimensional objects instead of single points and applying corrections when particles interact with channel walls, MOPSA can accurately simulate particle behaviors that particle tracers in existing commercial software tools cannot. Consequently, adopting MOPSA into the design process for microfluidic chips that contain cells, droplets, and other particles should be beneficial.

### A. Limitations of MOPSA and future directions

MOPSA makes several assumptions about the particles it simulates: the particles are circular or spherical, they have the same density as the fluid around them, they have smooth surfaces, they do not interact with other particles, and so on. These assumptions may limit the utility of MOPSA for some applications. For example, in its current form, MOPSA cannot simulate effects of inertia on either the particles or the fluid in a chip; these effects have been used as the basis for particle sorting and focusing in microfluidics.^{25} However, for at least some particle sizes and types, MOPSA's assumption of rigid spheres seems to still yield results that are consistent with experimental observations.

For particles that cannot be assumed to be rigid spheres, a number of modeling techniques have been proposed for simulating the behavior of these particles.^{26} For example, dissipative particle dynamics^{27–29} has been successfully applied to simulate the deformation of red blood cells in DLD chips.^{26,30} Zhu *et al.*^{31} modeled deformable cells as fluid-filled capsules enclosed by neo-Hookean membranes.^{32} Krüger *et al.*^{33} used the immersed-boundary method^{34} to model membrane dynamics during cell deformation. These and other models can be used in situations where MOPSA fails to accurately predict the behavior of a non-rigid, non-spherical particle; some of them could even be incorporated into future versions of MOPSA.

We also demonstrated that MOPSA can be easily extended to simulate other particle properties like particle density. This same approach can be used to enable MOPSA to support other important properties of particles. For example, by altering the *wallEffect* algorithm to allow cells to stick to channel walls, one could simulate the behavior of microfluidic devices that intentionally capture cells in this manner.^{17} Additionally, to include the effects of other forces such as electrostatics and acoustics on a particle's trajectory, additional equations can be added to Line 18 in Algorithm 1. To simulate particle-particle interactions that may be important in applications with high particle concentrations, it may be necessary to periodically recalculate the fluid velocity vector field while including particle positions (and particle-induced drag) in the finite element analysis calculation. In our experience, including particle positions in this manner has a negligible effect on simulation accuracy but significantly increases the required computational time, but it may be necessary for some applications.

Finally, MOPSA is currently a two-dimensional simulation technique and cannot be used to simulate three-dimensional microfluidic channel networks. However, nothing about MOPSA precludes extending the algorithm to three dimensions. This would enable MOPSA to simulate particle trajectories in emerging 3D-printed microfluidic chips that can have channels in all three spatial dimensions.^{35}

## V. SUPPLEMENTARY MATERIAL

See supplementary material for fluid velocity field for the simulations in Figs. 4, 5, and 7(a)–7(f); results from using COMSOL Multiphysics to simulate the deterministic lateral displacement (DLD) chips in Figs. 7(a)–7(f); closeups of the calculated fluids velocity fields for the DLD chips in Figs. 7(a)–7(f); calculated particle trajectories from MOPSA for the simulations in Figs. 4, 5, and 7(a)–7(f) in MATLAB *MAT* format; and CAD files for the chip designs in Figs. 4, 5, and 7(a)–7(f) in *DXF* format.

## ACKNOWLEDGMENTS

This work was supported in part by the National Science Foundation Division of Biological Infrastructure under Award No. DBI-1353974, the National Science Foundation Division of Computer and Communication Foundations program under Award No. CCF-1351115, and the National Science Foundation Division of Industrial Innovation and Partnerships under Award No. IIP-1640757.