A Criterion for the Complete Deposition of Magnetic Beads on the Walls of Microchannels

This paper analyzes numerical simulations of the trajectories of magnetic beads in a microchannel, with a nearby permanent cubical magnet, under different flow and magnetic conditions. Analytically derived local fluid velocities and local magnetic forces have been used to track the particles. A centered position and a lateral position of the magnet above the microchannel are considered. The computed fractions of deposited particles on the walls are compared successfully with a new theoretically derived criterion that imposes a relation between the sizes of the magnet and the microchannel and the particle Stokes and Alfvén numbers to obtain the complete deposition of the flowing particles on the wall. In the cases in which all the particles, initially distributed uniformly across the section of the microchannel, are deposited on the walls, the simulations predict the accumulation of the major part of particles on the wall closest to the magnet and near the first half of the streamwise length of the magnet.


Introduction
The application of magnetism in fluidic microsystems has been used since the boom of microfluidics in the early 2000 [1] with the concept of micro total analysis systems. Magnetic forces can be used to manipulate magnetic objects as magnetic particles or magnetically labelled cells. The movement of magnetic beads in microfluidic applications is usually controlled by permanent magnets or by magnetic fields generated by electric currents. Magnetohydrodynamic pumps, based on the generation of a magnetic field perpendicular to an electric field, can be used, as an alternative to pressure driven or electroosmotic flows, to produce flow in a conducting fluid [2,3]. The same principle can be used for mixing [4,5]. Plugs of ferrofluid can be moved with magnets to induce or to block the flow of an immiscible liquid [6,7]. More information about these few examples and other applications involving magnetic forces in microflows can be found, for example, in the reviews of Gijs et al. [8], Pamme [9], Berthier and Silberzan [10] and, recently, in Yang et al [11].
Microsized magnetic beads are constituted by nanoparticles of iron oxides embedded in a biologically-compatible polymer (latex or polystyrene) sphere. The external surface of the sphere can be coated with biological molecules such as enzymes or DNA fragments which can be easily transported, by applying a magnetic field, towards specific locations or biological targets. Magnetic beads are used mostly for in vitro applications (biodiagnostics and biorecognition) and recently for in vivo applications, such as cancer treatment. In this case functionalized magnetic particles can by transported by the blood flow and retained by a magnet implemented in the treatment zone. For the most common applications the sizes of the magnetic particles range from 5 nm to 6 μm. The smallest particles (5 nm-100 nm) can be used for applications in which the particles interact with proteins, viruses and genes [12]. Magnetic beads used for biotechnology are usually superparamagnetic because it is desired that the magnetic force vanish when the externally imposed magnetic field is switched off.
This study is focused on the prediction of the trajectories of magnetic beads in microchannels and the identification of the required magnetic and relevant physical properties needed for the complete deposition of the magnetic beads on the walls of the channel.

Physical and Mathematical Model
The physical model is shown in Fig 1. A permanent cubical magnet is located above a straight microchannel with rectangular cross section. The distance between the top wall of the microchannel and the bottom of the magnet is L zm . Two different positions of the magnet have been considered. A lateral position L xm = 0 (see Fig 1A) and a centered position for which the magnet is located symmetrically at the center of the microchannel (see Fig 1B). Note that the lateral position generates a weaker magnetic field within the channel than the centered position but, in experiments, it allows the visual inspection of the region close to the magnet where magnetic beads are expected to deposit. In both cases the gravity vector acts along the negative z direction.
The flow is assumed to be incompressible, laminar, steady and the physical properties of the fluid constant. The particles are assumed to be spherical and they do not affect the flow   according to the dilute approximation (i.e. one-way coupled Eulerian-Lagrangian approach). Under these hypotheses the flow can be described by the closed solution of the axial momentum equation for the pressure-driven fully developed flow. The axial velocity distribution, scaled with the averaged velocity, within the cross section of the microchannel is given by Eq 1 (see for example Shah and London [13]) The particles are tracked with a Lagrangian approach based on the numerical solution of the kinematic equation and the force balance equation for each particle. The non-dimensional forms of these equations can be written as, The length and velocity scales used to obtain the non-dimensional variables are the height of the microchannel (h) and the average velocity of the flow ( V ). The three terms on the right hand side of Eq 3 are the gravity force, the drag force and the magnetic force. For micron-sized particles suspended in liquid laminar flows the Stokes number (St) and the particle Reynolds number (Re p ) are much lower than unity and the drag coefficient can be computed as C D = 24/ Re p . Simulations carried out with the addition of the lift force to Eq 3 at the largest Reynolds number considered (Re = 10) indicate that this force [14] can be neglected because trajectories are not affected significantly by the inclusion of this term. As a first approximation, the nearwall drag modification [15] or the lubrication force [16] has been neglected.
The form of the magnetic term in Eq 3 is obtained from the dimensional form of the force per unit mass given in Eq 4 [17], in which the initial magnetization of the particles is neglected [18] Eq 4 assumes a linear dependence of magnetization of the particle with the applied magnetic field (see for example [18]). This approximation is valid for magnetic fields small enough to avoid the magnetic saturation of the particles. The magnetization of superparamagnetic particles, or their magnetization curve, can be modeled with the Langevin's law [10] that relates the degree of saturation as a function of the applied magnetic field. Typical magnetization curves show saturation of the particles between 3 A m 2 /kg [18] and 68 A m 2 /kg [19] with large variability depending on the size and specific composition of the particle.
The magnetic flux density has been computed using the analytical solution for permanent rectangular magnets reported by Furlani [20]. The components of the magnetic flux density are indicated in Eqs 5 to 7.
where F, H and G are defined as, The system of reference, indicated with the prime symbol, used for Eqs 5 to 10 has its origin in the center of mass of the magnet (see Fig 1) and ðx 0 1 ; y 0 1 ; z 0 1 Þ and ðx 0 2 ; y 0 2 ; z 0 2 Þ are the coordinates of two diagonally opposed vertices of the magnet which has the magnetization aligned with the z' direction. For example, for a cubical magnet of size L m , x 0 . Eqs 2 and 3 were numerically integrated using the second order Crank-Nicolson scheme. The Lagrangian tracking in-house code used here, developed in FORTRAN language and parallelized using OpenMP directives, has been previously applied for the simulation of the turbulent dispersion of particles in forced [21] and natural convection [22] flows.
The fluid velocity at the position of the particle was computed using Eq 1. In Eq 3, the term ðB Ã Ár Ã ÞB Ã , at the particle position, has been computed analytically. The exact expressions for the three components were derived using Mathematica [23] and the output obtained with the FortranForm command was directly copied and pasted in the simulation code to avoid errors. A typical time step for the simulations was Δt Ã = 10 −3 .
The computational domain for the lateral position of the magnet had dimensions while half of the microchannel was considered for the centered position because of symmetry with respect to the plane x Ã = 2.5 (see Fig 1B). In this case the dimensions were L Ã The computational domain was divided into 49x499x19 equal volumes for the lateral case and into 99x499x19 equal volumes for the centered case. A particle was placed at the center of each volume located at the inlet of the channel (y Ã = 0) (i.e. 49x19 = 931 particles for the lateral case and 99x19 = 1881 particles for the centered case). The time marching scheme was initialized setting the velocity of the particle equal to that of the fluid at the specific location of the particle. The positions of the particles were stored and the instantaneous concentration of particles was determined by computing the number of particles in each volume.
This information was used to calculate the joint conditional probability for a particle to be at location x, y, z at time t, given that the particle was released at the position x 0 , y 0 = 0, z 0 at time t 0 = 0. This probability can be used to extract information about the behavior of a continuous sources of particles located at the inlet of the channel. A similar approach is used in the simulation of scalar dispersion in turbulent flows at high Schmidt numbers [24]. This procedure based, on the Lagrangian tracking, overcomes the use of very fine grids, needed by the Eulerian approach, to capture the thin mixing interface of scalars (or clouds of particles) with very low molecular (or Brownian) diffusivity. For example, for 1 micron particles the Brownian diffusivity [10] at ambient temperature is 4Á10 −13 m 2 /s and the corresponding Schmidt number is 2Á10 6 . Even for the 5 nm particles the Schmidt number is about 10 4 , which makes the numerical solution of the transport equation for the concentration of particles computationally very expensive using the conventional Eulerian approach.

Results and Discussion
As suggested by Eq (3), the trajectories of the particles are dominated by the drag force, that is proportional to ðd Ã p Þ À2 and by the magnetic force, proportional to the non-dimensional group M, defined in Eq 11 and to the non-dimensional term ðB Ã Ár Ã ÞB Ã .
The first contribution to the magnetic force is the group M that depends on the magnetic characteristic of the particles (χ) and on the magnetization of the magnet. The second contribution, the term ðB Ã Ár Ã ÞB Ã , depends only on the size of the magnet and its relative position with respect to the microchannel.
As an example of the spatial distribution of the force generated by a magnet, Fig 2 shows the contours of the three components of the term ðB Ã Ár Ã ÞB Ã for the lateral (Fig 2A) and centered (Fig 2B) position of the magnet. The field corresponds to a cubical magnet (L m = 5 mm) located at L zm = 1 mm, above a microchannel (h = 0.1 mm, W = 0.5 mm). This distance is the typical width of the plastic cover of the microchips. The plots correspond to top views (x-y plane) of the microchannel. It can be seen that the lateral configuration (Fig 2A) generates significant vertical (z) and lateral (x) forces towards the magnet. The centered case (Fig 2B) generates intense vertical forces with maxima located close to the edges of the magnet and maximum axial (y) forces near the edges of the magnet.
To estimate the size and characteristics of the magnet needed to capture the magnetic beads which are flowing within a fluid we assume that the magnet is located at the top of the microchannel with a size larger or equal to the width of the microchannel.
The vertical velocity induced by the magnet can be estimated from the steady state vertical force balance on the particle. (see Eq 3) Rearranging Eq 12, the dimensional vertical velocity can be written as, The particles near the bottom wall of the channel need to travel vertically a distance h to reach the top wall of the channel, which is located below the magnet. If the vertical velocity is given by Eq 13, then the time for the vertical travel is To verify that all the particles are captured by the magnet we can consider that the particles located at the bottom wall of the channel should not travel along the streamwise direction a distance larger than the length of the magnet during the period of time given by Eq 14. If we take the average flow velocity, as an upper limit, an estimation of the length of the magnet (L m ¼ t z V ) can be written as, The magnitude of the term ðB Ã Ár Ã ÞB Ã z needs to be evaluated for the particular shape, size and location of the magnet. Fig 3 shows the magnitude of this term for different sizes of cubical magnets (L m ) located at the top of the microchannel of width W. The height of the microchannel is h = 100 μm and the bottom of the magnet is 1mm above the top wall of the microchannel. The value has been averaged in the volume of the microchannel below the magnet. It can be seen that the average value of this term for the lateral configuration depends strongly on the aspect ratio of the microchannel while for the centered case this dependence is smaller. In the case of the lateral position of the magnet and for comparison purposes, we include in Fig 3A  the  For cubical magnets with sizes between 4 mm and 10 mm the term jðB ÁrÞB z j is in the range 5000 m -1 -6000 m -1 . Additionally if we consider that the typical values of the magnetic saturation M s for permanent magnets are between 5Á10 5 -10 6 A/m (B s = 0.6-1.5 T) and that 0.1 χ 1 [18], ρ p % 1600 kg/m 3 and h%100 μm the contribution of the gravity force can be neglected in comparison with the magnetic force and Eq 15 can be written as Eq 16 indicates that the non-dimensional size of the magnet needed to capture the particles depends on the particle Stokes number and on the term MðB Ã Ár Ã ÞB Ã z that represents the ratio between the applied local magnetic force and the particle inertia (r p V 2 ). This ratio is known in magnetohydrodynamics as the Alfvén number, Al ¼ B 2 =m o r V 2 (see for example Lee and Choi [25]). Similarly we can introduce and define the particle Alfvén number as,  Fig 1A) (B) the centered position (see Fig   1B). using the volume-averaged magnetic force in the microchannel below the magnet (i.e.
MhðB Ã Ár Ã ÞB Ã z i. This quantity can be computed with Eqs 5 to 10 or obtained from Fig 3 for the specific conditions considered in this figure. Eq 16 can be rewritten as, Under physical conditions expressed above and at Re ¼ V h=n ¼ 1; ð V % 1 mm=sÞ, the lengths predicted by Eq 16 for M s = 5Á10 5 A/m and 10 6 A/m are 1.5 mm and 0.4 mm, respectively for χ = 1, and 15 mm and 3.8 mm, for χ = 0.1. This variability of the sizes, depending on the magnetic properties of the magnet (M s ) and of the particles (χ), suggests that the selection of the size of the magnet should be carried out with the knowledge of these properties.
In these particular examples the maximum degree of saturation of the particles depends on the maximum value of the magnetic field within the channel, which is located at the top wall. Specifically for M s = 5Á10 5 A/m (L m = 1.5 mm) the maximum values of the magnetic field for the centered position of the magnet and for the lateral position are 0.06 T and 0.04 T, respectively. The magnetization of the particle can be considered linear in the range of magnetic field from 0 to these maximum values depending on the particular magnetization curve of the particles, which some manufacturers of particles facilitate. For the present conditions, the linear dependence of the magnetization on the magnetic field for 0<B < 0.06 T is a reasonable approximation (as shown in Fig 7 of [19], Fig 5 of [26], Fig 4 of [27]) that correspond to particles of different materials and sizes.
The criterion expressed in Eq 18 is valid for magnetic particles with negligible Brownian motion. It is known that trajectories of nanoparticles can be affected by thermal random fluctuations [28][29][30] . Fig 4 illustrates a sketch of the possible paths of a particle initially located at the bottom of the channel (point A). The effect of the Brownian motion, which is indicated by the standard deviation, σ, of the displacements around a mean path (see Fig 4), can be handled by an extension of the size of the magnet of e m to capture the fraction of particles that would reach The standard deviation of the random displacements can be associated with the Brownian diffusion coefficient (D Br = k B T/3πμd p ) and the time for the vertical travel of the particle as, The substitution of Eq 14 into Eq 20 and using Eq 19 leads to, where the effect of gravity has been neglected (i.e. Fr −2 = 0). Eq 21 relates the extension of the magnet required, given the characteristics of the particles and the dimensions of the channel, to capture the particles that have been drifted along the streamwise direction by the action of the Brownian motion. If the dimension of the magnet is much larger than the extension (L m ) e m ) and the channel height (L m ) h), Eq 21 can be rewritten as, (1) The volume averaged non-dimensional magnetic force, hðB Ã Ár Ã ÞB Ã z i, is 0.37 for the lateral cases and 0.59 for the centered cases (L m = 5 mm and W = 0.5 mm) according to Fig 3. doi:10.1371/journal.pone.0151053.t002 As an example, considering hðB Ã Ár Ã ÞB Ã z i ¼ 0:5, T = 300 K, χμ o = 1.26Á10 −6 N A −2 and M s = 10 6 A m −1 , extensions of 2% and 20% are required for particles with diameters of 200 and 50 nanometers, respectively. In any case, if there are limitations for the increase of the size of the magnet, the selection of a magnet with a larger magnetization can be an alternative option to satisfy Eq 21.
Numerical simulations of the trajectories of the particles have been carried out for the lateral position of the magnet (cases L1 to L4) and for central position (cases C1 to C4). The physical parameters and the conditions considered for the simulations are summarized in Tables 1 and  2. Table 2 shows that in cases L1, L2, C1 and C2, L m /h is larger than 1/(St Al p ) and, in agreement with Eq 17, all the particles initially released at the inlet of the microchannel are deposited on the walls.  In Cases L3, L4, C3 and C4 only a fraction of the particles is captured by the magnet as indicated in Table 2. Fig 6 marks the initial positions of the particles that are not deposited on the walls and consequently they leave the computational domain through the outlet of the microchannel. It can be seen for the lateral cases L3 and L4 (Fig 6A) that the magnet deposits the particles initially located close to the top and to the right lateral walls. In Cases C3 and C4 ( Fig  6B), in which the magnet is centered with respect to the microchannel, the captured particles are initially located close to the top and to both lateral walls because of the relatively low velocity of the flow near these lateral walls. Fig 7 shows the regions with high rates of particle deposition for cases L1 (Fig 7A-1 and 7A-2) and C1 (Fig 7B-1 and 7B-2) for which all the particles are deposited on the wall. The time is larger than 525 non-dimensional units which is large enough to allow the major part of the particles to reach their final location of deposition (i.e. one throughflow corresponds to 100 nondimensional time units). To compute the deposition rate the computational domain has been divided into 49x499x19 equal volumes for the lateral case and into 99x499x19 equal volumes for the centered case. Note that initially a particle is located in the center of the volumes situated at the inlet of the computational domain. At each time step the number of new particles in each volume is stored and summed up during the simulation. Thus, this quantity represents the rate of increase of the number of particles in each volume per time step, or equivalently, the rate at which the concentration of particles increases with respect to the concentration of  Table 2).   Table 2). The vertical scale of the channel has been enlarged for better visualization and the projection of the magnet located between 25<x * <75 is shown in grey.
doi:10.1371/journal.pone.0151053.g007 particles at the inlet per time step. The isosurfaces plotted in 7Fig A and 7B correspond to a value of 16 and 6, respectively. In the case of the lateral position of the magnet ( Fig 7A) the particles are deposited near the lateral wall closest to the magnet, while for the centered position ( Fig 7B) the particles are deposited on the top wall. In both cases the simulations predict the location of the accumulations below leading edge of the magnet near the first half of the streamwise length of the magnet.

Conclusions
In this study we analyzed numerical simulations of the trajectories of magnetic beads in a straight microchannel under the influence of a nearby permanent cubical magnet for different flow and magnetic conditions. Analytically derived local fluid velocities and local magnetic forces have been used to track the particles. A centered position and a lateral position of the magnet above the microchannel are considered. It has been found theoretically that the fraction of the particles deposited on the walls of the microchannel depends on the Stokes and Alfvén particle numbers and that the size of the magnet should be larger than the height of the microchannel divided by the Stokes and Alfvén particle numbers to capture all the particles uniformly distributed across the section of the microchannel. The results of the numerical simulations are in agreement with this criterion. The Lagrangian tracking of the particles has shown that accumulations of deposited particles occur on the top wall of the microchannel for the centered position of the magnet, while for the lateral position the accumulation is located on the lateral wall closest to the magnet. In both cases the simulations predict the location of the accumulations below leading edge of the magnet near the first half of the length of the magnet.