Theory and computational study of electrophoretic ion separation and focusing in microfluidic channels

Abstract. In this work we describe the theory and 2D simulation of ion separation and focusing in a new concept of microfluidic separation device. The principle of the method of ion focusing is classical in the sense that it consists in opposing a hydrodynamic transport ensured by the solution flow to an electrophoretic driving force so that any ionic sample results poised within the microchannel at the point where the two forces equilibrate. The originality of the concept investigated here relies on the fact that thanks to the use of an ion-conducting membrane of variable thickness in electrical contact with the channel the electrophoretic force is varied continuously all along the channel length. Similarly, changing the geometric shape of the membrane allows a facile optimization of the device separation and focusing properties.

A second limitation, valid for any separative microdevice, is related to the fact that characterization of the targeted analyte generally requires that sophisticated detection means (e.g., mass spectrometry) [32][33][34][35][36][37] are connected to the exit of the separation microchannel, thus reproducing the same constraints as in classical separation methods but with a higher difficulty since the amount of analyte is minute.
To obtain a constant flow rate the microchannel must have a constant cross-section.Therefore, in order to create and adjust at will the gradient (dϕ/dx) one needs to divert part of the current flow along the microchannel length [59,60].This may be performed in a continuous mode by carrying out the separation in the presence of an excess of supporting electrolyte and implementing an ion-conducting membrane with nanopores (e.g., Nafion R ) with variable thickness into one of the channel walls as represented schematically in Fig. 1(a).Note that a conceptually analogous approach in which the electrical field experiences a continuous gradient within the microchannel has been recently proposed [61][62][63][64].However, both approaches differ essentially in the sense that in this previous work the variable electrophoretic driving force results from paired electrochemical reactions occurring at each end of a bipolar electrode inserted in the channel.This imposes a strong localization of the electrophoretic force gradient near each end of the bipolar electrode while in the present concept the gradient variations are spread at will over the whole channel length.
Under the conditions of Fig. 1(a), at any given position within the microchannel, the current flow faces two parallel resistances representative of each ionic transport mode, one through the solution and the other through the membrane as sketched in Fig. 1(b) for an element dx of the channel length.If the membrane nanopores are too small for allowing the analytes of interest to be carried through the membrane, they will be transported within the solution while experiencing a gradient of the electrical force which may then oppose more and more efficiently to the counter-transport by the solution flow.The ensuing decrease in transport rate along the electrical gradient proceeds up to the point when the selected analyte necessarily stops and may then be detected by any imaging technique such as fluorescence.
The main goal of this work is to investigate the theoretical aspects of this separation method and evidence by a series of 2D-simulations how such devices perform in particular when two analytes need to be focused at two different points along the microchannel length.For this purpose we will consider that the ion-conducting membrane has been deposited so that its thickness varies linearly along the channel length (Fig. 1a) but evi- dently other shapes may be implemented experimentally or in the simulations to test and predict if this would lead to better separation efficiencies under particular circumstances.Indeed, the present model involving a combination between analytical formulations and simulations is sufficiently versatile to test any configuration, hence to allow optimization of the device towards a particular case of interest.
A computer program written by the authors allows performing 2D simulations for different experimental situations in order to optimize the focusing properties.The simulation can readily account for both parabolic and constant (electroosmotic) flows and can be easily modified to cover more complicated situations.

Physico-chemical model
Consider a flow channel as depicted in Fig. 1(a) of height h and lateral width w (in the direction perpendicular to the cross-section in Fig. 1(a); w h).The solution within the channel is supposed to flow due to an external pressure gradient.In order to create conditions for separating charged particles according to their effective charge an electric field is generated within the channel using two electrodes located at a distance l from each other in the direction of the flow (l h and l w).The potential difference between these electrodes is U = V (x = L) − V (x = 0) and it is assumed that the electrical equipotential planes are vertical within the part of interest of the channel, viz. the field is planar and parallel to the channel length.Finally, for simplification we consider here that the electroosmotic drive due to the electrical field is negligible and that the main drive of electrochemical origin is by electrophoresis (note that in the experiments shown here the solution hydrodynamic flow is ensured only by electroosmotic forces due to the configuration of the experimental set-up).The theory developed below may easily be expanded to consider electroosmotic forces [65][66][67] and streaming potential [68] if this revealed of importance.
This electric field induces an electric current due to the movement of charged analytes and inert ions of supporting electrolyte with respect to the flowing solution.Additionally, the ionic current may also pass through a layer of ion-conducting membrane which is in electrical contact with one of the channel walls and has a varying thickness described by a function f (x).Thus the overall conductance of the solution and the membrane layer as a function of the coordinate x may be adjusted to achieve the required potential (field) distribution.Note that the exact positioning of the membrane vis-a-vis the microchannel is irrelevant within the level of approximation used in this work provided that at least one wall of the channel is electrically connected to the membrane.
Consider a vertical slice (of thickness dx) of the channel together with the ion-conducting wall (see Fig. 1(a)).The overall electric current I flowing through the system is distributed between the solution and ion-conducting film according to their relative resistances.The channel slice can thus be represented by an equivalent circuit shown in Fig. 1(b) where dR s (x) represents the resistance of the elementary solution slice and dR m (x) is the elementary resistance of the ion-conducting membrane slice at a lateral position x.
Denoting the ionic conductivities of the solution and membrane as γ s and γ m (which we assume for now to be constant), respectively, the elementary resistances may be expressed as [69]: The overall resistance of the channel slice is then: From the Ohm's law the potential drop at this elementary channel slice is dϕ = IdR, so that the potential distribution along the channel is described by the following integral: Noting that the potential drop along the whole distance between the electrodes is ϕ(l) = U , one can eliminate the current I and channel width w from (4) to obtain: For example, in the case of a wedge-shaped ion-conducting layer (i.e. with f (x) varying linearly with a wedge-slope k and from a zero value at x = l as shown in Fig. 1(a) the ensuing potential distribution along the channel is logarithmic: This is represented in Fig. 2 for different values of the wedge-slope k and shows that the potential variations and those of electrical field may be adjusted at will simply by changing the wedge-slope of the membrane.This evidences the great versatility of the present concept.Equation ( 5) was obtained under the assumption of constant conductivities of the solution and membrane.However, when this is not the case (viz., when there is a nonuniform distribution of ionic species throughout the solution and membrane due to ion exchange between them) this expression needs to be modified only slightly.Indeed, the resistances (1) and (2) of elementary vertical slices of the solution and membrane layer in this case become: where γ s (x, y) and γ m (x, y) are location-dependent conductivities of the solution and membrane layer, respectively.Therefore the integrand in both the numerator and denominator of (5) becomes: Nonlinear Anal.Model.Control, 2012, Vol.17, No. 4, 431-447 and the potential distribution (5) may be rewritten as: If the solution contains a target molecule (or particle) with the effective charge q the electric field within the channel will exert a force on it equal to: When this force is opposite in direction to the hydrodynamic force pushing the molecule in the direction of the flow, and provided that the channel is long enough, there is necessarily a position, x 0 , along the channel where these two forces equilibrate each other, which corresponds to the point where the analyte will accumulate.Equivalently, the analyte will stop when its hydrodynamic velocity is exactly compensated by the velocity imposed by the electric field.Then from the definition of mobility (u = v/E, where v is the speed of the molecule and E is the field) and the Einstein-Smoluchowski relation (u = |z|F D/RT , where z = q/|e| and D is the diffusion coefficient of the molecule) one readily obtains the relative speed of the charged analyte imposed within the flowing solution by the electrical field [70]: If for simplification we consider that the whole solution flows hydrodynamically at its average velocity, v avg , the charged analyte has an absolute velocity (i.e., vs. the channel walls) given by v avg − v(x) (note that this is essentially correct for an electroosmotic solution flow).This shows that any charged analyte will stop at the position x 0 such that A more exact determination of the position of the focused species sample and its width requires implementing the exact shape of the flow velocity.This will be performed below on the basis of the full analysis of mass transport conditions within the channel and of simulations for the solutions of the ensuing equations.

Mathematical model
The actual distribution of the molecules is described by the following mass-transport equation [13]: where v x (y) is the flow profile function (e.g., parabolic Poiseuille flow) and the formulation of Eq. ( 15) takes into account that the potential distribution is independent of the coordinate y.Note that the consideration of a constant flow profile (viz., corresponding to electroosmotic drive) simplifies the problem at hand by practically reducing it to a single spatial dimension along the channel axis.Therefore we do not consider this situation here and focus below on the parabolic hydrodynamic flow.
It is assumed that initially a plug sample of analyte is injected at some axial position x s (which lies either within the separation device or upstream of it) and has a width w s .Hence, at t = 0 the concentration c = c 0 for x s − w s /2 ≤ x ≤ x s + w s /2 and 0 ≤ y ≤ h and is zero elsewhere.
There is no interaction of the target species with the channel walls, so ∂c ∂y = 0 for y = 0 and y = h.
Since the sample has a finite size, at large distances from the separation device the concentration remains zero at all times: Let us introduce the following dimensionless variables and parameters: Then the mass transport equation ( 15) rewrites as and the initial condition becomes where the upper-case symbols represent the normalized parameters corresponding to lowercase ones.Boundary conditions ( 16), (17) retain their form in the normalized terms and therefore are not repeated here.Note that the species is theoretically "focused" only at large times, i.e., when τ → ∞, which corresponds to steady state conditions, though the plots shown in Fig. 3(e) establish that this is met asymptotically within reasonable times.In this case the initial position of the sample injection is not important and the resulting concentration distribution depends only on five dimensionless parameters: the charge of the species q = z|e|, the normalized channel length L, the normalized initial sample width W s at the point of injection, the Peclet number of the carrier flow, Pe, and the normalized potential distribution along the channel Φ(X) (compare with Fig. 2).In (e) are shown the variations of the maximum and half-width of the sample distribution at the centre of the channel or near its floor.
3 Simulation results

Time-evolution of a single species containing sample
Figures 3-6 show normalized concentration distributions for a sample containing a single charged molecule simulated for a linear ion-conducting layer with a wedge-slope (unless otherwise stated) k = 1 (see Eq. ( 7)) and the following dimensioned parameter values:  The plots in Figs.3-6 show that at short times the concentration distribution spreads diffusionally and moves downstream (i.e., to the right in all figures) being almost unaffected by the electrophoretic forces.This is evidenced in Fig. 3(e) where one observes that the sample concentration moves initially linearly with time.At long times (Fig. 3(e)) the sample approaches the point where the electrophoretic velocity becomes roughly equal to the average hydrodynamic velocity so that the sample displacement is gradually slowed down and eventually stopped.Note that since the solution hydrodynamic flow varies across y, the point at which the analyte concentration distribution in the center of the channel (viz., at y = h/2, i.e., Y = 0.5) is stopped and that at which it stops near the channel wall (viz., at y = 0 or y = h, i.e., Y = 0 or Y = 1) must differ.This is evidenced in Fig. 3(e) by comparing the X max positions at Y = 0.5 and Y = 0 or by examining the variations of the sample half-width variations (W s ) with time again at Y = 0.5 and Y = 0.These latest plots also establish that while the sample starts to enlarge its half-width due to diffusion over the range of time in which it progresses almost linearly along the channel length, as soon as the electrophoretic contribution becomes sufficiently important the sample is progressively focused, i.e. its concentration distribution becomes narrower than when it was traveling at intermediate times.Thus, not only the sample is ultimately poised at a given position within the channel, but it is also focused so as to eliminate part of the effect of diffusional spreading incurred during its initial travel within the channel.
More precisely this focusing property of the present device (viz., the width and shape of the focused concentration distribution obtained at steady state) depends on the intensities of flow and the electric field.This is demonstrated in Fig. 4 for three different combinations of the Peclet number (Pe = 10, 50 and 100) and the voltage applied between the two electrodes (U = 10 V, 50 V, 100 V, which correspond to dimensionless values Ψ = F U /RT = 389.2,1946.1 and 3892.2,respectively).One can see from Fig. 4 that although the sample focusing position is the same for all three situations considered the resulting width of the sample is the smaller the steeper the potential gradient, which obviously results in better focusing.In the case of Fig. 4 this had to be achieved by simultaneous increase in both the hydrodynamic flow rate and the electric potential applied along the channel.On the other hand, the same effect can be achieved by changing the geometry of the ion-conducting membrane, viz., by varying the electric field profile without changing the applied potential difference (see Fig. 2).This is exemplified in Fig. 5 for three wedge-shaped membranes with different slopes.It follows that the width of the focused sample depends on the diffusivity of the analyte as well as on the electric field strength.In order to achieve better focusing, the latter can be increased not only via increasing the potential difference applied along the separation channel but also by changing the shape of the ion-conducting membrane.This may allow achieving the same focusing quality using lower voltages than those used for illustration in Figs.4-5 if this is desired experimentally.On the other hand, it should be noted that employing higher voltages and higher solution flow rates allow the reduction of focusing times which may be advantageous in experimental practice.The plots in Figs.3-5 evidence the great versatility of the method since they demonstrate that one can adjust the electrical and hydrodynamic conditions to obtain a desired result.

Separation of two differently charged species
If the sample consists of molecules carrying different electric charges and/or different diffusion coefficients each component will be focused at different positions along the channel.Figure 6 illustrates the limiting distribution of two species with z = 1 and z = 2 (and an identical diffusion coefficient for simplicity) for the dimensionless parameters: Pe = 100, Ψ = 3892.2,k = 5.The initial samples of the species were superimposed and had the same dimensionless width W s = 1 as well as equal concentrations.However, as evidenced in Fig. 6, when focused, the two samples concentration distributions have different amplitudes and widths since they traveled different lengths and are ultimately trapped in areas with different electric field strength (see, e.g., Fig. 2 for k = 5).

Computational details
A computer program for the numerical solution of the model ( 16)- (20) was developed by the authors.The time-dependent PDE (19) was solved using the Alternating Direction Implicit finite difference method [71] which affords extremely rapid and accurate simulations.The computer program features a user-friendly interface allowing easy parameter input and graphical presentation of the computed concentration distributions as well as dynamic sample parameters (position and width).The program can be obtained from the authors upon request.

Conclusions
This theoretical work based on the concept recently introduced by some of us demonstrates that the focusing of selected analytes in microfluidic separation channels thanks to the opposition between a constant transport by the hydrodynamic flow and a gradient of the electrophoretic driving force provides extremely interesting opportunities for chemical and biochemical separations.Beyond its high experimental simplicity, the main advantage of this concept is that the key electrical gradient is directly provided by that of the thickness of an ion-conducting membrane so that the system may be optimized for specific analytes of interest by very simple operations.As proven by the theory presented in this work and the corresponding simulations the method is then highly versatile in the sense that these theoretical analyses showed that specific shapes of the ion-conduction membrane thickness profile may ensure the simultaneous focusing of several selected analytes at the same time, a property which may be extremely useful when searching for example for a specific protein distribution pattern in biological fluids.

Fig. 1 .
Fig. 1.(a) Schematic description of the focusing microchannel system geometry; (b) equivalent electrical circuit of a vertical slice of thickness dx of the microchannel shown in (a).

Fig. 3 .
Fig. 3. Electrophoretic focusing of a single species for Pe = 30 and Ψ = F U /RT = 1194.66(real parameter values: U = 30 V, vavg = 0.03 cm s −1 ).Simulated normalized concentrations shown in 3D after different durations: (a) τ = 0; (b) τ = 0.1; (c) τ = 1.5 (corresponds to maximum sample width); (d) τ = 7.5.See text for other parameter values.In (e) are shown the variations of the maximum and half-width of the sample distribution at the centre of the channel or near its floor.

Fig. 6 .
Fig. 6.(a) 3D-representation of the simultaneous focusing of 2 species with z = 1 (peak on the right) and z = 2 (peak on the left), respectively; (b,c) variations of the positions of the two sample maxima in the center ((b) Y = 0.5) and near the bottom ((c) Y = 0) of the channel as a function of time: z = 1 and z = 2.The shaded areas around each curve filled with the same color represent the widths of each sample during the focusing operation.See text for other parameter values.