This Article Abstract Full Text (PDF) Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Albrecht, D. R. Articles by Bhatia, S. N. Search for Related Content PubMed PubMed Citation Articles by Albrecht, D. R. Articles by Bhatia, S. N.

Geometric and Material Determinants of Patterning Efficiency by Dielectrophoresis

Dirk R. Albrecht ^*, Robert L. Sah ^*  and Sangeeta N. Bhatia ^* 

Departments of ^* Bioengineering, Medicine, and Whitaker Institute of Biomedical Engineering, University of California-San Diego, La Jolla, California 92093

Correspondence: Address reprint requests to Sangeeta N. Bhatia, MD, PhD, Microscale Tissue Engineering Laboratory, Dept. of Bioengineering, PFBH 217, University of California-San Diego, 9500 Gilman Dr., La Jolla, CA 92093-0412. Tel.: 858-822-3142; Fax: 858-822-4203; Email: sbhatia{at}ucsd.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION DIELECTROPHORETIC FORCE MODEL PATTERNING EFFICIENCY MODEL MODEL VALIDATION DISCUSSION CONCLUSIONS REFERENCES

 
Dielectrophoretic (DEP) forces have been used extensively to manipulate, separate, and localize biological cells and bioparticles via high-gradient electric fields. However, minimization of DEP exposure time is desirable, because of possible untoward effects on cell behavior. Toward this goal, this article investigates the geometric and material determinants of particle patterning kinetics and efficiency. In particular, the time required to achieve a steady-state pattern is theoretically modeled and experimentally validated for a planar, interdigitated bar electrode array energized in a standing-wave configuration. This measure of patterning efficiency is calculated from an improved Fourier series solution of DEP force, in which realistic boundary conditions and a finite chamber height are imposed to reflect typical microfluidic applications. The chamber height, electrode spacing, and fluid viscosity and conductivity are parameters that profoundly affect patterning efficiency, and optimization can reduce electric field exposure by orders of magnitude. Modeling strategies are generalizable to arbitrary electrode design as well as to conditions where DEP force may not act alone to cause particle motion. This improved understanding of DEP patterning kinetics provides a framework for new advances in the development of DEP-based biological devices and assays with minimal perturbation of cell behavior.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION DIELECTROPHORETIC FORCE MODEL PATTERNING EFFICIENCY MODEL MODEL VALIDATION DISCUSSION CONCLUSIONS REFERENCES

 
The physical manipulation of biological particles (from DNA to viruses to cells) is a vital component of miniaturized biotechnological platforms such as "lab-on-a-chip" devices and arrays for high-throughput assays (Huang et al., 2001; Kapur et al., 1999; Ozkan et al., 2003; Voldman et al., 2002). To date, cellular arrays have primarily been generated via selective surface modification or micropatterning techniques (Bhatia and Chen, 1999; Folch and Toner, 2000). However, these passive techniques rely on cell adhesion processes that occur over relatively long timescales (hours), are limited to planar surfaces, and cannot be generalized to nonadherent cell types. In the future, "active" patterning techniques that are independent of cell adhesion would be advantageous to decrease patterning time and extend capabilities to the manipulation of nonadherent cell types or cells that alter their behavior when spread against a rigid surface.

Active patterning can arise through the application of a variety of physical forces, including mechanical, fluidic, optical, acoustic, and electromagnetic (Ashkin, 1997; Iwasaka et al., 2001; Jager et al., 2000; Matsue et al., 1997; Ozkan et al., 2003; Wu, 1991). One particularly advantageous strategy utilizes dielectrophoresis (DEP), the translational motion of polarizable matter (neutral or charged) within a spatially nonuniform electric field (Hughes, 2003; Jones, 1995; Pohl, 1978). The DEP force moves particles toward regions of high field intensity (positive DEP or +DEP) or low field intensity (negative DEP or –DEP), depending on electrical properties of the particle and suspending medium. By varying electrode shape and excitation, DEP forces have been used for particle trapping, arraying, levitation, translation, fractionation, filtration, orientation, and characterization based on dielectric properties (Arnold and Zimmermann, 1988; Frenea et al., 2003; Gascoyne and Vykoukal, 2002; Matsue et al., 1997; Voldman et al., 2001). However, the general use of DEP forces for patterning cells to fixed locations has not been extensively pursued since its initial proposal (Matsue et al., 1997), in part because DEP forces are inherently transient and disappear when the field is removed. Because living cells appear to tolerate high electric fields for relatively short times (Archer et al., 1999; Docoslis et al., 1999; Glasser and Fuhr, 1998), long-term biological experimentation requires a method to stabilize cell position in the absence of DEP forces. Recent strategies for restraining cell migration after active DEP localization include the incorporation of cell-adhesive proteins (Gray et al., 2003) or etched microwells (Frenea et al., 2003). Alternatively, encapsulation of patterned cells within hydrogel biomaterials (Elisseeff et al., 2000) would further enable incorporation of cells that are nonadherent or require three-dimensional (3-D) microenvironments to maintain tissue-specific functions (Abbott, 2003; Benya and Shaffer, 1982; Cukierman et al., 2002; Schmeichel and Bissell, 2003). To achieve this goal of positioning cells by DEP and immobilization via hydrogel entrapment, it is desirable to identify experimental parameters that: 1), maximize DEP patterning efficiency, 2), support cell viability without perturbing cell function, and 3), preserve the chemistry of hydrogel formation. These requirements suggest that a model-guided, rather than empirical, determination of experimental parameters is advantageous.

Existing analytical and numerical models of DEP force illustrate spatial variations within a device but historically have been limited to specific geometries or approximate boundary conditions (Chang et al., 2003; Clague and Wheeler, 2001; Green et al., 2002; Heida et al., 2001; Masuda et al., 1987; Morgan et al., 2001; Schnelle et al., 1993; Wang et al., 1993, 1996). Practical applications of these models include predicting levitation height, particle velocity, and selected particle trajectories (Markarian et al., 2003; Morgan et al., 2001; Qiu et al., 2002), although none have integrated these models into a single parameter that describes the overall time for DEP motion or "patterning" to occur.

The purpose of the work, described herein, was to develop and validate a predictive model of DEP patterning efficiency. Toward this goal, an expression was derived for the time at which model particles localize to a parallel line pattern within a thin, rectangular fluid volume. The patterning chamber was bounded by a planar interdigitated bar electrode array that is common to many DEP devices. Patterning kinetics were based on a new analytical solution for electric potential that allowed continuous variation of electrode and chamber geometry. In contrast to previous models, this solution specified a finite, variable chamber height and an improved set of boundary conditions, guided by numerical model data, which resulted in superior accuracy. The variation in particle kinetics was compared for changes in chamber geometry (electrode spacing and width, chamber height), material parameters (particularly viscosity and medium conductivity), patterning configuration (+DEP and –DEP), and particle buoyancy. The model closely predicted experimental kinetics of model particles, even for the unexpectedly slow patterning at very small chamber height. Therefore, the model serves as a useful tool for optimization of materials, geometries, and field conditions for efficient patterning, particularly for the patterning of living cells where viability is crucial.

	DIELECTROPHORETIC FORCE MODEL

TOP ABSTRACT INTRODUCTION DIELECTROPHORETIC FORCE MODEL PATTERNING EFFICIENCY MODEL MODEL VALIDATION DISCUSSION CONCLUSIONS REFERENCES

 
The parallel, planar, interdigitated electrode array is depicted in Fig. 1 A. When a potential is applied across alternating electrodes, this array establishes a nonuniform electric field dependent on dimensional and material properties. Particles suspended in a fluid experience a spatially varying DEP force that governs particle kinetics. Therefore, we begin with a model of DEP force for any chamber geometry to investigate the geometric determinants of patterning kinetics.

View larger version (19K): [in this window] [in a new window]   FIGURE 1  (A) Schematic of the DEP patterning chamber with interdigitating bar electrodes (shaded) at bottom surface. Top and bottom surfaces (shaded) are nonconducting glass. Geometric variables include chamber height (h), electrode spacing (d), and electrode width (w). (B) The 2-D problem space for the numerical model, depicting volume dimensions and boundary conditions in the x-z plane. The rectangular liquid chamber containing water (w = 80 0; w = 10–4 S/m) is bounded at the top and bottom by glass (gl = 4.5 0; gl = 10–12 S/m; hgl = 1 mm). Potential is specified at the electrode surfaces (solid bars) only. (C) The simplified problem space used for the analytical solution. The boundary condition at the bottom plane is of Neumann type, where function f(x') is determined from the FEM model (see text). (D) Neumann BC at z' = 0. Normalized values from 13 FEM models spanning the geometry range of interest (points) are fit to a third-order polynomial f(x') for analytical solution BC (shaded line). Electrode (solid bar) edge is at x'/w' = 0.5.

The quasielectrostatic form of Maxwell's equations is appropriate for the low currents and high frequencies typical to DEP. For a homogeneous linear dielectric with conductivity and permittivity , the electric potential at a sinusoidal steady state is determined by:

(1)

where is oscillation frequency, and the tilde () indicates a complex variable (phasor). For a homogeneous medium, this reduces further to Laplace's equation for the real portion of the electric potential,

(2)

The time-averaged dipole approximation to the DEP force is given by (Jones, 1995):

(3)

where E_rms is the root-mean-square electric field vector, _m is the medium permittivity, R is the particle radius, and is the Clausius-Mossotti (CM) factor,

(4)

dependent on complex permittivities of the particle (p) and medium (m). Equation 3 is valid except where the spatial field nonuniformity is very large, such as at the electrode edges, and higher-order terms are necessary (Voldman et al., 2001).

Dimensionless equations
The electric field depends on the applied electrode voltage (V_rms), electrode width (w) and spacing (d), and chamber height (h). To reduce the number of model variables, dimensionless variables, denoted by prime ('), are defined as: potential displacement vector electrode width chamber height gradient operator and electric field The time-averaged DEP force becomes:

(5)

where the term in braces represents a characteristic force constant, . Dimensionless DEP force,

(6)

is obtained by solving Eq. 1 or 2 for potential throughout the solution space. As indicated by the subscript notation, this solution needs to be computed only once for a particular nondimensional electrode width w' and chamber height h'. DEP force is then scaled by the characteristic force constant, according to the applied voltage, electrode spacing, particle radius, and frequency-dependent dielectric properties.

Analytical solution for electric potential
DEP forces throughout the solution space can be predicted readily from the electric potential. Analytical solutions are advantageous because: 1), geometric parameters can be varied continuously and independently to provide physical insight; 2), analytical expressions can be obtained directly for electric field and DEP force; such that 3), errors are not introduced by discretizing the volume and numerically estimating gradients; 4), accurate field and force determinations can be made at volume boundaries; and 5), the approach is computationally more efficient than numerical techniques.

Even for a simple geometry such as the interdigitated planar electrode array, no exact analytical expression exists for the electric potential. A boundary condition (BC) of mixed type (Dirichlet and Neumann) best represents the electrode substrate-fluid interface plane at z' = 0, hereafter designated as the "electrode substrate." However, mixed boundary value problems are difficult to solve analytically and solution methods are limited and complex, involving Green's functions or integral transforms (Sneddon, 1966). In contrast, analytical solutions with a single type of BC across each distinct boundary are significantly simpler to obtain. For example, linear approximations to the electric potential at the electrode substrate between interdigitated electrodes were assumed in previous Fourier series (Masuda et al., 1987; Morgan et al., 2001) and Green's function (Clague and Wheeler, 2001; Wang et al., 1996) solutions. These approximate models, which we term "linear BC" solutions, have limited accuracy (Green et al., 2002), prompting the development of complex analytical solutions with improved BCs (Chang et al., 2003; Wang et al., 1996). However, improved solutions to date do not include the effects of a finite top boundary, and are therefore unsuitable for this geometric analysis.

Numerical methods, in contrast, are able to specify the correct mixed BC without difficulty. In this manner, electric field solutions were obtained for several electrode geometries, using point charge, charge density, finite difference, integral equation, and finite element methods (Green et al., 2002; Heida et al., 2001; Schnelle et al., 1993; Wang et al., 1993). However, each geometry requires a separate computation, and general design principles are difficult to extract.

In this article, we employ a hybrid numerical/analytical solution method to obtain Fourier series expressions of both electric potential and DEP force, in which geometric and material properties can be continuously varied. A significant advantage of these solutions is that they do not assume an unbounded half-space and are therefore appropriate for finite upper boundaries. Further, accuracy is greatly improved over linear approximations with only a moderate increase in model complexity. Although this hybrid approach remains an approximation, we demonstrate its suitability in this work and suggest that its simplicity may find further application in other areas.

Modeling approach
The model utilizes a combination of numerical and analytical methods. First, numerical solutions that allow the correct mixed BC definition at the electrode substrate were obtained for a range of chamber geometries. Then, a geometry-independent curve fit to either the electric potential or the electric field provides a single BC type at the entire boundary plane, such that a simple analytical expression of the electric potential can be derived. We found that the latter option, specifying the electric field normal to the wall (a Neumann condition) provides better accuracy, especially for narrow electrodes.

Boundary conditions
Because the electrodes are typically long relative to their width, the problem can be considered two-dimensional (Fig. 1 B). At volume boundaries and symmetry planes, the Neumann condition reflects the absence of current conduction across these planes, as described in detail elsewhere (Green et al., 2002). This approximate BC is valid at all frequencies for the potential in the fluid (water; _w = 80₀; _w = 10^–4 S/m) bounded by an insulating wall with reduced conductivity and permittivity (glass; _gl = 4.5₀; _gl = 10^–12 S/m). The electrodes are considered infinitely thin, with a constant potential equal to the applied voltage (' = –1 or +1 at z' = 0).

For the numerical model, upper and lower glass walls were included (Fig. 1 B). However, because the conductivity of glass is negligible, the solution space for the analytical model is simplified to contain only the liquid volume (Fig. 1 C). By the same argument as above, current conduction normal to the surface is zero at the liquid-glass boundaries, i.e., . The numerical model, described below, determines the unknown Neumann BC at the electrode surface.

Numerical finite element model
The commercially available CFD-ACE+ finite element software suite (CFDRC, Huntsville, AL) was used to solve Eq. 1 for the sinusoidal steady-state electric potential. This package utilized a finite volume method technique and a structured mesh with nodes concentrated near the electrode edges.

Fifteen separate solutions of potential and electric field at the electrode wall (z' = 0) were obtained for various chamber geometries: h' = 0.33 – 1.0, and w' = 0.1 – 0.5. Larger heights were not included because of negligible effects on electrode wall BCs above h' > 1. Numerical values for electric field normal to the electrode plane for all geometries converge to a single curve upon normalization in the x- and z-directions (Fig. 1 D). A third-order polynomial fit,

(7)

approximates the Neumann BC across the entire lower boundary for any electrode width and chamber height.

Improved analytical solution
The Fourier series solution for electric potential (Eq. 2) is defined by coefficients:

(8)

where at z' = 0. For these integral coefficients to have an analytical solution, is limited to linear combinations of polynomial, exponential, and sinusoidal functions. These functions may be discontinuous. Power, logarithmic, and inverse functions, for example, would not lead to analytical Fourier coefficients.

Using Eq. 7, Fourier coefficients depend on electrode width, w':

(9)

The electric potential is then described by an infinite series:

(10)

Here, is a scaling function such that potential at the electrode center:

(11)

In this article, the first 250 terms were computed for each series.

Solving Eq. 10, Fig. 2 A illustrates the potential and electric field boundary conditions at the electrode substrate (z' = 0), for h' = 2/3, w' = 0.2. The "improved BC" analytical solution closely matches the numerical model at the boundary, in contrast to the linear BC solution described elsewhere (Morgan et al., 2001).

View larger version (36K): [in this window] [in a new window]   FIGURE 2  Comparison of analytical solutions using the improved boundary condition (left, "Improved BC") or a linear approximation between electrodes (right, "Linear BC") and the numerical solution for w' = 0.2. (A) Electric potential and field components at the electrode boundary, z' = 0, demonstrate a closer match between the numerical model and the improved BC analytical model compared to the linear BC model. (B) Comparison between electric field magnitude from numerical solution and analytical solutions |E'|, where contour shading indicates relative error: The improved BC solution (left) deviates from the numerical solution by <3% for a majority of the solution space (dark red/orange), whereas the linear solution (right) deviates by >10% throughout most of the volume (white). Similar accuracy is achieved with models of different geometry.

By employing advantages of both numerical and analytical methods, the improved solution method allows quick computation of the electric potential, , with greater accuracy than previous analytical solutions. All variables are continuous, including particle position (x') and chamber geometries (w' and h'), an important advantage over discrete numerical solutions.

DEP force calculation
Nondimensional DEP force (Eq. 6) is conveniently recast in indicial notation:

(12)

where i,j represent vector indices. Because electric potential is given by an analytical expression (Eq. 10), gradients of potential are also analytical expressions that are readily computed using symbolic mathematics software, such as Mathematica (Wolfram Research, Champaign, IL). As a result, only a single series calculation is necessary for determining the DEP force vector at a particular position within the solution space.

	PATTERNING EFFICIENCY MODEL

TOP ABSTRACT INTRODUCTION DIELECTROPHORETIC FORCE MODEL PATTERNING EFFICIENCY MODEL MODEL VALIDATION DISCUSSION CONCLUSIONS REFERENCES

 
The velocity and position of a model particle can be tracked during simulated patterning using the continuous analytical solution of DEP force exerted upon it. The time required for the slowest particle to pattern is a measure of DEP patterning efficiency.

Theory and model
Force balance
Particle motion arises from a balance of applied DEP force, viscous drag, gravity, inertia, surface friction/adhesion, and interparticle attraction, and is further influenced by convection of the suspending fluid due to pressure, electrothermal, or electroosmotic flow. In this analysis, we assume negligible inertia due to low Reynolds number flow (typically, Re < 10^–5 in microfluidics), no convection, and no surface forces. Thus, a balance between DEP force (Eq. 5) and drag force defines kinetics of neutrally buoyant particles:

(13)

(14)

where µ is fluid viscosity and K(x') is a scaling factor equal to unity for motion of a rigid particle far from a wall (Stokes drag) but increases near a wall to reflect increased drag force (Goldman et al., 1967; Keh and Chen, 2001; Oseen, 1927).

The particle velocity field,

(15)

is nondimensionalized by defining a time constant for DEP,

(16)

and a dimensionless time, such that:

(17)

The position of any particle at any time, can be determined by solving this system of coupled differential equations with the initial condition Here, the initial position vector identifies which particle is tracked. A particle reaches its equilibrium patterned location, at a dimensionless time, that is dependent on its initial (and final) positions.

Patterning kinetics
Patterning kinetics for selected neutrally buoyant particles are demonstrated in Fig. 3. For comparison, the electrode plane is oriented on the bottom for the +DEP case and on top for the –DEP case, such that particles pattern downward to the lower plane in both illustrations. Path lines that track particle centroid over time curve away from the pattern location in the upper portion of the chamber and converge in the lower regions, for both +DEP and –DEP. Also, particle velocity is greatest near the electrodes, as indicated by increased separation between symbol points representing regular time intervals. To understand the evolution of pattern formation over time, consider the contours of dimensionless patterning time, in Fig. 3. At a particular time, all particles initially within the area where will have translated from their initial locations to the final pattern . As increases, this depleted space will increase until it encompasses the entire chamber area. By –DEP, the center column below the electrode patterns first, and expands over time at a nearly uniform rate. By +DEP, however, particles in a semicircular area surrounding the electrodes pattern initially, and this area grows at an exponentially decreasing rate. These distinct patterning kinetics occur because particles accelerate toward the field maxima by +DEP, whereas they slowly converge to the field minima by –DEP. Patterning time is shorter by +DEP than by –DEP for any initial particle position.

View larger version (37K): [in this window] [in a new window]   FIGURE 3  Patterning motion and kinetics by (A) –DEP and (B) +DEP. (Right) The DEP patterning motion varies with initial particle location ( gray open circles), as indicated by centroid pathlines. Symbol points are equally spaced in time, indicating higher particle velocity near the electrode (solid bar). Near the pattern location ( black open circles), motion is slow by –DEP but rapid by +DEP. (Left) Contours of dimensionless patterning time for various initial particle locations. This value is lowest near the pattern location and greatest at the opposite wall in between electrodes.

(18)

where initial and final positions are or and or for –DEP or +DEP, respectively, and is the dimensionless particle radius.

The predicted time to achieve complete patterning for any chamber geometry is:

(19)

Effects of geometry on patterning efficiency
Model results are presented below in two forms, based on nondimensionalization of geometric parameters to electrode spacing, d, or chamber height, h. The former case is convenient for choosing an optimal height given a uniform electrode pattern, whereas the latter formulation is ideal for varying pattern dimensions or when a particular height is required.

Constant electrode spacing
In many applications of DEP patterning, control over the dimensions and geometry of the resulting particle pattern is desired. For example, a parallel line pattern is specified by interdigitated electrode spacing, and variations in electrode width and chamber height could be made to accelerate patterning time.

To determine the effects of chamber geometry, the complete nondimensional patterning time, is computed for varying dimensionless chamber height and electrode width. For all geometries, patterning by +DEP (Fig. 4, B and D) is more rapid than by –DEP (Fig. 4, A and C). Increasing electrode width decreases patterning time for both +DEP and –DEP for all chamber heights, although this effect is more pronounced for thinner chambers.

View larger version (61K): [in this window] [in a new window]   FIGURE 4  Effects of geometry on patterning efficiency by –DEP (A, C) and +DEP (B, D) for constant electrode spacing. Dimensionless complete patterning time, T', varies with chamber height h' and electrode width w'. Model calculations are presented as contour plots (log contours at 1, 2, 5 gradations) above (A, B) and as families of curves below (C, D) for particular electrode widths. To explain the chamber height effect, patterning kinetics (E) and electric field lines (F) are shown for three different chamber heights (a < b < c), where b represents the optimal height with fastest patterning for w' = 0.4. Particle motion is tracked from x' = 0.5 to the field minimum at x' = 0. As chamber height increases beyond optimal (e.g., c), decreased electric field strength increases patterning time. However, at suboptimal chamber heights (a), patterning time increases dramatically due to confinement of the electric field nonuniformity near the electrode. Therefore, patterning is very slow in regions away from the electrodes (*), although it is very fast near the electrodes (**).

Constant chamber height
Many applications of DEP patterning require a specific chamber height, but allow variable electrode geometry as design parameters. Some examples include microfluidic devices that often contain uniform channel geometry, and the patterning of cells within a gel material whose thickness is defined by the chamber height (Albrecht et al., 2002). Other applications, such as particle depletion from the bulk of a microfluidic channel, would be less concerned with specific pattern dimensions than with chamber height. Finally, optimizing arbitrary particle patterns would require model efficiency information for a variety of relative heights (h') and widths (w').

The complete patterning time can be nondimensionalized in terms of chamber height h rather than electrode spacing d:

(20)

where dimensionless geometric variables are recast as electrode width, w* = w/h = w'/h' and electrode spacing, d* = d/h = h'^–1. Then,

(21)

Fig. 5 illustrates the improved analytical model solution for the nondimensional patterning time of a neutrally buoyant particle, for varying nondimensional electrode spacing (d*) and electrode width (w*). Patterning kinetics show similar trends with constant chamber height as with constant electrode spacing. Positive DEP patterning (Fig. 5, B and D) is more rapid than negative DEP (Fig. 5, A and C) for all geometries, and wider electrodes accelerate patterning for all chamber heights and for both +DEP and –DEP. Patterning time is also highly dependent on electrode spacing, increasing above and below optimal values. This variation is due to the balance between: 1), the direct dependence of on electrode spacing by (d*)⁴ (Eq. 20), and 2), the effects of concomitant inverse changes in relative chamber height, h' = (d*)^–1. The former effect dominates at electrode spacing d* > 1. However, as spacing is reduced below d* < 1, thereby increasing h' > 1, the latter effect dominates as patterning time increases exponentially with relative chamber height (Fig. 4). Optimal electrode spacing, which varies with electrode width, occurs when these dependencies are balanced.

View larger version (39K): [in this window] [in a new window]   FIGURE 5  Effects of geometry on patterning efficiency by –DEP (A, C) and +DEP (B, D) for constant chamber height. The complete patterning time, nondimensionalized to chamber height, , varies with electrode spacing d* = d/h and electrode width w* = w/h. Model calculations are presented as contour plots (log contours at 1, 2, 5 gradations) above (A, B) and as families of curves below (C, D) for particular electrode widths.

View larger version (40K): [in this window] [in a new window]   FIGURE 6  Optimal chamber geometry parameters for –DEP () and +DEP (•) determined by minimizing the complete patterning time T' or . (A) For a constant electrode spacing (d), the optimal chamber height (h' = h/d) decreases as electrode width (w' = w/d) increases. (B) Similarly, for a constant chamber height (h), the optimal electrode spacing (d* = d/h) varies with electrode width (w* = w/h). Shaded or hatched regions indicate parameters resulting in efficiency within 10% of optimal. Overall, patterning time decreases with wider electrodes (top panels). (C) An example of chamber height optimization for +DEP patterning of neutrally buoyant particles within a device containing three distinct electrode geometries: d = 4, 10, and 20 µm, w' = 0.5. Most rapid patterning for all three geometries occurs at chamber height h 7 µm (*).

A final consideration regarding the chamber height is the limitation imposed by the particle size, requiring for all cases h/2R > 1. Near this limiting case, an increase in drag force due to wall effects may occur, thus increasing the factor K in Eq. 14 and also the predicted patterning time. Remarkably, this effect is minimal for nearly all chamber geometries. For example, drag force on rigid spheres translating along a wall is increased <10% for a chamber height twice that of the particle diameters (h/2R = 2) and 80% for the limiting case (h/2R 1) (Keh and Chen, 2001). Thus, the model conclusions presented in this article are valid even for extremely thin chambers.

Gravitational effects
For particles that are not neutrally buoyant, gravitational force is included in the force balance:

(22)

where is the difference in density between the particle and the fluid, and g is the gravitational acceleration in the direction of unit vector . Then,

(23)

where the nondimensional gravitational term,

(24)

describes the relative influence of gravitational force to DEP force. Typical values of are of the order 0.01–1 for polymer beads and cells, patterned d = 100 µm apart at 1–10 V_rms. Therefore, DEP force is stronger than gravitational force for most conditions (hence the ability to levitate these particles). In contrast, metal particles may have 10–1000 under the same conditions.

To illustrate the influence of gravitational forces, dimensionless patterning time was computed for varying and chamber height in Fig. 7. Here, the chambers were oriented such that gravity would accelerate patterning. Thus, for particles denser than the fluid, the electrode array was oriented above for –DEP and below for +DEP patterning (as in Fig. 3, gravity acting downward). The greatest potential influence of gravity is illustrated by comparing neutral-bouyancy ( = 0) results with the asymptotic case, , representing maximum gravitation where all particles instantly settle to the bottom wall. This case also characterizes the situation where particles are allowed to settle completely in the absence of an electric field before patterning.

View larger version (22K): [in this window] [in a new window]   FIGURE 7  Effects of gravity on patterning efficiency by –DEP (A) and +DEP (B). Nondimensional complete patterning time is calculated for varying gravitational factor, (defined in text), and chamber height. Electrode width is w' = 0.2. Above, the relative patterning time compared to neutral buoyancy ( = 0) indicates that the potential for gravity to assist particle patterning is greater for +DEP than for –DEP.

To explain these phenomena, it is helpful to refer back to Fig. 3. The arcing pathlines indicate that downward particle motion occurs mainly in the upper half of the chamber, whereas particles converge laterally to the pattern in the lower half of the chamber. Because DEP force and particle velocity are greatest near the electrodes, downward motion is rapid for –DEP and slow for +DEP; conversely, the lateral convergence is slow for –DEP and rapid for +DEP. Because gravity accelerates the downward motion that is already rapid with –DEP but slow with +DEP, it has only a slight effect on patterning by –DEP but a profound effect by +DEP.

	MODEL VALIDATION

TOP ABSTRACT INTRODUCTION DIELECTROPHORETIC FORCE MODEL PATTERNING EFFICIENCY MODEL MODEL VALIDATION DISCUSSION CONCLUSIONS REFERENCES

 
The DEP patterning kinetics model was evaluated using polystyrene microspheres as model particles. Because the polystyrene beads in water experience negative DEP at high frequencies, the patterning chamber was oriented onto the microscope stage such that the upper chamber wall contained the electrode array (Fig. 8 A). The beads were then introduced into the chamber and allowed to settle randomly onto the bottom slide (shaded circles), such that patterning would occur in a single z-plane observable by conventional light microscopy. Upon electrode excitation, particles moved laterally toward field minima located opposite the electrodes (open circles).

View larger version (41K): [in this window] [in a new window]   FIGURE 8  (A) Schematic of video microscopy setup for validation of DEP patterning kinetics. Beads initially settle randomly (shaded circles). Application of AC voltage aligns beads (open circles) under the electrodes (solid and vertical shaded bars). Spacing between upper and lower glass slides (h) was adjusted between 34 and 110 µm. Electrode width varied from w = 20–50 µm, whereas electrode spacing was constant at d = 150 µm. Two bead diameters, 2R = 7.2 and 9.7 µm, were utilized. (B) Typical video images (h = 55 µm; w = 50 µm; 2R = 9.7 µm) of patterning beads at time t = 0, 1, and 6 s. Scale bar: 100 µm. (C) Typical patterning kinetics plotted as the distance away from the patterned location (mean ± SD, n 6) for h = 75 µm; w = 30 µm; 2R = 9.7 µm. (D–F) Summary of validation experiments measuring tpat,60 (or T'0.4), i.e., time for a bead to move from x = 60 µm (or x' = 0.4) to 1 µm away from the pattern location, for w' = 0.2. The improved analytical model (solid line) correctly predicts the optimal chamber height and the slower patterning kinetics at greater and smaller heights. Larger beads pattern faster (D; solid symbols and lower lines), and data points converge upon nondimensionalization (E). The improved BC solution results in closer prediction of patterning efficiency than the linear BC (dotted line). (F) Summary of validation experiments varying electrode width for height h = 55 µm and comparison to improved (solid line) and linear (dotted line) analytical models.

Interdigitated electrodes, spaced d = 150-µm apart and w = 20–50-µm wide, were microfabricated using conventional processes. Glass slides coated with transparent indium tin oxide (ITO) to a sheet resistance of 4–8 per square were purchased from Delta Technologies (Stillwater, MN). Slides were cleaned and photolithographically patterned using S1813 photoresist (Shipley, Marlborough, MA) and ultraviolet exposure through a custom emulsion mask. Next, exposed ITO was etched for 4 min with gentle agitation in a solution of 20% hydrochloric acid and 5% nitric acid, heated to 55°C. Etched slides were then sonicated in acetone to remove photoresist and cleaned. Adhesive copper tape provided electrical contact with the ITO film.

The patterning chamber was formed by sandwiching a thin silicone spacer between the interdigitated electrode array and a bare glass slide drilled with holes for fluidic entry and exit. The spacer thickness determines chamber height, h, and is composed of poly(dimethyl siloxane) (PDMS; Sylgard 184, Dow Corning, Midland, MI) to provide a watertight seal. The resulting fluid chambers were 25-mm long, 5–10-mm wide, and 34–110-µm high. To reduce particle adhesion, chamber walls were treated with poly(ethylene glycol)-disilane (molecular weight of 3400, Shearwater Polymers, Huntsville, AL).

The beads were introduced into the chamber via syringe and allowed to settle randomly onto the electrode array. Electrical excitation (500 kHz sine wave, 9.3 V_rms) generated by an Agilent 33120A signal generator (Agilent, Palo Alto, CA) caused lateral particle motion toward field minima located opposite the electrodes. An oscilloscope connected in parallel measured applied voltage. Particle motion was observed using a Diaphot 300 inverted microscope (Nikon, Melville, NY) with an analog video camera and processor (Dage-MTI VE-1000, Michigan City, IN). Video frames captured at 0.5–5 Hz were calibrated and segmented with Scion Image (Scion, Frederick, MD) to obtain particle centroid data. Particle kinetics were then analyzed with MATLAB (The MathWorks, Natick, MA). Data were nondimensionalized according to Eq. 16, using measured values for most parameters (Table 1) and Eq. 4 to estimate the Clausius-Mossotti factor.

Validation results
In Fig. 8 B, representative video images depict bead patterning behavior by –DEP. For actual validation experiments, fewer beads were introduced into the chamber such that nearby particles would not be encountered during patterning. Approximately six to eight beads per microscope field were selected for analysis based on: 1), their isolation from nearby beads, to prevent disturbances from particle-particle attraction or local fluid flow, and 2), their initial location far from the electrode pattern, to provide the most information on patterning kinetics. Thus, beads located originally at x > 60 µm (x' > 0.4) were accepted for tracking to obtain sufficient data points.

Typical patterning kinetics are described in Fig. 8 C. Because the electrodes were spaced d = 150-µm apart, no particle was ever more than d/2 = 75 µm away from the nearest field minimum (pattern location), located below the nearest electrode at x = 0. Particles that were initially far (40–75 µm) away accelerated toward this position, but near the pattern (<40 µm away), particle velocity gradually decreased. The "improved BC" analytical model predicts the sigmoidal position versus time curve closely, especially where particles are further away (larger x). Closer to the patterned location, the improved model predicts faster motion and shorter patterning time than experimental measurements. The "linear BC" analytical model predicts velocities about two times higher and patterning time about three times faster.

To demonstrate the effect of chamber geometry on patterning time, experiments were repeated with six chamber heights (h = 34–110 µm; h' = 0.23–0.73), two bead diameters (2R = 7.2 and 9.7 µm), and three electrode widths (w = 20–50 µm; w' = 0.13–0.33) (Fig. 8 D). The patterning time was recorded from an initial location at x = 60 µm (x' = 0.4) and designated t_pat,60 (or T'_0.4 in dimensionless form) to ensure that sufficient beads were tracked per measurement (n 6). Larger beads patterned faster for all geometries, and upon nondimensionalization, data points and model curves for different bead sizes converged to a single dimensionless patterning time for a given chamber height (Fig. 8 E). Notably, the model prediction of optimal chamber height (40 µm) is consistent with validation data. Heights >50 µm (h' > 0.33) resulted in exponentially longer patterning times, as exhibited in both the improved and linear analytic models. Below h 40 µm (h' < 0.27), patterning time also slowed, reflecting the situation predicted by the model wherein a diminished spatial field gradient slows patterning in the central regions between electrodes (Fig. 4, E and F; curve a). In fact, for h = 34 µm (h' = 0.23), particles >65-µm away from the pattern locations (x' > 0.43) were not observed to move at all during the experiment. This upturn in patterning time is not reflected in the linear BC analytical model (Fig. 8 E).

Electrode width effects are also consistent with the improved model (Fig. 8 F). Wider electrodes result in faster patterning, although no additional decrease in patterning time is seen for w' > 0.5. The linear BC model again underestimates patterning time, especially in the range used for validation, w' = 0.13–0.33.

	DISCUSSION

TOP ABSTRACT INTRODUCTION DIELECTROPHORETIC FORCE MODEL PATTERNING EFFICIENCY MODEL MODEL VALIDATION DISCUSSION CONCLUSIONS REFERENCES

 
Summary
This report examined the influence of various geometric and material parameters on the speed at which particles pattern by DEP forces. Faster patterning may be important to minimize damage to sensitive particles, especially living biological cells. An expression was derived for the time at which all particles travel to an equilibrium position, by +DEP or –DEP, within an interdigitated electrode array chamber of varying geometry. Particle kinetics were solved by a balance of forces, utilizing a new analytical expression of DEP force.

In contrast to previous solutions for electric potential, the solution presented here: 1), utilized a new method for specifying the BCs at the electrode array plane, and 2), specified a finite chamber height. First, numerical solutions for different geometries were combined into an analytical expression for the Neumann BC at the electrode substrate. The resulting solutions for electric potential, field, and DEP force have significantly greater accuracy than previous analytical solutions. In addition, they are continuous in position throughout the solution space, as well as in geometric parameters.

Variations in these geometric parameters were investigated to determine limitations of DEP patterning and for experiment optimization. The theoretical models were validated using inert microspheres, demonstrating good prediction of both particle kinetics and patterning time (within <50%) without the use of fitted or empirical scaling factors. Improved accuracy in the analytical DEP force solution translated to a twofold better prediction of patterning time compared to previous solutions (Fig. 8, C, E, and F).

Implications for cell patterning
The use of DEP forces for manipulation of living cells requires minimization of potentially cytotoxic physical, chemical, and electrical conditions. In particular, strong electric fields may have a significant physiological impact on a biological cell, due to current-induced Joule heating of the medium and direct field interactions inducing transmembrane potentials (Glasser and Fuhr, 1998). Although these effects can be reduced with low conductivity suspending media and high-frequency fields, respectively, even mild electric field conditions may cause subtle gene upregulation (Archer et al., 1999). Furthermore, electric field-induced damage is dependent on both duration and intensity of exposure (Glasser and Fuhr, 1998). Therefore, selecting the optimal patterning geometry (using Fig. 6) and materials (using Eq. 19) is important to minimize the exposure time (via maximal field gradient and DEP force) for a particular applied field intensity. Alternatively, the optimized chamber geometry may permit a decrease in the electric field strength while maintaining an acceptable patterning rate, again reducing possible cell damage. In general, patterning at the fastest overall velocity (i.e., at the greatest tolerated applied potential) appears to be advantageous, with the following rationale. Transmembrane voltage (V_tm) is proportional to applied voltage (V) but patterning time varies by V^–2, such that a linear measure of total exposure, , is minimized at greater applied voltage. However, biological systems rarely display linear responses to stimuli and this analysis offers only a general guideline. Rapid patterning further diminishes the exposure of cells to DEP buffers that often lack standard tissue culture medium components (ions, growth factors, proteins, and serum), and also reduces the relative influence of additional forces (bulk convection, Brownian motion) that may compromise pattern quality.

The parameters that lead to rapid patterning are summarized in Table 2. From Eq. 19, low viscosity and high permittivity of the suspending medium are desired. Although standard cell culture media meet these requirements, biomaterials for cell encapsulation (e.g., hydrogels) may have significantly greater viscosity and thereby hinder patterning speed. Patterning of polystyrene beads via –DEP and cells by +DEP is accelerated by low conductivity media that also minimize Joule heating. Reducing the characteristic feature length (electrode spacing) greatly speeds patterning, because field strength is increased and particles have less distance to travel, as do larger particles because DEP force is proportional to volume. Increasing the applied voltage accelerates patterning, but also raises potentially detrimental electric field effects. For mammalian cells, a limited range of electric field frequency is ideal, 1–10 MHz, that balances the maximal CM factor magnitude (Eq. 4) with minimal induced transmembrane potential at higher frequencies (Gray et al., 2003).

Because geometric variables influence patterning time so greatly, the principal advantage of the efficiency model is optimization of these parameters (Fig. 6). When a particular pattern is desired (such as cells located a uniform distance apart), an optimal height is determined using Eq. 19 and Figs. 4 and 6 A. However, in microfluidic devices, chamber (or channel) height is typically constant, and optimal electrode spacing is determined using Eq. 21 and Figs. 5 and 6 B. Optimal height prediction from Fig. 6 (h' = 0.28–0.40 for –DEP, w' = 0.2) was consistent with validation experiments, as the h' = 0.30 chamber showed the fastest patterning. It must be noted that the optimization figures do not consider particle size, although the model is valid for chamber heights greater than twice the particle diameter ("Chamber optimization using model"). When optimal parameters cannot be chosen exactly due to other system constraints, we found experimentally that it is advantageous to select a larger rather than smaller height. This is because suboptimal chamber heights show increased variation in particle velocity during patterning (curves a versus c in Fig. 4 E), resulting in a zone of poor patterning control with greater susceptibility to other forces (e.g., bulk fluid flow).

When several geometries are incorporated onto a single device, overall optimization is easily achieved by plotting dimensional patterning time and variables for each geometry. In "Gravitational effects", we chose as an example a set of three geometries previously used for analysis of +DEP-based particle aggregation, and showed that efficiency could be increased over 200-fold for neutrally buoyant particles by a simple reduction in chamber height from 30 µm to 7 µm (Fig. 6 C). The actual experiments (Markarian et al., 2003) reported the behavior of nonbuoyant particles ( = 0.01–0.16), and showed unsuccessful patterning as feature size decreases, despite the theoretical increase in DEP force as . We attribute these results in part to a chamber that was too tall, as when the reported experimental parameters were applied to the model (including gravity, using Eq. 19 and Fig. 7 B), efficiency improved 15-fold for all array geometries with the same fourfold height reduction. Thus, optimal chamber height can extend the reported lower limit on electrode dimension, although further limitations due to particle or electrical effects may still exist, as discussed below ("Model assumptions and additional forces").

Gravitational forces can evidently accelerate DEP patterning, although the improvements are more pronounced for +DEP than –DEP (Fig. 7). By –DEP, patterning time can be reduced by 50% with gravitational settling. In contrast, +DEP patterning time is significantly decreased because gravitational settling brings particles toward regions of high field nonuniformity where DEP force is strong. In effect, large chamber heights become feasible for +DEP patterning of particles with significantly different mass density. Also, electric field exposure can be minimized by allowing particles to settle in the absence of the field before DEP patterning. However, biological cells, biomolecules, and biopolymers typically have similar density to that of the suspending aqueous solution (1.0–1.1 g/cm³), such that << 1, settling is relatively slow, and chamber height should be carefully selected.

Using these tools, we have utilized DEP patterning in preliminary biological experiments with a variety of media and hydrogel formulations, several cell types (primary and cultured mammalian fibroblasts, chondrocytes, and hepatocytes), and a number of geometries (parallel lines, planar arrays, arbitrary shapes) (Albrecht et al., 2002), and demonstrated a reduction in patterning time consistent with the model (data not shown).

Model assumptions and additional forces
The patterning efficiency model assumes that particle motion is dominated by DEP forces. Several reports have demonstrated that electric field-derived fluid flow and interparticle effects may contribute significantly to particle motion. To aid in the use of the models for predicting patterning efficiency, the conditions where DEP is the dominating force are briefly discussed below.

Electrical effects
The use of high electric fields in DEP devices with microscale dimensions may give rise to significant temperature gradients, because power generation () is concentrated in a small volume. Electrohydrodynamic forces arise from the spatial variation in fluid permittivity, conductivity, density, and viscosity due to local temperature changes from Joule heating. However, fluid flow by natural convection is generally negligible for microelectrode devices and aqueous solutions (Ramos et al., 1998).

The electrical force on a fluid is composed of Coulomb and dielectric forces arising from gradients in conductivity and permittivity, respectively. Coulomb forces dominate at low frequencies and dielectric forces dominate at high frequencies. For aqueous media, the crossover frequency is approximately (Ramos et al., 1998), where is the charge relaxation time of the liquid. Electrical forces are proportional to medium conductivity and are reduced when dielectric forces dominate, i.e., at frequencies greater than (Ramos et al., 1998). The crossover frequency ranges from 100 kHz for water (_m 0.1 mS/m) to 1 GHz for standard cell culture medium (_m 1.5 S/m).

Electroosmotic forces may also contribute to fluid flow, due to the interaction of induced charges appearing at the electrodes and the electric field. Fluid velocity is proportional to surface charge density. The time to establish this charge is on the