INTRODUCTION
============

Current DNA microarray technology utilizes multiple oligonucleotide probes to detect the concentration of target molecules. These probes, even though interrogating the same target, often yield very different signal intensities. Without understanding the physicochemistry underlying this problem, the quantification of absolute gene abundance is unattainable and inter-probe comparison is unjustified, leaving DNA microarray technology severely compromised.

A number of physical models have been proposed to address this problem, mostly in the form of Langmuir derivatives (1--10). The Langmuir model is a generic mathematical form that also fits the description of first-order chemical reactions, which is frequently used for probe/target binding on DNA microarrays: where θ is the fraction of occupied probes, *T* free target concentration, *K* dissociation constant.

According to the Langmuir model, all probes should saturate at the same level, which is clearly not the case in microarray hybridizations. Various modifications were proposed to accommodate this difference in saturation levels. A generic version may be written as where χ is a probe-specific factor. While a physical meaning of χ is difficult to obtain, some ([@B7],[@B8]) tried to explain χ through the washing step in microarray experiments. That is, all probes reach the same saturation level in the end of hybridization, but they lose the bound targets to different extents during the washing step. This 'washing model' suggests a significant loss of signals upon each washing cycle. In experimental observations, the first washing cycle usually removes a considerable amount of partially bound targets, but it is clear that signal intensities do not decrease dramatically after extra washing cycles ([@B11]). This contradicts the above 'washing model'. Furthermore, the Langmuir derivatives predict that, in response to increasing target concentrations, probes with higher binding affinities saturate first. In experimental observations, on the contrary, low-affinity probes generally saturate first. Although Langmuir models seem to work well on simple surface hybridizations, no Langmuir derivative has adequately predicted probe signals in 'real' experimental settings, such as those in the Affymetrix Latin Square data with complex backgrounds ([@B12]).

The best prediction of probe signals to date was reported by Zhang *et al.* ([@B13]). They accounted for both specific binding and nonspecific binding in the form of , where *T̂* is total target concentration, while fitting 83 parameters to the data. Mei *et al.* ([@B14]) also sought a linear composition of binding energy, where the single base energy contribution alone used 75 parameters. Over-parameterization has been a concern in all these previous studies and invited criticism on their general applicability ([@B15]).

After all, a valid physical model of microarray hybridization will have to explain the probe difference through sequence-specific thermodynamics, as its oligonucleotide sequence is the defining property of a microarray probe. The free energy of polynucleotide hybridization in bulk solution has been successfully described by a nearest neighbor (NN) model ([@B16],[@B17]). However, this NN model is widely regarded as not applicable to high-density microarray hybridization, as it was either modified and re-parameterized ([@B5],[@B7],[@B9],[@B13],[@B18]), or abandoned ([@B1],[@B14],[@B19]).

We will first demonstrate that the NN model is applicable to probes free of secondary structures. With the thermodynamic component calculated from the NN model, we then propose a new competitive hybridization model to describe the kinetics. Our model, using only four global parameters that are fitted to Affymetrix Latin Square data, can successfully predict signal intensities of individual probes, and therefore, achieve the absolute quantification of target concentrations.

METHODS
=======

The Affymetrix Latin Square spike-in data U133A were retrieved from ([@B12]). They contained 14×3 hybridizations where spike-in targets were added at various concentrations from 0 pM to 512 pM. Probe information was obtained through ([@B20]), where only 30 of the 42 probesets were found. A total of 365 probes matched to target sequences. Among them, 10 probes with very low signal intensities (under 900 at highest target concentration) were removed. In total, 355 probes are included in this study. Background was taken as the signal intensity at zero spike-in concentration, and subtracted from data at other concentrations. No normalization was performed on these data. The probe self-folding energy, , was computed by RNAStructure \[version 4.5, function OligoWalk ([@B21])\]. Duplexing energy, , was computed by the current NN model with the parameters from Ref. ([@B17]). Compiled data and computational scripts used in this study are available upon request.

RESULTS AND DISCUSSION
======================

Thermodynamic predictability in DNA microarrays
-----------------------------------------------

Microarray probes and targets may form secondary structures by intramolecular self-folding. These structural effects are not accounted for in the NN model, posing a problem to the thermodynamic calculation. As a first step, we investigated the structural effects through the self-folding energy of probes, . In the Affymetrix Latin Square data (see Methods Section for details), about 45% of the probes can be selected by the criterion . For these probes, a clear correlation appears between log signal intensities () at the highest target concentration and the duplexing energy that are computed by the current NN model with existing parameters ([Figure 1](#F1){ref-type="fig"}, ). If the selection criterion is relaxed to , 75% probes are included and the correlation has (data not shown). However, the correlation diminishes at lower target concentrations (data not shown). These observations suggest that the current NN model offers a certain degree of predictability, but they cannot be accommodated by previous, Langmuir-like models. A new kinetic model is needed. Figure 1.Duplexing energy calculated by NN model correlates with signal intensity at 512 pM, the highest spike-in target concentration, for probes free of secondary structures (*ΔG*~0~ \> −1, 160 probes from Affymetrix U133A data). Each dot represents a probe.

A competitive hybridization model
---------------------------------

We treat DNA microarray hybridization as two subprocesses, the binding of targets to probes and the dissociation of target/probe duplexes. Assuming that equilibrium is reached at the end of hybridization and the binding rate is the same for all target molecules (see below), the dissociation rate is governed by the duplexing energy between paired target/probe. A kinetic equilibrium between binding and dissociation should be observed.

Two types of targets are explicitly modeled: 'specific targets' (perfect match) with probe-specific dissociation rate *k*~d~, and 'cross-hybridizing targets' with dissociation rate *k*~n~. These cross-hybridizing targets are present in large quantities because partially matching sequences are abundant in a transcriptome. For the moment, we simplify them as a uniform mixture with a probe-nonspecific *k*~n~.

The target/probe duplex formation is commonly believed to start with an initiation step, the base-pairing between a small number of nucleotide bases, and then extend to the rest of complementary regions ([@B22],[@B23]). If the initiation step sets the rate limit, the binding rate should be hardly specific to probe sequences. We therefore assume a single binding rate, *k*~b~, for all target molecules. How the specific factors ([@B24],[@B25]), including adsorption and electrostatics ([@B26]), steric and brush effects ([@B27]) and labeling ([@B19],[@B28]), come into play is not yet entirely clear. In this study, we postulate that the available area of probe spots is the limiting factor in adsorption, so that the binding is described as where is the number of target molecules going into the exposed probes over a unit of time, *N*~A~ the Avogadro constant, *V* the volume of hybridization solution. On the right side, α is the fraction of probes bound to specific targets, β the fraction of probes bound to cross-hybridizing targets. *p* is the total number of probes in unit of molar concentration (for simplicity, as if they were dissolved in the hybridization solution).

On the other hand, the dissociation is described by where is the number of target molecules leaving target/probe duplexes over a unit of time; *k*~d~ and *k*~n~ are dissociation rates for specific targets and cross-hybridizing targets, respectively.

At equilibrium between binding and dissociation, Equilibrium is established for both specific and cross-hybridizing targets. The proportions of specific targets and cross-hybridizing targets are determined by their concentrations: where is the concentration of free specific targets, the concentration of free cross-hybridizing targets.

Equations ([6](#M6){ref-type="disp-formula"}) and ([7](#M7){ref-type="disp-formula"}) can be combined to express β as: Then, Equations ([5](#M5){ref-type="disp-formula"}) and ([8](#M8){ref-type="disp-formula"}) give the fraction of specific binding Here , the concentration of free specific target molecules, is less than nominal spike-in concentration by the amount of probe binding: with *T̂* as the nominal spike-in concentration (total amount).

We assume the concentration of cross-hybridizing targets, , is large and can be treated as constant in this model. Let then Equation ([9](#M9){ref-type="disp-formula"}) becomes An analytical solution of Equation ([12](#M12){ref-type="disp-formula"}) is where

It can be shown that the other analytical solution of Equation ([12](#M12){ref-type="disp-formula"}), which bears a plus sign before the square root, has no valid physical meaning and merits no further discussion.

So α, the fraction of probes bound to specific targets, is described by three global parameters: *p*, *k*~b~ and γ, one probe-specific parameter *k*~d~ and one variable *T̂*. *k*~d~ can be expressed as: where *R* is the molar gas constant, *T* absolute temperature (318 K in the Affymetrix hybridization experiments), the energy computed from NN model, ξ as a scaling factor to account for binding to immobilized probes.

The physical meaning of our model is clear. Both specific binding α and cross-hybridization β compete for the same probe sites. As a result, high affinity probes (small *k*~d~) can achieve a higher fraction of specific binding, while low-affinity probes (large *k*~d~) saturate at a lower fraction. γ serves as a cross-hybridization factor. We made assumptions that are important to real experimental settings: a large quantity of cross-hybridizing targets are present; *k*~b~ is uniform for all targets and the adsorption is limited by the available area of probe spots. These assumptions make our model fundamentally different from previous competitive kinetic models ([@B29],[@B30]).

Experimentally, signal intensity is what is observed after washing, where most of cross-hybridized targets have been washed off: where is the observed signal intensity, τ the residual intensity from cross-hybridized targets, ι scanner bias, *A* the detection coefficient of fluorescence. As the unit of signal intensities is arbitrarily digitized, it only comes to a physical meaning through *A*.

Explanation to the correlation
------------------------------

First of all, we shall demonstrate that our model is capable of explaining the correlation at high target concentration in [Figure 1](#F1){ref-type="fig"}.

Equation ([12](#M12){ref-type="disp-formula"}) can be rearranged to a logarithmic form:

Note that the second item on the right side still contains the probe-dependent variable α. However, at high target concentration, the bound targets are minor comparing to free targets. This means, , and . Hence, Equation ([17](#M17){ref-type="disp-formula"}) at high target concentration is approximated as:

In a long range for , a linear approximation can be drawn between and . As in [Figure 2](#F2){ref-type="fig"}, Figure 2.Linear relationship between log (α/(1 − α)) and logα.

Combining Equations ([15](#M15){ref-type="disp-formula"}), ([18](#M18){ref-type="disp-formula"}) and ([19](#M19){ref-type="disp-formula"}), we get

At high target concentration, both cross-hybridization and scanner bias can be neglected. Therefore Equation ([16](#M16){ref-type="disp-formula"}) can be simplified to . We substitute the α in Equation ([20](#M20){ref-type="disp-formula"}) with : where , a constant for fixed *T̂*. Thus, is inversely correlated to . The observed correlation is explained by our competitive hybridization model. At low *T̂*, the premise is less valid; as a result, is less correlated to . A similar effect may be created by a very low , where a large fraction of targets is bound to probes and taken out of solution.

A bonus here is the determination of ξ. Since the coefficient for in Equation ([21](#M21){ref-type="disp-formula"}) should equate the slope in [Figure 1](#F1){ref-type="fig"}, we get .

Procedure of fitting model to Latin Square data
-----------------------------------------------

In DNA microarray experiments, signal intensities are measured in place of fluorescent densities of bound targets. However, common photomultiplier tube scanners usually carry a significant nonlinearity for low signal intensities ([@B31]). This means, the lower end of these Affymetrix data may deviate from the true fluorescent densities, a problem difficult to correct without knowledge of the specific instrument calibration data. And the signals from targets below 1 pM are hardly distinguishable from backgrounds, therefore, data from spike-in concentration 1 pM and above are used for our modeling.

The model fitting is to match the theoretical calculation of signal intensity, , to the experimentally observed counterpart is defined from Equation ([16](#M16){ref-type="disp-formula"}): Here, the background levels τ are observed values in these Affymetrix data (signal intensities at zero spike-in concentration). With background τ subtracted, the signal intensity should approach zero when the target concentration approaches zero. However, there is usually a deviation from zero that is known as scanner bias, ι, which can be therefore estimated by extrapolating the signal intensities at low target concentrations. For the data in this study, is taken. This value of ι is relatively small and has no significant effect on our model parameters. Though it is useful for stabilizing the small numbers in the fitting process.

With the theoretical value

Equation ([13](#M13){ref-type="disp-formula"}) can be written as where Γ is defined in Equation ([14](#M14){ref-type="disp-formula"}). In this equation, is known, the observed value for , and *k*~d~ can be calculated from Equation ([15](#M15){ref-type="disp-formula"}). So we only need to fit four global parameters: *A*, *p*, *k*~b~ and γ.

We use a fitness function of weighted squares \[similar to ([1](#M1){ref-type="disp-formula"})\]. For a probe *i*, the fitting error is calculated as where is observed signal intensity, the calculated value by Equation ([24](#M24){ref-type="disp-formula"}), *t* one of the nominal target concentrations (1--512 pM). Our model in Equation ([24](#M24){ref-type="disp-formula"}) is fitted to the training data by minimizing the sum of *E~i~* through brute-force searches as heuristic ranges of the four parameters can be obtained based on their physical meanings.

A useful constraint to the fitting is the value of . This is the signal intensity in Equation ([23](#M23){ref-type="disp-formula"}) when , often referred as the saturation level of hybridization. It is obvious that *P*~0~ should be larger but not infinitely larger than the highest signal intensity observed in the experiment. When varying *P*~0~ is used, as shown in [Figure 3](#F3){ref-type="fig"}, the overall fitting results are not very sensitive to *P*~0~ beyond a certain value. So we choose here. Figure 3.The fitting result is not sensitive to *P*~0~ beyond a certain value. All four parameters, *A*, *p*, *k*~b~ and γ, are fitted simultaneously. P~0~ = *A* · *p*; the fitting error is computed as in Equation ([25](#M25){ref-type="disp-formula"}).

Probe signal intensities can be successfully modeled
----------------------------------------------------

[Figure 1](#F1){ref-type="fig"} shows that (hence *k*~d~) can be reasonably approximated by the current NN model for the probes free of secondary structures. We use half of these probes to fit our competitive hybridization model, and determine the four global parameters, *A*, *p*, *k*~b~ and γ. The evaluation is then performed on the rest of probes.

The results indicate that our model captures the probe properties well. [Figure 4](#F4){ref-type="fig"}A shows the modeling of individual probe signals on both the training data and evaluation data. Overall, the prediction on training data has ([Figure 4](#F4){ref-type="fig"}B), and on evaluation data ([Figure 4](#F4){ref-type="fig"}C). If we relax the probe selection criterion to , about 75% of total probes are included, with prediction ([Figure 4](#F4){ref-type="fig"}D). The rest 25% of probes, which are presumably under stronger influence of secondary structures, can still be modeled with the same parameters but less accuracy at . Figure 4.Our model can successfully predict probe signal intensities. (**A**) The prediction on randomly chosen probes at random target concentrations. Top: the training data; Bottom: the evaluation data. (**B**) Scatter plot of all training data. (**C**) Scatter plot of the evaluation data. The training data are consisted of half of the probes from [Figure 1](#F1){ref-type="fig"}, and evaluation data from the other half. (**D**) Extended evaluation on 266 (75% of total) probes that satisfy Δ *G*~o~ \> −2.5. All signal intensities are in log scale. The parameters in this figure are *A* = 33.408 (pM)^−1^, *p* = 898 pM, *k*~b~ = 1.348*E*-3 *s^-1^* and γ = 245 500 pM s.

In the previous, heavily parameterized models, the best prediction on was correlation coefficient in Ref. ([@B14]) and in Ref. ([@B13]). In comparison, our model of four parameters produces for all probes, and for 75% probes after a preliminary selection by secondary structures (i.e. [Figure 4](#F4){ref-type="fig"}D). In conclusion, our competitive hybridization model can not only predict probe signals successfully, but also opens up paths to future improvements.

Prediction of target concentrations
-----------------------------------

With the four global parameters, target concentration can be calculated from Equation ([12](#M12){ref-type="disp-formula"}):

If we substitute ,

Since *k*~d~ calculation is more accurate for probes free of secondary structures, we focus on 19 out of the 30 probesets (transcripts) in this study that have five or more probes with . For these transcripts, Equation ([27](#M27){ref-type="disp-formula"}) is applied to calculate a target concentration from each probe. And the final concentration of a transcript is taken as the median of the data from its probes ([Figure 5](#F5){ref-type="fig"}A). [Figure 5](#F5){ref-type="fig"}B shows the prediction at gene level for all 19 transcripts. In fact, comparable results can be obtained by using the few probes with alone. At low concentrations, the predicted values in [Figure 5](#F5){ref-type="fig"}B tend to be higher than the nominal concentrations. We think this is likely to be a reflection of scanner nonlinearity in the low signal range, which can be corrected by an instrument calibration. Figure 5.Prediction of transcript concentration. (**A**) Example of the 11 probes for transcript 205267_at. Dots are probes with Δ*G*~o~ \> −1, other probes in crosses (slightly shifted horizontally for clarity). The transcript concentration (dashed line) is taken as the median value of all probes. (**B**) Prediction of 19 transcripts that have five or more probes with Δ*G*~o~ \> −1. Correlation coefficient between nominal concentrations and the predictions is *r* = 0.89. Error bars are standard deviations of the 19 transcripts. The predicted values bend away from the ideal line (dashed) at low concentrations probably because of scanner nonlinearity.

DISCUSSION
----------

In DNA microarray experiments, systematic variations stem from sample preparation and instrument operations. They are likely reflected in the global parameters of our model, *A*, *k*~b~ and γ. Therefore, batch variations can be expected in these parameters. The highest signal intensity in the Latin Square data was about 16 000. Comparing with the saturation level , this means about half of those probe sites are bound to specific targets at Thus, the fitted value of (as in [Figure 4](#F4){ref-type="fig"}) seems to be reasonable. Since our model has only four global parameters, they can be easily calibrated if control probes are built into array design. For instance, a set of targets complementary to the control probes can be spiked into the hybridization at various concentrations. Signal intensities of the control probes along with the known target concentrations can then be used to calibrate our model every time a hybridization experiment is performed.

We would like to emphasize that *k*~d~ is the only probe-specific factor in our model, and therefore plays a pivotal role in model accuracy. The accuracy of *k*~d~ or in this article is limited by the NN model, which is only a coarse approximation and affected by probe/target secondary structures. This can be improved but beyond the scope of this current study.

We assumed a constant cross-hybridization factor γ for all probes, which may not be the case. Further research on γ may improve the accuracy of our model. We did not deal with the background levels in this study, which are not important to signals at high target concentration but will affect signals at low concentrations. Background levels have a clear dependency on , and are well addressed in other studies ([@B32],[@B33]).

CONCLUSION
----------

Our study presents the first model of DNA microarray hybridization that explains probe signal intensities through sequence-based thermodynamic properties without excessive parameter fitting. This fills in the long standing knowledge gap in DNA microarray hybridization. Our model provides a mechanism of absolute quantification, and shall improve the quality control and reproducibility of the technology. With only four global parameters, this model can be easily calibrated through control features that are built into microarrays, and adopted in practice. We expect new design and quantification algorithms to take advantage of our results.

FUNDING
=======

National Oceanic and Atmospheric Administration (NA05NOS4261163 and NA06NOS42600117).

*Conflict of interest statement*: None declared.
