Flow topology during multiplexed particle manipulation using a Stokes Trap

TL;DR

This paper investigates the flow structures involved in multiplexed particle manipulation using a Stokes trap, revealing that optimal control depends on flow patterns with minimal stagnation points, which aids in precise particle control.

Contribution

It provides a detailed characterization of flow topologies during active particle control in a Stokes trap, highlighting the relationship between control algorithms and flow structures.

Findings

Optimal control relies on flow patterns with zero or one stagnation point.

System sensitivity analysis clarifies the influence of control parameters.

Results guide applications in colloidal particle manipulation and soft particle deformation.

Abstract

Trapping and manipulation of small particles underlies many scientific and technological applications. Recently, the precise manipulation of multiple small particles was demonstrated using a Stokes trap that relies only on fluid flow without the need for optical or electric fields. Active flow control generates complex flow topologies around suspended particles during the trapping process, yet the relationship between the control algorithm and flow structure is not well understood. In this work, we characterize the flow topology during active control of particle trajectories using a Stokes trap. Our results show that optimal control of two particles unexpectedly relies on flow patterns with zero or one stagnation points, as opposed to positioning two particles using two distinct stagnation points. We characterize the sensitivity of the system with respect to the parameters in the control objective function, thereby providing a systematic understanding of the trapping process. Overall, these results will be useful in guiding applications involving the controlled manipulation of multiple colloidal particles and the precise deformation of soft particles in defined flow fields.

Figures

Equations

$${\mathbf{R}}^i$$

$$\mathbf{u}=\frac{1}{\pi H}\sum _{i=1}^6\frac{(\mathbf{x}-{\mathbf{R}}^i)q_i}{\Vert \mathbf{x}-{\mathbf{R}}^i{\Vert }^2}$$

$$\sum _{i=1}^6{q}_i=0$$

$$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},\mathbf{q},\mathbf{R})\triangleq \frac{1}{\pi H}\sum _{i=1}^6\frac{(\mathbf{x}-{\mathbf{R}}^i)q_i}{\Vert \mathbf{x}-{\mathbf{R}}^i{\Vert }^2}$$

$$\mathcal{J}$$

$$+\gamma \left(\Vert {\tilde{\mathbf{X}}}_K-{\mathbf{X}}_F{\Vert }^2\right)$$

$$\frac{d\tilde{\mathbf{X}}}{dt}=[{\mathbf{f}}^T({\tilde{\mathbf{x}}}_1,{\tilde{\mathbf{q}}}_k,\mathbf{R}),{\mathbf{f}}^T({\tilde{\mathbf{x}}}_2,{\tilde{\mathbf{q}}}_k,\mathbf{R}){]}^T,{\tilde{\mathbf{X}}}_k={\mathbf{X}}_0$$

$$\sum _{i=1}^6{\tilde{q}}_{k,i}=0\forall k=1,\dots ,K$$

$$P_i=P_0+q_i{Z}_i$$

Full text

In recent years, particle trapping methods have been used to study physical and biological phenomena with exquisite levels of precision. Particle trapping has been used for diverse applications including parallel manipulation of single particles Chiou et al. (2005), determination of DNA physical properties Wang et al. (1997), and direct observation of viral DNA packaging by a molecular motor Chemla et al. (2005). Trapping methods rely on a variety of physical mechanisms and fundamental forces for manipulating particles, including optical fields Ashkin et al. (1986); Grier (2003), magnetic fields Gosse and Croquette (2002); Sarkar and Rybenkov (2016), electrical fields Cohen and Moerner (2006); Ropp et al. (2010), acoustic forces Hertz (1995); Guo et al. (2016), and fluidic forces Petit et al. (2011); Tanyeri et al. (2011); Tanyeri and Schroeder (2013); Shenoy et al. (2015, 2016); Kumar et al. (2019).

One of the most fundamental experiments in particle trapping focuses on measuring particle deformation in response to a controlled force Gosse and Croquette (2002); Wang et al. (1997). Flow-based trapping has long been used to study the deformation of suspended particles such as immiscible droplets. G.I. Taylor developed a four-roll mill for studying drop deformation in extensional flows Taylor (1934), though this method required manual operation for trapping drops in flow. Several years later, Bentley and Leal Bentley and Leal (1986) developed a computer controlled version of the four-roll mill that enables the trapping of freely suspended drops for long times. Several researchers further developed microfluidic analogs of the four-roll mill to generate well-defined shear, extensional, and mixed flows Hudson et al. (2004); Lee et al. (2007). The microfluidic four-roll mill was used to observe the deformation and dynamics of suspended particles, such as tank-treading, tumbling, and trembling dynamics of single vesicles in flow Deschamps et al. (2009). However, these microfluidic devices did not incorporate active flow control to confine particles over long time periods.

In 2011, Schroeder and coworkers developed an automated hydrodynamic trap to confine single particles in free solution Tanyeri et al. (2011); Tanyeri and Schroeder (2013); Shenoy et al. (2015). This method relies on active flow control using an on-chip metering valve to modulate flow rates in a 4-channel cross-slot device. The initial version of the hydrodynamic trap offers a simple design with a linear control algorithm Tanyeri et al. (2011), thereby producing a single stagnation point to confine individual particles in a planar extensional flow. Using this approach, the stagnation point is translated over micron-scale distances using an on-chip valve without affecting the principal axes of compression and extension, thereby generating a well-defined flow field that enables estimation of hydrodynamic forces on single suspended objects such as single polymer molecules Schroeder (2018). In 2016, the Stokes trap was developed as a general method to confine and manipulate multiple particles using the sole action of fluid flow Shenoy et al. (2016). This technique utilizes a model predictive control formulation to modulate flow rates in an $N$-channel intersecting cross-slot device to achieve simultaneous center-of-mass trapping and precise manipulation of multiple particles. The Stokes trap was initially demonstrated with an $N$=6 channel device to manipulate two distinct micron-scale particles along pre-determined paths and for the directed assembly of two sticky particles using scheduled fluid flows Shenoy et al. (2016).

Although the Stokes trap enables manipulation of multiple particles, the flow structure and evolution of flow topology surrounding the suspended particles during the trapping process is not well understood. In theory, a six-channel microdevice is capable of generating two stagnation points (for trapping two distinct particles), though it is not known if the flow topology during particle trapping admits two stagnation points. Moreover, the evolution and location of putative stagnation points and streamlines during fluidic trapping experiments is generally not known a priori. The relation between the control algorithm used to determine optimal manipulation pathways and particle dynamics is not obvious, and the resulting particle trajectories may be complex depending on the trapping objective. From this perspective, achieving a clear understanding of the underlying flow topology would be useful for characterizing optimal manipulation techniques. For example, it is known that flow structures can interact with and influence the topology of molecular orientation fields in nematic systems Giomi et al. (2017). Moreover, adhesive interactions between freely suspended vesicles in flow are known to depend on the stress profile induced on bilayer membrane by external flows. Consider the case of two vesicles that are on a trajectory for close approach or collision. Approaching vesicles experience compressional stresses, resulting in wrinkling of membrane surfaces, which reduces the effective contact area of adhesion, delays the film drainage process, and increases the time required for adhesion Kantsler et al. (2007); Turitsyn and Vergeles (2008); Narsimhan et al. (2015). From this view, achieving a systematic understanding of how the underlying flow topologies evolve during particle manipulation and interaction experiments will be essential for quantitative descriptions of such processes.

In this work, we systematically study flow topologies during particle manipulation experiments using a Stokes trap. In particular, we characterize flow topology during three canonical manipulation scenarios for two particles, including the movement of two particles towards each other, away from each other, and the interchange of particle positions. For each case, particle trajectories are tracked using active flow control experiments and directly compared to simulation results. Our results show that the experimentally determined flow rates are in good agreement with simulation results for several different manipulation scenarios. In addition, we systematically vary the control parameters in the objective function to understand the influence of the control algorithm on the time required for the particles to reach their target positions. Together, these studies elucidate the flow topology during particle manipulation and reveal the impact of the tunable control parameters on the trapping performance. Overall, these findings will help guide multiplexed flow-based trapping experiments involving particle collision, adhesion, and assembly in a systematic fashion.

I Controller formulation

The flow model and the model predictive control (MPC) scheme for multiplexed particle manipulation have been previously described Shenoy et al. (2016); Kumar et al. (2019). The following section summarizes the salient features of these models for understanding flow topology during trapping.

I.1 Flow model for the particle center-of-mass.

Consider the six-channel cross-slot device shown in Fig. 1a. The intersection of the six channels creates a hexagonal cross-slot as shown in Fig. 1b, with the vertices of the hexagon lying on a circle of radius $R$ which is equal to the width $W$ of a channel. The center of the circle is chosen to be the origin. Inlet/outlet channels in the cross-slot are assumed to be point sources/sinks Schneider et al. (2011), where $\bm{R}^{i}$ is the location of the ith source/sink, as shown schematically in Fig. 1b.

[TABLE]

Under low Reynolds number flow conditions, the height averaged 2-D velocity $\bm{u}\in\mathbb{R}^{2}$ within the cross-slot can be expressed as a linear superposition of the velocity fields caused by point sources using a 2-D Hele-Shaw approximation Leal (2007):

[TABLE]

where $H$ is the height of the device, $\bm{x}\in\mathbb{R}^{2}$ is a position vector, $\bm{R}^{i}\in\mathbb{R}^{2}$ is the position vector corresponding to the $i^{th}$ point source/sink, and $\bm{q}\in\mathbb{R}^{6}$ is a vector whose $i^{th}$ element $q_{i}$ represents the volumetric flow rate through the $i^{th}$ point source. Prior work used computational fluid dynamics (CFD) simulations to demonstrate that proper selection of the dimensions of the microfluidic device ($H$ and $W$) can reduce the relative error between the velocity predicted by the model in Eq. 2 and the CFD velocity to as low as 2% Shenoy et al. (2016). The microfluidic cross-slot is a closed system, and hence the flow rates $\bm{q}$ must satisfy mass conservation:

[TABLE]

Flow rates for a point source are taken to be positive when the fluid in the corresponding channel flows into the cross-slot and negative when the fluid in the channel leaves the cross-slot. Together, Eq. 2 and Eq. 3 determine the relationship between the imposed flow rates in the channels and the velocity field within the cross-slot. In the absence of external forces, a particle’s 2-D velocity ($\dot{\bm{x}}\in\mathbb{R}^{2}$) is identical to the fluid velocity at the particle center-of-mass (COM) ($\bm{x}\in\mathbb{R}^{2}$):

[TABLE]

where the $\triangleq$ symbol denotes a definition of the function.

I.2 Model predictive control for particle manipulation.

Using the linear velocity model for particle motion in Eq. 4, we specify a control scheme for manipulating two particles from an initial position $\bm{X}_{0}=[\bm{x}^{T}_{1,0},\bm{x}^{T}_{2,0}]^{T}$ to a final position $\bm{X}_{F}=[\bm{x}^{T}_{1,F},\bm{x}^{T}_{2,F}]^{T}$, where $\bm{x}_{i,0}$ indicates the $i^{th}$ particle’s initial position and $\bm{x}_{i,F}$ indicates the target final position. A model predictive control (MPC) scheme is used to achieve precise manipulation of multiple particles Mayne et al. (2000); Shenoy et al. (2016); Kumar et al. (2019). Using MPC, the flow rates required for particle manipulation are obtained by minimizing an objective function that determines an optimal balance between: (1) trajectories that move particles along the shortest path between the initial and final positions and (2) trajectories that require minimal flow rates but may require complicated paths. For the current sampling instant at time $t_{0}$, particle positions $\bm{X}_{0}$ are experimentally determined. Next, a future time horizon $T$ is defined and divided into $K$ equal intervals, such that $t_{k}=t_{0}+k\Delta,\ k=1\ldots K$, and $\bm{q}_{k}$ indicates the flow rates to be applied during $k^{th}$ interval $[t_{k-1},t_{k})$. Given a series of flow rates throughout the time horizon $[\tilde{\bm{q}}_{1},\tilde{\bm{q}}_{2},\ldots,\tilde{\bm{q}}_{K}]$, the locations of all particles can be predicted subject to the imposed flow rates at the future $K$ sampling points. Future particle positions are denoted as $\tilde{\bm{X}}_{k}=[\tilde{\bm{x}}^{T}_{1,k},\tilde{\bm{x}}^{T}_{2,k}]^{T}$.

Using the MPC formulation, it is desired to obtain a set of flow rates that simultaneously minimize the distance travelled by each particle and the flow rates during the trapping event. To accomplish this task, an objective function $\mathcal{J}$ is minimized with respect to the particle trajectories $\{\tilde{\bm{X}}_{k}\}^{K}_{k=1}$ and flow rates $\{\tilde{\bm{q}}_{k}\}^{K}_{k=1}$:

[TABLE]

where $\beta$ and $\gamma$ are control parameters or weights that penalize high flow rates and large deviations at the end of the time horizon. In general, $\beta$ and $\gamma$ are tuned by the user to achieve the desired speed of response during trapping experiments. The objective function in Eq. 5 is minimized at each sampling instant to obtain the set of flow rates for the horizon $[\tilde{\bm{q}}_{1},\tilde{\bm{q}}_{2},\ldots,\tilde{\bm{q}}_{K}]$, but only the flow rates corresponding to the first interval $\tilde{\bm{q}}_{1}$ are applied to the system. In the next sampling instant, the particle positions are resampled again and the process is repeated. The objective function defined by Eq. 5 relies on a non-linear model, which complicates the determination of a closed-form solution. For experimental implementation of the MPC scheme, we use the toolkit for Automatic Control and Dynamic Optimization (ACADO) Houska et al. (2011); Quirynen et al. (2015), which efficiently minimizes the objective function using numerical optimization methods and generates high performance C++ code to enable millisecond time solutions of Eq. 5. Together, this approach enables real-time implementation of the non-linear model in optical microscopy and imaging experiments.

II Methods

II.1 Experimental setup.

Standard techniques in soft-lithography are used to fabricate the microfluidic device shown in Fig. 1 (Supplementary Materials) Qin et al. (2010). Microfluidic devices are mounted on the stage of an inverted microscope as shown in Fig. 2. The flow setup consists of 6 pressure regulators (MPV Series, Proportion Air) connected to a nitrogen source with an output pressure of 30 psi. The pressure regulators accept an analog input voltage ranging from 0-10 V, provide a linearly dependent output of 0-5 psi, and can measure the output pressure through an analog output. The regulators are actuated through a computer-controlled LabVIEW data acquisition device with analog input and output capabilities (NI 9264, NI9205 and NI9174, National Instruments). Each output pressure line from the regulators is connected to a fluidic reservoir containing the sample (Elveflow XS fluid reservoir), which is connected to the microfluidic device using FEP and PEEK tubing (IDEX Health and Science). Tube dimensions (inner diameter and length) are chosen such the fluidic resistance within the tubing is significantly larger than the fluidic resistance of the microfluidic channels. The microfluidic device is then mounted on an inverted microscope (Olympus IX71) equipped with a mercury lamp for fluorescent imaging, a 10x magnification objective lens, and a charge-coupled device (CCD) camera (PointGrey GS3 23S6M). Unless otherwise noted, all experiments, were performed using fluorescent polystyrene beads (2.2 $\mu$m, Nile Red, Spherotech) suspended in an aqueous glycerol solution with a viscosity of 0.140 Pa-s measured at 25 oC. The microfluidic device is situated such that the center of the cross-slot is defined as the origin in the laboratory reference frame.

II.2 Controller implementation.

Particle manipulation experiments are performed using a custom code written in MATLAB (Mathworks) and LabVIEW (National Instruments) Shenoy et al. (2016); Kumar et al. (2019). Briefly, at each sampling instant, LabVIEW acquires an image (single frame) from the CCD camera and performs binary thresholding, followed by particle localization to determine the centroid positions of particles within a region of interest defined by the user. Within a region of interest, the particles nearest to the two set points are selected, and the center-of-mass (COM) coordinates of these particles are transformed from the camera reference frame to the laboratory reference frame and scaled to dimensionless quantities Kumar et al. (2019). Scaling is used to avoid numerical underflow; in brief, particle coordinates are on the order of few hundred microns and flow rates are on the order of a few $\mu L/hr$, and determining the minima for the objective function in Eq. 5 can be strongly influenced by numerical error in the optimization algorithm. In this way, non-dimensionalization scales all quantities to be approximately $\mathcal{O}(1)$. All lengths are scaled by the particle diameter $d=2.2\ \mu$m, the time is scaled by a characteristic inverse strain rate $\dot{\epsilon}^{-1}=t_{s}=1$. Thus, Eq. 5 is solved using the dimensionless particle coordinates to obtain the dimensionless flow rates to be imposed on each channel.

In order to calculate the corresponding pressures for each desired flow rate, the fluidic setup is converted to a fluidic circuit, assigning a reference pressure $P_{0}$ to the center of the cross-slot and a fluidic resistance $Z_{i}$ to the fluidic path from the $i^{th}$ reservoir to the cross-slot. The pressure $P_{i}$ to be applied to the $i^{th}$ fluid reservoir is calculated using the following ohmic relation between the pressure and the flow rate $q_{i}$:

[TABLE]

Identical tubing and fittings are used for each channel, and the fluidic resistance due to the tubing is much larger than the microfluidic channels. We therefore assume the same value of the fluidic resistance $Z_{i}$ for the fluidic path from each reservoir to the cross-slot. Based on the viscosity of the fluid and the dimensions of the tubing and channels, $Z_{i}$ was found to be $1.4\times 10^{15}$ Pa-s/m3. $P_{0}$ is generally set to be in the range of 2.5-3 psi, which is chosen to be in the intermediate range of the maximum of 5 psi for each pressure regulator, thereby maximizing the range of pressure variations in the positive and negative direction for controlling flow rates. Finally, actual pressure values are calculated using the process described above and converted to the appropriate analog voltage values, which is then used to drive fluid flow within the cross-slot device. The entire process repeats with a feedback time loop constant of 33 msec.

III Results And Discussion

III.1 Flow topologies in a six-channel device.

We began by experimentally characterizing permissible flow topologies in a six-channel device. An aqueous suspension of fluorescently labeled microspheres are introduced into the cross-slot channel, and the pressures are varied in the six-channel device to perform particle tracking experiments. In Fig. 3, the three permissible flow topologies are depicted for an $N$=6 channel device based on particle tracking velocimetry (PTV) experiments. The three flow topologies are characterized based on the number of stagnation points in the resulting flow pattern. In the first flow topology (Fig. 3a), five inlet streams and 1 outlet stream give rise to zero stagnation points within the cross-slot. Fig. 3b,c show two situations in which the pattern of inflow and outflow channels are the same. In the second flow topology (Fig. 3b), in clockwise direction beginning with the horizontal channel on the right, channels 1, 3, 4, and 5 are inlet channels while channels 2 and 6 are outlet channels. In this configuration, there is a single stagnation point created in the cross-slot. In the third flow topology (Fig. 3c), there are two stagnation points created with channels 1, 3, 4, and 5 being inlet channels, and channels 2, and 6 being outlet channels. Both Fig. 3b,c have the same inlet and outlet channels, but by changing the quantitative flow rates in each channel, flow topologies with one or two stagnation points are possible.

III.2 Flow topologies during particle manipulation.

The Stokes trap allows for different types of trapping or manipulation experiments. In one type of experiment, two particles can be trapped at pre-determined separate locations inside the cross-slot. In this way, two particles are essentially confined at the location of two distinct stagnation points. In a second type of experiment, two particles can be independently moved along pre-determined trajectories from initial to final positions, where particle trajectories are selected in an optimal manner as characterized by Eq. 5. At the end of the particle translation experiment, the flow field will typically vanish, however, trapping can be re-initiated to the new particle positions if desired.

In the following, we examine flow topology during the multiple particle manipulation experiment, wherein the initial and final positions for particle trajectories are intentionally chosen such that the emergent flow topologies are not a result of the symmetry of the device. In this way, we perform three different particle manipulation protocols: (1) moving two particles closer to each other, (2) moving two particles away from each other, and (3) interchanging the positions of the particles, shown schematically in Fig. 4.

Flow topologies and particle trajectories are characterized using experiments and numerical solution of the flow model. For experiments, the particle position, applied flow rates, applied pressures, and measured pressures for each regulator are recorded at each sampling instance (sampled at a rate of 30 Hz). The control weights are set to $\beta=0.0001$ and $\gamma=1000$. For simulations, the initial and target positions of the two particles and the weight parameters in the control model are set identical to experimental values. In this way, the MPC algorithm Eq. 5 was used to determine the simulation flow rates for the next sampling interval, and Eq. 4 was numerically integrated using a suitable scheme (e.g. Runge-Kutta 45) to estimate particle position in the next sampling instant given these flow rates. This process was repeated until the the total simulation time matched the experimental duration. Brownian motion was not included in the model because particle Peclet numbers were generally larger than unity for all experiments.

We began by studying particle trajectories for the case of two particles approaching each other during active control, which corresponds to the schematic in Fig. 4a. Fig. 5a shows particle trajectories from experiments and numerical simulations. Particle 1 was initially trapped at a coordinate of [-55, 10] and particle 2 was initially trapped at the coordinate [-34, -37], where positions are given in dimensionless units. At time $t$=0, particle 1’s target location was set to [-31, -6] and particle 2’s target location was set to [-16, -17]. The dashed lines indicate the simulation trajectory and the solid lines show the experimental trajectories. As shown in Fig. 5a, good agreement is generally observed between particle trajectories from experiments and simulations. The individual ‘x’ and ‘y’ components corresponding to experimental and simulation trajectories for each particle are plotted in Fig. 5b-e. In general, our results show that both components of the simulation trajectories agree well with the experimental trajectories, with both converging exponentially towards their target positions (Supplementary Material Video 1).

We further compare the experimentally applied flow rates, experimentally measured (actual) flow rates, and simulated flow rates for the manipulation process in Fig. 6. The experimentally applied flow rates correspond to flow rates calculated using Eq. 6, whereas the actual flow rates are determined by using the experimentally measured pressure values $P_{i}$ (determined in real-time for each regulator), calculating $P_{0}$ using the relation $P_{0}=\sum^{6}_{i=1}{P_{i}}/6$, and then back-calculating $q_{i}$ using Eq. 6. For completeness, the applied and measured pressures for each regulator are shown in Fig. 7. In general, we observe that deviations in the experimental and simulated particle trajectories (as observed in the $x$ and $y$ components) are correlated to differences in the applied and measured pressures. The pressure regulators generally show slightly different response times depending on whether a positive or negative pressure differential is applied. In particular, pressurizing tends to occur faster than depressurizing. Accordingly, for $P_{1}$ and $P_{2}$ in Fig. 7, the measured pressures tend to lag behind the applied pressures because initially a large negative pressure differential is required. However, the applied and measured pressures track fairly closely for the other regulators. This asymmetric difference in response rates leads to an initial transient deviation in the applied and actual flow rates, thereby resulting in the corresponding deviations in initial trajectories.

Given the generally good agreement between the experiments and simulations for particle trajectories and flow rates, we next used flow rates determined from simulations to analyze flow topologies. In this way, we plot streamlines during particle manipulation events to visualize the flow topologies (Fig. 8). During the particle manipulation event corresponding to the trajectories in Fig. 5, the particles first move in a roughly unidirectional flow in the same direction (Fig. 8a). Next, a single stagnation point is created between the two particles which moves closer and closer to the centroid of the two particle positions (Fig. 8b-c). Finally, both particles move towards their target positions driven by planar extensional flow, with a single stagnation point present between them (Fig. 8d). In this way, the particles approach each other at a very close distance along the compressional axis (Supplementary Material Video 2).

We next characterized flow topologies for the case of particles moving away from each other, which corresponds to the experiment shown schematically in Fig. 4b. Particle trajectories for this case are shown in Fig. 9, where particle 1 is initially trapped at [-31, -6] and particle 2 is initially trapped at [-16, -17], and at time $t$=0, the particle’s target locations are set to [-55, 10] and [-34, -37], respectively. Again we generally observe good agreement between experimental and simulated trajectories for both particles (Supplementary Material Video 3). The deviations of experimental trajectories from the simulated trajectories are present only during initial transients, after which there appears to be good agreement. As in the prior case of particle motion during close approach, deviations are again caused by the slow depressurization response of the pressure regulators (Supplementary Material Fig. S1 and Fig. S2).

We also characterized the evolution of the streamlines for this process (Fig. 10). Here, we observe that the particles are again initially transported in a unidirectional flow in one direction (Fig. 10a), after which a single stagnation point is created in the vicinity of the two particles, and both particles then move along planar extensional flow streamlines towards their final target (Fig. 10b-d). Interestingly, the trajectories of the particles are not reversed versions of the trajectories observed in Fig. 5a. As seen previously, the choice of the objective function causes the particles first to be moved by a roughly uniform flow, which then transitions to a single stagnation point flow, thus leading to different trajectories when the particles are moved closer together rather than further apart (Supplementary Material Video 4).

Finally, we characterized flow topologies corresponding to an interchange of particle positions, corresponding to the manipulation event schematically in Fig. 4c. Particle trajectories for experiments and simulations are shown in Fig. 11. In this experiment, particle 1 is initially at [-61, 15] and particle 2 is at [-34, 7], and at time $t$=0, the set points for both particles are interchanged. This manipulation event results in a large change in the particles’ locations, so it is accompanied by correspondingly large values of flow rates and applied pressures (Supplementary Material Fig. S3 and Fig. S4). As a result, we initially observe deviations in the experimental and simulation trajectories for each particle (Supplementary Material Video 5). However, after the initial transients decay, the experimental and simulation trajectories generally show good agreement. In this process, a stagnation point is created between the two particles (Fig. 11b), and as the particles approach each other, the principal axes of the stagnation point flow rotate as the particles attempt to move past each other slowly (Supplementary Material Video 6). Comparing the simulation and the experimental trajectories near the instant of close approach, we posit that hydrodynamic interactions (HI) between the two particles may cause the simulation trajectory to approach the set point faster than the experimental trajectory, as particle-induced HI is not included in the simulations. Finally, once the particles are able to move away from each other, the remainder of the motion towards the final target location is fairly rapid, with the stagnation point rotating further to facilitate movement. In this experiment, the particles tend to approach each other closely approximately midway through the trapping process, which may result in localization issues as the tracking algorithm could confuse the identities of the two particles. To avoid this issue, we added an additional constraint to the MPC formulation such that particles cannot approach closer than 4 particle diameters.

In summary, flow topologies corresponding to the three canonical two-particle manipulation experiments consistently involve flow fields with zero stagnation points (for rapid repositioning) and one stagnation point (for adjusting and approaching final positions), rather than two stagnation points located at the center-of-mass of each particle. We conjecture that this set of flow topologies allows for more efficient, non-trivial dynamical manipulation of particles rather than a simple quasi-static displacement of trapping points. In general, the control algorithm tries to impose a unidirectional flow which moves the particles closer to their target locations, followed by switching to a stagnation point flow, where the algorithm is capable of both translating the stagnation point as well as rotating the principal axes of the extensional flow.

III.3 Effect of the controller weights.

We now focus on the impact of the controller parameters on the speed of manipulation, including the effect of the controller weight $\beta$ associated with the flow rates and the controller weight $\gamma$ associated with the end penalty term $\gamma$. The optimal trajectories sensitively depend on the initial and final positions of the particles. In general, we can divide these trajectories into two categories, including one case wherein the distance of each particle to the set point smoothly decays monotonically and a second case wherein the distance of a particle initially decays to zero, then increases, and finally decays to 0 again. We consider the effect of the controller weights on both cases in order to fully characterize the sensitivity to speed of manipulation. In all cases, the default values for the controller weights are $\beta=10^{-4}$ and $\gamma=1000$. We characterize the speed of manipulation by measuring the time taken by each particle to first reach within 1 particle diameter of its target location.

We began by considering the case where both particles’ distance to the target position smoothly decays to 0 nearly monotonically. We first vary the control weight $\beta$ while holding $\gamma$ constant. Particle 1 is trapped at [10,10], and particle 2 is trapped at [-10,-10]. At time $t$=0, the set points for these particles are changed to [-30,10] and [-10,-30], respectively. Five different values of $\beta=1.0\times 10^{-5},3.2\times 10^{-4},1.0\times 10^{-4},3.2\times 10^{-3},1.0\times 10^{-3}$ are selected, which are spaced equally on a logarithmic scale between $10^{-5}$ and $10^{-3}$.

Fig. 13a shows the time required for particles to reach within 1 particle diameter for a given value of $\beta$. In this case, the distances of both particles to the target position decay smoothly to 0 roughly monotonically (Fig. S5a,b). Interestingly, a linear relationship between the value of $\beta$ and the time required for manipulation is observed. $\beta$ is a regularizer for the flow rate applied during manipulation, so small values of $\beta$ allow for larger values of the flow rates for the same value of the objective function. As a result, the time required to reach the target location increases as $\beta$ increases, because the flow rates are restricted from increasing without limit. In addition, even though both particles traverse different distances along their trajectories, they require roughly the same amount of time to reach their target location. This behavior occurs because the deviation of each particle from its target location is equally weighted in the objective function, so the flow rates obtained by minimizing the objective function are such that they cause both particles to arrive at the target in the same time. Moreover, we observe that as $\beta$ decreases, the path begins to approach a straight line between the initial and final position (Fig. S5c).

We next varied the control parameter $\gamma$ corresponding to the end penalty term. We select five different values of $\gamma=1.0\times 10^{2},3.2\times 10^{2},1.0\times 10^{3},3.2\times 10^{3},1.0\times 10^{4}$, which are again equally spaced on a logarithmic scale between $10^{2}$ and $10^{4}$. Fig. 13b shows the effect of $\gamma$ on the time taken by each particle to reach within 1 particle diameter of the target location. We again observe that both particles require nearly the same amount of time to reach their respective target locations. As previously discussed, $\gamma$ plays a critical role in removing offsets at the end of the prediction horizon. When the value of $\gamma$ is increased, the deviation of particle position at the end of the prediction horizon becomes more important compared to the flow rates and particle deviations over the prediction horizon. For large values of $\gamma$, the controller chooses larger flow rates to reduce the offset in particle position in an attempt to minimize the contribution of the end penalty term. In this way, increasing $\gamma$ leads to a reduction in the time required to reach the target location, which is shown in Fig. 13b. In addition, increasing the value of $\gamma$ beyond a certain value appears to have a negligible effect on the time required to reach the target location (Fig. 13b). Particle trajectories corresponding to these simulations are shown in Fig. S6c, while the individual particle distances over time are shown in Fig. S6a,b. We observe that the initial curved trajectory for both particles begins to straighten out with larger values of $\gamma$, which suggests that the particles are approaching the straight line trajectory between their initial and desired final locations. We characterized behavior for multiple initial and end point positions that likewise generated smooth trajectories. In all cases, we observe a similar dependence of $\beta$ and $\gamma$ on the time required for particles to reach final position (Figs. S7-S12).

Next, we consider the case where the distance of either particle to the target location does not smoothly decay to zero. To probe this case, we initially trap two particles at [20,30] and [-50,-40] and then change their set points to [30, 40] and [-30,-20], respectively. Fig. 14a shows the effect of variation of $\beta$ on the speed of manipulation. In this case, the trajectory for particle 2 was non-smooth, as the distance to the set point did not monotonically decay to zero (Fig. S13b), whereas the distance of particle 1 smoothly decays to zero (Fig. S13a). Accordingly, we find that the time required by particle 1 to reach the final position increases with $\beta$, but for particle 2, the time first increases and then saturates (Fig. 14a). Interestingly, both particles require different times to reach within a particle diameter of the target position, in contrast to the smooth trajectory case. Fig. 14b shows the effect of variation of $\gamma$ on the time required to reach 1 particle diameter. Particle 1 has a familiar trend with $\gamma$ as the time decreases with increasing gamma initially, and then has negligible changes for larger values of $\gamma$. On the other hand, particle 2’s time increases initially with $\gamma$, and then saturates. Again, these results suggest that particle 1 has a smoothly decaying distance to the target location, whereas particle 2 has a non-smooth trajectory (Fig. S14a,b).

IV Conclusions

In this work, we systematically analyze flow topologies during multiplexed particle manipulation using a Stokes trap. In MPC, the controller determines optimal flow rates to minimize an objective function. Compared to the initial version of the automated hydrodynamic trap Tanyeri and Schroeder (2013), which only has a simplistic design with a single stagnation point, the flow topologies generated in the Stokes trap are complex and time-dependent. Achieving an understanding of the dynamic flow fields during particle trapping will be essential for future applications in soft matter and flow-induced particle deformation experiments using this technique.

Our results show that a six-channel cross-slot is capable of three distinct flow topologies with zero, one, and two stagnation points, which can be accessed by changing the flow rates in each channel. Using a combination of experiments and simulations, we find that optimal particle manipulation involves a two-stage strategy: in the initial stage, there are no stagnation points in the vicinity of the two particles, however, in a later stage, a single stagnation point is generated between the two particles. The action of the controller is then to translate this single stagnation point and to rotate the orientation of the principal axes of extension and compression to enable a precise, exponential approach to the desired particle positions. In this way, a simpler flow field topology with zero or one stagnation point is superior in performing particle manipulation as opposed to a more complicated flow pattern, where two stagnation points could potentially control particle position.

The apparent advantage of zero or one stagnation points for manipulating two particles manifests itself in both shorter duration of the relocation process and in smaller flow rates needed to perform the task, both of which are crucial limitations in practical applications. Our findings have significant consequences for particle manipulation in general, but in particular for applications that involve particle-particle interactions. For instance, if two vesicles are trapped in this setup, and they are brought close to each other in a controlled fashion for adhesion experiments, our work establishes that both the particles essentially experience planar extensional flow, instead of the complicated two stagnation point flow (as shown in Fig. 3b).

We further explore the impact of the controller parameters on the speed of manipulation using simulations. In general, our results show that the time required to reach the target location is the same for both the trapped particles when they follow smooth trajectories. In addition, we found that increasing $\gamma$ leads to faster manipulation initially for smooth trajectories, after which the performance saturates. Broadly, this work provides a systematic characterization of the performance of the Stokes trap by examining the evolution of flow topologies during multiplexed particle manipulation, and further shows how the controller parameters can be tuned to obtain the desired level of trapping performance. Taken together, these results will help guide future experimental studies on deformable soft materials in well defined flow fields using a Stokes trap.

Author contributions

A.S., S.H. and C.M.S conceived the project and designed the experiments. A.S. performed the experiments, while A.S. and D.K. analyzed the experimental data. S.H. and C.M.S supervised the research and all authors wrote the manuscript.

Acknowledgements.

The authors thank Prof. Christopher V. Rao for insightful discussions. This work was supported by the National Science Foundation by Award # NSF CBET 1704668.