Two-Aperture Microfluidic Probes as Flow Dipole: Theory and Applications

A microfluidic probe (MFP) is a mobile channel-less microfluidic system under which a fluid is injected from an aperture into an open space, hydrodynamically confined by a surrounding fluid, and entirely re-aspirated into a second aperture. Various MFPs have been developed, and have been used for applications ranging from surface patterning of photoresists to local perfusion of organotypic tissue slices. However, the hydrodynamic and mass transfer properties of the flow under the MFP have not been analyzed, and the flow parameters are adjusted empirically. Here, we present an analytical model describing the key transport properties in MFP operation, including the dimensions of the hydrodynamic flow confinement (HFC) area, diffusion broadening, and shear stress as a function of: (i) probe geometry (ii) aspiration-to-injection flow rate ratio (iii) gap between MFP and substrate and (iv) reagent diffusivity. Analytical results and scaling laws were validated against numerical simulations and experimental results from published data. These results will be useful to guide future MFP design and operation, notably to control the MFP “brush stroke” while preserving shear-sensitive cells and tissues.

fashion 4 , maskless lithography by dissolving photoresist on a coated surface 10 , lysis of human breast cancer cells and collecting and analysis of mRNA released from the lysate 11 , pharmacokinetic studies on cells 12 , staining of tissue sections 13 ,analysis of neutrophil dynamics during chemotaxis 14,15 , and perfusion of organotypic brain slices 16 . In all these studies 17 , the MFP design, flow rates and gap height were optimized empirically and calibrated experimentally based on a trial and error method. Previous analyses in other probe geometries (microfluidic quadrupoles) have demonstrated that the steadiness and high reproducibility of MFP patterns are due to the presence of Hele-Shaw flow patterns at sufficiently small gap sizes. Under these conditions, the flow underneath the probe can be considered quasi two-dimensional and readily analyzed using potential flow theory drawing upon the strong analogy between 2D flow fields and electrostatic fields 9 . Christ and Turner have studied convection in various two-aperture MFP configurations 18 . Their hydrodynamic pressure measurements at several points underneath the probe, supported by numerical models and relevant pressure scales thoroughly confirm the validity of potential flow assumptions when modeling the flow behavior under the probe. Yet, their analysis does not provide an analytical framework to quantify neither advection nor diffusion in two-aperture MFPs. Advection and diffusion are essential parameters to operating MFPs and control diffusion gradient generation, shear stress at the surface, confinement area, and the effect of port geometry.
Here, we introduce the relevant physical theory describing the various design and operation parameters within a MFP. The two-aperture probe is analyzed as a flow dipole (or doublet) 19 using the same formalism as for 2D electrostatic dipoles. The model also goes further than the conventional Hele-Shaw formalism by also taking into account diffusive mass transfer. We extract several analytical results as well as scaling laws to control important parameters such as HFC zone dimensions, shear stress and diffusion gradients within the gap below the mesa. The resulting model was used to establish universal rules for systematic probe design and operation. The validity of the analysis was tested against numerical models and experimental data. Finally, the shear stress on the substrate was calculated for various operating conditions and a working range based on the tolerance of cells to shear stress was established.

Results
In their most general form, dipolar microfluidic probes comprise two apertures of dimension a separated by a distance d from center-to-center. The apertures are made into a flat mesa of side dimensions L 1 × L 2 , which is suspended above another flat surface, forming a gap G (Fig. 1). Fluid is injected through one of the apertures at a flow rate Q inj and aspirated back by the second aperture at a flow rate Q asp = αQ inj , where α is a parameter expressing the ratio between aspiration and injection flow rates (>1 for confinement to exist). Detailed experimental procedures for MFP operations are explained elsewhere 20,21 .
HFC data described here was obtained from previously published experimental results using a particular MFP geometry (d = 50 μ m, square apertures with a = 25 μ m, and mesa dimensions of (550 μ m × 447 μ m) 4 . The gap size and ratio of injected to aspirated flow rates were changed from one experiment to another.
The size of the HFC was determined experimentally by injecting a solute, namely fluorescently labelled IgG molecules under the MFP, which was atop a substrate where the IgGs readily adsorbed. The size of the adsorption area following a 10 s exposure was determined by fluorescent imaging as illustrated in Fig. 2.
Leakage-free HFC was achieved for flow ratios α ≥ 2.5 for the particular MFP described in Fig. 2. For lower values of α, the injected fluid flow pattern is sufficiently large to extend into the area outside the mesa and is no longer confined. The binding adsorption patterns of IgGs on the substrate were controlled by the flow ratio α and strongest for flow ratios between α = 2.5 to α = 4. For increasing flow ratios, the surface patterns gradually vanished, indicating that the injected IgGs were instead recaptured by the aspiration inlet without reaching the substrate surface.
The confinement as a function of the gap between the mesa of the MFP and the substrate was also assessed experimentally (Fig. 2B). For gaps varying from 4 μ m to 50 μ m and for flow ratios of 2.5 to 4, the HFC is relatively insensitive to the gap size. However, for smaller gap sizes and large flow ratios, a "wing effect" can occur, by which the diffusion gradient is located underneath the injection aperture (See Figs. 2B(7) and 2B (10)). For a flow ratio of 4, the size of the adsorbed pattern shrinks with increasing gap size. At large enough α values, the adsorption pattern disappeared completely when the gap became too large. Thus for a range of flow ratios, the surface pattern shrinks for increasing gap sizes, but at the same time the intensity of the fluorescence weakens reflecting that the local concentration of the molecules in solution near the substrate diminishes as the MFP is moved further away (Fig. 2B). Quantitative criteria have been derived, based on the theory presented below, to determine the critical values of α and G under which the adsorption pattern is stable (see Supplementary Information).
Theory: Formulation of the general MFP problem. The general configuration describing flow in a gap of height G right below a MFP is introduced schematically in Fig. 1E. When the gap G between the plates is sufficiently small, the flow profile between the plates can be described as a quasi two-dimensional flow, or Hele-Shaw flow 22,23 (See supplementary information), mathematically akin to Darcy's law 24 .  where η is the fluid dynamic viscosity and p(x, y) is the pressure distribution. In open microfluidics, the Hele-Shaw condition is generally met as it requires both creeping flow (Re ≪ 1) and that the flow profile varies smoothly over a distance L much larger than the gap G (specifically (L 2 ≫ G 2 )) where L is the characteristic length in (x, y) of the observed flow behavior. When Re ≪ 1, the superposition principle can be applied to add up the flow potential of the injection and aspiration apertures to compute a doublet 19 , a flow profile with streamlines mathematically identical to the field lines in a two-dimensional electrostatic dipole, or even more physically accurate, to the current density lines in a two-dimensional current dipole.
Analysis of the ideal MFP. Unlike electric charges, no flow apertures can be practically reduced to a point source of streamlines. Instead, they must be viewed as openings (round, square, or others) of characteristic dimension α. Yet, when the distance between the apertures is much larger than the apertures, the point source approximation is very effective in modeling multipolar flow behavior. We consider here the circular aperture as a general case. Since we are in the presence of a potential flow, Gauss's theorem 9 (Eq. (3)) can be readily used to readily calculate in polar coordinates the velocity profile outside a round aperture of thickness G and radius a/2: In polar coordinates, the mass balance equation provided above yields the simple result: The flow profile under the MFP can be calculated as a superposition of two monopoles of different flow rates located at ± d/2 from the origin on the x-axis as illustrated in Fig. 1E. When a ≪ d, the point source approximation is accurate, and the dimension a can be neglected. The full velocity profile under the probe is then given in Cartesian coordinates by: ( )

Stagnation point and hydrodynamic confinement limit
When α > 1, a point of zero velocity (stagnation point) occurs outside the injection aperture at a distance χ SP of the center of the probe as shown in Fig. 1F. The location of this stagnation point satisfies the condition: The latter condition on υ y is satisfied for all y due to the symmetry of the problem along the x-axis. The condition on υ x specifies the position of the stagnation point: This key result describes the maximum distance R from the probe center at which the injected analyte will travel before being recaptured by the aspiration inlet. In other words, injected analytes are confined within closed streamlines as represented in Fig. 1F. The definition of R is highly analogous to the confinement length defined earlier for microfluidic quadrupoles, in which α α = / ( + )/( − ) R d 2 1 1 9 . In both cases, when α = 1, injection and aspiration flow rates are equal and the stagnation point is found at infinity, just as in a regular dipoles and quadrupoles. In Fig. 3A,B, the numerical values of R are plotted against their theoretical value obtained in Eq. (7) as funtions of α and G. Coefficients of determination R 2 are listed in Table 1. As expected, the best match is obtained for the probes where the ratio of aperture size to interaperture distance (a/d) is minimum (R 2 = 0.99 for a/d = 1/6) and the worst for the largest a/d ratios (R 2 = 0.96 and 0.92 for a/d = ½). Thus the R 2 value provided in Table 1 can be viewed as a quality test of the point source analytical model when compared to the full 3D numerical solutions.
An exact expression for the width of the HFC underneath an ideal doublet can also be calculated using a similar method involving stream functions (see supplementary information).  Table 1 as well as the coefficient of determination R 2 for every model used. Scaling laws of diffusion broadening along the HFC area. In most practical situations, analytes within the probe will have a non-negligible diffusivity and will induce a concentration gradient on the outskirts of the HFC area. This broadening will yield significantly larger surface adsorption patterns than expected from the HFC alone. Furthermore, this diffusive broadening can be exploited, as described elsewhere, to generate rapidly tunable, highly reproducible sharp floating concentration gradients 9 . The general advection-diffusion equation of a diffusive species C(x, y) within under the MFP is described in cartesian coordinates by the equation: where  v is found in Eq. (5). While easily solved using finite element methods, the unwieldy Eq. (8) offers little physical insight in its current form. However, a scaling law describing the characteristic diffusion length of the HFC can be obtained by linearizing the velocity profile near the stagnation point using a Taylor series expansion of Eq. (5). Therefore: To extract the Péclet number near the stagnation point we make Eq. (8) and Eq. (9) dimensionless, by using the following scales: Thus, after simplifications, Eq. (8) can be expressed under its canonic form: Although Eq. (11) has separable variables, we leave to future work the determination of its full solution, which would need to take into account the curvature of the diffusion interface near the stagnation point, which itself depends on several operation variables. In a first approximation, analogous to the one performed on quadrupoles elsewhere 9 , we make Eq. (11) one dimensional by setting the y-components of the concentration derivatives to 0. In most cases, this approximation is highly justified since the interface radius of curvature is much larger than the diffusion length scale. Furthermore, the diffusion broadening region, in almost all cases in Fig. 2, display a uniform thickness even far from the SP, only decreasing in thickness when nearing the aspiration aperture. The approximation is further justified when using square apertures as they flatten the diffusion profile at the SP even at small curvature radii (e.g. Figs 2A (7) and 2B (11)). Finally, the symmetry plane along the y-axis implies that ∂ /∂ =  C y 0 near y = 0. When considering the steady state solution only and setting ∂ /∂ =  C y 0 2 2 , the simplified Eq. (11) with its boundary conditions becomes: The boundary conditions in Eq. (12) are valid only for high Peclet number, i.e. when the diffusion broadening area is thin an outside of the injection aperture, a condition typically met in probe operation. Solving Eq. (12) yields:  Fig. 3. Strong agreement was found between the simulation and theory for the stagnation point location, confinement area as well as diffusion length, as described by the R 2 values compiled in Table 1. Calculations of scaling laws for the confinement area has been explained in further details in the supplementary information.
Shear stress analysis. The MFP has been used to expose cells and tissues with various chemicals, and an important parameter that can influence the treatment outcome is shear stress. Shear stress plays a major role in many biological phenomena bringing about a number of cellular/intracellular events [26][27][28] . For example, endothelial cells and neutrophils transduce the applied shear stress stimulus into intracellular biochemical responses, by which they regulate the vessel structure 29,30 . Similarly, these stimuli also affect the orientation of osteoblasts and neurons resulting in cellular migration and matrix outgrowth 31,32 . An important application of the MFP is detaching shear sensitive cells without damaging the cell membrane, and thus calibration of shear stress is necessary.
In Hele-Shaw flows, the shear stress is directly proportional to the normalized velocity profile: The maximum value τ max is found where the maximum value of  ν  lies, i.e. at (x = −(1 − a/d), y = 0) near the inner edge of the aspiration aperture (Fig. 4A). To reflect this observation, we rescale τ max using the aspiration flow rate α Q inj as the new characteristic flow rate and using the aperture size a instead of probe length d as the new characteristic length. It ensues that when = / a a d 1 or α ≫ 1, the contribution of the injection aperture to the maximum velocity becomes insignificant, giving the simple scales τ ατ ≈ / a max 0 and τ ατ ≈ .
/ a 0 866 max 0 , respectively for the round and square aperture probes for the maximum shear rate under the probe (see supplementary information). Thus maximum shear stress is inversely proportional to gap G (~1/G 2 , Fig. 4B) and varies linearly with the aspiration flow rate Q asp . Hence, we expect the maximum shear stress to vary linearly with the flow ratio, as shown in Fig. 4C. Maximal shear stress is an important experimental parameter since in biological systems cells are affected Scientific RepoRts | 5:11943 | DOi: 10.1038/srep11943 by shear, and for example the growth rate of endothelial cell was shown to be susceptible to shear as low as 1 Pa 33 .
High shear stresses can cause damage to the cell membrane, and consequently, lead to cell death 34 . To account for this issue in MFP design and operation, the maximal shear stress at the bottom of the substrate is characterized first by varying the flow ratio and then by varying the gap size while keeping the other parameter constant. Combining the knowledge from former studies of cell detachment using low shear stress (around 0.5 Pa) 32 to detach PC12 neural cells and maximum shear stress that can be tolerated by neurons (around 1.5 Pa) 32 , we performed a set of simulations with different gaps and flow rate ratios. By considering both limits described above, we outlined the safe operational parameters of the MFP to work with neurons (Fig. 4D).
HFC contour and width. Another important application of the current MFP model is to allow precise control over the shape and dimensions of the probe writing head based on user-controlled parameters. When used in surface patterning mode, the dimensions of the HFC area will dictate the shape of the "brush stroke" that the MFP leaves when it writes on a surface. One can proceed to visualize the shape of the hydrodynamic flow confinement (HFC) area by determining the associated stream function 22 of the velocity profile. In general terms, the stream function is defined as: where the integration constant can be set to zero. As the streamlines define the trajectory of an infinitesimal particle in the flow field, the streamline passing by the stagnation point (2x SP /d = (α + 1)/(α − 1) , 0) will describe the envelope of the HFC (neglecting diffusivity). In Eq. (11), the streamline corresponding to this condition is x y x y arctan 1 arctan 1 1 2 18 Hence, Eq. (18) is the function describing the upper half (y > 0) of the HFC contour. While the variables cannot be separated, the shape of the implicit stream function is plotted in Fig. 5A forvarious values of α.
The HFC area or simply confinement area (C. A.), i.e. the area of the spot left by the probe in surface patterning mode, can be found by plotting Eq. (18) for a specific value of α and numerically integrating under the curve. One exception occurs when α = 2 where Eq. (18), once simplified, yields (2x/d − 1) 2 + (2y/d) 2 = 4, which is the equation of a circle of radius x = d centered at the injection aperture (x = d/2) (Fig. 5A). The confinement area of the HFC in this particular case is C.A. = πd 2 . Another exception occurs when α = 3, which also has a simple solution: C.A. =  6). In all cases, confinement area will be proportional to the square of the inter-aperture dimension d and will decrease with increasing flow ratio α. Diffusion is set to 0 everywhere in this section to highlight the contribution of the flow profile alone to the C.A.
Theoretical HFC Width. Perhaps the simplest and most generally applicable result of this paper is the positions of the probe's stagnation point R as described in Eq. (7). Another useful measurement is the HFC maximum width W, the maximum extent of the HFC in the dimension perpendicular to the probe's main axis. The HFC width as a function of flow rate has been studied numerically and experimentally elsewhere 18 , but without the insight provided by a complete mathematical derivation, which is provided here for the ideal case of point sources under the Hele-Shaw flow assumption.
Since the contour of the HFC can be described theoretically by the stream function described in Eq. (18), we can compute the local maximum of the HFC width in the y dimension by differentiating the stream function with respect to x such that Physically, this represents the point at the edge of the HFC at which the y component of the velocity is identically zero. Using the normalized coordinates, we get: This is the condition that must be respected when v y = 0 and y = y max = W/2, where W is the HFC width. Substituting Eq. (20) in Eq. (18) yields: This is the transcendental equation that analytically describes the HFC's normalized full width ). The width function α ( ) ∼ W is plotted in Fig. 5B. Note that Eq. (21) is a normalized function and that it is therefore valid for any MFP geometry where the point source approximation can be made.

Discussion
Applications to brush stroke analysis and probe operation control. One of the most important applications of the current model is to be able to control the shape and dimensions of the probe writing head by specifying user-defined parameters such as the injection and aspiration flow rates or the probe gap G. The contour function described in Eq. (18), combined with the HFC maximum radius R (Eq. (7)) and the HFC width W (Eq. (21)) completely defines the area and maximum dimensions of the MFPs writing area. To account for the diffusive broadening around the contour, Eq. (14) can also be applied. However, all these parameters specify the shape of the probe's writing under static conditions. When the probe is moving at constant velocity U in any given direction with respect to a fixed substrate, viscous forces will tend to deform the static flow profile. The fully developed Hele-Shaw profile will be a superposition of a parabolic flow profile stemming from the probe apertures and a Couette flow profile of mean velocity U/2 induced in the direction of the probe displacement. This flow behavior can be reduced to a purely parabolic flow profile as in Eq. (1) by placing oneself in a referential moving in the direction of the probe at the mean fluid velocity U/2. In doing so, the height-averaged flow profile found in Eq. (2) is kept unchanged and the flow in this reference frame is identical to that of Eq. (5) (see proof in supplementary information). Once solved, the velocity profile in the probe's moving frame, using the dimensionless scales described in Eq. where θ is the angle of the probe velocity vector with the probe axis and β = πGUd/2Q inj is the dimensionless probe velocity. When β > 1, the probe velocity disturbs the HFC contour significantly. Yet when β ≪ 1, the contour can be assumed not to be perturbed by the probe movement. Using the typical Scientific RepoRts | 5:11943 | DOi: 10.1038/srep11943 operation parameters previously described (Q inj = 0.44 nL/s, G = 10 μ m, d = 50 μ m), we find that the probe perturbation due to its velocity is negligible when U ≪ 2Q inj /πGd = 0.6 mm/s. This critical velocity can be made arbitrarily large without changing the HFC contour by increasing Q inj and Q asp proportionally while keeping the ratio α constant. Within this approximation of small β, we now define a "brush stroke" for the probe operation where the width of the stroke will depend on the user-controlled parameters, including the probe's velocity and displacement angle (Figs 5(C-E)). The brush width can be decomposed in its two principal axes: the width of the probe when moving parallel and perpendicular to the probe's axis, respectively W and ⊥ W , where W is calculated using Eq. (21) and The linear combination describing the brush stroke at low velocities for an arbitrary xy probe displacement is generally described as: It is important to point out that the diffusion broadening has not been taken into account when deriving Eq. (24) and to do so would require needlessly cumbersome mathematical derivations. However, the effect of diffusion can be simply taken into account by making the approximation that the diffusive broadening is small compared to the HFC radius and increases it uniformly by a distance Δ x, obtained from Eq. (14), on each side except near the aspiration aperture. Under this approximation, x D and Eq. (24) can still be used.

Conclusion
In this paper, we have provided a theoretical analysis of ideal dipolar probes (d ≫ a ≫ G) and tested the validity of our transport model under probes with finite aperture sizes and arbitrary geometries. The analysis is then validated with data obtained from simulations and previously published experimental data. The general set of analytical results and scaling laws are summarized in Table 2. The simple results established will allow precise control of HFC dimensions and confinement area, diffusion length and gradient properties at the HFC's edges and shear stress control with respect to flow rate ratio α and gap size G. They will provide precise design criteria (d,a, and mesa dimensions, during fabrication process) and control the probe's brush stroke and velocity during operation (α, G, Q inj , D). For example, if the user seeks a circular HFC of area A at the tip, the inter aperture distance d to obtain this exact area can readily be calculated (in this special case, when α = 2, the HFC area is precisely A = πd 2 , see supplementary information). Once the value of d is found, the dimensions of the mesa need only to be a few times larger than the maximal R value of the HFC (with R = 3d/2 when α = 2). Then, using the formula for the HFC outer edge (Table 2), the precise shape of the HFC could then be calculated for any value of α for  Table 2. Design criteria, formula and validity domain of two apertures MFP. This parameters and related formula provide substantial guide to select best experimental conditions for cell study applications.
this probe. Finally, the outer edge diffusion layer can be adjusted during the experiment by modulating the gap size G while maintaining the HFC's circular shape constant. In the previous example, we can also determine in advance whether we need square/round, or large/small apertures. Small apertures will be preferred in situations where α is large to prevent the "wing effect" (Fig. 2B-10) that occurs when the stagnation point is found inside the aperture. Large apertures are, on the other hand, preferable when the shear stress needs to be minimal (Fig. 4B,C).
We also provide a complete characterization of the MFP's "brush stroke" when operating in surface patterning mode. The analysis suggests that the simplest way to precisely control the width of the stroke is to move the MFP either parallel or perpendicular to its axis since the width of the stroke ( , ⊥ / / W W ) is then completely independent of the flow velocity and given a simple formula ( Table 2) that can be easily modified to account for diffusion by adding the diffusion broadening scale (DL, Table 2) to the width.
The results presented in Table 2 can also be used to indirectly measure and calibrate probe height above a surface by looking at diffusion length with respect to G (∝ √G) or shear stress as a function of the applied flow rates. Finally, the general theoretical framework provided here, and the set of analytical results that stem from it, enables a more precise control for MFPs of any shape or form in a variety of applications including neuron migration 35 , stem cell differentiation 36,37 and shear-dependent cell poration 38,39 without using costly, complicated trial and error experimental method.

Methods
Numerical Modeling. The numerical models were established using the finite element method with commercial software (COMSOL Multiphysics 3.5a, USA) on a computer with 8 cores, 64-bit CPU and 26 GB RAM. The computational model was built by combining Stokes flow and steady state convection-diffusion mass transfer equations. Boundary conditions were assumed to be "no-slip" conditions at solid surfaces and "open boundaries" at atmospheric pressure at the perimeter of the MFP. To compute the mass transfer model, the MFP's rigid walls and substrate surface are considered as insulators and the flux was set to zero. The concentration at the perimeter was set to zero reflecting the assumption that under HFC no solute reaches the edge, and that the large volume constitutes a sink. Concentration at the injected aperture was normalized and set to 1 to provide a relative scale. At the aspiration aperture, diffusion was neglected owing to the local dominance of convective mass transport.
Numerical calculations of the HFC were performed and the calculated concentration as a function of position at the substrate surface was reported as function of flow ratio in Fig. 2A, and as function of the gap size, Fig. 2B and Fig. 2B. There is a good visual agreement between experimental and numerical results, which was also quantified by measuring the surface area for each case, and comparing it for different flow conditions (Fig. S1). The observed discrepancy between numerical and experimental results can be mainly attributed to the indirect measurement of protein adsorption at the surface rather than the direct measurement of the flow profile. Uncertainty on the experimental gap, exposure time, as well as vibrations of the MFP while running the experiment will cause the protein footprint on the surface to spread. As a consequence, the pattern on the surface is not a record of the average pattern, but of the maximal extent of the flow profile that may occur due to various fluctuations during the 10 s of streaming across the surface. Nevertheless, the numerical models are in good agreement with the adsorbed protein experiments. To extend the validity of our models to any probe geometry and to extract useful design criteria based on scaling laws (see section "Theory"), numerical simulations of the MFP with various geometries and gaps, and under different flow rates were performed using full 3D finite-element models.
To test the limits of our advection/diffusion models, we have proceeded to perform the following series of simulations emulating probe geometries already used in the literature, Table S1. Water was used as fluid (solvent) with a density of 998.2 kg/m 3 and dynamic viscosity η = 0.001 N·m/s (at 20 °C). Immunoglobulin G (IgG) was used as solute and a diffusion coefficient in water of D = 4 × 10 −11 m 2 /s was used 7 . The flow regime was estimated prior to the simulations to help select the appropriate models. The Reynolds number in all regimes was on the order of Re ≈ O(10 −2 ) ≪ 1 indicating creeping flow conditions to be treated using Stokes formalism 22 . Similarly, the Peclet number representing the ratio between convective and diffusive mass transport was found to be high (Pe ≫ 1) 40 for all MFP operation regimes, indicating that convective mass transport dominates diffusion except very close to the stagnation point.