An Approach To Extract Rate Constants from Reaction−Diffusion

Anal. Chem. 2005, 77, 3417-3424

An Approach To Extract Rate Constants from Reaction-Diffusion Dynamics in a Microchannel Jean-Baptiste Salmon,* Claire Dubrocq, and Patrick Tabeling

Microfluidique, Mems et Nanostructures, ESPCI, 10 rue Vauquelin, 75005 Paris, France Sandrine Charier, Damien Alcor, and Ludovic Jullien

CNRS UMR 8640 De´ partement de Chimie, EÄ cole Normale Supe´ rieure, 24 rue Lhomond, 75231 Paris Cedex 5, France Fabien Ferrage

New York Structural Biology Center, 89 Convent Avenue, New York, New York 10027

A theoretical model is proposed to extract rate constants of second-order chemical reactions down to the millisecond time scale from the observation of reactiondiffusion processes in a microchannel. We validate this theoretical approach by examining an appropriate model reaction. The measured rate constant is in excellent agreement with this obtained from nuclear magnetic resonance experiments.

The analysis of the mechanisms of chemical reactions as well as the determination of dynamical parameters, such as rate constants or diffusion coefficients, is essential in chemical engineering, chemistry, and biology. Several strategies have been developed to address the corresponding issues. In relaxational approaches, the kinetics of the reaction is extracted from investigating the relaxation of the chemical system toward the thermodynamic equilibrium, after small induced perturbations1,2 or spontaneous fluctuations.3,4 Excellent temporal resolution is achieved using such methods, since they do not require mixing the chemical reactants. However, exploring the dynamics of the system only close to the chemical equilibrium * To whom correspondence should be addressed. Current address: LOF, CNRS-Rhodia, 178 avenue du docteur Schweitzer, 33608 Pessac cedex, France. Fax: 0033556464790. E-mail: [email protected]. (1) Eigen, M., de Mayer, L., Eds. Relaxation Methods in Techniques of Organic Chemistry, 2nd ed.; Interscience Publishers; John Wiley and Sons: New York, 1963. (2) Sandstro ¨m, J. Dynamic NMR Spectroscopy; Academic Press: London, 1982. (3) Magde, D.; Elson, E.; Webb, W. W. Phys. Rev. Lett. 1972, 29, 705. (4) Krichevsky, O.; Bonnet, G. Rep. Prog. Phys. 2002, 65, 251. 10.1021/ac0500838 CCC: $30.25 Published on Web 04/26/2005

© 2005 American Chemical Society

may hinder mechanism discrimination. The approaches that rely on mixing reactants do not suffer from the latter drawback: the relaxation from far-from-equilibrium conditions can be investigated. Unfortunately, the temporal resolution is now strongly limited by the mixing time of the reactants. The mixing step can easily be achieved down to 1 s with classical means, and down to a few milliseconds in conventional stopped-flow devices, where the mixing occurs via turbulence. However, these approaches may be limited either by the sensitivity of the detection method or by the large amounts of samples required to reach sufficient flow rates (mL/s) to induce turbulence.5,6 Recently, some groups have developed continuous approaches using microfluidics to overcome these drawbacks7,8 (see ref 9 for a review on microfluidics). More specifically, the reactants are continuously injected in a microchannel, then mixed, and the concentrations of the products/reactants are measured at different positions downstream. Since the chemical species are convected in the microchannel, one can perform stationary measurements of the chemical kinetics along the flow. Microfluidics is essential to develop such techniques, because large velocities are required to access fast chemical rates, and one has to avoid turbulence (laminar flows with velocities as large as 10 (5) Pilling, M. J.; Seakins, P. W. Reaction Kinetics; Oxford University Press: Oxford, U.K., 2001. (6) Steinfeld, J. I.; Francisco, J. S.; Hase, W. L. Chemical Kinetics and Dynamics, 2nd ed.; Prentice Hall: Upper Saddle River, NJ, 1999. (7) Pabit, A.; Hagen, S. Biophys. J. 2002, 83, 2872. (8) Song, H.; Ismagilov, R. F. J. Am. Chem. Soc. 2003, 125, 14613. (9) Stone, H. A.; Stroock, A. D.; Ajdari, A. Annu. Rev. Fluid Mech. 2004, 36, 381.

Analytical Chemistry, Vol. 77, No. 11, June 1, 2005 3417

Figure 1. (a) Schematic setup: reactants A and B respectively injected at a constant flow rate into Y-shaped microdevice. The product of the reaction appears in the diffusion cone of A and B downstream. The typical heights and widths of the main channel used in our experiments are d ) 20 µm and w ) 230 µm. (b) Top view of an experimental illustration using fluorescence microscopy. Concentration map of the fluorescent acidic state (C) resulting from the reaction and diffusion of the nonfluorescent solutions of protons (A) and of the basic state (B) of a pH indicator both diluted in acetonitrile (arbitrary units). The flow rate inside the main channel is Q ) 4 µL/min, and the initial concentrations are a0 ) 1.02 × 10-1 M for the proton (region Y < 0) and b0 ) 1.1 × 10-3 M for the basic state of the pH indicator (region Y > 0). The origin of the Y-axis and of the X-axis is set at the Y-junction.

cm/s are easily obtained in microdevices). Moreover, small amounts of samples are consumed in such experiments since the typical flow rates are on the order of a few microliters per minute. However, quickly mixing the reactants at these small length scales is a challenge since it occurs only by molecular diffusion.9 Various ingeneous techniques have been developed to overcome this difficulty,10-13 and some proved to be successful in measuring chemical characteristic times as short as a few milliseconds.7,8 Unfortunately, the mixing time in such devices cannot be decreased indefinitely, and a complete description of the flow during the mixing step is needed to access faster kinetics. This point has encouraged other groups to follow a rather different approach.14,15 This approach is illustrated in Figure 1 in the specific case of the following reaction that is frequently encountered in chemistry and biology:

A+BhC

(1)

and the product of the reaction, C (the fluorescent acidic state of a pH indicator in our case), appears in the diffusion cone of A and B. The evolution along the flow of the concentrations of the reactants and the product contains all the dynamic information related to the diffusion and reaction processes. The main problem is to extract this information in a relevent and easy way. In the absence of chemical reaction, appropriate models are available for extracting the diffusion coefficients from the concentration maps, as demonstrated by Kamholz et al. in the case of various fluorophores and fluorescently labeled proteins.14,16 However, the treatment has to be revisited and extended in the presence of a chemical reaction such as (1) (vide infra).14,15 The present work presents a convergent method to extract rate constants of second-order reactions down to the millisecond time scale from the analysis of the reaction-diffusion process in a microchannel, without any computational fluid dynamics (CFD) tools. The paper is organized as follows. The first part revisits a reaction-diffusion model of a second-order reaction in a microchannel. After certain assumptions and using dimensionless variables, one obtains equations that only depend on the thermodynamics of the reaction. We then reveal two regimes of the solutions of these equations. In particular, we show that the rate constant can only be extracted from a limited zone of the concentration maps, whose extension is given by the characteristic time of the reaction and the mean velocity. Eventually, we briefly discuss the range of validity of this approach. The second part presents the experimental data obtained in the case of a proton exchange reaction. We analyze these data with our model and show that the measured rate constants are in line with those obtained by nuclear magnetic resonance (NMR) relaxation experiments. THEORY Two-Dimensional Model. In the present work, we only consider diluted solutions. The relevant dynamical parameters, such as the diffusion coefficients and the rate constants, are thus constant and there is no coupling between the flow and the chemical reaction. Therefore, the following equations can be used to describe the reaction-diffusion process occurring in the setup drawn in Figure 1, in the case of reaction 1:

(v‚∇)A ) - kAB + k′C + Da∆A (v‚∇)B ) - kAB + k′C + Db∆B (v‚∇)C ) + kAB - k′C + Dc∆C

The reactants A and B (in our experiments, the proton and the nonfluorescent basic state of a pH indicator) are separatly injected at the same constant flow rate in the two arms of a Y-shaped microdevice. They mix by diffusion in the main channel (10) Knight, J. B.; Vishwanath, A.; Brody, J. P.; Austin, R. H. Phys. Rev. Lett. 1998, 80, 3863. (11) Johnson, T. J.; Ross, D.; Locascio, L. E. Anal. Chem. 2002, 74, 45. (12) Stroock, A. D.; Dertinger, S. K. W.; Ajdari, A.; Meziæ, I.; Stone, H. A.; Whitesides, G. M. Science 2002, 295, 647. (13) Song, H.; Tice, J. D.; Ismagilov, R. F. Angew. Chem., Int. Ed. 2003, 42, 767. (14) Kamholz, A. E.; Weigl, B. H.; Finlayson, B. A.; Yager, P. Anal. Chem. 1999, 71, 5340. (15) Baroud, C. N.; Okkels, F.; Me´ne´trier, L.; Tabeling, P. Phys. Rev. E 2003, 67, 060104.

3418

Analytical Chemistry, Vol. 77, No. 11, June 1, 2005

(2)

Da, Db, and Dc, and A, B, C respectively designate the diffusion coefficients and the concentrations of the three species A, B, and C. v is the velocity field in the main channel displayed in Figure 1, and k and k′ are the rate constants associated to the forward and backward reactions (1). In the most general case, it is rather difficult to extract the rate constants from the comparison of the solutions of eqs 2 to experimental data such as those displayed in Figure 1b. Indeed, one has first to compute the velocity profile v(X, Y, Z) for the geometry displayed in Figure 1a; second to compute the concen(16) Kamholz, A. E.; Schilling, E. A.; Yager, P. Biophys. J. 2001, 80, 1967.

tration fields by solving eqs 2 using CFD tools for a given set of parameters (k, k′, Da, Db, Dc); and finally to compare the results of the simulation to the experimental concentration maps. Such a fitting procedure is laborious to perform since one simulation is required to evaluate each value of the rate constants. The fitting process may be advantageously simplified by taking into account the specific geometry displayed in Figure 1a. Indeed, in a microchannel with a high aspect ratio (w/d . 1), the flow is mainly uniform along the Y-direction and Poiseuillelike along the Z-direction.17 Moreover, if the height d of the channel is small enough, the concentrations are homogeneous along the Z-direction, and a two-dimensional (2d) model may correctly describe the reaction-diffusion process. These assumptions are discussed in the following. In this case, eqs 2 become

vj∂XA ) - kAB + k′C + Da∂Y2A vj∂XB ) - kAB + k′C + Db∂Y2B vj∂XC ) + kAB - k′C + Db∂Y2C

(3)

where vj is the mean velocity in the channel. To derive eqs 3, we assume that the convective transport is more efficient that the diffusive transport along the X-direction. In microfluidics, this approximation is frequently valid because the Peclet number defined by Pe ) vjd/D is generally greater than 1 (D is the typical diffusion coefficient of the reactants).14-20 In the present series of experiments, Pe ranges between 100 and 1000. Moreover, we choose Db ≈ Dc for the sake of simplicity. This assumption is valid in the experiments discussed later because the protonation of a base does not considerably affect its diffusion coefficient. Moreover, such an assumption does not fundamentally modify the results described below.21 Using the following nondimensionalized variables21

a ) A/a0, x)

c ) C/xa0b0

b ) B/b0,

kxa0b0 X, vj

y)

xx k

a0b0 Y DaDb

(4)

∂xa ) - (ab - βΓc)/β + χ∂y2a ∂xb ) - (ab - βΓc) β + ∂y2b/χ (6)

where

β ) xa0/b0,

χ ) xDa/Db,

and Γ ) Ke/a0

(7)

Ke designates the dimensionless equilibrium constant of the reaction that numerically identifies to k′/k upon assuming the solution to be ideal and by taking 1 M as the standard concentra(17) Kamholz, A. E.; Yager, P. Biophys. J. 2001, 80, 155.

tion (the reference system is an infinitely diluted solution). With such nondimensionalized variables, eqs 6 involve only three parameters (β, χ, Γ) and do not depend anymore independently on k and k′, but only on the ratio k′/k. Equations 6 explicitely involve the equilibrium properties of the reaction only. Actually, the kinetics of the process is now only contained in the dimensionalized solutions through eqs 5. As we will show later, the fitting process is considerably simplified in contrast to the general case described by eqs 2. Indeed, only a unique simulation of eqs 6 is ultimately needed to extract the kinetics of the reaction for a given set of parameters β, χ, and Γ. Note that there should be no length scales present to use these nondimensionalized equations. Therefore, it is essential that the concentration fields at x ) 0 are given by

a(0, y) ) 1 (b(0, y) ) 0 resp)

if

y0

(8)

(5)

where a0 (b0, respectively) is the initial concentration at the entrance of the main channel of A (B, respectively), one considerably simplifies the model since eqs 3 read

∂xc ) (ab - βΓc) + ∂y2c/χ

Figure 2. (a) Concentration map of c inferred from eqs 6 with χ ) 1.5, β ) 10, and Γ ) 0.001. The initial conditions at x ) 0 are given by eqs 8. (b) Map of ab - βΓc. Null values correspond to the local chemical equilibrium. (c) Concentration map of c vs yr ) (y - x1/2)/x1/2.

and that the velocity profile is fully developed at x ) 0 (the origin of the y-axis and of the x-axis is set at the Y-junction). Such assumptions will be addressed later, when entrance effects will be discussed. In the following, we briefly describe the solutions of eqs 6 to evidence the regime where the kinetic information can be extracted. Figure 2a displays the calculated concentration field of the product of the reaction with the initial conditions described by eqs 8 and for χ ) 1.5, β ) 10, and Γ ) 0.001. These parameters correspond roughly to our experiments. As shown in Figure 2a, the reaction-diffusion process is not symmetric since the maximum of the production of c occurs in the region of the less diffusive and less concentrated species. The asymmetry depends on the values of χ and β (see refs 21 and 22 for more details). (18) Beard, D. A.; J. Appl. Phys. 2001, 89, 4667. (19) Dorfman, K. D.; H. Brenner J. Appl. Phys. 2001, 90, 6553. (20) Ismagilov, R. F. Stroock, A. D.; Kenis, P. J. A.; Whitesides, G.; Stone, H. A. Appl. Phys. Lett. 2000, 76, 2376. (21) Ga´lfi, L.; Ra´cz, Z. Phys. Rev. A 1988, 38, 3151. (22) Taitelbaum, H.; Y.-E. Koo, L.; Havlin, S.; Kopelman, R.; Weiss, G. H. Phys. Rev. A 1992, 46, 2151.

Analytical Chemistry, Vol. 77, No. 11, June 1, 2005

3419

Figure 2b displays the map of ab - βΓc that measures the distance to the condition of the local chemical equilibrium: ab ) βΓc (see eqs 4 and 7). The local chemical equilibrium is almost reached downstream x ≈ 1 that corresponds in real units to X ) vj/(k(a0b0)1/2) ) vjτr, where τr is a characteristic time of the reaction. This map thus evidences two regimes: (i) for x j 1, the concentration field evolves toward the chemical equilibrium and is therefore controlled by the reaction and the diffusion processes; (ii) for x J 1, the system is close to the equilibrium and thus only evolves through diffusive processes. To evidence this last point, Figure 2c shows the concentration map of c versus yr ) (y - x1/2)/x1/2. Using this normalized unit, one easily observes that c(x, y) converges for x J1 toward a function F that only depends on (y - x1/2)/x1/2, i.e., c(x, y) ≈ F[(y - x1/2)/x1/2] for x > 1. The function F thus depends on Γ, χ, and β and does not depend anymore on k as can be seen using eq 5. This feature indicates that, above x J 1, the concentration field only depends on diffusive processes, since both the width and the position of the maximal concentration of the product evolves along the flow as x1/2. As a conclusion, this model reveals two regimes: if x j 1 (i.e., X j vjτr), the concentration map of the product contains both the reactant diffusion and the kinetics of the reaction 1. In contrast, it is only governed by diffusion beyond x ≈ 1 (i.e., X ≈ vjτr): no kinetic information can be extracted in this spatial domain. Examination of the Assumptions of the Model. Two different kinds of assumptions were made to derive the previous model. We first assume that the velocity profile at the entrance of the microchannel is fully developed and that eqs 8 correspond therefore to the initial concentrations at the entrance of the microdevice. In fact, the velocity profile has a complex shape at the entrance of the main channel and is fully developed only beyond an entrance length Le.16,18,19 At low Reynolds numbers (as is the case in microfluidics), Le only depends on the geometry and Le ≈ w. Therefore, eqs 6 and 8 are only appropriate to investigate chemical reactions associated with characteristic times τr larger than the residence time in the entrance zone; i.e., vjτr > Le. For d ) 20 µm, w ) 230 µm, and a flow rate of 15 µL/min, vjτr > Le ≈ w gives τr > 4 ms: chemical kinetics faster than 4 ms cannot be studied within the framework of the model described by eqs 6 with such geometry and flow rate. One could increase the flow rate or decrease the channel width w to overcome this limitation. One could also envisage using CFD tools to study the dependence of Le on the geometry of the junction to decrease the prefactor of the relation Le ≈ w. To our knowledge, such work has not been done for Y-shaped microdevices. The second assumption relies on the 2d approach to describe the process of reaction-diffusion in a microchannel. As recently shown by Ismagilov et al., 3d structures of the concentration fields along the Z-direction of the microchannel may appear due to the Poiseuille velocity profile.20 Actually, the concentration field of the diffusing molecules is homogeneous over the height d only for distances X J vjd2/D ) Ped, where D is their diffusion coefficient. As a consequence, the 2d assumption is only valid for X J Ped. Thus, one should have vjτr > Ped, or τr > d2/D to adopt a 2d model. For typical values d ≈ 20 µm and D ≈ 10-9 m2/s, this condition yields τr > 400 ms, which should restrict this technique to rather slow kinetics. However, most imaging techniques (e.g., fluores3420 Analytical Chemistry, Vol. 77, No. 11, June 1, 2005

Figure 3. Protonation reaction of the basic state B of the fluorescent pH indicator that was used in the present series of experiments. The asterisk indicates the carbon atom bearing the proton, the signal of which was used to measure the rate constants by 1H NMR.

cence microscopy in the present work) provide concentration maps averaged over the height of the microchannel.14,15 We consequently compared the solutions of the 2d model and the solutions of a 3d model, when averaged over the height of the microchannel. The details of this work are given in the Supporting Information. The main result is that the difference between the two solutions is relatively small and that it is valid to use a 2d model to describe the experimental concentration maps of the reaction-diffusion process if they correspond to averaged values over the height of the channel. Quantitative measurements of the error originating from such an assumption are given in Supporting Information. EXPERIMENTAL SECTION Presentation of the Studied Chemical System. A proton exchange reaction was chosen to evaluate the present theoretical model since its dynamics was extensively covered at the relevant time scale.23 In a nonaqueous solvent such as acetonitrile (spectroscopic grade, 99.5%; Aldrich), where solvent ionization is negligible 28 this category of reactions obeys a simple one-step mechanism (1) for most systems encountered in chemistry. Photophysical and kinetic constraints were taken into account to design the investigated system. One state of the fluorescent pH indicator (acidic or basic) has to be fluorescent whereas the other is not. In addition, the reaction of proton exchange has to be rather slow to examine its kinetics with the present setup. We consequently adjusted the steric hindrance around the basic site to get appropriate rate constants. In fact, proton exchange reactions are often limited only by diffusion,23 and this often leads to very fast relaxation times of the reaction down to the nanosecond range. Figure 3 displays the coumarin (retained as a pH indicator) and its proton exchange reaction equivalent to (1), after identification of the proton to A and the basic and acidic states to B and C, respectively. The coumarin B [6,7-dimethoxy-4-(4-hydroxy-2,2,6,6tetramethylpiperidinomethyl)coumarin] was synthesized according to the following protocol. 4-Bromomethyl-6,7-dimethoxycoumarin (239 mg, 0.8 mmol, Aldrich) and 4-hydroxy-2,2,6,6tetramethylpiperidine (628 mg, 4.0 mmol, 5 equiv, Aldrich) in acetonitrile (20 mL) were stirred at 80 °C for 3 days. After evaporation of the solvent, water was added to the residue to eliminate the excess of starting amine, and the remaining solid was dissolved in methylene chloride. The organic phase was washed twice with water, dried over anhydrous sodium sulfate, and evaporated. Thus, the coumarin B was obtained as a yellow powder (275 mg; 92%) that was recrystallized in toluene (175 mg; 58%): mp 188 °C; 1H NMR (250 MHz, CDCl3) δ (ppm) 6.99 (s, (23) Eigen, M.; Angew. Chem. 1963, 75, 489.

Figure 4. Photophysical properties of the coumarin at 298 K in acetonitrile. Thick lines: single-photon absorption. Thin lines: fluorescence emission after one-photon excitation at λexc ) 348 nm (normalized by the respective quantum yields of fluorescence emission). The dashed lines correspond to the acidic state C, and the solid lines correspond to the basic state B. The figure also displays the band-passes of the filters used to record the images by fluorescence microscopy.

1H), 6.89 (s, 1H), 6.87 (s, 1H), 4.15-4.00 (m, 1H), 3.99 (s, 3H), 3.96 (s, 3H), 3.81 (s, 2H), 1.95 (dd, 2J ) 12.2 Hz, 3J ) 3.9 Hz, 2H), 1.51 (dd, 2J ) 3J ) 11.7 Hz, 2H), 1.18 (s, 6H), 1.01 (s, 6H); 13C NMR (62.9 MHz, CDCl ) δ (ppm) 162.4, 157.8, 152.6, 149.4, 3 146.1, 111.8, 110.8, 103.8, 100.4, 63.7, 56.8, 56.4, 56.3, 49.9, 43.9, 32.8, 22.8; MS (CI, NH3) m/z (%) ) 377 (49), 376 (100) [M + H]+, 374 (17), 221 (26), 158 (16); elemental analysis calcd (%) for C21H29NO5 (375.46 g‚mol-1): C, 67.18; H, 7.78; N, 3.73, Found: C, 67.07; H, 7.87; N, 3.61. Photophysical Properties. UV/visible absorption spectra were recorded on a Kontron Uvikon-940 spectrophotometer. Corrected fluorescence spectra were obtained from a Photon Technology International LPS 220 spectrofluorometer. The concentrations of the solutions for fluorescence measurements were adjusted so that the absorption was below 0.15 at the excitation wavelength. Figure 4 displays the absorption and emission spectra of the acidic and basic states of coumarin in acetonitrile. The present coumarin displays a large one-photon absorption at λmax ) 352 nm for the acidic state ( = 1.1 × 104 M-1 cm-1) and at λmax ) 339 nm for the basic state ( = 1.1 × 104 M-1 cm-1). The basic state (B) of the coumarin is not fluorescent (quantum yield of fluorescence ΦF < 0.01 with a maximum emission at 443 nm). In contrast, the acidic state (C) is fluorescent (quantum yield of fluorescence ΦF ) 0.06 ( 0.01 with a maximum emission at 450 nm). The overall fluorescence quantum yields ΦF were calculated using

ΦF ) Φref

( )

Aref(λexc) I n A(λexc) Iref nref

2

(9)

where the subscript ref stands for standard samples, A(λexc) is the absorbance at the excitation wavelength, I is the integrated emission spectrum, and n is the refractive index for the solvent. The uncertainty in the experimental value of ΦF was estimated

to be (10%. The standard fluorophore for the quantum yields measurements was quinine sulfate in 0.1 M H2SO4 with Φref ) 0.55.24 Measurement of the Protonation Constant Ke of the Coumarin. The protonation constant Ke was obtained by following the evolution of the photophysical properties of a 1.1 × 10-4 M solution of the basic coumarin in acetonitrile, upon addition of an acetonitrile solution containing the basic coumarin (total concentration 1.1 × 10-4 M) and trifluoromethanesulfonic acid (CF3SO3H, triflic acid, total concentration 2.2 × 10-3 M) (see Figure S-4, Supporting Information). The emission spectra IFλexc([H+]) were analyzed following ref 25 to determine Ke. We found Ke ) (2.0 ( 0.1) × 10-5 at 298 K (reference system, infinitely diluted solution; standard state, 1 M). Determination of the Rate Constants k and k′ for the Protonation Reaction of the Coumarin Using NMR. In a comparison with the results obtained by the present approach, we also measured the rate constants of reaction 1 by 1H NMR. The measurements of the rate constants at 298 K were carried out on an Avance Bruker DRX 600 spectrometer equipped with a TXI probe with three orthogonal pulsed field gradients (PFG). Sample temperature was controlled within 0.1 K. Three samples 6.4 × 10-3 M in coumarin were prepared in perdeuterated acetonitrile. In the first, no other solute was added so that the basic state of the coumarin was predominant (acidic state not detectable on a proton spectrum). In the second sample, a large excess of triflic acid was added, so that the acidic state was predominant (basic state beyond detection). The third sample was prepared with a concentration of triflic acid of 2.5 × 10-3 M. In the latter sample, the ratio between the respective signals of the basic and acidic states was 1.53 ( 0.03. The pulse sequence employed is a variation of the double PFG spin-echo nuclear Overhauser effect26 in which only one echo was performed. The selective 180° pulse has a Gaussian-shaped amplitude and a duration of 35 ms. All gradient pulses have a 1-ms duration and a sine shape to minimize eddy currents. Each spectrum was acquired with 512 scans. The aromatic proton indicated by an asterisk in Figure 3 was used for the analysis due to the relative isolation of each of these with respect to all the other protons: in the absence of any vicinal proton, no scalar-coupling effect may interfere with the excitation scheme and no nuclear Overhauser effect is expected during the mixing time. Two series of experiments were performed, each starting with a selection of one of the signals of the chosen proton. For each series, the following mixing times (in ms) were employed: 3, 50, 100, 200, 300, 400, 600, 800, and 1000, with a repeat of the 3, 200, 400, and 1000 time points. The longitudinal relaxation rates of each species were measured in the first two samples. The same sequence was employed, but the selective pulse had a duration of 10 ms. Each spectrum was the sum of 64 scans. For both samples, the relaxation delays employed were as follows: 3, 354, 747, 1192, 1706, 2313, 3057, and 4016 ms, with a repeat of the 3and 747-ms time points. For each spectrum, the intensities of each (24) Demas, J. N.; Crosby, G. A. J. Phys. Chem. 1971, 75, 991. (25) Jullien, L.; Canceill, J.; Valeur, B.; Bardez, F.; Lefe`vre, J.-P.; Lehn, J.-M.; Marchi-Artzner, V.; Pansu, R. J. Am. Chem. Soc. 1996, 118, 5432. (26) Stott, K.; Stonehouse, J.; Keeler, J.; T-.L. Hwang, Shaka, A. J. J. Am. Chem. Soc. 1995, 117, 4199.

Analytical Chemistry, Vol. 77, No. 11, June 1, 2005

3421

signal were sampled over three points around the maximum of the peaks. The relaxation decays were fitted using the program CurveFit.29 The longitudinal relaxation rates R1 of the investigated proton is R1(ac) ) 0.350 ( 0.008 s-1 in the acidic state, and R1(bs) ) 0.248 ( 0.008 s-1 in the basic state. The difference between the rates may be interpreted by the increased number of protons or additional impurities introduced with the triflic acid solution. Even if this difference is intrinsic, it is negligible with respect to the expected rate constants. For the exchange experiments, the value obtained for the signal appearing through chemical exchange was divided by the one of the initially selected signal. The kinetic rate data were then processed so that the ratio between the two peaks obtained in each spectrum was the only value manipulated in the data fitting procedure. This permitted us to make the method particularly robust with respect to long-term variations of the experimental conditions, such as the magnetic field inhomogeneity, as well as artifacts, such as imperfections in the selection part of the pulse sequence. The quantity obtained does not depend on longitudinal relaxation rates, as long as the difference between the rates of the two species is negligible with respect to the kinetic rates. The normalized values were plotted against time and fitted by the function

1 - e-(ka+kb)t f (t) ) ka kb + kae-(ka+kb)t

(10)

where ka and kb are the effective rate constants in the apparent process: ka

Ps y\ z Pg k b

(11)

where Ps and Pg are the polarizations of the selected and growing signal, respectively. Equation 11 was obtained by the direct integration of the differential equation system for Ps and Pg without taking into account the longitudinal relaxation. When the signal of the basic form was initially selected, the rate constants measured were as follows (see Figure S-5a, Supporting Information): ka ) 2.29 ( 0.06 s-1 and kb ) 4.08 ( 0.14 s-1; kb/ka ) 1.78 ( 0.03. Taking into account that ka ) k [H+] and kb ) k′, we found k ) 2.1 ( 0.1 × 105 (M s)-1 and k′ ) 4.1 ( 0.1 s-1. When the signal from the acidic form was selected, the measured rate constants were ka ) 3.95 ( 0.18 s-1 and kb ) 2.62 ( 0.15 s-1; ka/kb ) 1.51 ( 0.16 (see Figure S-5b, Supporting Information). With ka ) k′ and kb ) k[H+]: k ) 2.0 ( 0.1 × 105 (M s)-1 and k′ ) 4.0 ( 0.1 s-1. The two sets of data were in good agreement, and the ratios of the rate constants agreed with the integral ratio between the peaks of the basic and acidic forms in a regular proton spectrum: 1.53 ( 0.03. Microfluidics Experiments. Microfluidic chips were fabricated in poly(dimethylsiloxane) (PDMS) using standard soft photolithography techniques.27 At the time scale of our experiments, we checked that acetonitrile did not make our PDMS swell (27) Mcdonald, J. C.; Whitesides, G. M. Acc. Chem. Res. 2002, 35, 491. (28) When solvent ionization is significant (for instance water), proton exchange may be catalyzed by the solvent. This may considerably complicate the analysis of the dynamics of proton exchange.23

3422

Analytical Chemistry, Vol. 77, No. 11, June 1, 2005

and that the reactants did not diffuse into the PDMS. Several different channel geometries were tested during the experiments, but we only present here results concerning one specific geometry. The systems consist of Y-shaped microdevices, as depicted in Figure 1, with d ) 20 and w ) 230 µm. The angle between the two arms of the chip ranges between 15 and 45°. To investigate the kinetics of the studied reaction, the two reactants were injected at the same constant flow rate into the microdevice using a single syringe pump (Kd Scientific, KDS101). In the present series of experiments, the concentration of coumarin in acetonitrile was set at b0 ) 1.1 × 10-3 and 1 × 10-3 M, whereas the concentration of triflic acid in acetonitrile was varied between a0 ) 5.6 × 10-3 and 1.12 × 10-1 M. The parameters β ) (a0/b0)1/2 and Γ ) Ke/a0 thus range between 2.26 and 10.59, and between 1.8 × 10-4 and 3.6 × 10-3, respectively. The experiments were performed at room temperature (∼25 °C). The investigated flow rate was between Q ) 2 and 15 µL/min corresponding to mean velocities ranging from vj ) Q/(dw) ≈ 5 mm/s to 5 cm/s. Such values correspond to Reynolds numbers Re ≈ 0.5-5 well below the threshold above which hydrodynamic instabilities appear (Re ≈ 2000). The flow is thus always laminar in the present series of experiments, and the velocity profile is fully developed after some entrance length Le of the order of w. If one considers the following approximation D ≈ 10-9 m2/s for the diffusion coefficients of the studied molecules, the Peclet numbers in our experiments range between 100 and 1000. Therefore, we can neglect the transport by diffusion along the flow, compared to the transport by convection since Pe . 1. Moreover, the 2d assumption is strictly valid only for distances in the main channel downstream X ≈ dPe. It corresponds to distances ranging between a few millimeters and a few centimeters. All our measurements were performed before such values. However, we numerically showed in the Supporting Information that one can use a 2d description of the reaction-diffusion process if one measures the concentration map averaged over the height of the microchannel. More precisely, we quantified the error originating from such an approximation, with the help of the nondimensionalized number Da ) ka0d2/D. In our experiments, NMR measurements showed that k ≈ 2.1 × 105 (M s)-1, thus Da ranges roughly between 102 and 104 depending on the initial concentration a0. According to the results of Figure S-3 (see Supporting Information), the difference between the 2d model and the 3d model concerning the averaged maps of the concentration can be neglected. Optical Detection. As exposed above, the acidic state of the coumarin is strongly fluorescent. Its formation was observed by fluorescence microscopy using a filter cube with the following specifications: excitation filter λexc ) 375 ( 15 nm, dichroic mirror λd ) 420 ( 20 nm, and emission filter λem ) 460 ( 25 nm (see Figure 4). The images of the main channel were collected using a 5× magnification objective (NA ) 0.12, Leica) and were digitalized with a charge-coupled device (CCD) camera (Cohu). After subtraction of a black image to account for the noise of the optical device, the images were corrected from inhomogeneous illumination by dividing with a reference image of the microchannel filled with a uniform concentration of acidic coumarin. For all the experiments, the depth of field of the objective is larger than the channel height. Thus, the recorded images average the signal

arising from the acidic coumarin over the height of the microchannel. At a millimolar concentration in fluorophores, dynamical self-quenching phenomena can be neglected. In addition, reabsorption of emitted fluorescence is not significant in such thin microchannels. Under such conditions, the intensity recorded by the camera is directly proportional to the concentration of acidic coumarin averaged over the height of the microchannel. Numerical Calculations and Fitting Procedures. The solutions of eqs 6 and eqs S-2 are computed using a custom-made solver in C for partial derivative equations. The number of nodes typically ranges between 50 and 200. No flux boundaries were used. Typical simulation times range between 10 and 100 s. Each recorded image from the CCD camera is analyzed using Matlab (The MathWorks, Inc.). After correction of the illumination (vide supra), the measured fluorescence intensity is directly proportional to the concentration in the acidic state C, but we do not know a priori the constant of proportionality. To extract rate constants from the experimental data, we used the following fitting procedure of the experimental images. First, we extracted from each image the experimental curve Cmax(X), which corresponds to the maximal concentration along each line Y as a function of X (the origin of the Y-axis and of the X-axis is set at the Y-junction). Second, we calculated the solution of eqs 6 and, therefore the theoretical curve cmax(x), for the corresponding experiment. To extract k, we dimensionalized the theoretical curve cmax(x) using eq 5 for several values of the rate constant k: one thus obtains k several functions cmax (X). We divided these dimensionalized theoretical curves with the experimental curve Cmax(X). We kept the value of the rate constant that corresponds to the value of k k for which cmax (X)/Cmax(X) is constant over the whole dynamical range. To avoid artifacts due to the diffusive regime where no kinetic information can be extracted, the fit was circumvented to the kinetic regime, i.e., when Cmax(X) significantly varies. This procedure was performed for different flow rates and for different set of initial concentrations. The uncertainty of our fits is mainly of experimental origin (position of the Y-junction, flow rate, quality of the recorded fluorescence images). RESULTS AND DISCUSSION Figure 5a displays the concentration map of the fluorescent acidic coumarin when the two reactants (the protons and the basic coumarin both diluted in acetonitrile) were injected at the same constant flow rate (Q ) 4 µL/min) in the two arms of a Y-shaped microdevice. The respective concentrations are a0 ) 1.02 × 10-1 M and b0 ) 1.1 × 10-3 M, yielding β ) 9.63. This map is recorded using standard fluorescence microscopy. In line with the larger concentration in protons than in basic coumarin (β > 1), and with the anticipated larger diffusion coefficient of the proton with regard to that of the coumarin (χ > 1), the diffusion of the produced fluorescent acidic state is directed into the region of the basic coumarin (region Y > 0, the origin of the Y-axis and of the X-axis is set at the Y-junction) (vide supra). As shown in Figure 5b, the maximal concentration in C along the flow, i.e., Cmax(X), only depends on the time τ ) X/vj elapsed from the Y-junction: the curves obtained for different values of the flow rate collapse. Cmax(X) increases during ∼10 ms and then saturates along the direction of the flow. As exposed above, rate constants can only be extracted from the first kinetic regime. The

Figure 5. (a) Concentration map of the fluorescent acidic state (C) resulting from the reaction and diffusion of the nonfluorescent proton solution (A) injected from the region Y < 0, and of the basic state (B) injected from the region Y > 0 (arbitrary units). The flow rate inside the main channel is Q ) 4 µL/min, and the initial concentrations are a0 ) 1.02 × 10-1 M and b0 ) 1.1 × 10-3 M for the proton and the basic state of the pH indicator, respectively. The parameters are thus β ) 9.63 and Γ ) 1.96 × 10-4. (b) Maximal concentration in C, Cmax, vs τ ) X/vj for several flow rates: 4 (×), 6 (0), and 10 µL/min (O). For the same field of view (the scanned channel length is of the order of 1500 µm), the scanned field of time is reduced from ∼0.1 s for Q ) 4 µL/min to ∼0.04 s for Q ) 10 µL/min. (c) Concentration map of the acidic state in the main channel vs τ ) X/vj and Yr ) (Y - X1/2)/ X1/2.

second regime only contains the information related to the diffusion process that we first address in the following. Figure 5c displays the concentration map plotted against τ ) X/vj and Yr ) (Y - X1/2)/X1/2 to evidence the diffusive regime as explained in the previous section (note that the same holds when plotted against Y/X1/2). The concentration of C evolves according to diffusion processes only, since it can be written as C(X, Y) ≈ F[(Y - X1/2)/X1/2] for τ J 10 ms. In this regime, the particular shape of the concentration profile along Y only depends on β, χ, and Γ. To obtain the function F, we averaged the map displayed in Figure 5c along the X-direction beyond the kinetic regime. Figure 6 displays such functions F for the cases β ) 3.05 and 9.63. Note that the measured functions F do not vanish for large Y as expected theoretically. This discrepancy arises from the nonzero brightness of the basic coumarin. The functions F are fitted in the space of the reduced variable yr ) (y - x1/2)/x1/2. The parameters governing the shape of F are β, χ, and Γ: a0, b0, and Ke are known; χ is therefore the only unknown experimental parameter that remains involved in the fit. Since there is an affine relation between yr and Yr, the theoretical functions F can easily be compared to the experimental ones if we impose these functions to have the same half-width, as well as the same position of their maximum (see Figure 6). Using such a procedure, one only fits the shape of the function F, which only depends on χ. we found χ ) 1.5 ( 0.2: the diffusion coefficient of protons in Analytical Chemistry, Vol. 77, No. 11, June 1, 2005

3423

Figure 6. Experimental functions F measured in the diffusive regime for β ) 3.05 (a) and β ) 9.63 (b), by averaging the corresponding concentration maps plotted against the dimensionless variable Yr (see Figure 5c for instance). Q ) 4 µL/min was chosen to have a large diffusive regime in the field of view of the microscope and thus to get well-defined functions F. The error bars are the standard deviations around the mean values. The continuous lines are the best theoretical fits using eqs 6 and with χ ) 1.5.

Figure 7. Maximal concentrations along the Y-direction Cmax vs τ for β ) 2.26 (a), 3.05 (b), and 9.63 (c). In each case, the flow rate in the main channel is 14 µL/min and b0 ) 1.1 × 10-3 M. The continuous lines are the theoretical curves cmax(τ) for χ ) 1.5 and the corresponding β with k ) 2.1 × 105 (M s)-1. The ticks on the τ-axis indicate τe ) 4.5 ms corresponding to the estimated entrance length w ) 230 µm at this flow rate.

acetonitrile is thus ∼2.25 larger than the diffusion coefficient of the basic state of coumarin. Let us now turn to the kinetic regime. In this regime, the value of k explicitely affects the concentration maps. To perform the fit of the experimental data, one generates dimensionalized theoretical functions cmax(X) for different values of k by using the solutions of eqs 6 and using eq 5 (see the Experimental Section for more details). We retain the value of k that provides the best agreement between the experimental Cmax(X) and the theoretical cmax(X) curves. Figure 7 displays three experimental curves Cmax against τ for different experimental conditions. The flow rate Q ) 14 µL/min, was chosen to zoom on the 0-30 ms dynamical range.29

3424 Analytical Chemistry, Vol. 77, No. 11, June 1, 2005

The results of the fits obtained for the smallest β values perfectly agree with the results of the NMR measurements: we found k ) 2.1 ( 0.2 × 105 (M s)-1 for β ) 2.26 and 3.05. In contrast, significantly lower values (∼5.5 ( 1 104 (M s)-1) were derived for the largest values of β (7.14, 9.63, and 10.59). This discrepancy arises from the absence of complete development of the velocity profile in the microchannel at the time scale of the chemical reaction.16 Indeed, the approximative entrance length Le ≈ w ) 230 µm yields τe ) w/vj ≈ 4.5 ms at Q ) 14 µL/min that is larger than the chemical reaction characteristic times τr at large β. A complete description of the coupling between the reactiondiffusion process and the flow development at the entrance of the channel may be important to access faster kinetics. CONCLUSION In the present work, we introduced the solutions of the simplest model describing the processes of reaction and diffusion of a second-order chemical reaction in a microchannel. With appropriate assumptions, we evidenced two regimes: (i) beyond a distance of the order of vjτr, where τr is a typical time associated to the reaction, and vj is the mean velocity of the flow, the concentration maps only depend on diffusive processes; (ii) for X < vjτr, the concentration maps depend on both the kinetics and the diffusion of the reactants. From an experimental point of view, we successfully applied the theoretical model to investigate the kinetics of a proton exchange reaction. We evidenced the two expected regimes of the reaction-diffusion dynamics. In the kinetic regime, we showed the influence of the entrance effects on the measured rate k, and we successfully measured the rate constant of this reaction for some experimental conditions. Under appropriate conditions, the extracted rate constants are in excellent agreement with those measured by NMR. ACKNOWLEDGMENT The authors are deeply grateful to the members of the microfluidics group at ESPCI, and more specifically to Armand Ajdari and Laure Me´ne´trier-Deremble for many fruitful discussions. J.-B.S. also acknowledges the members of the LOF, RhodiaCNRS laboratory, for many discussions. SUPPORTING INFORMATION AVAILABLE The comparison between the solutions of the 2d model and the averaged 3d model, Figure S-4 concerning the measurement of the protonation constant Ke of the coumarin, and Figure S-5 concerning the NMR relaxation experiments. This material is available free of charge via the Internet at http://pubs.acs.org. Received for review January 14, 2005. Accepted April 6, 2005. AC0500838 (29) Available from A. G. Palmer at: http://cpmcnet.columbia.edu/dept/gsas/ biochem/labs/palmer/software/curvefit.html.