2. Materials and Methods A three-dimensional mathematical model was implemented to describe the hybridization of a DNA target to a surface-immobilized probe in a glass wafer microfluidic biochip (obtained from Invitrogen, formerly Xeotron, [13]). The microfluidic biochip had several independent hybridization wells (15 μm deep and 56 μm in diameter) where probes were attached and connected through several channels according to light-directed synthesis technology previously reported [13,21]. The utility of this biochip for microarrays including the related validation experiments (e.g., materials used, probe attachment parameters, etc.) was previously studied and published [13,21]. In this study, to evaluate the hybridization kinetics, only a single hybridization well was considered (Figure 1). In physical experiments, the target DNA is captured in many micro-reactors available for hybridization. Though the effect of ion concentration in the hybridization buffer is known to significantly impact hybridization, we have left the hybridization buffer constant according to published methods as this has been previously optimized [22]. The target DNA solution flows into the micro channel at a constant flow rate. The flow of an aqueous solution of DNA inside a small channel is assumed as uniform. The target DNA molecules are transported by convection along the flow direction where they are free to diffuse. Fluid flow is circulated through the chip multiple times, so the fluid flow rate evaluated in this study is, in fact, the recirculation flow rate. When the target DNA molecule reaches the hybridization well, the well-known association–dissociation reaction occurs between the target DNA and the immobilized probes. Both the association and dissociation events are modeled here using a differential kinetic equation. 2.1. Mathematical Model Only one of the many hybridization wells, connected with a microchannel, is considered in the three-dimensional model (Figure 1). The fluid containing target DNA flows from left to right where the incoming flow profile is characterized by fully developed laminar flow, i.e., it is parabolic with zero velocity at the channel walls. Fluid flow in the channel follows the Navier–Stokes equation, (1) ρ∂u∂t−η∇2u+ρ(u⋅∇)u+∇p=F (2) ∇⋅u=0 where u denotes the velocity field vector, ρ is density, η is the dynamic viscosity, and p is pressure. At steady-state, the first time-dependent term disappears. The external force, F, such as gravity can be neglected in such miniaturized systems. The incoming flow has a small concentration of the target DNA which is described by a convection–diffusion equation of the form, (3) ∂c∂t+u⋅∇c=D∇2c where c denotes the concentration of solution-phase targets, u is the identical velocity field vector given in Equation (1). The diffusion coefficient, D, was estimated using Table 1 of Chan et al. [23], where it was reported that the diffusion coefficient of DNA depends on its length and decreases with increasing base pairs [23]. Once the target DNA reaches the bottom of the hybridization well, the association and dissociation kinetics are described by using an ordinary differential equation, (4) dBdt=konc(Rt−B)−koffB where kon represents the association rate constant, koff is the dissociation rate constant, Rt is the total surface concentration of probes, and B is the surface concentration of bound targets at time, t. The association rate constant, kon, depends on the DNA sequence involved, the temperature, and the ionic strength of the medium but not to the point that it would drastically affect the order of magnitude. Based on the experimental results reported in literature [11,24,25], it was assumed that the association rate constant does not change appreciably over the temperature range used in our simulation (25–55 °C). Therefore, a value of 106 M−1s−1 for kon was used in these simulations. The dissociation rate constant, koff, was calculated using the thermodynamic model presented in literature [12]. In this model, standard Gibbs free energy and hybridization energy for different DNA sequences was used to calculate the dissociation rate constant at a specific temperature. The model also assumes that the DNA targets do not diffuse on the surface and that there is no leakage of the molecules at the edges of the surface. 2.2. Numerical Simulations Numerical simulations of a single hybridization well were performed with COMSOL Multiphysics (FEMLAB), a finite element software package, which uses two geometries and three application modes to solve the formula. The first geometry is three-dimensional and represents the channel connecting the micro well. Within this geometry, the model used two application modes: incompressible Navier–Stokes, and convection and diffusion to model the transport of the target DNA. The second geometry is two-dimensional and simulates the reaction surface with the diffusion application mode. It is assumed that the hybridization reactions take place on a 2D surface (bottom surface of the micro-wells) whereas bulk flow of incoming target DNA is resolved on a 3D geometry. To correctly model the multi-scale physics, extrusion-coupling variables were defined. These variables couple the mass transport from 3D bulk flow using a surface boundary condition to the hybridization reactions on a 2D surface. In COMSOL Version 4.2 and later, instead of using extrusion-coupling variables, a Surface Reaction interface can be used with a single 3D geometry. This interface resolves surface coverage of the adsorbed species, with a mass source or sink automatically coupled to the surface boundary condition in the bulk mass transport equation. The meshing around the channel was not greater than 10 μm and around the hybridization well it was smaller than 0.2 μm. In all simulations, an “insulation” boundary condition was chosen for the non-reactive surfaces at the top and the bottom of the microchannel except where the hybridization occurred. The model is solved in two steps using different solvers. First, it solves the incompressible Navier–Stokes application mode with a non-linear solver followed by the time-dependent solver to simultaneously solve the convection and diffusion and the diffusion application modes. The model predicts the concentration of duplexes formed due to interactions between the target DNA and immobilized probes inside the hybridization well. In order to correlate the duplex concentration obtained through our model with the hybridization signal intensities, we used Equation (5) [26], (5) I=L1+e−d/θ where I is log fluorescence intensity, d is log dye concentration, θ defines the spread and slope of the linear range of the curve and the “background” level (set to 1.0 for Cy3 dye), L is the upper limit of the dynamic range (set to 3.0). Preset values were taken from the study [26]. For our simulation, we assumed that there was only one dye molecule attached per molecule of target DNA because it was end-labeled. The complete algorithm and the data flow used for our model are illustrated in Figure 2. 2.3. Experimental Validation The mathematical model was validated using the real-time hybridization experiments performed on Xeotron hybridization station (Xeotron Corp., Houston, TX, USA—now part of Invitrogen). Eighteen targets were designed to complement the 18 probes for the hybridization experiments (Table 1). However, this is not a large enough sample size to make predictions about the effect of individual melting temperatures on hybridization kinetics due to other confounding factors such as mismatches [27]. The targets (45–52 mers) were synthesized and amine modified on the 3′ end by Operon Biotechnologies Inc. (Huntsville, AL, USA - part of MWG-Biotech AG, Ebersberg, Germany). The targets were coupled with Cy3 (Amersham Biosciences, Little Chalfont, United Kingdom, Cat. No. PA23001) and then dissolved in dimethyl sulphoxide (Sigma-Aldrich, St. Louis, MO, USA, Cat. No. D2438 in 75 mM sodium carbonate buffer at pH 9.0). A hybridization buffer containing 35% deionized formamide (Thermo-Fisher Scientific, Waltham, MA, USA), 6x SSPE (pH 6.6; Invitrogen; 1x SSPE contains 0.18 M NaCl, 10 mM NaH2PO4, and 1 mM EDTA pH 7.7), and 0.4% Triton X-100 (Sigma-Aldrich) [22]. Several hybridization experiments were performed to assess the effects of temperature, concentration and flow rate. The real-time data images of nucleic acid hybridization were obtained using OSA Reader v1.91 software (IMSTAR, Paris, France) in conjunction with a CCD camera (resolution 10 μm per pixel) attached to a Nikon Eclipse 50i fluorescent microscope (100× magnification; Nikon Instruments, Melville, NY, USA).