Microfluid Nanofluid (2008) 4:273–285 DOI 10.1007/s10404-007-0169-0

RESEARCH PAPER

Characterization of transport in microfluidic gradient generators Bryan R. Gorman Æ John P. Wikswo

Received: 14 December 2006 / Accepted: 1 March 2007 / Published online: 30 May 2007
Springer-Verlag 2007

Abstract We present a two-dimensional model that describes the concentration profile of a class of previously reported microfluidic devices which are of particular interest in cellular taxis research. The devices generate stable concentration gradients by mixing and dividing two or more external inputs into a large number of discrete streams. This study focuses specifically on modeling the confluence of the discrete streams in a long chamber. We derive a closed-form solution for gradient generators with any arbitrary number of sampling streams. By relating the physical dimensions to the Pe´clet number, we create a model independent of flow rate and therefore dependent only on the specific nature of the boundary condition provided by the upstream network. As a result, the modeling method we propose may help evaluate the effectiveness of competing gradient generation schemes. Finally, our analytical work introduces a framework for developing simple design rules of interest to experimentalists working with these devices.

B. R. Gorman
J. P. Wikswo Vanderbilt Institute for Integrative Biosystems Research and Education, Vanderbilt University, Nashville, TN 37235, USA B. R. Gorman (&)
J. P. Wikswo Department of Biomedical Engineering, Vanderbilt University, Nashville, TN 37235, USA e-mail: [email protected] J. P. Wikswo Department of Molecular Physiology and Biophysics, Vanderbilt University, Nashville, TN 37235, USA J. P. Wikswo Department of Physics and Astronomy, Vanderbilt University, Nashville, TN 37235, USA

Keywords Microfluidics
Gradient generation
Chemotaxis
Convection–diffusion mass transport
Pe´clet number

1 Introduction Microfluidics has recently emerged as a powerful technology enabling chemical and biological assays that manipulate the chemical environment at extremely small spatial scales, occur rapidly, and use minimal amounts of reagent. Among the unique physical properties of microscale fluids is the laminar nature of flow (Squires and Quake 2005). Laminarity requires that mixing occur primarily by molecular diffusion, as opposed to the convective turbulence that dominates mixing in macroscale fluids. As a consequence, mixing in microchannels is inherently inefficient and has been the focus of a vast body of mass transport research (Nguyen and Wu 2005). Numerous micromixers have been devised implementing active and passive mixing strategies. Active mixing, which requires the addition of external energy, is generally more efficient, while passive mixing can be implemented on disposable chips. However, a divergent path of microfluidic mass transport research has recently emerged, with the goal not to mix streams fully, but instead to maintain a heterogenous concentration profile of desired shape and magnitude at a given location on-chip. These ‘‘gradient generators,’’ originally developed by Jeon et al. (2000), utilize the uniform mixing properties of laminar flows to pattern solution and surface properties at the microscale. The fabrication and usage of these devices are described in detail by Rhoads et al. (2005). A typical gradient generator is illustrated in Fig. 1a. All gradient generators are composed primarily of two

123

274

(a)

C* = 0

(b)

C* = 1

Concentration Inlets Gradient Generation Network

H z

y

0 x W

(c) Sampling Streams

Gradient Chamber

Normalized Concentration, C*

Fig. 1 a A pyramidal gradient generation scheme, after Jeon et al. (2000). b The geometry and coordinate system of the model. The location y = 0 is at the entrance to the gradient chamber. The flow is in the positive y direction. c Examples of staircase input boundary conditions (solid) and associated impulse response functions (dashed), after Sager et al. (2006)

Microfluid Nanofluid (2008) 4:273–285

1

0.8

0.6

0.4

0.2

0.2 Waste Outlet

components: a network that divides a small number of concentration inlets from a source external to the chip into a larger number of discrete ‘‘sampling’’ streams, and a chamber that draws each sampling stream in confluence to achieve a continuous concentration profile of desired shape inside the chamber. In Fig. 1, the network depicted is a pyramidal scheme (Jeon et al. 2000) with two inlets and five generated streams. In general, such a network may consist of any number of inlets of arbitrary concentration dividing into a greater number of discrete streams. At a given stage of the pyramidal network, each of P vertical streams is divided into two separate horizontal streams and then mixed with its adjacent neighbors to form P + 1 vertical streams at the next stage. In the depicted network, flowing buffer in one inlet and solution of concentration C0 in the other will result in streams of concentration 0, 0.25 C0, 0.5 C0, 0.75 C0, and C0 at the third and final stage of the network. These streams are then brought together at the input of the gradient chamber to create a linear profile. Alternatively, one could forgo the network and simply combine streams with different prepared concentrations, as was done in the gravity-driven gradient device presented by Zhu et al. (2004). Most development of this technology has come from investigators utilizing microfluidic gradient generation to study cellular chemotaxis mechanisms. Chemotaxis, directional motility induced by concentration gradients in solution, is a ubiquitous physiological mechanism underlying such diverse phenomena as inflammation (Harris 1954), embryogenesis (Behar et al. 1994), and cancer metastasis (Murphy 2001). While the most common response is for cells to move up the gradient toward a

123

0.4

0.6

0.8

1

x* = x / W

chemokinetic agent, chemorepulsion, or ‘‘fugetaxis,’’ has recently been demonstrated in leukocytes (Vianello et al. 2006). In contrast to traditional chemotaxis assays, such as Boyden (1962), Zigmond chambers (Zigmond and Sullivan 1979), and Dunn (Zicha et al. 1991) chambers, microfluidic gradient generators can achieve stable concentration profiles of arbitrary shape and maintain them indefinitely. To deliver cells to the gradient chamber, a cellular injection port is typically added to the device (Jeon et al. 2002). Cellular migration is then easily observed through a thin glass substrate. Microfluidics technology has thus presented a superior method for investigating chemotaxis mechanisms. Design improvements have led to gradient generators with complex gradient shapes and precise dynamic control. Dertinger et al. (2001) extended the pyramidal scheme to greater than two device inlets and demonstrated the creation of non-monotonic concentration functions by running multiple networks in parallel. Lin et al. (2004) modified the pyramidal scheme to create nonpolynomial gradients that were adjusted dynamically through changes in inlet flow rates. Irimia et al. (2006a) devised a rapid switching system that allows the gradient direction to be switched in less than 5 seconds. Irimia et al. (2006b) developed an entirely new gradient generation network architecture that allows the systematic creation of non-polynomial concentration profiles via a series of dividers that split and recombine flow. A property of the sampling streams generated by this network is that they have unequal flow rates. Campbell and Groisman (2007) recently developed an alternative network architecture to systematically generate non-polynomial concentration profiles, but which maintains the

Microfluid Nanofluid (2008) 4:273–285

property that each sampling stream entering the gradient chamber has equal flow rates. Recent theoretical work has focused on mathematically characterizing the network schemes. It was observed originally by Dertinger et al. (2001) that pyramidal schemes generate a profile that is polynomial in nature, and the degree of the resulting polynomial profile is one less than the number of concentration inlets. This was later proved mathematically by Sager et al. (2006), who further demonstrated that the discrete streams generated by the network can be thought of as sampling a specific function, termed the ‘‘impulse response polynomial.’’ They developed a general method for calculating the concentration of each sampling stream and the coefficients of the impulse response polynomial, given the concentrations of the inlets and the properties of the network. Figure 1c illustrates input boundary conditions and associated impulse response functions generated by a pyramidal network with three inputs and nine sampling streams. However, the authors are aware of only one rigorous model of gradient generator transport to date. Wang et al. (2006) recently presented a complete systems model of gradient generators, taking into account the transport throughout the gradient generation network as well as in the gradient chamber. Their modeling is based upon a Fourier series expansion of the concentration profile, as is the model to be described presently. The systems approach is particularly useful in the design of partial-mixing gradient generators, i.e., gradient generation schemes based upon the juxtaposition of partially mixed streams. However, while partial-mixing devices are generally more compact and can often create complex gradient shapes more naturally, complete-mixing gradient generators hold two key advantages that will ensure their continued relevance. First, the design of complete-mixing devices is more intuitive. That is, there is an obvious mathematical framework relating the design of the network to a specific concentration function one wishes to sample. Second, as we demonstrate in the present study, the behavior of complete-mixing gradient generators is largely independent of the flow rate and channel dimensions. In contrast, changing the flow rate in a partial-mixing device would not likely yield the same concentration profile in the gradient chamber and hence would require the design and fabrication of a new device. For complete-mixing gradient generators, the systems approach of modeling every single channel is unnecessarily complex, especially if complete mixing can be ensured at each stage of the network. Fortunately, it is quite simple to determine whether the serpentine mixers of a classic gradient generator sufficiently ensure complete mixing: simply calculate whether the length L  UW 2 =D: In this study, we model only the confluence of the sampling

275

streams in the gradient chamber. That is, we assume everything upstream of the gradient chamber is well known and completely mixed. As a result, we are able to derive extremely simple expressions describing the concentration profile of any arbitrary gradient generators that require only a few lines of Mathematica (Wolfram Research, Champaign, IL) code to implement. It is through this simple analytical framework that we seek to extend the quantitative understanding of transport within microfluidic gradient generators. As microfluidic gradient generators become increasingly used by both biologists and engineers, research utilizing them would be greatly aided by the development of simple design rules based upon a rigorous mass transport model. The analytical model we propose will provide such design rules for working with these devices. Finally, we utilize numerical and analytical results to determine the conditions under which secondary gradients arise naturally in these devices and may cause the net gradient direction to deviate significantly from that intended.

2 Mathematical formulations 2.1 Background theory Mass transport in gradient generators is modeled most generally by the classic species conservation equation (Deen 1998), @c þ u
rc ¼ Dr2 c; @t

ð1Þ

where c is the concentration of the analyte, u is the velocity of the fluid, and D is the molecular diffusivity. Equation (1) is based on the assumptions that the fluid is incompressible, the diffusivity is constant, and the analyte to be modeled is dilute. The geometry of the gradient chamber is defined in Fig. 1b using a Cartesian coordinate system, where x is the cross-stream coordinate (width), y is the axial coordinate (length), and z is the depthwise coordinate (height). We restrict our analysis in this study to the steady-state concentration profile (¶c/¶t = 0). Furthermore, we exclude entry effects and assume fully developed flow throughout the channel. This is done primarily to draw theoretical comparisons between different gradient generation network configurations, so that device-specific features (e.g., the angle at which incoming streams enter the channel) do not factor into the analysis. The exclusion of entry effects is justified by Kamholz et al. (2001), who implemented full 3D numerical simulations of T-mixer devices to estimate the error introduced by assuming fully developed flow at less than 5%.

123

276

Microfluid Nanofluid (2008) 4:273–285

With the assumptions of steady-state transport and fully developed flow, the only non-zero component of the velocity is the axial component; hence, u ¼ uðx; zÞ^j; where ^j is the unit vector in y. Equation (1) then reduces to
2  @c @ c @2c @2c u ¼D þ þ : @y @x2 @y2 @z2

ð2Þ

The velocity profile obtained by analytically solving the Navier–Stokes equation with no-slip boundary conditions (u = 0) at the channel walls is (Chatwin and Sullivan 1982) G ðHz  z2 Þ 2l     4GH 2 X cosh np x  12 W =H npz ;  sin 3 3 H lp n odd n cosh ½npW=2H 

uðx; zÞ ¼

ð3Þ

where –G is the axial pressure gradient (constant in magnitude), and l is the viscosity. We will make use of the depth-averaged velocity, u
,     GH 2 8GH 2 X cosh np x  12 W =H u
ðxÞ ¼  ; 12l lp4 n odd n4 cosh ½npW=2H 

ð4Þ

and the expression for average velocity, U, (Chatwin and Sullivan 1982) ( ) GH 2 192 H X tanh ½npW=2H  U¼ 1 5 : p W n odd n5 12l

ð5Þ

As the aspect ratio, e = H/W, approaches zero, the velocity profile of Eq. (3) asymptotically approaches the classic ‘‘parallel-plate’’ profile, uðzÞ ¼

6U ðHz  z2 Þ; H2

we simplify Eq. (2) to obtain the partial differential equation for depth-averaged concentration, C(x,y),
2  @C @ C @2C ¼D U þ 2 : @y @x2 @y

Analysis of Eq. (7) is simplified by scaling each term as follows: x ¼

x y y C ; y ¼ 2 ; C ¼ ¼ ; W W U=D WPe C0

@C @C 1 @C ¼ þ : @y @x2 Pe2 @y2

ð9Þ

The solution to Eq. (9) was presented previously by Wu and Nguyen (2005) using a Fourier series expansion, with no-flux boundary conditions (¶C*/¶x* = 0) imposed at the lateral walls in x and a fully-mixed boundary condition (¶C*/¶y* = 0) imposed at y* fi ¥. Our analysis differs in the input boundary condition. Instead of assuming a specific form for the concentration profile at the entrance of the gradient chamber, we begin by specifying the general Dirichlet condition at y* = 0, C  jy ¼0 ¼ f ðx Þ:

ð10Þ

As a result, we obtain the concentration profile for any general input function to the gradient chamber, C  ðx ; y Þ ¼ A0 þ

1 pffiffiffi X 2 An cosðnpx Þ n¼1

2.2 2D analytical model

ð8Þ

where Pe = U W/D is the Pe´clet number, a measure of the ratio of characteristic velocities of convection to diffusion, and C0 is the maximum input concentration. Thus, Eq. (7) is transformed to

! 2 n2 p2  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiy ; exp  1 þ 1 þ 4n2 p2 =Pe2

ð6Þ

where the velocity profile is parabolic in the shallow depthwise (z) direction and constant in x. For e  1; the simplified model of Eq. (6) is valid throughout the crosssection, except near the lateral cross-stream sidewalls, where the velocity decays to zero over a distance of approximately one channel thickness, H (Deen 1998).

ð7Þ

ð11Þ

where

A0 ¼

Z1

f ðx Þdx

ð12Þ

0

and pffiffiffi Z 2 f ðx Þ cosðnpx Þdx : 1

Equation (2) may be further simplified in the small aspect ratio limit, e  1: Given the small length scale in z, we assume that depthwise variation in concentration is sufficiently small to justify a reduction in dimensionality. Moreover, we assume that the solute is convected with the average velocity of the fluid. Based on these assumptions

123

An ¼

ð13Þ

0

The Pe´clet number under most microfluidic flow conditions is in the range of 10–105 (Stone et al. 2004). For flows in the high Pe´clet regime (large width, high velocity,

Microfluid Nanofluid (2008) 4:273–285

277

low diffusivity), the axial diffusion term can be neglected and Eq. (11) becomes C  ðx ; y Þ ¼ A0 þ

1 pffiffiffi X   2 An cosðnpx Þ exp n2 p2 y : n¼1

ð14Þ Equation (14) demonstrates the generalization of a given gradient generator configuration under all flow conditions Pe  1. Since Eq. (14) does not contain the Pe´clet number, the concentration profile depends only on the properties of the generation scheme upstream of the gradient chamber given by f(x*). 2.3 Coefficients for gradient generators Having derived the general form of the concentration profile, we will now present simple formulas for calculating the A0 and An coefficients that do not require symbolic integration. Gradient generators most generally consist of a scheme for creating an arbitrary number, N, of discrete inputs of homogeneous concentration, which are brought together at the entrance of the gradient chamber. For the assumed velocity profile, Eq. (6), to be valid, it is required that each input stream have the same viscosity (Wu and Nguyen 2005), a reasonable assumption for streams varying only by the concentration of a dilute species. We denote C*i as the concentration of the ith input stream, which can be calculated generally using methods described elsewhere [e.g., Sager et al. (2006)]. For gradient generators where all inputs have equal flow rates, as is the case where the sampling streams are generated by the classic pyramidal network, the integrals in Eq. (12) and (13) simplify to the sums A0 ¼

N 1X C N i¼1 i

ð15Þ

A0 ¼

N 1X C  qi Q i¼1 i

ð17Þ

and An ¼

pffiffiffi N 1 2X np

ðCi



 Ciþ1 Þ sin

i¼1

np

Pi

j¼1

qj

Q

! ð18Þ

:

2.4 Discrete Fourier transform solution As a final step in the generalization of solutions to the 2D model, we use the Fourier transform in the analysis of Eq. (9) to get   @ C^ 1 @ C^ 2 ^ C ¼ x þ ; Pe2 @y2 @y

ð19Þ

where C^ is the Fourier transform of C* and x is the spatial frequency. Applying the same boundary conditions in y as before, we obtain the solution ! 2 2x pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiy ; C^ ðx; y Þ ¼ f^ðxÞ exp  1 þ 1 þ 4x2 =Pe2

and pffiffiffi N 1
 2X  npi  An ¼ ðCi  Ciþ1 Þ sin : np i¼1 N

unequal flow rates to generate non-polynomial impulse functions. Another use of unequal flow rates is employed by the addition of ‘‘shield flows’’ on one or both sides of the input channel to focus and compress a gradient dynamically. For streams with unequal flow rates, a slightly more general form of the above equations can be expressed. Following Wu and Nguyen (2005), we assume that each stream has the same fluid density so that the ‘‘size’’ of each stream in the input function (ignoring hydrodynamic focusing effects occurring at the entrance of the channel) will be determined by the mass fraction of the solvent, which is given by the relative flow rate. Therefore, let qi be the flow rate of the ith stream and Q the total flow rate. The coefficients A0 and An are then



ð16Þ

Under this paradigm, the concentration profile of a T-mixer (i.e., a long mixing channel with two streams leading into it) can be described with parameters N = 2, C*1 = 0, and C*2 = 1. Thus, our model treats T-mixers as a special case of gradient generator. A number of recently developed gradient schemes employ unequal flow rates. The gradient generation scheme recently developed by Irimia et al. (2006b) uses

ð20Þ

where f^ðxÞ is the Fourier transform of the input function, f(x*). Thus, the evolution of the cross-stream concentration profile along the channel axis may be thought of as the convolution of the channel input function, f(x*), and a blurring function dependent upon the axial distance, g(y*), which has a form very similar to the classic Gaussian kernel. Expressed in Fourier space, this becomes, by the convolution theorem,  C^ ðx; y Þ ¼ f^ðxÞ^ gðx; y Þ:

ð21Þ

123

278

(a)

(d)

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0 0

0.2

0.4

0.6

0.8

1

0

0.4

0.6

0.8

1

x* add mean remove padding FFT -1

subtract mean pad symmetrically FFT

(b)

(c)

800

800

600

600 ^ g(ω,y*)

400 200

400 200

0

0 0

50

100

150

200

250

0

We have thus reduced the entire gradient generator mass transport problem into a simple multiplication of two functions. However, we have so far neglected to mention that the Fourier transform is defined on an infinite interval. Therefore, it does not strictly apply to solution of Eq. (9). Nevertheless, a very convenient approximation to the solution can be achieved using the discrete Fourier transform (DFT). The basic method of solution using the DFT is illustrated in Fig. 2 and involves four major steps. First, a discrete set of points representing the concentration profile on the domain is obtained and the mean is subtracted. These points could either come from experimental measurements or from sampling the model input boundary condition. The no-flux boundary conditions in x are imposed by symmetrically padding the input values in both directions. Second, the fast Fourier transform (FFT) is applied to obtain a discrete approximation to f^ðxÞ: Third, the input function is multiplied by the blurring function g^ðx; y Þ given in Eq. (20). Finally, the inverse FFT is applied, the padding is removed, and the mean is added back, to obtain an approximation to the solution.1 This approach is especially useful under certain circumstances. For instance, consider the case of trying to determine the concentration profile downstream based on a set of experimental measurements recorded at the entrance to the channel. With the analytical Fourier series solution, one would have to interpolate between points to achieve a continuous function and then symbolically integrate that 1 We implemented the DFT solution algorithm in MATLAB (Mathworks, Inc., Natick, MA). Further details are available from the authors upon request.

50

100

150

200

250

ω

ω

123

0.2

x*

·

Fig. 2 The discrete Fourier transform (DFT) method for calculating concentration profiles. a An input function, f(x*), represented by 80 samples. b DFT of input function, f^ðxÞ: c Multiplication of b by the blurring function, g^ðx; y Þ: d Concentration profile at a distance y* downstream. For comparison, the exact solution is plotted as a solid line

Microfluid Nanofluid (2008) 4:273–285

function. However, it is much more efficient to apply the blurring convolution directly to the discrete values, as we have just demonstrated. Another instance where the DFT approach would be useful is if the input function one wishes to model cannot be easily integrated symbolically. With the DFT approach, one would simply sample the function in order to obtain sufficient information to model its transport. 2.5 3D numerical model To test the validity of the assumptions utilized in reducing the transport equations to a two-dimensional model that could be solved analytically, we constructed a simple and stable numerical algorithm to solve the full 3D transport represented by Eq. (2) (assuming steady-state and fully developed flow). For high Pe flows, we can again assume that axial diffusion is negligible. Thus, Eq. (2) becomes u ðx ; z ; eÞ

@c @ 2 c @ 2 c ¼ þ ; @y @x2 @z2

ð22Þ

where z* = z/W and u* = u/U is the normalized velocity at a given point, computed from Eqs. (3) and (5). There are computational and conceptual advantages to solving Eq. (22) as opposed to Eq. (2). The advantage computationally is that y* becomes a ‘‘one-way’’ space coordinate. That is, the concentration at a point on the domain is only influenced by the concentration at previous values of y*. In that sense, it is not unlike a time coordinate. Therefore, we only need to solve for the concentrations on a single x*–z* crosssectional plane at a time and then integrate to move on to solving the next plane downstream. The conceptual

Microfluid Nanofluid (2008) 4:273–285

advantage of Eq. (22) is that we must consider only one parameter in analyzing the behavior of the solution, e. Extensive analysis of Eq. (22) was previously done by Kamholz and Yager (2002) in analyzing the effects of channel sidewalls in T-mixer mass transport. We implemented no-flux boundary conditions at the walls of the channel and specified the input condition at y* = 0 as before with Eq. (10), except imposed uniformly in z. The solution to Eq. (22) was obtained via a finite volume approximation (Pantankar 1980). We used a highresolution quadrilateral mesh aligned with the direction of the flow so that numerical diffusion effects were nonexistent. We employed a systematic procedure for ensuring grid-converged solutions at a desired cross-sectional plane. First, we set u* = 1 and increased the number of elements in x and y until the numerical solution equaled Eq. (14), the exact solution in the uniform velocity case, to less than 10–4 (no variation in z considered at this point). We then added the effect of the cross-stream sidewalls by setting u* equal to Eq. (4) divided by Eq. (5) (still considering no variation in z). We then continued to increase the number of elements until the L2 norm error with respect to Eq. (14) had converged to within two significant figures of its gridconverged value. Finally, we moved to the full 3D problem by setting u* equal to the full velocity profile of Eq. (3) normalized by Eq. (5). We then increased the number of elements in z until the L2 norm error of the depth-averaged concentration with respect to Eq. (14) had converged to within two significant figures of its grid-converged value. We found that the most computationally efficient results were obtained by using a graded mesh that became very fine near the transverse sidewalls and points of discontinuity in the input boundary condition (e.g., at x* = 0, 0.2, 0.4, 0.6, 0.8, and 1 for an N = 5 linear gradient generator). Furthermore, we found that increasing the spacing between cross-sectional planes from fine to coarse along the channel axis in y* was more computationally efficient than uniform spacings. As a final cost-saving measure, we note that since u is symmetric about z = H/2, one only needs to solve Eq. (22) on half of the domain. Mesh sizes and solution times varied somewhat depending on the y* location and the input boundary condition. For instance, the results presented in Sect. 3.4.1 for an N = 5 sampling stream linear gradient generator, with e = 0.25 and y* = 0.00272, were obtained using a 796 · 862 · 199 (x* by y* by z*) element mesh. The resulting linear systems were iteratively solved using the line-by-line TDMA method (Pantankar 1980). The algorithm was custom-coded in Fortran 902 and typically took several hours on a 2.0 GHz Intel Xeon processor.

2

Further details are available from the authors upon request.

279

3 Results and discussion 3.1 Analytical concentration profiles In Fig. 3, we display concentration profiles described by Eq. (11) using different entrance boundary conditions which illustrate the versatility of the analytical method we propose for solving a wide variety of problems. Figure 3a is a plot of an N = 9 quadratic gradient, whose entrance boundary condition was derived by a standard pyramidal scheme. Thus, each sampling stream had the same flow rate, and the Fourier coefficients were given by Eqs. (15), (16). Figure 3b is an exponential gradient whose boundary condition was derived by the recently proposed ‘‘Universal Gradient Generator’’ (Irimia et al. 2006b). Each sampling stream had a different flow rate; thus, the Fourier coefficients were given by Eqs. (17), (18). In our final example, Fig. 3d, we used the non-monotonic boundary condition composed piecewise of linear and power-law functions shown in Fig. 3c. The Fourier coefficients for this gradient were obtained by symbolically integrating this function in Eqs. (12) and (13). It is clear that regardless of the entrance boundary condition, the axial evolution of concentration follows a general pattern. At y = 0, the profile is the piecewise staircase function generated by confluence of the discrete streams. As we move farther downstream, diffusion gradually smooths out this function. Eventually, the concentration profile approaches the average concentration at its fully mixed state. 3.2 Error analysis of 2D model Figure 4, a direct comparison of concentration values generated by the 2D analytical model of Eq. (14) and full 3D numerical model of Eq. (22), illustrates excellent agreement at all axial positions. Figure 4 also illustrates several key principles of the spatial distribution of approximation errors on the domain. The approximation error is greatest near the sidewalls and minimal at the center. As the axial position increases, this error increases and becomes significant throughout a greater portion of the channel, until a maximum value of error is reached after which both numerical and analytical solutions approach the same fully mixed state as y fi ¥. Therefore, the primary contribution to the observed error is a breakdown in the constant velocity assumption near the lateral sidewalls in x. The specific nature of this spatial distribution is illustrated in Fig. 5, which displays approximation error as a function of channel position. Figure 5 was generated by solving Eq. (22) with a linear boundary condition (i.e., f(x*) = x*) so that there was no bias towards one particular type of gradient configuration. Figure 5 shows that the maxi-

123

280

Microfluid Nanofluid (2008) 4:273–285

Concentration, C*

(a)

(b) 1

0.5 0.4 0.3 0.2 0.4

0.1 0 0

Concentration, C*

Fig. 3 Concentration profiles predicted by Eq. (11). a The concentration profile of a ninestream (N = 9) quadratic gradient device. b The concentration profile of a sixstream (N = 6) exponential gradient device generated by the ‘‘Universal Gradient Generator’’ proposed by Irimia et al. (2006b). The size of each stream is proportional to its relative flow rate. c A nonmonotonic boundary condition composed piecewise of linear and power-law functions. d The concentration profile generated by the entrance boundary condition in c. For all results, Pe = 10

0.8 0.6 0.4 0.2 0.4

0.3 0.2

0.4

0.3

0 0

0.2

0.2

0.1

0.6

0.8

x/W

1

0.2 0.4

y/W

0

0.6

0.1 0.8

x/W

1

y/W

0

0.10

y / W 0.075

(d) Concentration, C*

Concentration, C*

(c) 1 0.8 0.6 0.4 0.2

0.05 0.025 0

1 0.8 0.6 0.4 0 0.25 0.5

0.2

0.4

0.6

0.8

1

x/W

x/W

y* = 0.00005

1

0.75 1

3.3 Gradient optimization analysis

y* = 0.01500 0.9 y* = 0.05000

Concentration, C*

0.8 0.7

y* = 0.15000 0.6 where y* = y / (W Pe)

0.5 0.4 0.3 0.2 0.1 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x* = x / W

Fig. 4 Comparison of 2D analytical model (solid) to numerically calculated full 3D model (dashed). A five-stream (N = 5) linear gradient generator was modeled at high Pe and e = 0.1 at several axial locations. The analytical results are from Eq. (14), and the numerical results are from Eq. (22), averaged through the channel depth

mum (L¥) error over the entire domain is only approximately 1.2% for e = 0.05, 2.2% for e = 0.1, and 3.8% for e = 0.2. Thus, the approximation errors decrease as the aspect ratio approaches zero. However, even for moderate aspect ratios that do not meet the required e  1 limit, such as e = 0.25, we find that the analytical model is extremely accurate over the length scales where the gradient is optimal (i.e., within 1% of W Pe). For e = 0.25, we calculated approximation errors of 0.383% (L2) and 0.828% (L¥) for an N = 5 stream linear gradient generator at y*opt = 0.00272 (the optimal axial location as determined using methods in Sect. 3.3).

123

Our analytical work immediately introduces a method of calculating simple design rules governing the ideal axial location where the transverse concentration profile is optimized with respect to a desired ‘‘target function,’’ g(x*)3. These design rules follow the simple form y* = y*opt, where y*opt is a constant depending solely on the target function and the properties of the mixing channel input (e.g., input stream concentrations, number of input streams). The method we describe allows experimentalists to quickly determine the ideal axial location of the gradient under virtually any flow conditions and gives designers a rule-of-thumb estimate for placing device features such as cellular injection ports. The idea of treating the concentration profile as an optimization problem for finding the ideal axial position for use in experiments was originally developed by Walker et al. (2005). They demonstrated how adjusting the flow rate changed the position of the optimal axial position and the length scale over which the gradient was sufficiently good (higher flow rates meant farther downstream optimal positions and larger length scales). By scaling with the Pe´clet number, however, we have shown that the concentration profile folds neatly into the single description given by Eq. (14), regardless of flow rate, channel size, or diffusivity. Therefore, it is redundant to consider the optimization problem at different flow rates. Instead, determining 3

Typically, the target function, g(x*), and the impulse response function, h(x*), are equivalent. However, it may be necessary to approximate a complicated target function with a simpler impulse response function, and thus it would be appropriate to minimize the RMS deviation with respect to the target function.

Microfluid Nanofluid (2008) 4:273–285

281

0.5

3.5 3

y*

Approximation Error (Infinity Norm)

H/W = 0.2

H/W = 0.1

2.5

2.2%

0.45

2.0

0.4

1.8

0.35

1.6 1.4

0.3

1.2

0.25

1.0 0.2

2

0.8

0.15

1.5

0.6

0.1

H/W = 0.05

0.4

0.05

1

0

0.2 0

0.2 0.4 0.6 0.8 1

x* = x / W

0.5 0 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

y* = y / (W Pe)

Fig. 5 Spatial distribution of approximation errors introduced by simplifying the model to two dimensions. Solid L¥ approximation error of Eq. (14) with respect to the depth-averaged 3D profile of Eq. (22). Inset Absolute value of difference between 2D model (Eq. 14) and depth-averaged 3D model (Eq. 22) for each position in the channel. The results were obtained using the linear input boundary condition f(x*) = x*. The horizontal bars show the channel crosssection for each curve

a single optimal scaled axial position that is always true for a given input boundary condition provides us with a convenient platform for comparing different types of gradient generator configurations. We first assume that the Pe´clet number is high, so that the concentration profile is approximated by Eq. (14). The criterion we use to determine the optimal location is the root-mean-squared deviation4 (E) as a function of axial position,

Eðy Þ ¼

8 1
½C ðx ; y Þ  gðx Þ2 dx

0

912 = ;

:

ð23Þ

Consider as an example an N = 5 linear gradient generated using the standard pyramidal scheme (Jeon et al. 2000; Dertinger et al. 2001). The impulse response polynomial for the device is h(x*) = (5/4) x* – 1/8. We take this polynomial to be the target function, but restrict the contribution of the error to only the interval where the function h is positive, [1/10,9/10]. The equation for the error is then 91 8 2 9=10  >
 2 > = <5 Z 5 1 Eðy Þ ¼ x  C  ðx ; y Þ  dx : > > 4 8 ; :4

We numerically minimized this function to determine an optimal axial distance of y*opt = 0.00272. The optimized cross-stream concentration profile is plotted in Fig. 6, along with the target function. In a similar manner, we calculated y*opt coefficients for pyramidal linear gradient generators with sampling streams of 2– 12. The coefficients and the associated E(y*opt) are listed in Table 1. To illustrate the properties of the modeling error functions, we plotted E(y*) for linear gradient generators of sampling streams 2–6 in Fig. 7. The error starts at a local maximum at y* = 0 where the transverse concentration profile is a staircase function. Then, the error decreases steeply until a minimum is attained at the optimal axial distance, y* = y*opt. Farther along the channel axis, the error increases and asymptotically approaches a maximum as the transverse concentration profile becomes fully mixed as y* fi ¥. Figure 7 also demonstrates the effect of N, the number of discrete streams generated by the pyramidal mixing scheme, on the error profile. As N increases, the optimal axial position, y*opt moves closer to the channel entrance. The corresponding minimized RMS deviation at that position decreases as well. These results provide theoretical support to the insight of Sager et al. (2006) that the N streams may be thought of as sampling a desired function, where the granularity of the sampling decreases with N. 3.4 Secondary concentration gradients In this section, we consider the gradient of the concentration profile, rc ¼

@c ^ @c ^ @c ^ i þ j þ k; @x @y @z

ð25Þ

1

0.8

Concentration, C*

4%

0.6

0.4

0.2

ð24Þ

1=10

4 Other norms may be considered. However, the general principle of calculating a single constant representing all flow conditions Pe  1 remains.

0.2

0.4

0.6

0.8

1

x* = x / W

Fig. 6 The optimal transverse concentration profile for a five-stream (N = 5) linear gradient generator. The concentration profile (solid) is given by Eq. (14), at y*opt = 0.00272. Also shown is the impulse response function, h(x*) = (5/4) x* – 1/8 (dashed), and the staircase input function (solid)

123

282

Microfluid Nanofluid (2008) 4:273–285

Table 1 Design rules governing the optimal axial location of linear gradient generators. Optimal locations were calculated by minimizing the root-mean-square deviation, E(y*), with respect to the associated impulse response polynomial, h(x*) = mx* + b N

m

b

y*opt

E(y*opt)

2

2

–1/2

0.01413374

0.0245309

3

3/2

–1/4

0.00688257

0.0115719

4 5

4/3 5/4

–1/6 –1/8

0.00409252 0.00272487

0.0072030 0.0050923

6

6/5

–1/10

0.00195109

0.0038722

7

7/6

–1/12

0.00146961

0.0030871

8

8/7

–1/14

0.00114899

0.0025442

9

9/8

–1/16

0.00092438

0.0021491

10

10/9

–1/18

0.00076070

0.0018502

11

11/10

–1/20

0.00063760

0.0016171

12

12/11

–1/22

0.00054260

0.0014309

^ are the unit vectors in x, y, and z, where ^i; ^j; and k respectively. The purpose of the microfluidic devices we have hereto discussed is to develop and maintain a uniform concentration gradient of given magnitude in a single direction, the cross-stream (x) direction. We shall therefore refer to gradient component ¶c/¶x as the ‘‘primary’’ concentration gradient. However, as we will demonstrate, variation in concentration in the y and z coordinates may lead to ‘‘secondary’’ concentration gradient components of comparable magnitude. We shall refer to gradient components ¶c/¶y and ¶c/ ¶z as the axial and depthwise secondary gradients, respectively. If these secondary gradients exist at a given location, then the true gradient is not in the expected direction, but at some angle to it. For chemotaxis studies, secondary gradients may be unacceptable, since the infor-

Modeling Error, E(y*)

25%

20

15 Increasing N sampling streams 10

5

N=2 N=3

0.01

0.02

0.03

0.04

0.05

y* = y / (W Pe)

Fig. 7 Error as a function of scaled axial position for linear gradient generators with N = 2 to 6 sampling streams. The minimum of each curve represents the optimal axial position

123

mation delivered to cells by the device will not be exactly that intended. We use our analytical and numerical modeling tools to investigate secondary concentration gradients. 3.4.1 Depthwise gradients Since the velocity profile is nonhomogeneous in the depthwise (z) direction, the solute is convected at different speeds corresponding to which streamline it is carried by, and thus concentration also varies as a function of z (a contour plot of the velocity profile of Eq. (3) is shown in Fig. 8a). The development of z gradients is simultaneously opposed by diffusion between streamlines, the characteristic time of which decreases with the square of the channel height. As suggested by Beard (2001), we can expect that concentration variations in z will be fully relaxed at an axial distance y once the solute residence time in the channel, y/U, is much greater than the characteristic time for depthwise diffusion, H2/D. In our scaled notation, the inequality y/U  H2/D becomes y  e2 :

ð26Þ

We showed in Table 1 that the optimal value of y* for linear gradient generators is approximately 0.01 for N = 2 sampling streams and is even less for higher N. Thus, the inequality of Eq. (26) does not appear to be satisfied for normal-sized values of e (i.e., e = 0.05–0.25) even at the most important downstream locations in the channel. An exhaustive study of depthwise concentration variation as a function of e and y* in T-mixers was presented by Chen et al. (2006). To investigate depthwise gradients, we used the numerical solution of Eq. (22) for an N = 5 stream linear gradient generator with aspect ratio e = 0.25 at y*opt = 0.00272. We then calculated the derivatives in x* and z* and constructed a map of the net gradient throughout the channel cross-section (assuming no axial secondary gradient). The gradient of the entire cross-section is shown in Fig. 8c, and a smaller portion of the cross-section is shown in Fig. 8d to illustrate the gradient at higher resolution. We observe from Fig. 8 that depthwise secondary gradients exist throughout the channel cross-section and are especially significant near the sidewalls of the channel where x* < 1/N and x* > (N–1)/N for gradient generators with streams of equal flow rates. For chemotaxis studies, these results lead us to recommend that one not place cells within a distance of either H or W/N from the cross-stream sidewalls of the gradient chamber in order to ensure consistency of information presentation to the cells. Within a distance of H, the velocity of solvent is slowed by the proximity to the sidewalls, and hence the shear stress on

Microfluid Nanofluid (2008) 4:273–285

(a) 0.25 z* = z / W

0.2 0.15 0.1 0.05 0

(c)

0.25

z* = z / W

(b)

0.15

0.2

0.1 0.05 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x* = x / W

(d) 0.125 z* = z / W

Fig. 8 Depthwise concentration gradients in an N = 5 stream linear gradient generator for aspect ratio e = H/W = 0.25 at y*opt = 0.00272. a Contour plot of the velocity profile, Eq. (3). b For comparison purposes, the input boundary condition at y* = 0. c A map of the net gradient throughout the channel cross-section. The maximum arrow length corresponds to a gradient magnitude of 1.7. d A portion of the cross-section magnified to illustrate the gradient in higher resolution. The maximum arrow length corresponds to a gradient magnitude of 1.6

283

0.1 0.08

Direction of primary gradient

0.06 0.04 0.02 0

0

0.05

0.1

0.15

0.2

0.25

0.3

x* = x / W

the cell is different. Within a distance of W/N, the primary gradient component is smaller, and there are also significant depthwise variations that cause the net direction of the gradient to be at an angle in the cross-sectional plane. 3.4.2 Axial gradients We noted in Sect. 2.2 that W Pe is a scale for y and W is a scale for x. Thus, for Pe ~ 1, the concentration changes over the same length scale in y as it does in the primary direction, x. That is, ¶C/¶y ~ ¶C/¶x, and the net direction of the gradient will be at an angle to the primary direction, pointing either partially upstream or downstream. We investigated axial secondary gradients using the 2D analytical model defined by Eq. (11). The analytical model is particularly convenient for this purpose, since expressions for the partial derivatives are easily obtained by differentiating the Fourier series term-by-term. For this analysis, we again considered N = 5 linear gradient generators. The values of gradient components ¶C/¶x and ¶C/ ¶y were calculated for Pe´clet numbers of Pe = 15 and Pe = 50 at the respective optimal axial distances y*opt = 0.00512 and 0.00332 (determined using the methods described in Sect. 3.3) and plotted in Fig. 9. Figure 9 illustrates the significance of the axial gradient component at different Pe´clet. At Pe = 15, the magnitude of the normalized axial gradient component, ¶C*/¶(y/W), reaches a maximum value of 0.324 on the working interval

(x* = 0.2 through 0.8), which is approximately 26% of the normalized primary gradient value, ¶C*/¶x* = 1.25. The result is a maximum deviation in gradient direction of 14.5 from the primary direction. At Pe = 50, the magnitude of the normalized axial gradient component reaches a maximum value of 0.105, which is 8.4% of the primary gradient value, for a maximum deviation in gradient direction of 4.81. As noted in Sect. 2.2, almost all microfluidic flow conditions involve high Pe´clet. Nevertheless, one could envision the need to present to cells gradients of solutes with high diffusivity, such as oxygen or hydrogen, or in environments with slow-moving fluid (low shear stress). Such conditions may result in low Pe´clet numbers. Our results suggest that there is a lower bound on the Pe´clet number that may be successfully used in the type of microfluidic gradient generators considered in this study. For low Pe´clet, the net gradient will not point in the intended direction. We suggest that groups planning gradient generator experiments involving Pe < 100 model the transport with Eq. (11) to determine whether the axial gradients are sufficiently small for their purposes. Walker et al. (2005) and Sai et al. (2006) observed that cells in gradient generators migrate at an angle downstream dependent upon the relative values of the flow rate and the magnitude of the induced gradient. They hypothesize that shear stress and concentration gradients exert relative tractile forces on chemotaxing cells, and therefore the net

123

284

Microfluid Nanofluid (2008) 4:273–285

Concentration Gradient Magnitude

(a)

(b)

Pe = 15

1.5

Pe = 50

1.5 Primary

1

1 0.5

0.5 Secondary (Axial)

0.2

0.4

0.6

0.8

1

0.2

0.4

0.6

0.8

1

-0.5

-0.5

x* = x / W

x* = x / W

Fig. 9 Comparison of the primary gradient component, ¶C/¶x (solid), and the secondary axial gradient component ¶C/¶y (dashed) at different Pe. a Pe = 15. Components are plotted at y*opt = 0.00512.

b Pe = 50. Components are plotted at y*opt = 0.00332. Each gradient component was normalized by multiplying by W and dividing by C0

migration direction is a combination of the two forces. In both studies, the minimum flow rate was 1 lL/min, which corresponded to a minimum Pe´clet number of 667. At Pe = 667, the maximum deviation in gradient angle pointing upstream or downstream is less than 1. Thus, axial secondary gradients, as we have defined them, do not account for the observed angular migration pattern. Finally, we note that the modeling methods we have proposed do not take into account disruption of the flow field and gradient by the introduction of an obstruction (e.g., biological organism) to the gradient chamber. Undoubtedly, the magnitude of such a disruption would depend on the relative size and shape of the obstruction and its location in the channel, along with potentially several other factors. An investigation elucidating these factors and a model detailing the concentrations seen at the surface of the obstruction would be a valuable addition to the study of gradient generators. If the size of the obstruction is much smaller than the size of the channel, then it may not be necessary to model the entire channel, but simply the local space around the obstruction. The modeling tools we have presented may be useful in providing boundary conditions for such analysis.

with a set of discrete points from experimental measurements. By numerically solving the full 3D conservation equation, we find that the analytical 2D model is accurate to within a worst-case error of 4% with respect to the depthaveraged concentration for aspect ratios of e = 0.2 or less. Additionally, we used our numerical and analytical results to investigate the existence of secondary gradient components that cause the net direction of the gradient to deviate from the intended direction, and we developed rule-ofthumb estimates for conditions under which such deviations may be significant. Moreover, the high Pe´clet assumption (which is valid for typical microfluidic flow conditions) leads to an extremely convenient simplification to the analytical model. By relating the physical dimensions to the Pe´clet number, we find that for Pe  1 the model folds into a single description encompassing all flow rates, solute diffusivities, and channel dimensions. The only parameter influencing the behavior of the model is the input (y = 0) boundary condition provided by the specific gradient generation network upstream of the gradient chamber. Thus, the modeling tools we have presented lead to direct comparisons of different gradient generation networks. Despite recent interest in developing network architectures that can generate concentration profiles representing any function, there currently exists no method known to us for evaluating the performance of competing schemes that is based upon a model of the concentration profile in the gradient chamber itself. Our modeling tools may provide such a method. Finally, we investigated simple design rules for gradient generators. We extended the work of Walker et al. (2005) in developing a method for calculating the scaled axial position where the concentration profile is optimized with respect to a given target function. The advantage of our approach is that only one coefficient representing all

4 Conclusions We have developed an analytical 2D model describing transport in microfluidic gradient generators. This model allows the implementation of any piecewise continuous boundary condition at the entrance to the channel. We derived equations for implementing the boundary conditions provided by gradient generators with equal flow rate inputs (Eqs. 15, 16) and unequal flow rate inputs (Eqs. 17, 18). Additionally, by applying the discrete Fourier transform, we present an alternative method to using Fourier series that is particularly convenient when one is working

123

Microfluid Nanofluid (2008) 4:273–285

flow conditions is necessary to describe a given gradient generator. We anticipate that these results will assist the design and experimental usage of microfluidic gradient generators. Acknowledgements We thank Michael Hwang for programming assistance. We are further indebted to Drs. G. Kane Jennings, Dmitry Markov, and Robert Roselli for reviewing early drafts of this manuscript. The Advanced Computing Center for Research and Education (ACCRE) at Vanderbilt graciously provided computational resources. We acknowledge support from the NSF-sponsored VaNTH ERC, the Systems Biology/Bioengineering Undergraduate Research Experience (SyBBURE), the Vanderbilt Institute for Integrative Biosystems Research and Education (VIIBRE), NIH Grant 5U01AI061223, and a Whitaker Foundation Special Opportunity Award.

References Beard DA (2001) Taylor dispersion of a solute in a microfluidic channel. J Appl Phys 89:4667–4669 Behar T, Schaffner A, Colton C, Somogyi R, Olah Z, Lehel C, Barker J (1994) GABA-induced chemokinesis and NGF-induced chemotaxis of embryonic spinal cord neurons. J Neurosci 14:29–38 Boyden S (1962) The chemotactic effect of mixtures of antibody and antigen on polymorphonuclear leucocytes. J Exp Med 115:453– 466 Campbell K, Groisman A (2007) Generation of complex concentration profiles in microchannels in a logarithmically small number of steps. Lab Chip (in print). doi:10.1039/b610011b Chatwin PC, Sullivan PJ (1982) The effect of aspect ratio on longitudinal diffusivity in rectangular channels. J Fluid Mech 120:347–358 Chen JM, Horng TL, Tan WY (2006) Analysis and measurements of mixing in pressure-driven microchannel flow. Microfluid Nanofluid 2:455–469 Deen WM (1998) Analysis of transport phenomena. Oxford University Press, New York Dertinger SKW, Chiu DT, Jeon NL, Whitesides GM (2001) Generation of gradients having complex shapes using microfluidic networks. Anal Chem 73:1240–1246 Harris H (1954) Role of chemotaxis in inflammation. Physiol Rev 34:529–562 Irimia D, Liu SY, Tharp WG, Samadani A, Toner M, Poznansky MC (2006a) Microfluidic system for measuring neutrophil migratory responses to fast switches of chemical gradients. Lab Chip 6:191–198 Irimia D, Geba DA, Toner M (2006b) Universal microfluidic gradient generator. Anal Chem 78:3472–3477 Jeon NL, Dertinger SKW, Chiu DT, Choi IS, Stroock AD, Whitesides GM (2000) Generation of solution and surface gradients using microfluidic systems. Langmuir 16:8311–8316 Jeon NL, Baskaran H, Dertinger SKW, Whitesides GM, de Water LV, Toner M (2002) Neutrophil chemotaxis in linear and complex

285 gradients of Interleukin-8 formed in a microfabricated device. Nat Biotechnol 20:826–830 Kamholz AE, Yager P (2002) Molecular diffusive scaling laws in pressure-driven microfluidic channels: Deviation from onedimensional Einstein approximations. Sens Actuator B Chem 82:117–121 Kamholz AE, Schilling EA, Yager P (2001) Optical measurement of transverse molecular diffusion in a microchannel. Biophys J 80:1967–1972 Lin F, Saadi W, Rhee SW, Wang SJ, Mittal S, Jeon NL (2004) Generation of dynamic temporal and spatial concentration gradients using microfluidic devices. Lab Chip 4:164–167 Murphy P (2001) Chemokines and the molecular basis of cancer metastasis. N Engl J Med 345:833–835 Nguyen NT, Wu ZG (2005) Micromixers—a review. J Micromech Microeng 15:R1–R16 Pantankar SV (1980) Numerical heat transfer and fluid flow. Hemisphere, New York Rhoads DS, Nadkarni SM, Song L, Voeltz C, Bodenschatz E, Guan JL (2005) Using microfluidic channel networks to generate gradients for studying cell migration. Methods Mol Biol 294:347–357 Sager J, Young M, Stefanovic D (2006) Characterization of transverse channel concentration profiles obtainable with a class of microfluidic networks. Langmuir 22:4452–4455 Sai J, Walker G, Wikswo J, Richmond A (2006) The IL sequence in the LLKIL motif in CXCR2 is required for full ligand-induced activation of Erk, Akt, and chemotaxis in HL60 cells. J Biol Chem 281:35931–35941 Squires TM, Quake SR (2005) Microfluidics: Fluid physics at the nanoliter scale. Rev Mod Phys 77:977–1026 Stone HA, Stroock AD, Ajdari A (2004) Engineering flows in small devices: Microfluidics toward a lab-on-a-chip. Annu Rev Fluid Mech 36:381–411 Vianello F, Papeta N, Chen T, Kraft P, White N, Hart W, Kircher M, Swart E, Rhee S, Palu G, Irimia D, Toner M, Weissleder R, Poznansky M (2006) Murine B16 melanomas expressing high levels of the chemokine stromal-derived factor-1/CXCL12 induce tumor-specific T cell chemorepulsion and escape from immune control. J Immunol 176:2902–2914 Walker GM, Sai JQ, Richmond A, Stremler M, Chung CY, Wikswo JP (2005) Effects of flow and diffusion on chemotaxis studies in a microfabricated gradient generator. Lab Chip 5:611–618 Wang Y, Mukherjee T, Lin Q (2006) Systematic modeling of microfluidic concentration gradient generators. J Micromech Microeng 16:2128–2137 Wu ZG, Nguyen NT (2005) Convective-diffusive transport in parallel lamination micromixers. Microfluid Nanofluid 1:208–217 Zhu X, Chu L, Chueh B, Shen M, Hazarika B, Phadke N, Takayama S (2004) Arrays of horizontally-oriented mini-reservoirs generate steady microfluidic flows for continuous perfusion cell culture and gradient generation. Analyst 129:1026–1031 Zicha D, Dunn GA, Brown AF (1991) A new direct-viewing chemotaxis chamber. J Cell Sci 99:769–775 Zigmond SH, Sullivan SJ (1979) Sensory adaptation of leukocytes to chemotactic peptides. J Cell Biol 82:517–527

123