Simulation Tool Coupling Nonlinear Electrophoresis and Reaction

Article pubs.acs.org/ac

Simulation Tool Coupling Nonlinear Electrophoresis and Reaction Kinetics for Design and Optimization of Biosensors Ofer Dagan and Moran Bercovici* Faculty of Mechanical Engineering, TechnionIsrael Institute of Technology, Haifa 32000, Israel S Supporting Information *

ABSTRACT: We present the development, formulation, validation, and demonstration of a fast, generic, and open source simulation tool, which integrates nonlinear electromigration with multispecies nonequilibrium kinetic reactions. The code is particularly useful for the design and optimization of new electrophoresis-based bioanlaytical assays, in which electrophoretic transport, separation, or focusing control analyte spatial concentration and subsequent reactions. By decoupling the kinetics solver from the electric field solver, we demonstrate an order of magnitude improvement in total simulation time for a series of 100 reaction simulations using a shared background electric field. The code can efficiently handle complex electrophoretic setups coupling sharp electric field gradients with bulk reactions, surface reactions, and competing reactions. For example, we demonstrate the use of the code for investigating accelerated reactions using isotachophoresis (ITP), revealing new regimes of operation which in turn enable significant improvement of the signal-to-noise ratio of ITP-based genotypic assays. The user can define arbitrary initial conditions and reaction rules, and we believe it will be a valuable tool for the design of novel bioanalytical assays. We will offer the code as open source, and it will be available for free download at http://microfluidics.technion.ac.il.

E

interest in paper-based assays. Ge et al.8 used electrophoresis to separate three electro-inactive amino acids (serine, aspartic acid, and lysine) on a microfluidic paper-based analytical device (μPAD). Moghadam et al.9 reported the integration of ITP onto a nitrocellulose-based paper microfluidic device to get a 900 fold increase in initial sample concentration. Fast computational tools are an essential part of assay development, as they can be used to provide insight into complex and coupled physics, as well as for assay design and optimization, often significantly reducing experimental time and cost. A rapid simulation tool coupling electrophoresis physics with multispecies reaction kinetics would be highly beneficial for the development of novel bioassays and biosensing strategies. Generic simulation tools such as Comsol Multiphysics could, in principle, be used for such simulations. However, those are not optimized for nonlinear electromigration problems and may require hours or even days to complete.10 Several simulation tools dedicated to the analysis of kinetic reactions also exist. KinTek, Copasi, and Dynafit11−13 are examples of programs for fitting kinetic data based on numerical integration of rate equations. These programs are aimed at understanding complex and multistage chemical

lectrophoretic separation and concentration techniques are extensively used in a wide range of chemical and biochemical applications, including drug discovery, genetics, and food analysis. Recently, there has been growing interest in the use of on-chip electrophoretic techniques for rapid biosensing and point-of-care diagnostics, requiring not only separation and focusing of analytes but also reaction, binding, and hybridization of participating species.1,2 For example, Kawabata et al.3 combined electrophoretically mediated microanalysis (EMMA) and isotachophoresis (ITP) in a mixing and reaction assay for enhanced sensitivity of immune reaction. Persat et al.4 demonstrated simultaneous extraction and detection of nucleic acids by combining ITP with molecular beacons hybridization for a sequence specific detection and quantification of miRNA in human liver from prepurified total RNA. Bercovici et al.5 applied a similar method for rapid detection and identification of bacterial urinary tract infections. Eid et al.6 demonstrated a 12-fold improvement in limit of detection over previous ITP-based hybridization by using a sieving matrix and an ionic spacer to achieve a second step of separation after reaction. Hughes et al. 7 introduced a microfluidic immunoblotting strategy that leverages electrokinetic transport to perform multistep immunoassays, which include separation, immobilization, and binding with probing antibodies. Using this strategy, the authors have demonstrated detection of prostate specific antigen (PSA) in biological samples, without the use of pumps or valves. Recently, due to the need for low cost point-of-care biosensors, there is growing © 2014 American Chemical Society

Received: May 10, 2014 Accepted: July 10, 2014 Published: July 10, 2014 7835

dx.doi.org/10.1021/ac5018953 | Anal. Chem. 2014, 86, 7835−7842

Analytical Chemistry

Article

In the Governing Equations and Numerical Methods section, we present the governing equations and numerical methods used in the code, and in the Results and Discussion present its validation against existing experimental data. We then present its applicability in analyzing acceleration of surface-based reactions using counterflow ITP, where target molecules are continuously focused on top of a surface functionalized with capture probes. The code provides predictions for the rates of reactions and signal gains and reveals two distinct reaction regimes in which the ITP assay changes its efficiency. We also present the use of the code to study the reaction kinetics of competing reactions (e.g., a wildtype and mutant competing for a probe), revealing an optimal detection time for achieving maximal signal-to-noise ratios (SNR). We believe that the code’s capabilities together with its online availability will make it valuable for the design and optimization of electrophoresis driven bioassays.

processes and reactions. However, they all lack coupling with transport physics. Similarly, several rapid simulation tools for prediction of nonlinear electrophoretic transport have been developed. Simul,14 GENTRANS,15 and SPRESSO16,17 are the most common simulation tools reported. These solvers differ from each other in their numerical schemes, but all use an area averaged 1D formulation of the Nernst−Planck equations and solve for coupled multispecies electromigration−diffusion− advection problems. In their original form, these simulations accounted only for equilibrium acid base reactions for determining local pH. Recently, GENTRANS18 was extended for the simulation of chiral separation accounting for analyte− selector complexes in solution. However, these complexations are restricted to chemical equilibria between two monovalent acids and bases. Dubsky et al.19 introduced an extension to Simul, SimulChir, which includes a model for complexation of chiral selectors with enantiomers. More recently, Hruška et al. and Svobodová et al. reported a new version of SIMUL (SIMUL 5 Complex)14,20,21 which accounts for equilibrium complexation of individual analytes and ligands. Despite significant progress in development of electrophoresis simulation tools and the growing importance of on-chip electrophoretic bioassays, no simulation tool currently exists which allows rapid prediction and evaluation of electrophoresis-driven bioassays combining nonequilibrium reaction kinetics. Detailed reviews of existing simulation tools can be found in Thormann et al.,22,23 Mosher et al.,24 and Smejkal et al.25 We here present a new, fast, generic, and open source simulation tool, which integrates nonlinear electromigration with multispecies nonequilibrium kinetic reactions, in the bulk (homogeneous) and on surfaces (heterogeneous). The code is built as an extension to SPRESSO,16,17 inherits its existing numerical methods (SLIP scheme, and adaptive grid) and physical capabilities (ionic strength dependence, nonuniform electroosmotic flow), and adds a new layer of reaction kinetics computations. As illustrated in Figure 1, this enables simulation of electrophoresis-driven assays involving nonequilibrium reaction kinetics of multiple hybridization, binding, and reaction steps. The user can define arbitrary initial conditions and reaction rules, allowing direct design and optimization of novel biosensors and biosensing schemes.

■

GOVERNING EQUATIONS AND NUMERICAL METHODS Formulation and Governing Equations. Accounting for electromigration, diffusion, and reactions, the general threedimensional mass transport relation is given by the Nernst− Planck equation ∂ci , z ∂t

+ ∇·(μi , z Ec⃗ i , z + uc⃗ i , z) = ∇·(Di , z ∇ci , z) + R i , z

i = 1...N , z = ni...pi

(1)

where E⃗ is the electric field vector; u⃗ is the fluid velocity vector; and ci,z, μi,z, Di,z, and Ri,z are respectively the concentration, ionic mobility, molecular diffusivity, and production or consumption rate of an ionic species of the family i, having a valence z. Here, we follow the notation of Štědrý et al.,26 where “family” refers to all ionization states, z, of a species i. ni and pi are the minimum and maximum valences of that family, and N is the total number of species families. When assuming only acid−base reactions, the set of equations can be significantly simplified by summing, within each family i, over all valence states z. This approach, first suggested by Saville and Palusinksi,27 yields N equations for the total concentrations ci = ∑piz=nici,z, ∂(A(x) ci̅ ) ∂ (A(x)μi E ̅ ci̅ + A(x)u ̅ ci̅ ) + ∂t ∂x ∂c ⎞ ∂ ⎛ ⎜A(x)Di i̅ ⎟ i = 1...N = ⎝ ∂x ∂x ⎠

(2)

Here, we formulate the equation in an area averaged form.17 A(x) is the cross-sectional area, and the overbars denote areaaveraged quantities. In such formulation, the reaction terms completely vanish, and instead effective mobilities and diffusivities are used to describe the weighted mean of mobilities and molecular diffusivities of the family ci, pi

μi =

Figure 1. Schematic illustration of the main capabilities included in the code. We combine nonlinear electromigration with multispecies nonequilibrium reactions, in the bulk and on surfaces. This enables simulation of electrophoresis-driven assays involving multiple hybridization, binding, or reacting steps. The user can define arbitrary initial conditions and reaction rules, allowing direct design and optimization of novel biosensors and biosensing schemes.

∑ gi ,zμi ,z , z = ni

pi

Di =

∑ gi ,zDi ,z z = ni

(3)

where gi,z, which is a function of pH and pKa, is the fraction of dissociation of species i, ionized in valence state z. We refer the reader to Bercovici et al.16 for more details about the derivation of gi,z. 7836



Article

this requires computation of the effective mobilities of all such species, which relies on a computationally expensive solution of an algebraic equation for the local pH and (ii) a relatively low cost computation which uses the electric field and pH obtained in stage i and explicitly integrates eqs 4 and 5 in time for the analytes or reacting species of interest. For a single isolated simulation, this approach does not yield any improvement in computation time. However, it is highly advantageous for optimization problems where different reactant properties such mobilities, pKa’s, reaction rates, and concentration can be simulated without recomputing the electric field. To test the improvement in computation time, we simulated a 5 min reaction of a target DNA and a complementary probe within an ITP interface and compared the computational time between a full computation approach and our decoupled approach. Results show a significant improvement in computation time, which increases with the number of simulations. For 30 simulations, the gain in computational time reaches 8-fold. Section S1 in the SI provides full details on this comparison, as numerical results of the speed ratio versus number of simulations. Numerical Implementation. Our code is built as an extension to SPRESSO, an open source, nonlinear electrophoresis solver developed by the Santiago group at Stanford University.16,17 We use the latest version,17 which makes use of a SLIP scheme solver,28 together with an adaptive grid algorithm. The algorithm uses a weighting function that clusters grid points in regions where the numerical scheme adds numerical dissipation. This guarantees unconditional stability and enables robust high-resolution simulation in a variable cross sectional channel. The solver implements a solution of the Nernst−Planck equation and provides the electric field, hydronium ions concentration, cH, and the grid point distribution as a function of time and space. To enable decoupling of the solution, we modified SPRESSO such that these parameters, for each time step, are saved to an output file. Our reactions module, SpressoRXN, then uses the same grid point distribution provided to it by SPRESSO at each time point. Since the length scales associated with electric field gradients are smaller than those associated with analyte distributions,29 the grid distribution inherited from SPRESSO is sufficient for maintaining a stable solution in the simulation of analyte reactions. As further demonstrated in the Results section, we have confirmed this for a large number of challenging cases, including ITP focusing. Surface Reactions. Immobilized species play an important role in a large number of biosensors as they provide inherent separation of reactants from complexes. Their simulation however presents additional challenges due to the analytes’ abrupt transition from moving to stationary species and vice versa. Stationary species are defined the same way as standard species, with the mobility and diffusion properties defined as μi = Di = 0. The concentrations ci in eqs 4 and 5 are related to the surface densities, bi, by ci = (bi/H) with H being the channel height, assuming a low Damkohler number. For more details justifying this assumption, see section S2 in the SI. To simulate a finite surface containing immobilized probes of surface density b0, confined between the boundaries XL and XR, a top-hat distribution of the immobilized species’ initial concentration is a natural choice:

Our interest here extends beyond acid base reactions. In order to capture hybridization or binding processes, reactions between different families must be accounted for. As this is most appropriate for biomolecules, we choose to model second j,k order reactions of the total concentrations; Ri = kj,k oncjck − koffci, where i marks the product of two different species j and k (j ≠ j,k k) and kj,k on, koff are respectively the on- and off-rates of their reaction. Accounting for such reactions, the equation for a product ci formed by a pair of reactants cj, ck, takes the form ⎡ ∂(A(x) ci̅ ) ∂c ⎤ ∂ A(x)⎢μi E ̅ ci̅ + u ̅ ci̅ − Di i̅ ⎥ = R iP + R iR + ⎣ ∂t ∂x ∂x ⎦ (4)

i = 1...N

RPi

RRi

where and are the production or consumption rates of the species in its role as a product or reactant, respectively, and are defined by the following equations: j,k j,k R iP = A(x)kon cj ck − A(x)koff ci

j , k = 1...N

N

R iR =

∑

i,j i,j −A(x)kon ci cj + A(x)koff ck

j , k = 1...N

j = 1, i ≠ j

(5)

It is important to note that each species can serve both as a reactant and a product in the same system. This formulation enables simulation of competitive reactions, such as in the case of a wildtype target DNA and a mutant both competing on the same probe. Furthermore, it allows for multistep reactions, where each species can serve both as a reactant and as a product in the system. To close the set of equations given in eq 4, the area averaged electric field, E̅ , must also be determined. Invoking current conservation and net neutrality27 and assuming that the electric field is locally axial, the electric field can be expressed as17 E̅ =

1 ⎛ I (t ) ∂S ̅ ⎞ + ⎜ ⎟ σ ̅ ⎝ A (x ) ∂x ⎠

(6)

where σ̅ and S̅ are respectively the area-averaged conductivity and diffusive current, given by N

σ̅ =

pi

∑ ∑ zci̅ ,zμi ,z F , i = 1 z = ni

N

S̅ =

pi

∑ ∑ zci̅ ,zDi ,zF i = 1 z = ni

(7)

with F denoting Faraday’s constant. Decoupling of Electric Field and Reaction Kinetics Computations. In principle, the sets of eqs 4 and eqs 6 and 7 could be directly integrated in time to yield the solution for the concentration field at any time. However, an important simplification can be made to achieve a more efficient solver: for the majority of assays of interest, such as capillary electrophoresis, isoelectric focusing, or ITP, the concentrations of reactants are significantly lower than those of the background electrolyte. For example, in peak mode ITP, buffer concentrations are on the order of 10−100 mM, while biomolecules of interest are in the range of 1 fM to 1 μM. Hence, the contribution of the reacting species to the conductive and diffusive current in eqs 6 and 7, as well as to the local pH values, can be neglected. The numerical solution can thus be decoupled into two sequential computations: (i) a relatively high cost computation of the electric field distribution at each time point from eqs 2 and 6 using all high concentration (typically buffer) species 7837



c0, i

⎧ b0 ⎪ XL ≤ x ≤ XR , = ⎨H ⎪ ⎩ 0 x < XL or x > XR

Article

(8)

However, the high concentration gradients at the edges of the surface result in high wave numbers and lead to numerical oscillations. To overcome this fundamental problem, we instead set the initial distribution of the immobilized species to be uniform along the entire computational domain and only enforce a zero reaction term in areas outside the immobilization boundaries, thus the right-hand side of eq 4 take the form: R iP = j,k j,k ⎧ XL ≤ x ≤ XR j , k = 1...N ⎪ A(x)kon cj ck − A(x)koff ci ⎨ ⎪ x < XL or x > XR ⎩0 R iR = ⎧ N i ,j ⎪ ∑ − A(x)kon ci cj XL ≤ x ≤ XR j , k = 1...N ⎪ j = 1, i ≠ j ⎨ i ,j ⎪ + A(x)koff ck ⎪ x < XL or x > XR ⎩0

(9)

Figure 2. Validation of the SpressoRXN code against experimental ITP-aided hybridization data from ref 30. Species A is mixed in the TE at various concentrations as noted in the figure, while species B is mixed with the LE at a fixed concentration of 1 nM. The two species concentrate and focus at the ITP interface. Shown is the fraction of B hybridized versus time. We obtain good agreement with experimental results, with slight underprediction at short times and overprediction at long times. Full details of the experimental and simulation conditions are provided in section S3.1 of the SI.

(10)

■

RESULTS AND DISCUSSION Code Validation. ITP is an electrophoretic technique which uses a discontinuous buffer to focus species of interest at a high electric field gradient formed at the interface between high (LE) and low (TE) mobility ions. Recent work demonstrated the use of ITP for accelerating the hybridization reaction between two DNA strands.30 Such ITP focusing, coupled with reactions, serves as an excellent validation case for our code, as it involves sharp electric field gradients, discontinuous buffers, and kinetic rates which vary with space and time. Bercovici et al.30 performed controlled experiments in which ITP-based hybridization of DNA molecular beacons (MBs) with target DNA was monitored in time. We here use the experimental data for validation of the code. The experimental conditions and the full list of parameters used in the corresponding simulations are provided in section S3.1 of the SI. Figure 2 presents the fraction of hybrid concentration, f ITP, versus time for the three different concentrations of species A and a fixed concentration of B. We obtain good agreement with experimental results, thus validating the underlying assumptions and implementation of the simulation. Prediction of ITP-Aided Reactions under Realistic Conditions. ITP-aided reactions have now been employed in several assays including rapid and highly sensitive detection of miRNA cancer markers4 and detection of urinary tract infections.5 While prediction of reaction kinetics under ITP, and in particular the results presented in Figure 2, can be obtained from an analytical model,30 this model is based on several restrictive assumptions including the existence of one species in excess concentration, constant interface width, overlapping Gaussian-shaped concentration distributions, and pH-independent analyte mobility. Using our code, we are able to provide prediction and optimization for more realistic and general assays in which reacting species have different mobilities (which are pH dependent), reactant concentrations are on the same order of magnitude, where dissociation rates

play an important role, and where nonuniform electroosmotic flow increases dispersion. We demonstrate this with the simulation of two different cases presented in Figure 3. Case I presents a comparison between the analytical model30 and our simulation results for the case in which all model assumptions are maintained. As shown in the inset for case I, one species is present at excess concentration, and the spatial distribution of the concentrations of the two species perfectly overlaps. In this case, as expected, the analytical model predicts well the hybridization process and matches the numerical solution. In case II, we deliberately violate the model assumptions by taking equal concentrations of reactants, different and pH-dependent mobilities, and a relatively large reaction constant (Kd = 10−6 M) and account for dispersion due to nonuniform electroosmotic flow (EOF). In this case, the analytical model significantly overpredicts both the hybridization rate and the fraction of hybridization at the steady state. This is primarily due to the dispersed shape of the analyte concentration which reduces the overall concentration of each of the species, as well as their region of overlap. This example elucidates the need for a simulation tool in accurate prediction of realistic cases. Simulations of Surface Based Biosensors. Surface biosensors are based on reaction of analytes with receptors immobilized to a surface and are used in a variety of biological and chemical applications and in particular in genotyping assays (e.g., DNA arrays) and immunoassays (e.g., ELISA). Such assays traditionally make use of pressure driven flow,31 and in some cases electroosmotic flow,32,33 to introduce the sample, as well as any secondary probes or wash buffers. In several recent works, electrophoresis has been used to deliver targets to patterned surfaces. For example, Hoettges et al.34 combined dielectrophoretic effects with electrohydrodynamic fluid flow to concentrate particles on active sensor surfaces and improve the readout signal. More recently, the Herr group35 demonstrated 7838



Article

Figure 3. Prediction of rapid DNA hybridization using ITP. Comparison of the analytical model30 and our numerical simulation for (I) a case in which all analytical model assumptions are maintained. Results are in good agreement, providing additional validation for our numerical simulation. (II) A more general hybridization case in which the analytical model assumptions are deliberately violated. Shown is the fraction of reactant B hybridized vs time. In case I, we simulated the hybridization of 100 nM of species A premixed with TE and 10 nM of species B premixed with LE. Furthermore, to meet model assumptions, we used equal mobilities μA = μB = μAB, and a low equilibrium constant koff /kon = 3 × 10−11 M. In case II, we deliberately deviated from model’s assumptions, taking equal initial concentrations of A0 = B0 = 10 nM, different mobility values, μA = μAB ≠ μB, and a relatively large dissociation constant koff/kon = 10−6 M. The simulation enables divergence from the restrictive analytical assumptions and investigation of the more interesting and realistic cases of ITP-aided hybridization. The simulation used 400 grid points on an 8-mm-long channel, representing a 50-mm-long channel, and was solved in a frame of reference moving with the interface. LE and TE composition was identical to that described in Figure 2, and the current density was 2830 Am−2.

was applied. See section S3.2 in the SI for complete details on the simulation parameters. We demonstrate the difference between the two cases by running the simulation with different initial target concentrations and observe the ratio of hybridized surface molecules (b(t)) to total surface molecules (bm) at different times. Figure 4 presents a comparison between counterflow ITP surface reactions (based on simulation results) and standard surface reactions in which target continuously flows over the surface (based on eq 11). While standard flow-based reaction results in a uniform distribution of bound targets on the surface, the spatial distribution of the target under ITP results in a nonuniform distribution on the surface. Figure 4a shows the typical distribution of the product (bound targets) on the surface. To obtain a scalar parameter which quantifies the extent of the reaction, we define the average bound surface concentration as

the use of polyacrylamide gel photopatterned with antibodies to achieve automatic immunoblotting of proteins. The same group also demonstrated measurement of kinetic rates using electrophoresis of finite length sample zones through a functionalized polyacrylamide gel.36 Recently, Karsenty et al.37 demonstrated the use of ITP to focus a sample of interest and deliver a high concentration target to a prefunctionalized surface, thus enabling rapid reaction at the sensor site. Due to growing interest in electrophoresis-driven surfacebased assays, we here demonstrate the use of the code in studying the use of ITP for improving the sensitivity of surface biosensors. Following the notation of Squires et al.,38 let bm denote the total surface density [sites/m2] of receptors on the sensor and b(t) denote the surface density of receptors that are bound to target molecules, such that the unbound sites are given by bm − b. For a reaction limited process (i.e., assuming perfect mixing), in which a target concentration is c0, the fraction of bound probes is given by38 c0/Kd b(t ) (1 − e−(1 + c0/ Kd)koff ·t ) = bm 1 + c0/Kd

b ̅ (t ) =

−δ /2

∫+δ/2

b(t ) dx

(13)

where δ = (RμT)(μLEμTE)/[FμLEELE(μLE − μTE)] is the characteristic ITP interface width. Figure 4b presents the average fraction of hybridization, b̅/ bm, for initial nondimensionalized concentration, c0/Kd, ranging from 10−9 to 103 (corresponding to 1 aM to 1 μM), and at three different times: 10 s, 1 min, and 1 h. At all times, the signal from ITP is higher as compared to standard hybridization; however at short times and low initial concentrations, the gain achieved from ITP is the most substantial. At short times (10 s), the gain from ITP is roughly constant across all tested concentrations and improves the signal by approximately 2 orders of magnitude. However, for longer times (1 h), the results reveal two regimes of hybridization. In the first regime, c0/Kd > 10−2, the initial concentrations are sufficiently high such that both the ITP and the standard reactions are dependent mostly on the on rate, kon. In this regime (marked as I), the ITP sensor approaches saturation and the gain from ITP is smaller. The second regime (marked as II) corresponds to c0/Kd ≤ 10−2, where in the standard case the reaction is determined by the off rate, koff, while in the ITP case, the reaction is still kon

(11)

where kd = koff/kon is the equilibrium dissociation constant and kon and koff are the on and off rates of the reaction. In a reaction limited case, the time scale for the sensor to equilibrate is thus given by −1 τR = koff (1 + c0/Kd)−1

1 δ

(12)

Clearly, both the steady state fraction of hybridization and the reaction rate depend directly on concentration. This suggests that increasing the local concentration in the vicinity of the surface would directly benefit the signal and reaction time of such assays. Using the code, we explore the case of an ITP interface which is balanced by counterflow,39−41 while continuously accumulating sample on top of the surface. We performed simulations using the SLIP scheme in an 8 mm long channel discretized with 800 grid points in a uniform and constant grid. Current density was set to 8488 Am−2, resulting in an ITP velocity of VITP = 2.2 × 10−4 m s−1. To balance electromigration, an equal and opposite bulk velocity 7839



Article

Figure 4. Prediction of ITP-based surface reactions and comparison with standard flow based reacdtions. (a) Schematic showing definitions used for analyzing the ITP hybridization case: the solid blue line describes the product concentration under ITP reactions. The maximum concentration is denoted bpeak, and δ is the characteristic ITP interface width. The shaded blue region indicates the amount of product within this region. The dashed line denotes the average product concentration, b̅(t), in this region. (b) Fraction of bound receptors (logarithmic scale) for ITP-aided reactions (O symbols) and standard reactions (X symbols) against nondimensional initial target concentration, at fixed times of t = 10 s, 1 min, and 1 h. Results reveal two reaction regimes, which depend on the initial concentration of the sample, with a maximum gain of approximately three orders of magintude at short times. (c) Ratio of bound receptors (logarithmic scale) against time (linear scale) for several initial concentrations. The two regimes are clealry visible as well, with maximum gain obtained at low concentrations. (d) Normalized peak target (cT, solid line) and product (bpeak, dashed line) concentrations against time. At low initial concentration, target concentration reaches a constant value in which the rate of target molecules entering the ITP interface equals the rate of targets binding to the surface.

every target molecule that enters the ITP interface is immediately captured by the surface. Multispecies Reactions. The quality of a biosensor, designed for identification of targets within a sample, relies not only on its sensitivity but also on its specificity. The challenge arises when other species in the sample compete with the target analyte for the same reaction. This competition mostly depends on the reaction affinity and concentration of the participating species. In this section, we demonstrate the use of the simulation for the optimization of an ITP biosensor used for identification of a target in a competitive environment. As illustrated in Figure 5a, we consider again the case of DNA detection using molecular beacons. However, in addition to the probe (P)−target (T) reaction, we also introduce a mutant (M) which can also hybridize to the probe, though with lower affinity. The kinetic constants we use here are based on experimental measurements by Gotoh et al.42 and are detailed in Table S1 of the SI. The signal-to-noise ratio obtained from a mixture of probes, targets, and mutants can be defined (see complete derivation in section S3.3 of the SI) as43

dependent due to the high concentration at the interface. As a result, the gain from ITP is approximately 3 orders of magnitude. Figure 4c shows the ratio of reaction fractions as a function of time for the same set of cases. The two regimes can be easily seen here as well; at short times the gain across all concentrations is nearly identical. At high concentrations (regime I), this gain quickly diminishes as the number of available sites decreases and the sensor approaches saturation. In contrast, at low concentrations (regime II), we obtain nearly constant gain with a slight increase in long times as most of the target molecules, reaching the surface, hybridize. Figure 4d presents the normalized target (solid line) and product (dashed line) peak concentration against time, for three different initial concentrations. At high concentrations (belonging to regime I), the surface approaches saturation, and thus the concentration of targets at the ITP interface continues increasing. In contrast, at low concentration (belonging to regime II), the concentration of bound probes, b, continuously increases while the concentration of focused targets at the interface, cT, reaches a constant. This could be explained by the fact that as the target concentration increases, so does the reaction rate with the surface. As long as the surface is far from saturation, a steady state is obtained in which the flux of targets entering the ITP zone equals the reaction flux with the surface. Under the conditions simulated here, the normalized target curve asymptotically approaches cT/c0 = 5500. Here, ITP creates highly efficient reaction conditions in which essentially

γ

SNR =

c T− P + cM− P + α cP γ

control control cM −P + α cP

(14)

where cT−P, cM−P, and cP are respectively the concentrations of the hybridized target probe, mutant probe, and free probe and 7840



Article

used in our simulation and reaction rates of the target and mutant to the probe are defined in Table S1 of the SI. Figures 5b presents the signal-to-noise ratio as a function of time for three different initial probe concentrations. Figure 5c presents the absolute signal vs time for these conditions. As expected, results show a clear tradeoff between signal and SNR; i.e. a high probe concentration contributes to higher signal but also increases the noise signal and thus reduces the SNR. However, while at long times and at low concentrations, the SNR asymptotically declines to 1 and it is essentially impossible to detect the target, the kinetic simulation reveals an optimal time (in this case, approximately at 10 s), in which the SNR achieves a maximum value and the target can be detected. The peak in SNR is formed due to the competition of characteristic hybridization time scales associated with the target−probe and mutant−probe hybridizations; at short times, the higher on rate of the target−probe reaction results in faster hybridization of the target compared to the probe, and thus a steady increase in SNR. However, as the free targets deplete, the reaction rate of the mutant becomes dominant, resulting in an increase in noise and decrease in SNR. In contrast to standard hybridization, here the reaction is strongly coupled to focusing dynamics which govern the concentration at all times. Since the ITP interface continuously electromigrates along the channel, the existence of an optimal detection time translates in a practical application into an optimal detector location. This important consideration must be weighed against the standard consideration for a signal, whose optimization places the detector as far downstream as possible.

■

CONCLUSIONS We have developed and validated a new simulation tool which couples nonlinear electromigration with reaction kinetics. Our 1-D simulation tool is aimed at assisting the development, analysis, and optimization of new electrophoresis based biosensors and assays involving nonequilibrium chemical reactions, such as DNA hybridizations and immunoassays. Importantly, the simulation can be applied to both homogeneous reactions (in the bulk) and surface base reactions (e.g., with immobilized probes). We have shown that in addition to serving as an engineering tool for bioassay design, the code can also be used to investigate and obtain deeper understanding of complex electrophoretic setups which include coupled physical mechanisms such as competitive reactions. For example, we studied the use of counterflow ITP for acceleration of surface-based binding reactions and identified two distinct concentration regimes, which dictate the gain compared to standard hybridization by flow. In addition, we demonstrated that the detection time in ITP hybridization assays could be optimized to maximize SNR. For example, a wildtype target that could not be detected at the steady state in the presence of a competing mutant could be easily detected in a relatively early time window. By decoupling the kinetics solver from the electric field solver, we are able to obtain fast simulation times. As demonstrated, this method is in particular beneficial in optimizations where multiple reaction cases are performed on the same background electric field. While only a finite and specific set of examples could be shown within the scope of this paper, the simulation can be used to design and optimize a wide range of new linear and nonlinear electromigration-reaction assays. We will offer the code as an open source so that other researchers may be able to

Figure 5. Simulation of competing hybridization reactions under ITP. (a) We simulate the case of a probe sequence, P, which can hybridize to either its target DNA strand, T, or to a competing mutant DNA strand, M. (b) SNR as a function of time, for constant target and mutant concentrations of 10 nM and probe concenrations varying from 1 nM to 100 nM. The normalized contributions to fluorescence intensity of the hybridized T−P, M−P, and the free probe are α = β = 1 and γ = 0.1. Results reveal that while in low concentrations detection of the target cannot be achieved at a steady state due to low SNR, reaction kinetics results in a short time window in which SNR is significant. (c) As expected, the absolute signal increases monotonically in time. There is thus a trade off between obtaining sufficient SNR and a sufficient signal.

α and γ are constant factors relating the concentration of each state to a fluorescent intensity. This SNR value should be interpreted as a quantitative measure for the ability to detect the target of interest in the presence of a competing mutant. We consider the case of a sample containing a target concentration T0 and mutant concentration M0. The sample is mixed with the TE, which also contains a fixed concentration of molecular beacon probes. Under ITP, all three species focus and competitively hybridize at the interface. To compute the SNR, we performed two sets of simulations: one using only 10 nM of mutant sequences (for computing the noise) and one containing an equal amount of target sequences (for computing the signal). This is a relatively extreme case, in which the competing species is at an equal concentration. For each case, we varied the probe concentration in the well between 1 nM and 100 nM. Parameter values 7841



Article

use and modify it to fit their specific needs. It will be available for free download at http://microfluidics.technion.ac.il.

■

(24) Mosher, R. A.; Breadmore, M. C.; Thormann, W. Electrophoresis 2011, 32, 532−541. (25) Smejkal, P.; Bottenus, D.; Breadmore, M. C.; Guijt, R. M.; Ivory, C. F.; Foret, F.; Macka, M. Electrophoresis 2013, 1493−1509. (26) Štědrý, M.; Jaroš, M.; Hruška, V.; Gaš, B. Electrophoresis 2004, 25, 3071−3079. (27) Saville, D. A.; Palusinski, O. A. AIChE J. 1986, 32, 207−214. (28) Jameson, A. Int. J. Comput. Fluid Dyn. 1995, 4, 171−218. (29) Rubin, S.; Schwartz, O.; Bercovici, M. Phys. Fluids 2014, 26, 012001. (30) Bercovici, M.; Han, C. M.; Liao, J. C.; Santiago, J. G. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 11127−11132. (31) Kartalov, E. P.; Zhong, J. F.; Scherer, A.; Quake, S. R.; Taylor, C. R.; French Anderson, W. BioTechniques 2006, 40, 85−90. (32) Jacobson, S. C.; Hergenroder, R.; Moore, A. W. J.; Ramsey, J. M. Anal. Chem. 1994, 66, 4127−4132. (33) Burke, B. J.; Regnier, F. E. Anal. Chem. 2003, 75, 1786−1791. (34) Hoettges, K. F.; Hughes, M. P.; Cotton, A.; Hopkins, N. A. E.; McDonnell, M. B. IEEE Eng. Med. Biol. Mag. 2003, 22, 68−74. (35) He, M.; Herr, A. E. J. Am. Chem. Soc. 2010, 132, 2512−2513. (36) Kapil, M. A.; Herr, A. E. Anal. Chem. 2014, 86, 2601−2609. (37) Karsenty, M.; Rubin, S.; Bercovici, M. Anal. Chem. 2014, 86, 3028−3036. (38) Squires, T. M.; Messinger, R. J.; Manalis, S. R. Nat. Biotechnol. 2008, 26, 417−426. (39) Everaerts, F. M.; Vacík, J.; Verheggen, T. P. E. M.; Zuska, J. J. Chromatogr. A 1971, 60, 397−405. (40) Reinhoud, N. J.; Tjaden, U. R.; van der Greef, J. J. Chromatogr. A 1993, 641, 155−162. (41) Bhattacharyya, S.; Gopmandal, P. P.; Baier, T.; Hardt, S. Phys. Fluids 2013, 25, 022001. (42) Gotoh, M.; Hasegawa, Y.; Shinohara, Y.; Shimizu, M.; Tosu, M. DNA Res. 1995, 2, 285−293. (43) Bonnet, G.; Tyagi, S.; Libchaber, A.; Kramer, F. R. Proc. Natl. Acad. Sci. U. S. A. 1999, 96, 6171−6176.

ASSOCIATED CONTENT

S Supporting Information *

Additional information as noted in text. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was supported by the Israel Science Foundation (grant nos. 512/12 and 1698/12). We also gratefully acknowledge support from FP7Marie Curie Career Integration Grant No. PCIG09-GA-2011-293576.

■

REFERENCES

(1) Dittrich, P. S.; Tachikawa, K.; Manz, A. Anal. Chem. 2006, 78, 3887−3908. (2) Trietsch, S. J.; Hankemeier, T.; van der Linden, H. J. Chemom. Intell. Lab. Syst. 2011, 108, 64−75. (3) Kawabata, T.; Wada, H. G.; Watanabe, M.; Satomura, S. Electrophoresis 2008, 29, 1399−1406. (4) Persat, A.; Santiago, J. G. Anal. Chem. 2011, 83, 2310−2316. (5) Bercovici, M.; Kaigala, G. V.; Mach, K. E.; Han, C. M.; Liao, J. C.; Santiago, J. G. Anal. Chem. 2011, 83, 4110−4117. (6) Eid, C.; Garcia-Schwarz, G.; Santiago, J. G. Analyst 2013, 138, 3117. (7) Hughes, A. J.; Lin, R. K. C.; Peehl, D. M.; Herr, A. E. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 5972−5977. (8) Ge, L.; Wang, S.; Ge, S.; Yu, J.; Yan, M.; Li, N.; Huang, J. Chem. Commun. 2014, 50, 5699−5702. (9) Moghadam, B. Y.; Connelly, K. T.; Posner, J. D. Anal. Chem. 2014, 86, 5829−5837. (10) Dixon, D. R.; Clark, S. B.; Ivory, C. F. Electrophoresis 2012, 33, 880−888. (11) Johnson, K. A.; Simpson, Z. B.; Blom, T. Anal. Biochem. 2009, 387, 20−29. (12) Hoops, S.; Sahle, S.; Gauges, R.; Lee, C.; Pahle, J.; Simus, N.; Singhal, M.; Xu, L.; Mendes, P.; Kummer, U. Bioinformatics 2006, 22, 3067−3074. (13) Kuzmič, P. In Methods in Enzymology; Johnson, M. L., Brand, L., Ed.; Academic Press, 2009; Vol. 467, pp 247−280. (14) Hruška, V.; Jaroš, M.; Gaš, B. Electrophoresis 2006, 27, 984−991. (15) Breadmore, M. C.; Mosher, R. A.; Thormann, W. Anal. Chem. 2006, 78, 538−546. (16) Bercovici, M.; Lele, S. K.; Santiago, J. G. J. Chromatogr. A 2009, 1216, 1008−1018. (17) Bahga, S. S.; Bercovici, M.; Santiago, J. G. Electrophoresis 2012, 33, 3036−3051. (18) Breadmore, M. C.; Kwan, H. Y.; Caslavska, J.; Thormann, W. Electrophoresis 2012, 33, 958−969. (19) Dubský, P.; Tesařová, E.; Gaš, B. Electrophoresis 2004, 25, 733− 742. (20) Hruška, V.; Beneš, M.; Svobodová, J.; Zusková, I.; Gaš, B. Electrophoresis 2012, 33, 938−947. (21) Svobodová, J.; Beneš, M.; Hruška, V.; Ušelová, K.; Gaš, B. Electrophoresis 2012, 33, 948−957. (22) Thormann, W.; Caslavska, J.; Breadmore, M. C.; Mosher, R. A. Electrophoresis 2009, 30, S16−S26. (23) Thormann, W.; Breadmore, M. C.; Caslavska, J.; Mosher, R. A. Electrophoresis 2010, 31, 726−754. 7842


Simulation Tool Coupling Nonlinear Electrophoresis and Reaction

Recommend Documents