Modeling

It has previously been demonstrated that for target-limited conditions and hybridization times up to 12–24 h the hybridization reaction timescale is markedly shorter than the diffusion timescale

[35]. Consequently, reaction terms can be neglected in describing the spatiotemporal evolution of the target concentration and a Fickian diffusion model accurately captures the dominant dynamics of the full diffusion-reaction model

[35]. This simplification is supported below by a rigorous scaling analysis, which we hope might be of interest also to future hybridization studies.

Thermodynamic models, such as those based on Langmuir adsorption theory, are a popular approach to modeling microarray hybridization

[36],

[37],

[38],

[39]. Equilibrium Langmuir-based models, however, cannot adequately describe the spatial bias observed in this study, because all spots in our simplified array consist of the same probe molecules with identical thermodynamic properties (e.g. hybridization free energy) and would therefore be predicted to have identical signal intensity by thermodynamic models. Extensions of the Langmuir model that take into account the effect of probe depletion

[40] are also inadequate for this purpose because even in the case of global probe depletion, in which the behavior of each spot affects other spots, an assumption of perfect mixing is levied that disallows the possibility of location-dependent phenomena.

Extensions of Langmuir-based models to diffusion-reaction models are necessary to capture spatial variations because they incorporate transport processes

[35],

[41]. Gadgil and co-workers proposed and solved such a model

[35]. In the bulk fluid, reaction is irrelevant and target dynamics are described by the diffusion equation

[35]. Close to the spot surface, both diffusion and reaction could in principle be important. For target-limited conditions and hybridization times up to 12–24 h, Gadgil et al.

[35] found that the hybridization reaction timescale is markedly shorter than the diffusion timescale. Diffusion is thus negligible near the spot, where dynamics are dominated by reaction. This justifies the use of a perfect absorption boundary condition (

*C*=

0) at the spot surface. To show this explicitly, we present a formal scaling analysis that demonstrates the smallness of a dimensionless parameter,

*P*, measuring the ratio of the hybridization timescale to the diffusion timescale.

The diffusion-reaction equation for the concentration of target,

*C*, is

where

*D* is the diffusion coefficient of the target,

*r* and

*z* are the radial and vertical coordinates,

*k* is the hybridization rate constant, and

*B* is the concentration of probe. In equation (1), we have not included the dissociation of bound target strands because an assumption of irreversible hybridization kinetics is justified for perfectly complementary hybridization

[6], as is the case in our experiments. It is useful to derive a dimensionless form of equation (1) by introducing dimensionless variables (denoted by a hat) as follows:

The characteristic scales used to make variables dimensionless are the spot radius

*R*, the chamber height

*h*, a reference probe concentration

*B*_{0}, an (arbitrary) reference target concentration

*C*_{0}, and the diffusion timescale

*R*^{2}/D. Equation (1) can be written in terms of the dimensionless variables as,

where the dimensionless parameter

represents the ratio of the hybridization timescale

to the diffusion timescale

. Incidentally, we note that the limiting timescale for diffusion is the one pertaining to diffusion in the radial direction, rather than in the vertical direction, since in our case (and in

[6])

*R*<

*h* (if instead one had

*R*>

*h*, one would have to redefine

*P* as

.

The value of

*P* determines whether it is justifiable to neglect either reaction or diffusion in the vicinity of the spot. For

*P*1, reaction is too slow to matter and can be neglected, whereas for

*P*1 reaction is much faster than diffusion, and the latter can be neglected. Most hybridizations fall in the latter category. Importantly,

*P*1 in equation (3) yields the condition,

, which corresponds to a boundary condition of perfect absorption at the spot surface (

*C*=

0).

For example, for the case considered by Gadgil et al.

[35] (

*R*=

50 µm,

,

.,

,

=

78 µM), we find a reaction timescale

and a diffusion timescale, (

*D*/,

*R-2.*)-

*−1*.

*=250 s*, thus obtaining

*P*=

5×,10-

*−5*.<

*<1*. The smallness of this value supports the conclusion – obtained by Gadgil et al.

[35] through full integration of equation (1) – that at the spot surface reaction dominates over diffusion and thus

*C*=

0.

The same result applies to the system under consideration here. We performed a careful analysis of the parameters involved by resorting to the relevant literature (we note that our values are not the same as those of Gadgil et al.

[6], though these differences have little effect on the end conclusion). The radius of our spots was

*R*=

75 µm. We used a diffusion coefficient of

*D*=

1×,10-

*−9*., m-

*2*.,

*s*-

*−1*., based on published estimates for single-stranded RNA

[42]. We estimated

*k* according to the Wetmur-Davidson equation

[43] to be 5.0×,10-

*4*., M-

*−1*.,

*s*-

*−1*. Although

*k* is known to be sensitive to several factors

[44], this value agrees well with experiments

[45]. To estimate

*B*_{0}, we followed the method of Cheung et al.

[46] and assumed that the spot is originally hemispherical with radius

*R* and has a concentration

*B*_{S}, then dries to a disk of the same radius and height

*h*_{D}, while conserving the total amount of probe. Assuming

*h*_{D}=

2 µm

[35] and using the value

*B*_{S},=

50 µM from our experiments, this yields,

*B*_{0}=

(2/3)(

*R*/,

*h*_{D}.),

*B*_{S}=

1.25 mM. Using these values, we find a reaction timescale

. and a diffusion timescale

, thus obtaining

*P*=

2.8×,10-

*−3*.

*1*. It is thus justified to assume that the spot acts as a perfect sink, where

*C*=

0.

We solved the diffusion equation describing the target concentration under conditions representing the essential features of the chamber's geometry and operation, with the goals of (i) understanding the dominant process responsible for hybridization heterogeneity, and (ii) exploring the role of geometrical parameters in determining the magnitude of this heterogeneity. For this purpose, we solved the three-dimensional unsteady diffusion equation

where

*C* is the concentration of the target,

*D* is the target's diffusion coefficient,

*t* is time, and

^{2} is the Laplacian operator.

Equation (4) was solved numerically using a finite-element code (COMSOL Multiphysics) for the dimensions of the simplified microarray (), with the exception that the chamber height assumed in the numerical simulation was 150 µm or 750 µm. In detail, the numerical model represents a hybridization chamber of half-width

*L* and height

*H*, containing a square block of 30×30 spots that covers an area of

*M*×

*M* mm

^{2}. Each spot had a diameter of 150 µm and the inter-spot distance is 150 µm. Hence, the 30×30 spots covered 8.5×8.5 mm (

*M*=

8.5 mm). Both the height of the chamber and the relation between

*M* and

*L* will be seen to be important in determining hybridization heterogeneity.

Each simulation was carried out over 24 hours, with a diffusion coefficient of

. Importantly, a rescaling of the diffusion coefficient is simply equivalent to a rescaling of time, because the diffusion equation is linear. For example, if the diffusion coefficient is halved (larger target molecules), results will apply unchanged, except that all times will be doubled. The simulation proceeds by solving the three-dimensional diffusion equation in the entire chamber at subsequent instants in time. The initial condition was

*C*=

1 in the entire chamber, i.e. a constant concentration of target (well-mixed). Because the diffusion equation is linear in the concentration, the solution is independent of the particular initial target concentration (allowing us to set

*C*=

1 without loss of generality), unless the concentration is so high that individual spots become saturated (a case we do not consider here).

For each spot, the flux of target strands to the spot was calculated at each instant in time and subsequently integrated over the total time (24 h). For this, the boundary condition imposed at the spot was one of perfect absorption (

*C*=

0), consistent with reaction being considerably faster than diffusion, as detailed above. The flux was thus computed based on Fick's law as –

*DA*(

*C*/

*z*), where

*C*/

*z* is the gradient of target concentration normal to the spot and

*A* is the surface area of the spot. This yielded the total amount of target that had been bound to the spot over the simulation time. At all other surfaces (top, bottom and side walls of the chamber), which were impermeable to the target, no-flux boundary conditions were imposed.

Three different cases were considered. The first is the closest to the simplified array described above, with a large and shallow hybridization chamber (

*L*=

10.8 mm,

*H*=

150 µm), differing from the real one only in the chamber height (the actual chamber has

*H*=

250 µm). The second models a large and deep hybridization chamber (

*L*=

10.8 mm,

*H*=

150 µm), to explore the effect of chamber height. The third case is that of a small and shallow chamber (

*L*=

4.5 mm,

*H*=

150 µm), to understand the effect of the amount of surface area not covered with spots. In all cases, the number of spots and the area covered by the spots are the same.