Sorting cells by their dynamical properties

Recent advances in cell sorting aim at the development of novel methods that are sensitive to various mechanical properties of cells. Microfluidic technologies have a great potential for cell sorting; however, the design of many micro-devices is based on theories developed for rigid spherical particles with size as a separation parameter. Clearly, most bioparticles are non-spherical and deformable and therefore exhibit a much more intricate behavior in fluid flow than rigid spheres. Here, we demonstrate the use of cells’ mechanical and dynamical properties as biomarkers for separation by employing a combination of mesoscale hydrodynamic simulations and microfluidic experiments. The dynamic behavior of red blood cells (RBCs) within deterministic lateral displacement (DLD) devices is investigated for different device geometries and viscosity contrasts between the intra-cellular fluid and suspending medium. We find that the viscosity contrast and associated cell dynamics clearly determine the RBC trajectory through a DLD device. Simulation results compare well to experiments and provide new insights into the physical mechanisms which govern the sorting of non-spherical and deformable cells in DLD devices. Finally, we discuss the implications of cell dynamics for sorting schemes based on properties other than cell size, such as mechanics and morphology.

1 Experimental materials and methods 1

.1 Device fabrication and design
Devices were fabricated by replica molding, using the same method and equipment used in previous work by Beech et al. [1]. Adhesion of blood cells to the inner surface of our devices was reduced by formation of a polymer brush on the polydimethylsiloxane (PDMS). Immediately after O 2 plasma treatment and bonding, devices were filled with 0.2% PLL(20)-g[3.5]-PEG(2) (SuSoS AG, Dübendorf, Switzerland) and rinsed after 20 minutes with autoMACS ™ .
The experimental device has been designed to sort particles of different sizes, with successive sections within the device corresponding to different critical radii [2]. The device geometry [1,2] is defined by the post diameter D = 20 µm, lateral center-to-center spacing between posts λ = 32 µm, lateral shift in successive pillar rows ∆λ, and the height between the two enclosing plates H, as shown in Supplementary Figure S1. The lateral shift can conveniently be defined as a fraction of the post spacing ∆λ = λ. The experimental device consists of 13 consecutive sections of obstacle arrays, which are differentiated by the different lateral shifts ∆λ between successive rows as summarized in Supplementary Table S1. We note that some of the sections have row-shift fractions that are equivalent to M/N , where M and N are integers and M = 1, i.e. we have a "non-integer" row-shift fraction. This is due to limitations in the photomask fabricating process (a manufacturing grid of 200nm), which means that in many cases the desired integer row-shift fractions are not possible, in which case the nearest non-integer row-shift fraction is used. While non-integer row-shift fractions have been shown to give rise to two additional zig-zag modes with accompanying critical sizes [3,4], in the present study these modes are not observed. Two sequential devices have been used, distinguished by the distance between the top and bottom plates covering the obstacle arrays. One device had a height of H = 11 µm, which is larger than a RBC diameter of about 8 µm, while the other device had H = 4 µm, which is smaller than the RBC diameter but larger than the RBC thickness of about 2 − 3 µm. We will refer to these two devices as 'thick' and 'thin' DLDs, respectively.
Several inlets and outlets were employed to control flow and sample input/output, see Supplementary Figure S1. Both devices had one large fluid buffer inlet and one small buffer inlet to minimize the effects of walls and to create a well defined starting position to which total displacement can be compared. The sample inlet channel has been placed between the two buffer inlets and included a filter to remove any large particles which could cause clogging. A pressure gradient was used to drive the flow. Outlets were kept at atmospheric pressure and the overpressure at the inlets was maintained using a MFCS-4C pressure controller (Fluigent, Paris, France). Flow in the devices was driven by a pressure drop of 22 mbar between the buffer inlets and outlets.

Sample preparation
Blood was extracted from healthy volunteers via finger pricking. As mentioned in the paper, measurements were conducted both at a viscosity contrast of 5 and at a viscosity contrast of 1 between the cytosol and the surrounding  Table S1: Parameters defining the obstacle array geometry for each section of the DLD device. ∆λ is the lateral shift between successive rows for each section, = ∆λ/λ is the section's shift fraction, and Rc is the critical particle size [5,6]. The fourth column in the table specifies the number of rows in each section.
buffer. For high viscosity contrast measurements, the cells were suspended in autoMACS ™ running buffer (Miltenyi Biotech, Auburn, CA). This solution (pH 7.2) contains phosphate buffered saline (PBS) supplemented with 0.5% bovine serum albumin (BSA), 2 mM EDTA, and 0.09% azide. It acts to suppress blood clotting and is isotonic so it does not affect the cell shape due to a changed osmotic pressure. For low viscosity contrast measurements the autoMACS ™ running buffer was supplemented with Dextran-500 (# 700013-096, VWR International LLC, PA, USA).
To measure the viscosity (see Supplementary Figure S2), Ubbelohde viscometers were used (UBBEL Visco,Paragon Scientific Ltd, UK.). The temperature ranged between 21.8-22.°C, which is similar to the measurements conducted with cells in the DLD device. Dextran is widely used to change the viscosity of various samples. It is a neutral polysaccharide, i.e. the pH or salt concentration of the solution does not affect the resulting viscosity. In order to minimize the change to the osmotic pressure across the cell membrane, this larger Dextran molecule with an average molecular weight of 500 kDa was chosen. At a concentration of 11% the change in osmotic pressure is 0.22mM which is negligible compared to isotonic blood osmolarity of 300mM. Furthermore, microscopic examinations also confirmed that the low viscosity contrast samples did not visibly differ compared to the high viscosity contrast samples.

Data acquisition and analysis
All images were taken through an inverted Nikon Eclipse TE2000-U microscope (Nikon Corporation, Tokyo, Japan). High-speed images were taken using an EoSens mini MC-1370 camera (Mikrotron GmbH, Unterschleissheim, Germany). In all other cases an Andor iXon EMCCD camera (Andor Technology, Belfast, Northern Ireland) or Hamamatsu Orca Flash 4.0 (Hamamatsu, Shizuoka Pref., Japan) were used. Throughout the course of experiments, the dynamics and displacement of single RBCs in different device sections have been monitored.
The simultaneous transit of many RBCs through each section of the DLD device has been recorded to video using the microscope and camera set-ups previously described. Trajectories of individual RBCs were extracted from these recordings using the particle-tracking application MOSAIC in the image processing suite ImageJ [7]. Particle detection has been improved by subtracting a median average of all frames from the entire video and removing the background.
After extracting the sets of trajectories for different device sections from the video recordings, a pre-screening procedure has been performed. In the pre-screening process, trajectories which are too short or where direct interactions between RBCs have been detected, were discarded. As a rule of thumb, trajectories were considered to be too short if they were much shorter than the video domain, corresponding to RBC traversal over only a few pillars in the device. The short trajectories might be present at the beginning or the end of videos and also due to a detection failure if at some point video contrast becomes insufficient for tracking. The trajectories with direct interactions between RBCs Supplementary Figure S2: Viscosity measurement of different dextran-500 concentrations in autoMACS. In order to achieve a viscosity contrasts of C = 2, 1 and 0.25 between the medium and the cytosol, the suspending medium viscosity needed to be increased to two, five and twenty five times that of normal autoMACS buffer viscosity.
correspond to those where collisions between RBCs were detected. These trajectories were also excluded from further analysis, since such situation has not been considered in simulations.
In order to establish a common pattern followed by the RBCs in a given section of the device, all pre-screened trajectories were superimposed on top of one another and averaged into a single trajectory. Optimal positioning of trajectories for superposition has been found by the minimization of the Hausdorff distance between each member in the set of trajectories. The Hausdorff distance d H (X, Y ) between two sets (X, Y ) of points in 2D Cartesian coordinates is defined as where sup denotes the supremum and inf the infimum. d H (X, Y ) can be thought of conceptually as the largest of all distances from a point in one set to the closest point in a second set. This method for the alignment of trajectories is very similar to techniques used in image-recognition algorithms [8]. In our algorithm, we have fixed a base trajectory in place, while a second trajectory traverses an overlaid grid to find the position with a minimum Hausdorff distance to the former trajectory. Additional refinement of the RBC trajectories has been done during this alignment stage, where some outliers were removed if computed d H of a trajectory was larger than one standard deviation from the mean d H value of all trajectories in the given section. These outliers can arise for multiple reasons, such as the presence of slightly different transit modes within one section (see Supplementary Figure S4), possible blockage between two posts (see Supplementary Figure S5), or even due to missing posts which might be due to some problems in device fabrication.
A summary of the RBC trajectory refinement for each device is shown in Supplementary Tables S2, S3, S4, and S5, where the numbers of total, accepted, and rejected trajectories are given for all experiments. Supplementary Figure S3 presents visually the fraction of accepted trajectories for the thick and thin devices at C = 5. Notice that the sections with lowest trajectory acceptance are those where the transition between displacement and zigzag modes is occurring, or where extremely negative zigzag modes with slight variations in row-swapping frequency result in large differences in lateral displacement per pillar. The possibility of different transit modes within a single section is illustrated in Supplementary Figure S4 for the section 11 of the thick device. 69 aligned trajectories yield a main average trajectory in this section with a displacement-zigzag pattern {5,4,5,4,4,5,4}, while 17 rejected trajectories seem to correspond to two other displacement-zigzag patterns {5,4,4,4,5} or {5,4,5,5,4}. The differences can arise from natural variations in RBC properties indicating that this section might be sensitive to them. However, it is also possible that the recorded trajectories are too short to properly capture the period of transit mode in section 11 and the displacement-zigzag patterns above correspond to just portions of a full transit mode in this section. Currently, we cannot rule out one or another possibility. In addition, Supplementary Figure S5 illustrates a rejected trajectory, where a blockage between pillars has been detected. This example corresponds to a displacement mode and the sudden jump in the trajectory is due to the blockage of a single inter-pillar space identified visually from the video.

Fluid simulation -mesoscale hydrodynamics
To represent the suspending fluid, we employ a mesoscale hydrodynamic simulation approach which is a variation of the smoothed dissipative particle dynamics (SDPD) method [9], adapted to conserve angular momentum [10]. SDPD is a particle-based fluid-dynamics method [9] well suited to mesoscopic length scales. It improves upon the popular dissipative particle dynamic (DPD) method [11,12], by incorporating the force scheme used in smoothed particle hydrodynamics [13,14], which has been derived directly through the discretization of the Navier-Stokes (NS) equations. A caveat of the original SDPD method [9] is its violation of conservation of angular momentum caused by force components which act perpendicular to the inter-particle axis. Recent multi-particle collision dynamics simulations [15] and SDPD simulations [10] have demonstrated that conservation of angular momentum is necessary in order to properly describe the dynamics of two different fluid phases (e.g., extra-and intra-cellular fluids separated by a RBC membrane) in flow. Consequently, the original SDPD formulation has been extended to obtain a new SDPD+a method, which satisfies angular momentum conservation [10]. This method is necessary for accurate simulation of RBC dynamics in the DLD device.  which include conservative (C), dissipative (D), rotational (R), and stochastic (S) terms as follows where p i and p j are local particle pressures given by the equation of state p = p 0 (ρ/ρ 0 ) α − b, with p 0 , ρ 0 , α, and b being selected model parameters [10,14]. Particle density ρ is calculated locally as ρ i = j W L (r ij ), where W L (r) = 105 is the Lucy function [13] and r c is the cut-off radius. The weight function w(r) is calculated by ∇W L (r) = −rw(r), such that w ij = w(r ij ). The coefficients γ ij and σ ij determine the strength of the dissipative and random forces, where the friction coefficient γ ij is defined as with η 0 being the desired dynamic viscosity. The random force coefficient is σ ij = 2 k B T γ ij , while tr[dW ij ] and dW s ij refer to the trace of a random matrix of independent Wiener increments dW ij and its traceless symmetric part, respectively. For a formal derivation of these equations, the reader is referred to Ref. [10].
Newton's second law of motion governs the time evolution of the i-th particle's position, translational velocity, and rotational velocity as followsṙ where N ij = 1 2 r ij × F ij is the torque exerted on particle i by particle j. The equations of motion are integrated using the velocity-Verlet algorithm [16].
Section Total # Rejected Accepted Accepted Average # of traj. fraction of posts

Red blood cell model -triangulated surfaces
The RBC membrane is modeled as a triangulated network of springs [17][18][19][20][21], whose vertices are coupled to the fluid via frictional forces. A total of N v particles constitute the mesh vertices and N s = 3(N v − 2) springs follow the edges of the mesh, reproducing the elasticity of the membrane (U sp ). A total of N t = 2N v − 4 triangles make up the entire membrane surface and incident triangles have an associated potential energy (U bend ) given by the angle between them, associated with membrane bending rigidity. Furthermore, local and global area constraints (U area ) are enforced along with a global volume constraint (U vol ). Formally, the total energy of a RBC is given as The total contribution of springs is given by where l j is the length of spring j, l m is the maximum permitted spring extension, x j = l j /l m is the fractional extension towards maximum length, ζ is the persistence length, and k p is the spring constant. The equilibrium length of each spring l 0j is set in accordance with an initial triangulated mesh of a stress-free biconcave RBC shape [20,21]. The membrane bending rigidity in absence of spontaneous curvature is described by a bending energy where k b is the bending constant and θ j is the instantaneous angle between the two triangles incident on edge j. Finally, the area and volume constraints are accounted for by two potentials: where k a , k d , and k v are the global area, local area, and volume constraint coefficients, respectively. A is the instantaneous surface area of the membrane, A j is the instantaneous area of the j-th triangle in the network, and V is the   instantaneous RBC volume. The desired total surface area A r , individual triangle area A 0 j , and interior volume V r are set in accordance with the initial triangulation [20,21].
We relate the RBC model's variables to physical macroscopic properties of the RBC membrane by linear analysis for a regular hexagonal network [20,21]. The membrane shear modulus is related to the spring variables by with x 0 = l 0 /l m . The area-compression K and Young's Y moduli are found as K = 2µ 0 + k a + k d and Y = 4Kµ0 K+µ0 . The Helfrich model is employed to describe the bending coefficient k b in terms of macroscopic bending rigidity κ [22], The value of x 0 is set to 2.2 for all springs [21].

Simulation set-up and boundary conditions
The floor, ceiling, and pillar walls of the DLD device are modeled by a layer of frozen particles with a thickness of r c which share the same equilibrium structure as the suspending fluid. This ensures that all particle interactions occurring in the locality of boundaries do not display artifacts (e.g., particle density variations) which could arise in the absence of wall particles due to an improper distribution of conservative forces. Furthermore, the wall particles are included in calculation of fluid particle densities ρ near the wall. In order to prevent particles from penetrating the walls, RBC vertices and fluid particles are subject to bounceback reflections at walls. Bounce-back reflections are preferred over specular reflections as they achieve a better approximation of no-slip boundary conditions (BCs) at the wall. Finally, to fully guarantee no-slip BCs at the walls, an adaptive tangential force is applied to fluid particles within a distance r c from the walls [23].
To prevent mixing between intra-and extra-cellular fluids, bounce-back reflections for the solvent particles are also introduced at a RBC membrane. Furthermore, the RBC is coupled to fluid flow via viscous friction between the N v mesh vertices and local fluid particles. The dissipative and random force components of the DPD method are used to achieve these interactions [20]. No-slip BCs at membrane vertices are enforced by careful selection of the friction parameter γ in the dissipative force F D . A fluid sheared over the effective surface of a membrane vertex exerts a friction force on the membrane given as F v = V h ng(r)F D dV , where n is the fluid number density, g(r) is the radial distribution function of fluid particles about the membrane vertices, and V h is the hemisphere volume of fluid situated above the vertex. Equating this integral to the total force required by a continuum hydrodynamical description leads to an expression for the calculation of γ [20].
The bumper array is simulated using a 3D domain enclosing a single-column obstacle with its axis in the z direction, perpendicular to the roof and floor of the domain box, as shown in Supplementary Figure. S6. Efficient representation of an infinite bumper array environment is achieved with periodic BCs in the x and y directions and a shift in the y direction for each boundary-crossing event in the x direction. The assumption of an infinite obstacle array neglects any effects induced by the walls at either edge of the device. The missing wall effects are revealed by solving the NS equation for a bumper array in 2D with explicit wall-edge boundaries. We find that the edge walls enforce a zero net flow in the y direction. This condition is recreated in the 3D infinite array using an adaptive force applied to solvent particles such that the net flow in the y direction vanishes.
Simulations are performed for every section of the device, altering the lateral shift ∆λ induced at the periodic boundary to correspond with a specific section of the experimental device, as shown in Supplementary Table S1. This gives us a total of 13 trajectories for a single RBC transit through the 13 sections of each thin (H = 4 µm) and thick (H = 11 µm) device.
Establishing standard length and energy scales allows us to relate most simulation parameters to physical RBC properties. We use the RBC membrane area to define an effective RBC diameter D r = A r /π and an average bond length for a given number of vertices 3Nt . The experimental value for an average healthy RBC's surface area A r = 133.5 × 10 −12 m 2 [24] suggests a D r = 6.5 µm. For time scaling, we define a characteristic RBC relaxation time τ = η o D 3 r /κ r , where η o is the suspending fluid viscosity and κ r is the RBC membrane bending rigidity. The parameters used for the RBC model are summarized within Supplementary Table S6 in units of D r and the thermal energy k B T ; along side we present the corresponding average values for a healthy RBC in physical units.
The SDPD+a fluid parameters are given in Supplementary Table S7. The average distance between RBC vertices l 0 is chosen as the length scale, which has a value of l 0 = 0.4 in simulations. The exponent α in the equation of state takes the value α = 7 and the SDPD equilibrium density ρ 0 = n, where n = 9 is the fluid's number density. The speed of sound for the equation of state is given as c 2 = p 0 α/ρ 0 and the Mach number has been kept below 0.1 in all simulations, ensuring a good approximation for incompressible fluid flow.

Supplementary videos
In order to best demonstrate the differences in RBC dynamics for cells with viscosity contrasts C = 5 and C = 1, we present several videos of simulated and experimental RBC transit in sections 4 and 11 of the thick device.

Viscosity contrast C = 5
The three videos "SIM Sec4Cont5.avi", "EXP Sec4Cont5.avi", and "SIM Sec11Cont5.avi" show the behavior of simulated and experimental RBCs traveling through sections 4 and 11 of the thick device in a suspending medium which is five times less viscous than the intracellular fluid. Notice that row-swapping events occur; the RBCs are in a zig-zag mode. The RBCs display tumbling dynamics, which is especially prevalent before and after row-swapping events.

Viscosity contrast C = 1
The other two videos "SIM Sec4Cont1.avi" and "EXP Sec4Cont1.avi" also show the simulated and experimental RBCs traveling through section 4 of the thick device. However, the inner and outer fluids have the same viscosity and as a result, we see RBCs traveling in displacement modes without zig-zagging. The RBCs tumble to a lesser degree and there are clear local deformations as the cell membrane rearranges itself in a tank treading motion.  Table S7: Parameters defining the model values used for the RBC properties and their physical equivalents. Nv is the number of membrane vertices, Ar is the RBC membrane area, l 0 is the average bond length, Vr is the RBC volume, T is the temperature, Yr is the membrane Young's modulus, κr is the membrane bending rigidity, and k d , ka, and kv are the local area, global area, and volume constraint coefficients, respectively. In all simulations, we have chosen Ar = 133.5 and k B T = 0.4, which implies that Dr = 6.5 and l 0 = 0.4.