Computational modeling and fluorescence microscopy characterization of a two-phase magnetophoretic microsystem for continuous-flow blood detoxification

Magnetic beads can be functionalized to capture and separate target pathogens from blood for extracorporeal detoxification. The beads can be magnetically separated from a blood stream and collected into a coflowing buffer solution using a two-phase liquid-liquid continuous-flow microfluidic device in the presence of an external field. However, device design and process optimization, i.e. high bead recovery with minimum blood loss or dilution remain a substantial technological challenge. We introduce a CFD-based Eulerian-Lagrangian computational model that enables the rational design and optimization of such systems. The model takes into account dominant magnetic and hydrodynamic forces on the beads as well as coupled bead-fluid interactions. Fluid flow (Navier-Stokes equations) and mass transfer (Fick’s Law) between the coflowing fluids are solved numerically, while the magnetic force on the beads is predicted using analytical methods. The model is demonstrated via application to a prototype device and used to predict key performance metrics; degree of bead separation, flow patterns, and mass transfer, i.e. blood diffusion to the buffer phase. The impact of different process variables and parameters – flowrates, bead and magnet dimensions and fluid viscosities – on both bead recovery and blood loss or dilution is quantified for the first time. The performance of the prototype device is characterized using fluorescence microscopy and the experimental results are found to match theoretical predictions within an absolute error of 15%. While the model is demonstrated here for analysis of a detoxification device, it can be readily adapted to a broad range of magnetically-enabled microfluidic applications, e.g. bioseparation, sorting and sensing.


Introduction
The separation and detection of biomaterials from complex media such as biological fluids is fundamental and often challenging in multiple areas within biology, biotechnology and medicine. 1,2 The removal of toxic substances from blood (biotoxins, microorganisms, heavy metals, radioactive toxins or drugs) is of particular interest since their presence in the bloodstream can cause deleterious metabolic alterations and even death. [3][4][5] The removal of the diseasecausing agents from blood can be considered as the most direct treatment for different clinical conditions such as intoxication, bacteraemia or autoimmune diseases. 6,7 Among the different extracorporeal methods, the use of functionalized magnetic beads as toxin sequestering agents has been considered one of the best alternatives since it can improve the effectiveness of the treatment while reducing adverse side effects. 4 The interest in extracorporeal detoxification with magnetic beads stems from the limitations of the traditional methods (hemodialysis, hemofiltration, plasmapheresis, extracorporeal imunoadsorption, direct injection of chelators and antibodies, etc.). 8,9 Furthermore, advances in magnetic bead synthesis and functionalization enable multiple surface decorations that can be tailored to bind with a broad variety of blood borne pathogens. [10][11][12][13][14][15] Extracorporeal blood detoxification using functionalized magnetic beads is a two-stage process as seen in Fig. 1. First, the beads are mixed with blood taken from the patient and selectively bind with target toxins thereby forming toxin-bead complexes within the blood. In the second stage, the toxin-bead complexes are magnetically separated from the blood and the resulting toxin-free blood solution is returned to the patient's circulatory system. The ultimate goal of this process is to selectively remove pathogens from a patient's blood without altering blood constituents and function. Numerous magnetophoretic microseparators have been designed for the recovery of magnetic beads from biological fluids. [16][17][18][19][20] Recently, some studies have demonstrated the advantages of twophase continuous-flow separators, where the beads are deflected from the blood stream and collected into a separate coflowing buffer solution by a magnetic force applied perpendicular to the flow direction. 21 However, the complexity of these systems is high and many design challenges exist that need to be addressed. For example, Lee and coworkers 22 employed a dual inlet microfluidic system for the removal of magnetic nanoparticles (attached to E. Coli) from bovine blood, and although 25% of the bacteria was collected in a saline buffer, they found difficulties in maintaining two symmetric laminar flows because the difference between the saline and blood viscosities was very large (1 cP for saline and 10 cP for blood). To address this problem, they employed a single-inlet, dualoutlet device, where blood was infused and the particles remained trapped at the walls (i.e. a batch device). Yung et al. 23 employed a similar design for the separation of fungi from blood, and although they cleared 80% of bead-bound-fungi, significant blood loss (50-60%) was observed when the flow rates of both solutions (aqueous PBS and whole blood) were the same. The use of a more viscous solution as the buffer did not improve the phase separation, but instead, by modifying the flowrate ratio between the two fluids the blood loss was reduced.
On the other hand, Kang et al. 24 employed a two-phase microfluidic channel to separate bacteria attached to magnetic nanoparticles (128 nm) and found serious difficulties in recovering all the particles at relatively high velocities (80%). Although these authors did not observe blood dilution or loss, their particles were very small to obtain an adequate separation and thus, they had to add to their system micron-sized beads (1 µm) that act as local magnetic field gradient concentrators to magnetize and attract the smaller particles when exposed to an external magnetic field. Finally, none of the previous studies have considered the diffusion of blood components to the buffer solution as they flow side-by-side inside the channel, which degrades the quality of the biofluid. Therefore, achieving high bead recovery while minimizing blood loss or dilution and avoiding diffusion between fluid phases remains a technological challenge.
It should be noted that these previous studies, including the current work, focus on positive magnetophoresis. This method exploits the different magnetic properties (susceptibility values) between the particles and the surrounding fluid, including the cells present in blood. Positive magnetophoresis has been extensively applied for separating paramagnetic, ferromagnetic or superparamagnetic particles suspended in non-magnetic fluids, including biofluids, and especially, in batch separators. [25][26][27] In this case, the particles are attracted to the magnetic field as the surrounding media is diamagnetic.
However, there is another magnetic separation technique called negative magnetophoresis in which non-magnetic materials (including beads, bacteria, cells, etc.) are manipulated from magnetic fluids (ferrofluids or solutions of paramagnetic salts such as manganese or gadolinium solutions). 27 In this case, the target components can be considered as "magnetic holes" and are repelled from the magnetic field towards an area of field minima. 28 This technique is very interesting from a clinical point of view, since harmful bacteria could be separated from blood in a label-free manner (no incubation would be required prior to the magnetophoretic processing). Nonetheless, its application for blood detoxification is potentially limited for several reasons.
First, compatible magnetic solutions should be synthesized (for example, MgCl2 solutions commonly used in diamagnetophoresis are non-biocompatible at the concentration levels required for cell separation). In addition, blood is complex and comprised of billions of cells with sizes that range between 5-15 µm. 28 Some of these cells might separate from the blood sample and migrate along with the target biomolecules towards their outlet, as the separation is based mostly on differences in size. Finally, the magnetic fluid cannot be reintroduced into the patient and therefore, post-processing stages would need to be applied in order to restore blood properties and components. However, this method has proven to be very efficient for diagnostics and clinical research (i.e. detection and isolation of circulating tumor cells, HeLa cells, sickle red blood cells, bacteria, etc.). [28][29][30][31][32] More information about the fundamental differences between both techniques can be found in Ref. 27.
The progress of positive magnetophoresis technology as applied in blood detoxification requires the following; (a) the complete capture of the magnetic bead-toxin complexes using appropriate nonuniform magnetic fields (i.e. those that exhibit a gradient in the field distribution) and (b) the control of the neighboring coflowing fluids (blood and buffer) without intermixing the two inside the device.
In this article we introduce a CFD-based Eulerian-Lagrangian computational model for predicting the overall performance of a two-phase liquid-liquid microfluidic separator. The model predicts fluid flow within a fixed-mesh Eulerian framework, magnetic particle dynamics using a Lagrangian formulation, and mass transfer, i.e. blood diffusion to the buffer phase, using Fick's law. It takes into account the dominant magnetic and hydrodynamic forces on the particles as well as fully-coupled particle-fluid interactions that give rise to interparticle hydrodynamic coupling and particle motioninduced mixing of the coflowing fluids. We use the model to study the performance of a prototype microfluidic device with a Y-Y (inputoutput) flow configuration that has been employed for the recovery of magnetic beads from a flowing blood solution, as depicted in Fig.  1. Specifically, we investigate the recovery of micron-sized magnetic beads from human whole blood and their collection into an aqueous buffer. This work includes an analysis of bead dynamics, flow patterns, and mixing of the co-flowing blood and buffer fluids. The impact of different process variables -flowrates, bead and magnet dimensions and fluid viscosities -on bead recovery is quantified for the first time. We also characterize the device performance using fluorescence microscopy and use this data to validate the model. The experimental results are compared with corresponding theoretical predictions and found to be within an absolute error of 15%. The modeling approach is summarized below and a more detailed discussion of the theory can be found in our previous works. 25,[33][34][35] Finally, the model is well-suited for parametric analysis and optimization, thereby facilitating the development of novel microfluidic systems, not only for blood detoxification, but also for many other applications that utilize magnetophoresis such as bioseparation, sorting, sensing and drug discovery, etc.

Modeling Approach
In this section we summarize key details of the computational model. The model was developed by customizing a commercial multiphysics CFD software program, FLOW-3D from Flow Science Inc. (ver11.2, www.flow3d.com). Specifically, custom magnetic field and force algorithms were integrated into the program to predict the field distribution from various field sources and the corresponding force on magnetic particles. FLOW-3D has two preprogrammed models, the "Particle" and "General Moving Object" (GMO) models that predict particle dynamics taking into account various forces (e.g. hydrodynamic and gravitational) and coupled particle-fluid interactions (i.e. two-way momentum transfer) wherein the moving particles perturb the flow field. However, as of this writing, the program had no magnetic analysis capability, which is what we added. FLOW-3D solves for fluid flow (Navier-Stokes equations) using a Volume of Fluid (VOF), finite-difference based, Eulerian fixed grid approach, while a Lagrangian framework is used to model particle dynamics within the fluid. The model is based on the following assumptions: (a) all fluids are Newtonian and incompressible (blood is often assumed to be Newtonian with a viscosity value between 3 and 4 cP for medium to high shear rates), 25,36 (b) the magnetic particles have a linear magnetization curve with saturation, (c) interparticle magnetic dipole-dipole coupling is negligible because of a low particle concentration, (d) the field sources are ideal 3D rare-earth permanent magnets, and (e) there are no other magnetic materials present in the computational domain that would otherwise perturb the magnetic field. This last assumption is important as it allows for the use of analytical expressions for the magnetic field distribution and force. This eliminates the need for numerical field predictions, which greatly simplifies the analysis and reduces simulation time.

Magnetic analysis
We use a Lagrangian approach to predict the motion of the magnetic beads. In this analysis the beads are modeled as discrete objects and their trajectories are predicted in accordance with classical Newtonian dynamics: 27 Here, mp and vp are the mass and velocity of a bead and Fext represents all external force vectors exerted on it.
As previously noted, FLOW-3D has various forces preprogrammed, but not the magnetic force. We modeled the magnetic force Fm acting on a bead using the "effective" dipole moment method in which the bead is replaced by an "equivalent" point dipole with a moment mp,eff at its center as described by Furlani. 18,25,27,37 Fm is given by: where represents the permeability of the free space (4π·10 -7 H·m -1 ) and Ha is the applied magnetic field intensity at the center of the bead. The moment mp,eff can be determined using a magnetization model that takes into account self-demagnetization and magnetic saturation of the beads: where the function f(Ha) is calculated as follows: where Vp is the volume of a spherical bead of radius rp, , is the effective susceptibility of the bead that can be related to the intrinsic susceptibility 38,39 and Ms,p represents the saturation magnetization of the beads.
Thus, the expression for Fm takes into account conditions both below and above the saturation of the bead as a function of the external applied field and can be written as: In order to solve Eq. (6), an expression is needed for the field distribution generated by an external rare-earth magnet. The field and gradient produced by the magnet depend on its material properties, shape and dimensions. Three different commercially available rectangular-shaped rare-earth magnets were used for this study. These have dimensions LxHxW equal to 5x5x3, 8x8x4 and 10x5x3 mm 3 , respectfully.
To determine the 3D field distribution of the magnets, we employed an analytical model developed by Furlani. 40 It should be noted that we can use an analytical field expression because we have assumed that there are no other magnetic materials in the computational domain, which is true for the prototype system that we study. A closed-form expression for calculating the field of the magnet geometry shown in Fig. 2a, with its magnetization in the z direction, is: Here, Ms is the saturation magnetization of the magnet, i.e. 10 6 A·m -1 , which was estimated based on experimental magnetic field measurements using a custom MatLab program. In the field equations, x1=L-/2, x2=L/2, y1=-H/2, y2=H/2, z1=-W/2 and z2=W/2 represent one half of the magnet's dimensions (L, H and W) in the x, y and z directions, respectively, when the coordinate system centered with respect to the center of the magnet. It should be noted that in our subsequent device performance analysis we use a coordinate system that is centered at the center of the microchannel, as presented in Fig. 2b. Given the expressions for the field components, we can determine the magnetic force on a particle at any point in the domain.

Hydrodynamic force
The hydrodynamic force Fhd on the beads is predicted numerically using the following expression: where P is the pressure, Madded is the added mass (additional resistance for an accelerating or decelerating body), v is the fluid velocity at the location of the particle and Fdrag is the drag force acting on the beads. The amount of added mass can be calculated using the following equation: where the coefficient CM is theoretically predicted to be equal to 0.5 and ρ is the fluid density. 41 Finally, the drag force Fdrag can be obtained using a modified form of Stokes' approximation for the drag on a sphere: where AP is the bead cross sectional area, which can be written as AP = πrp 2 . CD is the drag coefficient for steady-state flow around a sphere and can be calculated from: where Rep is the particle Reynolds number for which the expression is valid. 42 More specifically, we have: where dh is the hydraulic diameter and η is the fluid viscosity.
We investigate the effect of the fluid properties on magnetic separation by employing viscous aqueous solutions (η values of approximately 0.001 and 0.004 Pa·s) and human whole blood (η=0.0035 Pa·s). The buffer solution used as the collecting phase in the bioseparator is modeled as water at 20ºC with a viscosity of 0.001 Pa·s. Finally, to calculate Fhd we need an expression for the fluid velocity, which is developed in the next subsection.

Flow and mass transfer
The fluid velocity is predicted numerically in FLOW-3D by solving modified Navier-Stokes and continuity equations (Eqs. 15 and 16) that account for the interactions between fluid and magnetic beads (see FLOW-3D user's manual for details): Please do not adjust margins Please do not adjust margins where the div(τ) term represents the viscous accelerations and V is the volume of fluid in a computational cell. The last term in Eq. 15 represents particle induced fluid accelerations and FP is given by: These equations are solved numerically using a finite differencebased Volume of Fluid (VOF) method (see FLOW-3D user's guide). As previously noted, all fluids are assumed to be incompressible.
In addition, the model also solves the mass transport of a species between both phases by solving the continuity equation and Fick's law. The presence of a solute is defined as a property of both fluids (advectable scalar) and it is treated as an interstitial solute. The advectable scalar in this work is defined as the concentration C, i.e. the mass per fluid volume in a cell and the transport equation that is solved is given as follows: where Vf and A are the fractional volume and area open to flow for each mesh cell, respectively, and D is the diffusion coefficient of the solute. In this work, we analyze the component diffusion between phases through their interface as they flow within the channel according to the component's diffusion time. This time represents the time required for partitioning of solutes between the phases and is calculated as: [43][44][45] = (19) where Hi is the height of the channel that phase "i" is flowing through and "j" represents the number of dimensions in which diffusion takes place. Thus, if the residence time exceeds the diffusion time of the solute inside the device, the component will have enough time to diffuse to the co-flowing phase. Therefore, we study the diffusion at different values of the ratio γ, described as follows: where Lc represents the total length of the chip (Lc=2 mm in our prototype) and � the average velocity of the phase 1 (inlet velocity), whereas H1 is the height of the branch where this phase flows (H1=150 µm). Fig. 2b shows the x-z plane of the microchannel used in the simulations. Although the force balance for each bead is solved within a 3D analysis (the magnetic field distribution and magnetic force is calculated in x, y, and z directions), the governing equations for the flow along the width of the channel (y axis) are not considered. The channel height (Hc) and length (Lc), in z and x directions, are 300 µm and 2 mm, respectively. Hc was modified for the simulations of particle recovery from blood as the glass microdevice employed for that experiments is 280 µm wide.

Methodology
The magnets are positioned at a distance d=200 µm from the lower channel wall with the left corner of the magnet placed at the beginning of the chip's straight section, at x=-Lc/2, as shown in Fig.  2b. This location in x direction was selected due to the higher magnetic gradients that occur at that position, and also, because of the easier visual positioning under the microscope. Beads enter the upper inlet (equally distributed along the diameter), at a constant concentration of 0.1 g·L -1 (0.2 g·L -1 for bead recovery from blood), and are tracked at every time step. Interbead effects (e.g. magnetic dipole-dipole coupling) are neglected due to the low concentration, i.e. values in the typical range of relevant detoxification studies. 4,10,46,47 The performance of the separator is evaluated through a dimensionless number J 48 that describes the relationship between the magnetic and fluidic drag forces. Thus, the J number scales both the magnetic force (in the z direction) and drag force (in the x direction, obtained using Stokes' approximation for the drag on a sphere 27,49 ) on the beads located in the middle plane of the channel. This location was chosen because of the average value of the magnetic force acting at that point. The J number can be written as follows: To calculate J, the magnetic force was calculated at the mid-plane of the channel (z=0) for every simulation using a separate MatLab program. For computing the denominator of the number, the product η1 � for phase 1 was chosen because the beads are initially located in the phase 1 flow, and therefore, when they cross the interface and move to the collecting phase, the separation is successful. Thus, the numerator of this number provides the average Fm,z that pulls the beads downwards, whereas the denominator represents the average Fdrag,x acting on the beads during separation.
As previously noted, the flow field and equations that govern bead motion were solved using FLOW-3D. We integrated the magnetic analysis, the magnetic field and force, via the development of an external Fortran subroutine, which was compiled in Visual Studio v.2013 (Microsoft). The subroutine passed the components of the magnetic force to FLOW-3D at each time step, which were incorporated into the force balance within the bead analysis model. Furthermore, MatLab v.2015 (The MathWorks, Inc.) was also used for rapid analysis of the magnetic field and force at different points of the channel.

Experimental methodology
The experiments for validating the model with aqueous solutions were performed using a prototype SU-8 Y-Y shaped device with a rectangular cross-section of 300 µm x 200 µm (Microliquid) as shown in Fig. 3a. This chip was placed in a methacrylate holder, and tygon tubing (ø=0.8 mm) was employed for the fluidic connections (Fig. 3b).
For the experiments with human whole blood, a glass microdevice with a U-shaped cross-section was employed (Fig. 3c). This device was fabricated with conventional photolithography and wet etching techniques. 50 The glass chip was mounted into an aluminum holder and connected to the pumps with fused silica capillaries (150 µm i.d., 375 µm o.d., Polymicro Technologies). Two syringe pumps (Legato 210 infuse/withdraw syringe pump, Kd Scientific) and single syringes (Omnifix 5mL luer, BRAUN) were used to control the flow.
For the fluidic analysis, fluids with different characteristics, i.e. two different fluid phases, were pumped into the two inlets, respectively: a) Phase 1 -Fluorescein sodium salt solution (employed as a coloring agent to distinguish the phases), (Scharlau) or human whole blood (explained below), and b) Phase 2 -deionized water or a mixture of ethylene glycol in water 50% v/v (PanReac AppliChem), that exhibits a density and viscosity values similar to human blood (1050 kg·m 3 and 4.1 cP at 20ºC). 51 For the mass transfer analysis, we employed a specific Phase 1 fluorescein solution (60 mg·L -1 , diffusion coefficient for fluorescein in water: 4.25·10 -6 cm 2 ·s -1 ), 52 whereas deionized water was pumped into the phase 2 inlet. To study the possible migration of blood cells between the fluid phases, we prepared aqueous solution of nonmagnetic 6-8 µm polystyrene beads (Spherotech, Inc., catalog number HISP-60-2) at a concentration of 1 g·L -1 . This solution was pumped into the upper inlet, whereas deionized water was introduced into the device through the lower outlet. Bead concentration was measured under the microscope after sample collection at the outlets.
For the magnetic bead recovery analysis, we prepared aqueous solutions of fluorescent magnetic beads of two different sizes, i.e.   Finally, a set of experiments with human whole blood was performed in order to evaluate the device performance under a more realistic scenario. Human whole blood from healthy subjects was collected with informed consents according to a protocol approved by the Advanced Separation Processes Group following the University of Cantabria guidelines. The samples were collected into 3.5 mL tubes with anticoagulant (sodium citrate 3.2%) and was processed within 5 hours of blood draw. Before the experiments, saline buffer (NaCl 9 mg·mL -1 ; Fresenius, Kabi) with bemiparin sodium (Hibor®, Rovi Pharmaceuticals Laboratories) was pumped into the device through both inlets for 10 minutes. For the co-flow experiments, blood was loaded into a 2-mL syringe and injected into the glass device at different flow rates. For the magnetic bead recovery, blood was spiked with the 4.9 µm magnetic beads at a concentration of 0.2 g·L -1 . Saline buffer was introduced into the device through the lower outlet for both studies. The syringes were agitated every 5-10 minutes in order to prevent blood cells from settling.
Although the size of the magnetic beads employed in this work was specified from the manufacturer, the susceptibility and saturation ( , , Ms,p) were not available. We performed a parametric estimation of these values from an experimental test and found that the saturation magnetization has more impact on the device performance. The values obtained were Ms,p = 1.5·10 4 A·m -1 and , = 0.25, which fall within the range of beads employed for similar bioapplications. 53,54 The magnetic field was provided by rectangular magnets (5x5x3, 8x8x4, 10x5x3 mm 3 ; Supermagnete, Germany). They were made from grades 45N-52N neodymium iron boron (NdFeB), with remnant magnetizations of Br=1.32-1.47 T. Although the company specifies the range of the saturation magnetization of the magnets, we nevertheless estimated the actual value from field measurements using an analytically based MatLab model. More specifically, the magnetic field produced by one of our magnets was characterized using a 3D magnetic field mapping instrument, the MMS-1-R from SENIS GmbH (www.senis.ch). A three-dimensional probe was used to scan the magnetic field at different z values above the center of the upper surface of magnet (Bx, By =0 at that points), and these values were employed for the estimation of the saturation (Fig. 4a). A value of 1013.1 kA·m -1 was obtained and used for all the simulations. The theoretical predictions are in excellent agreement with the measured data as seen in Figs. 4b and 4c where the comparison between them is shown for an x-y plane at z=1.5 mm from the magnet surface.
Fluorescence intensity caused by the presence of fluorescein or beads in the fluid phases was detected using a microscope (Nikon SMZ18) equipped with a green fluorescence filter (light wavelengths of around 550 nm) and a Jenoptik ProgRes C5 camera. Images were taken using the ProgRes® CapturePro software (CapturePro V2.10.0.0). Whereas the analysis of the flow patterns and diffusion phenomena was easily performed by analyzing the photographs taken with the camera (fluorescein diffusion is even visible to the naked eye), it was not possible to quantitatively calculate the bead separation effectiveness from the photographs. Therefore, we developed a custom code for image analysis using MatLab. The image processing is shown step-by-step in Fig. 5. Firstly, the user selects the vertice points delimiting the input and output channels using an image of the microchip taken under visible light (Fig. 5a). These points are used to delimit the channels in the following images taken under fluorescence light, and the channel limits are superimposed for visual assessment (Fig. 5b). From the fluorescent image, the code automatically splits the color components (Red (R), Green (G), and Blue (B)). Afterwards, the program enhances the G component (as fluorescence is emitted in green) by removing the R component that generates noise in the images (Fig. 5c). At this point, an intensity threshold is established and manually tuned until the pixels detected in the clean G component correspond to the beads inside the channel (Fig. 5d). Finally, the code accounts for the number and area of beads that exit the channel through each outlet in the Y junction by counting the areas of green pixels detected in the outlets. Thus, the experimental bead recovery can be estimated as follows: It should be noted that we have checked the accuracy of the developed method by comparing the calculated recovery of the beads (measured with the software) and the concentration measured at the outlets after sample collection and processing for one of the experimental tests of this work. We found that the results obtained with both methods differ by less than 10%, which validates the accuracy of our customized image processing program. Since this method has been automatized, once the picture of the device is taken, it takes only 10 seconds to obtain the recovery percentage of the device, in comparison to the sample collection, bead reconcentration and measurement process, which is much more tedious and time consuming. Furthermore, it can be used for realtime process monitoring and even integrated into a control loop for optimization, as a consequence of the continuous mode of the process.
For each experimental condition, different tests were performed and an average of 5 different micrographs were taken and analyzed. The value reported in the results section represents the average recovery. Nonetheless, for the experiments with human whole blood, the bead recovery was calculated by measuring bead concentration at the outlets after sample collection, since the high concentration of blood cells might interfere with the illumination of the beads in the upper branch and lead to different fluorescence intensity values in comparison with the co-flowing aqueous solution.

Fluid phase flow patterns and mass transfer
In this section, the flow patterns of phases with the same and different fluid properties are analyzed in order to examine the flow conditions that lead to the complete separation of the fluids at the bioseparator outlets. Furthermore, we study the flow conditions that maximize and minimize the interdiffusion between streams.
In Figs. 6a and 6b, the volumetric fraction occupied by aqueous phases (fluorescein solution -deionized water) inside the device when the flow rates are the same is presented. When both phases flow through the channel at the same mean velocity (set in this case at ≈ 0.5 cm·s -1 ), independent co-flow of the phases is achieved and mixing at the outlets is not observed. Furthermore the experimental images are in agreement with our theoretical predictions as seen in Fig. 6b. When a more viscous solution was pumped into the channel, i.e. fluorescein solution -ethylene glycol 50% v/v mixture or human blood -saline buffer, such behavior is not observed. However, as demonstrated in our previous work, 48 when the two fluids have a different viscosity there is a negative impact on the co-flow and phase separation.

Please do not adjust margins
Please do not adjust margins This is because the difference in viscosities between the fluids creates different pressure gradients in the respective branches of the channel, which then translated into the mixing at the outlets. Nevertheless, this can be minimized by ensuring the same pressure drop in each phase throughout the length of the channel. With the Hagen Poiseuille law the separation can be improved when the product (η �)i for each phase i is the same. Figs. 6c and 6d show the fluorescein -ethylene glycol separation at the outlets when this criterion was applied (average fluid velocities at inlets are 0.5 cm·s -1 and 2.1 cm·s -1 for ethylene glycol and water solutions, respectively). Furthermore, Fig. 6e shows the blood-buffer co-flow when the saline buffer is pumped at a flow rate almost 3 times higher than blood, and Fig. 6f demonstrates how the separation is substantially improved using this criterion in comparison with the use of the same flow rates for both fluids. This behavior was also observed in other studies. 23,55,56 When working with complex and multicomponent solutions (biological fluids such as blood), not only the solution loss or dilution at the outlet should be analyzed, but also, the possible diffusion of its constituents to the buffer solution. We analyzed the diffusion of fluorescein from phase 1 to phase 2 as a function of γ (ratio between residence and diffusion times as presented in Eq. (20)). In Fig. 7a we show the simulation of fluorescein diffusion for γ values close to 1. In this case, the residence time is approximately equal to the diffusion time (≈ 40 s) and the fluorescein has time to completely diffuse to phase 2. Thus, the mass transfer under this velocity condition is complete as the concentration reported at both outlets is the same (30 mg·L -1 ), which is half of the initial value at the phase 1 inlet (60 mg·L -1 ). As can be seen from the experimental image (Fig. 7b), the fluorescence intensity at the lower outlet is approximately the same as that reported at the upper outlet. We calculated the change in the fluorescence intensity across the channel width for different downstream positions x=-1, 0 and 1 mm (Fig. c). The plot demonstrates the diffusion phenomena for γ≈1 as the fluorescence intensity across the width drastically varies at the interface for x=-1 mm (inlet), but a small reduction in the value is found for the positions x=0 and x=1 mm (the intensity fluorescence across the width almost remains constant at x=1 mm). The different intensity fluorescence at channel width equal to 0 for the three x positions might be due to the inhomogeneous illumination conditions inside the device. In Figs. 7d and 7e we show the diffusion of fluorescein for γ values close to 0.01. In this case, the residence time is much lower than the diffusion time (≈ 0.4 s) and the fluorescein diffusion from phase 1 to phase 2 is neglected. As seen from the simulation (Fig. 7d), the fluorescein concentration at the lower outlet is around 6 mg·L -1 , 80% lower than the value reported when γ was approximately 1. Furthermore, the fluorescence intensity at the lower outlet is not perceived from the experimental image (Fig. 7e). For this γ value, a sharp change in the fluorescence intensity is reported at the interface between both fluid phases for all downstream positions due to the low diffusion reported between streams as perceived from Fig. 7f.
Finally, the migration of big macromolecules was investigated by using non-magnetic polystyrene beads. We employed beads with sizes of approximately 7 µm and observed their trajectories along the device. This is important as it provides some insight as to the transport of blood components, such as red blood cells (6-9 µm) or white blood cells (8-14 µm). 28 We simulated the trajectory followed by these 7 µm polystyrene beads (a flow rate range of 0.005-0.5 µL·s -1 was studied and no difference was observed over this range). The beads followed the fluid streams depicting a straight line from the inlet to the outlet, i.e. no deflection was found for any of the flow rates tested. This result is in good agreement with our experimental data and also expected under these flow conditions, since the 10 | J. Name., 2012, 00, 1-3 This journal is © The Royal Society of Chemistry 20xx Please do not adjust margins Please do not adjust margins diffusion coefficient of macromolecules of that size (blood cells) is around 10 -14 m 2 ·s -1 as calculated by the Stokes-Einstein equation.
Thus, working with γ values of approximately 0.01, the diffusion of relatively small molecules (such as fluorescein) and big macromolecules (blood proteins or cells) can be considered negligible due to their low diffusivity value, and thus, high diffusion time (the diffusion time increases up to several minutes for these macromolecules, 48,57 whereas the residence time of our fluid phases inside the bioseparator is less than 4 seconds for all the experiments performed for studying the magnetophoresis). Therefore, it can be assured that diffusion of blood components to the buffer solution is negligible at the time scales of interest.

Bead magnetophoresis
In this section the recovery of magnetic beads is analyzed under different magnetic field and fluid flow conditions. Different variables (flow rates) and parameters (magnet and bead dimensions, as well as the rheological properties of the phases involved) are studied in order to optimize the force balance acting on the beads.

Influence of magnet dimensions
The influence of the magnet dimensions on the magnetic field and force in the microchannel was studied using three different-sized magnets: 5x5x3 mm 3 ; 8x8x4 mm 3 and 10x5x3 mm 3 . The magnetic field in z direction provided by the three magnets inside the microchannel is analyzed in Table 1, for a x-y plane (limited by the chip contours) located at the midline of the channel, i.e. z=0. The magnetic force in z direction generated on 2.29 µm and 4.9 µm beads with the three magnets at the same plane is also reported in Table 1.
Only the z component of the field and force is presented since it is the component that causes bead deflection from phase 1 to phase 2.
As reported in Table 1, the magnetic field generated in the channel is nonuniform, with values that range from -10 to 40 mT (5x5x3 mm 3 ), from -20 to 40 mT (8x8x4 mm 3 ) and from -25 to 25 mT (10x10x5 mm 3 ). Although these values are not equal, the difference between the positive and the negative values is practically the same, so the magnetic gradient is very similar for the three structures. This is verified by examining the magnetic force generated by those magnets on the beads. In this case, the magnetic force acting on the 4.9 µm beads ranges from -0.1 to -0.22 nN for the 5x5x3 mm 3 magnet, and from -0.1 to -0.24 nN for the 8x8x4 mm 3 and 10x10x5 mm 3 magnets. Therefore, the magnetic bead recovery is expected to increase slightly as the magnet size increases due to the higher magnetic force.
Following the magnetic field and force analysis, bead separation under conditions of variable hydrodynamic forces with the magnetic force held fixed is discussed. In Figs. 8a, 8b and 8c, the theoretical and experimental bead recovery (4.9 µm beads) from aqueous solutions as a function of the applied flow rate is presented for the 5x5x3 mm 3 , 8x8x4 mm 3 and 10x5x3 mm 3 magnets, respectively. As seen in the figures, and due to the similar magnetic force generated by the 3 magnets, the bead recovery slightly increases with the size of the magnet, i.e. the theoretical bead separation increases an average of 13% when the 5x5x3 mm 3 magnet is replaced by the 8x8x4 mm 3 or the 10x5x3 mm 3 magnets.  For each figure, the magnetic force is kept constant for all the flow rates. However, the hydrodynamic force increases with the flow rates and when this value is higher than 1.1 µL·s -1 , bead recovery decreases by as much as 25%. Note that the experimental recovery is around 13%, 15% and 26% for the 5x5x3 mm 3 , 8x8x4 mm 3 and 10x5x3 mm 3 magnets, respectively. For these flow rates, the ratio J between magnetic and fluidic drag forces is around 0.1. On the other hand, as the flow rate is reduced, the magnitude of the fluidic drag force decreases as well, and for values equal to 0.25 µL·s -1 , the J ratio ranges from 0.64 to 0.81 (depending on the magnet). Under these low velocity conditions, bead recovery increases up to 89% for the 5x5x3 mm 3 magnet and up to 100% for the 8x8x4 mm 3 and 10x5x3 mm 3 magnets.
Furthermore, the simulated results follow the same trend as the experimental ones although there is a deviation between the two at some points of the curves for the 5x5x3 mm 3 and 8x8x4 mm 3 magnets. The difference between the experimental and simulated results might be caused by the inhomogeneous illumination conditions inside the device for each experiment or by an error in the image analysis. However, in almost all cases there is a satisfactory global value with an absolute error around 15%.

Influence of bead size
As presented in Eqs. (3) and (12)(13)(14), the bead size affects the fluidic drag force (proportional to the bead radius) and the magnetic force (proportional to the bead volume). Therefore, we analyzed the influence of bead size (i.e. the 2.29 and 4.9 μm beads) on magnetic recovery. Firstly, we simulated the magnetic force field acting on the beads (Table 1) inside the SU-8 chip. From this table, it can be seen that the magnetic force (in the z direction) is reduced by one order of magnitude when the bead size is reduced to half. This reduction in the magnetic force greatly affects the magnetic recovery.
In Fig. 9, the magnetic recovery of both bead types as a function of the applied flow rate is presented. As plotted in the figure, only 20% of the smaller beads are recovered from the solution for flow rates as small as 0.25 µL·s -1 . This is because for these conditions, the ratio J takes values smaller than 0.1. It should be noted that this J value for the bigger beads is obtained for flow rates higher than 1.3 µL·s -1 . Whereas the magnetic force is diminished by one order of magnitude for the 2.29 μm beads, the drag force is approximately reduced by half, but remains at the same order of magnitude as compared to 4.9 μm beads. The different force balance acting on both beads explains the different recovery ratios. For the 2.29 μm beads, the separation decreases between 20% and 77% for the flow rate range under analysis, and only by using flow rates smaller than 0.1 µL·s -1 , full recovery can be achieved. Finally, there is a good agreement between the theoretical predictions obtained with the CFD model and the experimental data, with an absolute error lower than 10%.

Influence of viscosity
In this subsection we analyze the effect of the rheology of the biofluids by deflecting the beads from a blood mimicking fluid (ethylene glycol mixture 50% v/v) and from human whole blood, named in this subsection as "BMF" and "HWB", respectively.
In Fig. 10 the separation of 4.9 μm beads from both aqueous and BMF solutions is presented for the SU-8 microdevice. By comparing both scenarios, we find that the separation is negatively affected when a more viscous solution is employed as phase 1. This is due to the drag force, proportional to the viscosity of the solution, which implies that this force is approximately 4 times higher when BMF solutions are treated. For this scenario, bead recoveries of around 37% are obtained for flow rates of 0.3 µL·s -1 (J ratios smaller than 0.15), whereas 100% were recovered from water for this flow rate value. Higher flow rates of the BMF resulted in bead recoveries lower than 30%, values between 23% and 74% lower than the base case in which water was employed as phase 1.
Finally, in Fig. 11 the separation of 4.9 μm beads from both aqueous and HWB solutions is presented. It should be first noticed that for these experiments, a glass microdevice with a different cross-section shape and area was employed, and thus, the flow rates are not comparable to the ones presented in Fig. 10.  Magnetic bead recovery (theoretical and experimental results) from an aqueous solution and from human whole blood (HWB) as a function of the applied flow rate for the 4.9 µm beads using the 10x5x3 mm3 magnet and the glass microdevice.
As happened with BMF, the separation of these beads from HWB is lower than the base case, i.e. in this case, the separation is reduced between 22% and 72% when HWB is employed. These values are very similar to the ones in which a BMF was employed as phase 1.
Therefore, to obtain high recovery values from blood solutions, flow rates of approximately 0.005 µL·s -1 are required (around six times lower than the flow rate required to achieve full separation from water in this device as seen in Fig. 11), which correspond to values of the J ratio of 1.2. This way the dominant forces acting on the beads would remain of the same order of magnitude, which favors the recovery. This flow rate results in an average velocity value of 0.65 mm·s -1 (for this device) with residence times of around 3 s, which is still a small value to allow diffusion of blood macromolecules to the buffer stream (γ ≈ 0.08 for fluorescein). Therefore, for the separation of beads from viscous biofluids such as blood, slower flowrates could be safely employed since the integrity of the biofluid will not be compromised.
Finally, it can be seen from Fig. 11 that there is a good agreement between the theoretical results and the experimental data (error less than 10%), showing the capability of our model to accurately describe magnetophoresis from complex fluids as blood.

Conclusions
In this work, we have introduced a CFD-based computational model for analyzing and optimizing magnetic bead separation from biological fluids in a liquid-liquid continuous-flow microdevice for blood detoxification.
Our main findings have been systematized in order to meet the three technological challenges that should be addressed in real detoxification processes. Firstly, the process should be designed so that the ratio between the flow rates ensures the same pressure drop for both phases along the device length. This way co-flowing of streams can be achieved thereby maintaining two completely separated phases at the channel outlets. Secondly, the residence and diffusion times (γ) should be taken into account in order to evaluate possible interdiffusion between streams. For avoiding such phenomena, γ value should preferably remain below 0.1.
Lastly, for the magnetic bead recovery analysis, the effect of competing magnetic and hydrodynamic forces has been investigated. We showed that these are affected by multiple variables and parameters (flow rates, particle size, rheological properties of the fluids, etc.) that are highly coupled with each other. Consequently, a systematic approach is needed for optimization. We showed that key parameters can be integrated into a dimensionless number J that scales both forces acting on the beads. The obtained data suggests that when the magnetic force is approximately 1.2 times the fluidic drag force, adequate bead recovery efficacies can be reached from human whole blood. Once the conditions that satisfy the three requirements have been identified, process design (parallelization of the devices as a function of the flow rate of blood that needs to be treated) should be carried out.
Finally, we have validated the model via the analysis of prototype devices and have used fluorescence microscopy to characterize key performance metrics: flow patterns, mass transfer between coflowing fluids and bead separation. The experimental data are consistent with the theoretical predictions within an absolute error of 15%. Therefore, we conclude that the computational model enables understanding of the fundamental underlying mechanisms of device performance and is useful for the development of novel magnetophoretic applications. In this case, it contributes to the development and optimization of magnetic bead separation from biological fluids in a multiphase continuous-flow microdevice that can potentially be integrated in a two-step system for extracorporeal detoxification of biofluids, which will be the focus of our future work.

Conflicts of interest
There are no conflicts of interest to declare.