BACKGROUND
presently, there are several high throughput methods of quantifying changes in gene expression including oligonucleotide microarrays, quantitative realtime pcr  and "next generation sequencing" .
of these, high density oligonucleotide microarrays are arguably the most important tools for genomic investigation.
although next generation sequencing is a promising alternative to microarrays for genome-scale expression profiling, and exhibit more sensitivity in the low-expression limit  <cit> , microarray technology is substantially less expensive and the resulting data sets require much less information processing.
moreover, microarrays have substantially higher throughput than qpcr.
microarrays consist of dna probes  which are orderly arranged on a glass slide.
probes may be attached to  the slide surface via  mask-dependent  <cit>  or maskless photolithographic dna synthesis technology  <cit> , or  robotic printing of pcr products or synthetic oligomers  <cit> .
the first two of these methods yield oligonucleotide arrays , which are more reliable than pcr product-based arrays; they comprise the majority of commercial arrays by market share.
thus, we restrict our considerations to microarrays of these types.
each probe is designed to hybridize with a specific cell-derived and fluorolabeled dna species such as cdna or crna  <cit> .
if these targets originate from mrna, then the fluorescence associated with each microarray feature is assumed to be proportional to the amount of each transcript  <cit> .
however, in recent years, this assumption has been called into question.
several studies have shown that different microarray assays yield different results when used to quantify differences in expression .
moreover, significant differences have also been reported among the results of microarray and qpcr assays  <cit> .
although the microarray quality control consortium confirmed a "high level of interplatform concordance in terms of genes identified as differentially expressed" for several commercial microarrays targeting the human transcriptome in  <dig>  <cit> , the reliability of many other array designs and experimental protocols have not been characterized systematically.
theoretical analyses of smaller systems have suggested that problems with probe reliability are thermodynamic in origin  <cit> .
sequence and gc content heuristics are commonly employed in the design of probes  <cit> .
however, these heuristics do not preclude the possibility of finite lengths of complementary sequence between probe candidates and non-target ssdna.
consequently, all probe-ssdna interactions will exhibit favorable energetics of interaction to some extent, and there will always be a finite amount of cross-hybridization between probes and non-target cdna.
however, experimental identification and quantification of the relative amounts of correctly-hybridized and cross-hybridized probes is impractical, since the only measurement of an array reader is the fluorescence associated with each microarray feature.
chemical dynamics of dna hybridization
however, the relative amounts of these species may be quantified using population-balances combined with an appropriate thermodynamic model of hybridization.
consider a microarray with np oligonucleotide probes that target many if not all of the nt solution-phase ssdnas  the hybridization network may be written as
where ℓ ∈ .
most arrays have multiple reporters per target  <cit> .
the deterministic time evolution of the process is completely specified by the following chemical population balance equations
in eq.
<dig>  v is the volume of cdna solution added to the array,  is the population of unhybridized probes of type m ,  is the population of unhybridized transcripts of type ℓ , and  is the population of hybrids composed of probe m and transcript ℓ. before the addition of cdna to the array, there are n unhybridized probes per feature and  molecules of each cdna ℓ. thus, the initial conditions of this system of np + nt + np nt ordinary differential equations are   = n,  =  and   =  <dig>  numerical solution of eq.
<dig> until t ~∞ yields the equilibrium populations of desired hybrids  as well as cross-hybrids.
unfortunately, this deterministic approach exhibits certain practical pitfalls.
since the reaction rate constants and cdna populations vary over many orders of magnitude, eqs.
<dig> are "stiff".
the size of the hybridization reaction network compounds the problem; typical genomic assays are designed to measure the expression levels of thousands of transcripts.
for example, baker's yeast  possesses approximately  <dig>  genes   <cit>  and humans possess approximately  <dig>   <cit> .
since most microarrays feature one or more distinct reporters  for each target, the size of the hybridization network  will be enormous.
for these genomes, eq.
<dig> represents millions to billions of stiff ordinary differential equations.
stochastic simulation
alternately, one may utilize the stochastic approach to chemical kinetics, which underlies the aforementioned rate equations  <cit> .
to begin, let us consider the general case where a volume of solution containing nt cdnas at populations   is added to an array with np surface-bound probes at populations  = n , where again, the superscript "0" denotes the initial state.
upon addition of the solution to the slide, unhybridized probes and targets will randomly hybridize in accordance with eq.
<dig>  assuming perfect mixing and isothermal hybridization - both ostensibly achieved with most assays - the state of this system may be defined in terms of the populations of unhybridized probes , unhybridized target dna , hybrids , and the hybridization volume, v. the probability of a transition from one state to another is defined by the stoichiometry of eq.
<dig>  since ssdna molecules directly interact to form a dsdna hybrid,
is the probability that probe m will hybridize with cdna ℓ within the imminent time interval δt, and
is the probability that any of the  hybrids composed of probe m and cdna ℓ will dehybridize within the imminent time interval δt.
we recognize  and  as the rates of the forward and reverse reactions of eq.
<dig>  respectively.
eqs.
<dig> and  <dig> are microphysically valid  <cit>  and are in fact the basis of the validity of the aforementioned rate equations .
we refer interested readers to gillespie's paper on this subject  <cit> .
eqs.
<dig> and  <dig> are the bases of the stochastic simulation algorithms   <cit> .
ssas simulate the time evolution of a chemical process by repetitively  selecting a reaction μ among a set m of potential reactions,  selecting the time τ until that reaction occurs,  updating the state of the system to reflect the occurrence of the selected reaction, and  updating the time.
exact ssas differ only in how the first two steps are implemented.
each methodology features its own memory and speed enhancements and restrictions.
direct method
the direct method samples two exact density functions to obtain the quiescence time and imminent reaction event.
the selection rule for the quiescence time is
where a <dig> is the sum of all reaction rates and r <dig> ~u  .
the selection rule for the imminent event μ is
in a microarary hybridization network, these rates are defined by eqs.
<dig> and  <dig>  as discussed, such that μ ∈ .
although the direct method is the faster and more memory efficient of the two ssas first developed by gillespie  <cit> , it is numerically intensive when applied to large chemical networks.
since o operations are required for eq.
<dig>  simulations with m reaction types require o calculations per time step.
since each probe may potentially bind every target, m = 2np nt for the reaction system described by eq.
<dig>  thus, there are o operations per time step with the direct method.
considering typical values for np and nt, there will be hundreds of millions of operations per time step in dm simulations of conventional microarrays.
next reaction method
in  <dig>  gibson & bruck proposed a new exact ssa called next reaction method.
this approach purportedly reduces the number of required calculations per time step from the o of the direct method to o, where k is tantamount to the number of chemical species with which the average chemical species will react  <cit>  .
the speed enhancement of this algorithm is most prominent when the reaction network is sparse , however, it should be much faster than the direct method for reaction networks as coupled as eq.
<dig>  where k ~np and m = 2npnt .
the next reaction method is significantly different than the direct method in both its data handling and mc selection rules.
first, the absolute times at which all reactions might occur  ~ exponential, ν ∈ ) are selected by mc, by contrast to the mc selection of just one quiescence time in the direct method.
one may show that the smallest of these times  is the time at which the next reaction  occurs, and that this selection rule is an exact mc selection from the reaction probability density function, like eqs.
<dig> and  <dig>  <cit> .
another noteworthy difference between the direct method and the next reaction method is the employment of a "dependency graph" by the latter.
this data structure reduces the number of calculations per reaction event, both for event selection and updating the reaction times.
however, the dependency graph requires o objects to store the dependencies of reaction rates upon the populations of reactants and products shared with other reactions.
for microarray hybridization, this translates to a storage requirement of o() reaction objects, which is prohibitively large for genome-sized microarray simulations; even current supercomputers with terabytes of memory cannot meet these requirements.
RESULTS
algorithm
in simulations of microarray hybridization, the extreme computational burden of the direct method and the memory burden of the next reaction method may be alleviated by judicious storage and summation of the terms in eq.
<dig>  our algorithm employs a data structure called a hybridization table  that stores partial sums of the terms in eq.
<dig>  for this reason, we call our approach the "method of partial sums" .
initialize
calculate the auxiliary variables , ϕj, j =  <dig>   <dig> ... np , α and ϕ
repeat
calculate the total reaction rate α <dig>  and quiescence interval τ
if the next reaction is a hybridization  then
select the hybridizing probe m and cdna ℓ
else
select the dehybridizing probe m and cdna ℓ
end if
update populations
update hybridization rates: αm and αj ≠ m, j =  <dig>   <dig> ... np
update dehybridization rates: ϕm, ϕ
t ← t + τ
until p > <dig>  for  =
algorithm 1: method of partial sums for the simulation of the coupled reaction network composed of the hybridizations of np oligonucleotide probes and nt cdna species.
hybridization table
the structure of the hybridization table is outlined in fig.
<dig>  in addition to storing the populations of targets {}, probes {}, hybrids {} and rate constants  and  , the table also contains the rates at which probes of type j will hybridize with any target,
and the rates at which all hybrids composed of probes j will dissociate,
two additional quantities are stored in the table: the total rate of hybridization
and the total rate of dehybridization
the total rate of reaction may then be calculated from
eqs.
<dig>   <dig>   <dig>  and  <dig> effectively subdivide the running sum in eq.
<dig> into independently manageable groups of information.
this in turn reduces the number of operations required for event selection reaction rate updates.
reaction event selection
the selection of the imminent hybridization or dehybridization event is accomplished by summing the terms of eq.
<dig> such that the total number of operations in eq.
<dig>  is minimized.
we begin by ordering the reactions according to column and row in the ht  such that .
eq.
<dig> is then implemented by summing the reaction rates corresponding to each entry in the ht, by column and then by row: first, by hybridization events, and then by dehybridization events.
one begins by identifying whether or not the imminent event will be a hybridization.
since the imminent event is defined by the reaction whose rate causes the sum of αν to exceed the quantity r2α <dig>  and all hybridization reactions precede the dehybridization reactions in the ht, it follows that
if eq.
<dig> results in the selection of a hybridization event, α >r2a <dig>  and one needn't consider the dehybridization rates in the selection of the event to come.
the next step is the selection of the probe  that will hybridize in the imminent event.
using the notation of the direct method and the order of summation, the index of the event to come  must be less than or equal to np nt, where  is the last of the hybridization rates in the table .
eq.
<dig> becomes
in essence, the event is selected by summing the rates corresponding to each column of the "hybridization" section of the ht, column by column, followed by row, until the quantity r2a <dig> is exceeded.
noting that the quantity is equal to αj , eq.
<dig> may be simplified to
which defines the probe that will hybridize in the imminent event.
the equation defining the selection of the target may be derived similarly first, we express the sums on the right and left sides of eq.
<dig> as
and
simplifying these expressions, we obtain
thus, the target  in the imminent event may be obtained by summing the quantities  until the quantity  is exceeded.
if the event to come is a dehybridization, then α <r2a <dig>  hence, the sum of all hybridization rates may be subtracted from each term in eq.
<dig>  leaving
where {aν } are the rates of the dehybridization reactions and we have explicitly noted the fact that all np nt hybridization rates are contained within α. expressing this in terms of the hybrid indices, we obtain
again, the selection of the the event may be simplified using the auxiliary variables in the ht.
if the probe m is released in the imminent event, the value of m may be defined by sequential addition of the values of {Φj} until the quantity  is exceeded
note that this also defines the identity of the probe that will dissociate from the target.
the cdna to be released in the dehybridization event  may be calculated by subtraction of eq.
<dig> from eq.
19:
the selection rules for hybridization  and dehybridization  substantially reduce the number of operations required by eq.
<dig>  whereas the direct method may require 2np nt operations to select a reaction, our selection rules require at most  operations.
for a genome-sized microarray simulations, with np ~nt ~  <dig>  our reaction selection approach will be several orders of magnitude faster than those of other algorithms.
system state accounting
upon selection of the imminent event, the populations of the species involved  must be updated in accordance with the stoichiometry of the reaction.
the structure of the ht facilitates this procedure: the row and column of the selected event designate the identities of both the reactants and products as well as the stoichiometric changes in population.
however, other quantities must be updated, most notably the partial sums of reaction rates that facilitate event selection.
the first quantity to be updated is the hybridization rate of the probe m
where i = - <dig> for hybridization and + <dig> for dehybridization.
since each partial sum αj is a function of the population of the ℓth cdna , these quantities must also be updated:
subsequently, a must be recalculated using eq.
<dig>  collectively, 2np operations are required for the update of the hybridization section of the ht.
after a hybridization or dehybridization event, only one of the dehybridization rates must be updated.
thus, only one partial sum requires modification:
moreover, inasmuch as Φm is the only affected partial sum, Φ may also be updated in one operation,
thus, only two operations are required to update the dehybridization section of the ht.
in summary, after selection of a reaction, the ht may be updated in o operations, a substantial improvement over the direct method.
this is largely a result of the fact that reaction rates  are not stored in the ht, precluding the need to update np rates of events featuring probe m and nt rates featuring target ℓ. the same argument applies to the dehybrization reaction rates {Φj}.
ultimately, the method of partial sums is efficient so long as the partial sums can be updated easily.
determination of equilibrium
microarray analysis is predicated upon the assumption that the probes and solution-phase cdna have equilibrated prior to scanning the slide and measuring the fluorescence associated with each feature.
we follow this experimental convention in silico.
the hybridization process is at equilibrium when the rates of change of all chemical populations in eq.
<dig> are zero, implying
as straightforward as this criterion appears, it is difficult to employ in practice.
eq.
<dig> represents millions of comparisons for genome-scale microarrays, requiring o operations per time step.
this many operations would also be required if one determined the steady states of all molecular species, which would additionally require storage of the current state as well as previous states.
clearly, this is memory-prohibitive.
furthermore, eq.
<dig> is exact only on average  <cit> .
hence, it will never be exactly satisfied at any point in time within a single simulation.
to circumvent these issues, we propose an alternate approach that employs the average total rates of hybridization  and dehybridization .
considering the definitions of these quantities  summation of eq.
<dig> over all ℓ and m yields the result that α = Φ at equilibrium.
this criterion may be established using student's t test for two populations with unknown means and standard deviations  <cit> .
strictly speaking, this is necessary but insufficient for the specification of thermodynamic equilibrium.
however, it is remarkably effective as a heuristic.
we implement it as follows:
repeat
save α and Φ to disk every o reaction steps
maintain the last ten saved values of α and Φ in memory as vectors α and Φ
keep running averages of the numbers in these vectors,  and .
if  <then
if p <  <dig>  for h0:  =   then
reinitialize α and Φ
end if
end if
until p > <dig>  for h0:  =
testing
heuristic for the determination of equilibrium
to evaluate the heuristic approach to the determination of the equilibrium state, we performed simulations of the hybridization of modified agilent probes with yeast cdna .
a typical example of the time evolution of a and f is presented in fig.
<dig>  where the time at which equilibrium is attained  is marked by the red circle.
at and after this point, the populations of all hybrid species were at steady state, fluctuating in accordance with the predictions of equilibrium statistical mechanics.
subsequent comparison of hybrid populations with eq.
<dig> confirmed the findings of our "equilibrium heuristic" for all simulations.
in standard hybridizations, imperfect mixing causes the transport of cdna or crna to be diffusion limited  <cit> .
this in turn presents an obstacle to hybridization, as the time required for a  target to diffuse to it's probe is large.
this, in turn affects the kinetics.
the most common solution to this problem is to increase the concentration of target cdna or crna, which results in an abundance of these molecules at equilibrium.
as we have discussed, our simulations feature perfect mixing.
thus, we may use substantially less solution-phase cdna.
as a result, free cdna is sparse at equilibrium, which introduces fluctuation into a. by contrast, f exhibits little fluctuation because hybrid populations are fairly large compared to the change in population accompanying a  reaction.
effect of reaction rate constants upon the steady state
as discussed, the values of the kinetic rate constants should not affect the equilibrium state of our simulations provided that their ratios are constant .
this assumption may be formulated as a testable hypothesis as follows: if the values of  and  affect the equilibrium state of simulation, then they will affect the fractional occupancy of each probe m ∈  defined by
in this expression
is the total population of probes m hybridized with any type of cdna at equilibrium.
we test these hypotheses at the system level by performing simulations using rate constants
that is, we perturb the rate constants generated via eqs.
<dig> and  <dig> by multiplying their results by a uniform random number rℓm.
since a unique random number is generated for each probe m and cdna ℓ, one may quantify the effect of kinetic perturbations on the equilibrium state of a simulation via the hypothesis
if kinetic rate constants significantly affect any of the fractional occupancies at the equilibrium state , then this hypothesis will fail.
simulations of the hybridization of all  <dig> modified agilent probes with  <dig> full-length yeast cdna molecules were conducted for two types of perturbations.
in the first study, r ~u , where u  is a uniform random number on the interval .
this allowed us to evaluate the effect of the widest variation of the rate constants.
in the second study, r ~u .
rate constants for all  <dig> , <dig> hybridization events were independently perturbed.
the results of our hypothesis tests are illustrated in fig.
<dig>
in both cases, our results clearly indicate that the equilibrium state is unaffected by the values of the rate constants, as expected.
interestingly, the time required to reach equilibrium correlates with the heterogeneity of the dehybridization rate constants: deviations of the results for the two cases at  <dig> hours of cpu time  indicate that the hybridization of many probes with solution-phase cdna was incomplete in cases where their rate constants were perturbed by factors less than 10- <dig>
analyses of the timeseries of the overall rate of reaction  revealed that the progression to equilibrium is considerably slower if r ~u  than if r ~u .
this conforms with the experimental observations of dai and coworkers  <cit> , which demonstrated differences between the kinetics of specific and nonspecific hybridizations.
comparison of stochastic simulation algorithms
results of all ssas applied to the process described by eq.
<dig> should yield statistically indistinguishable results since they share a common stochastic process.
in his seminal work  <cit> , gillespie showed that the "first reaction method" and "direct method" were equivalent.
subsequently, gibson and bruck demonstrated that their "next reaction method" is equivalent to gillespie's algorithms.
since the method of partial sums shares the mathematical underpinnings of the direct method, its results should be indistinguishable from those generated by the next reaction method or gillespie's algorithms, as well as the law of mass action.
we initially tested this hypothesis by conducting simulations of the hybridization of ten full length cdna species from yeast to ten probes for those species.
five simulations were conducted via both the method of partial sums and the next reaction method.
additionally, we solved the corresponding population balance equations  for this illustrative hybridization process.
our results  clearly show that all methods yield equivalent results.
average populations for all  <dig> hybrid species could not be distinguished by simulation or calculation method via the t-test .
interestingly, several common differential equation solvers could not integrate eqs.
<dig> for this ten by ten system from t =  <dig> until steady state, including the matlab packages ode <dig> and ode <dig>  the two hundred equations necessary to model this small array required the use of the ode15s, which is designed for stiff sets of odes.
by contrast, the stochastic simulation algorithms were unimpeded in their numerical progress.
although the stochastic simulation algorithms give indistinguishable results for the population states both in time and at equilibrium, their computational performances are significantly different.
the differences in the computational speeds of the ssas were evaluated by conducting simulations of the hybridization of full length yeast cdnas to microarrays featuring  <dig>   <dig>   <dig>   <dig>  and  <dig> probes targeting those cdnas .
additionally, a simulation with  <dig> probes and  <dig> targets was conducted with the method of partial sums.
the same initial populations of targets and the same set of probes  were used in both the next reaction and mps simulations.
estimation of the time at which equilibrium is attained was determined as discussed, and then used in simulations conducted with the next reaction method.
our results are illustrated in fig.
<dig>  and the contrast is stark.
the method of partial sums outperforms the next reaction method for simulations of microarrays of all sizes.
it requires one hour to complete a simulation for an array large enough to be used for genomic characterization.
by contrast, the next reaction method is not capable of performing such simulations due to its memory requirements.
if np ~nt ~  <dig>  the pointers required by the next reaction method will consume approximately three terabytes of  <dig> bit computer memory by themselves; the ht requires no such pointers.
next reaction simulations with as few as  <dig> probes required  <dig> hrs of cpu time to reach equilibrium, and for a computer with extraordinary memory, we forecast that simulations using the next reaction method would require a month to simulate the hybridization to the agilent yeast array.
in addition, the scaling of the cpu time with respect to the number of probes differs among the two stochastic simulation algorithms.
both algorithms feature power-law scaling, however, the method of partial sums scales as , whereas the next reaction method scales as  .
these differences arise due to the differences in the data structures underlying the two methods, which in turn affect the number of calculations required per time step.
our algorithm permits two concurrent simulations of an  <dig> million-reaction system  to reach equilibrium within two hours using  <dig>  ghz dual g <dig> apple servers with  <dig> gb memory.
<dig> gb of memory are sufficient for simulations featuring  <dig>  probes, e.g.
simulations of yeast arrays with one perfect match probe and one mismatch probe for each transcribed gene.
simulations of an array designed for the human genome with one probe per gene - at approximately  <dig>  genes - possess hybridization reaction spaces approximately  <dig> times larger, and an estimated run time of  <dig> -  <dig> hours based on cpu time scaling of t ~.
such simulations can be achieved on most university shared-resource machines , which commonly feature hundreds of gigabytes of ram.
illustrative example: characterization of cross hybridization
the in silico characterization of cross hybridization is only as sensitive as the number of probe molecules per feature.
agilent arrays have  <dig>  ×  <dig> probe molecules per feature, facilitating the hybridization and measurement of as many target cdnas and giving, in principle, as much resolution.
however, the time required for stochastic simulations is proportional to the number of molecules therein  <cit> .
balancing these resolution and population considerations, we selected a hybridization volume of  <dig>  nl corresponding to initial probe populations of  <dig> molecules/feature and concomitant scaling of the feature diameter to maintain the surface concentration.
we also employed a total concentration of  <dig> ng per  <dig> μl hybridization volume.
at this concentration, only a few if any probes become saturated.
this concentration is a tenfold dilution of that recommended by agilent's protocol, however, it is is within the range of commercial oligonucleotide microarray protocols  <cit> .
moreover, the agilent array effectively has a probe population four times higher than the one we used, as it has four redundant features per probe .
in fig.
<dig> we illustrate a small subset of the hybrid populations predicted from our simulations.
the average populations of hybrids were calculated from five replicate simulations, each of which required  <dig>  hours of cpu time as discussed .
the complete sets of results are provided in additional files  <dig> and  <dig>
many of the cross-hybridizing cdna species are exceptionally homologous.
for instance, a_06_p <dig> - the probe designed to target cdna for ypl281c  also hybridizes with cdna for yor393w ; in this case, the two orfs are identical.
the third hybrid, ymr323w , shares all but twelve bases of the other two.
other cross hybridizing cdna species have less overlap in sequence, but share the probe sequences completely or in part.
all cross-hybrids listed in fig.
<dig> share at least 90% of the target sequence.
a probe's proclivity to cross hybridize in the presence of thousands of potentially competitive cdnas may be expressed in terms of its selectivity, s, defined as the fraction of fluorescence intensity associated with a microarray feature that originates from the corresponding target.
mathematically,
where  is the population of hybrids composed of probe m and its target ), and  is calculated via eq.
<dig>  for example, the selectivity of the aforementioned agilent probe a_06_p <dig> is  <dig> %.
in fig.
<dig>  we present a summary of the selectivities for all probes in the agilent set.
as these results show, the vast majority of the probes for this commercial microarray are selective , and do not exhibit cross hybridization when they bind yeast cdna at the concentrations specified by the expression state.
for other microarrays and probe sets designed by various publicly available software packages, the cross hybridization may be more extreme.
since the populations of all hybrids depend upon the populations of potentially cross-hybridizing cdna molecules, the value of s for each probe depends upon the expression state of a cell.
interestingly, the distribution of s among the probes is independent of the expression state.
simulations with different initial conditions corresponding to four additional mc-generated expression states yielded selectivity distributions that were indistinguishable from the aforementioned distribution .
this result suggests a method of characterizing the overall reliability of a proposed probe set.
since the distribution of selectivities is invariant with respect to the initial target populations , a statistic that characterizes that distribution should be robust.
we propose that the average selectivity can serve as such a metric.
the use of such metrics should be used with care, however, inasmuch as they lend themselves to "ecological fallacy".
implementation
our software is implemented in c++ and executed on apple g <dig> processors running macos x .
unafold  <dig>  is available for download at the dinamelt server http://dinamelt.bioinfo.rpi.edu/ and compiles on a variety of operating systems.
additionally, we have implemented the next reaction method as described by gibson and bruck  <cit>  in c++ for the purposes of comparing its performance with that of our algorithm.
discussion
in recent years, advances in microarray manufacturing have opened the doors to custom microarray design.
individual researchers can now upload their own microarray probe sequences to one of many sites  and have a custom array manufactured.
in response, many probe design algorithms have been arisen to fill the needs of researchers intent on performing global investigations of gene expression  <cit> .
however, none of the probe sets generated by these algorithms have been evaluated in terms of their selectivity or proclivity to give a linear relationship between feature fluorescence and target concentration.
studies of this type carried out by the microarray quality control consortium  <cit>  were costly, involving dozens of participants.
such laborious quality control is not feasible for each and every probe set designed by a novel probe tool.
however, robust population-based simulations of the hybridization process may be employed to evaluate candidate probe sets, given robust estimates of the thermodynamic free energies of hybridization.
in the tests of our simulation algorithm, we have utilized equilibrium constants that were calculated via the nn model of allawi and santalucia.
this model predicts free-energies to within three significant digits of experimentally measured values for solution-phase hybridization  <cit>  and has been validated for oligonucleotide arrays in which the oligonucleotides are connected to the surface by linkers, making them more "solution like"  <cit> .
however, many studies have demonstrated that nn models do not accurately predict the thermodynamic properties of hybridization between solution-phase and surface-bound oligonucleotides .
electrostatic interference in conventional microarrays - which do not feature 3d linkers of the type employed by weckx and coworkers  <cit> , is a major reason for this discrepancy.
solution-phase na+ may shield the phosphate groups of both hybridized and single-stranded probes and targets  <cit> .
cations will also shield the negatively-charged glass surface to which the probes are bound by organizing the formation of a layer of counter ions.
the surface density the oligonucleotide probes also plays a role, partly via steric interference and partly via changes in the charge distribution at the glass-water interface  <cit> .
ultimately, surface electrostatic effects cause a location-dependent effect of mismatches  <cit> .
therefore, extra care must be taken when applying our algorithm to characterize cross hybridization in real microarray assays.
if equilibrium constants are calculated using nn models , corrections for the effect of the debye layer should be introduced at a bare minimum.
theoretical predictions of the effect of the resulting "debye layer" upon the melting temperature tm  <cit>  have been confirmed via experimental measurements  <cit> .
this correction may be applied to nn models for dna-crna duplexes as well , if simulations of crna-dna assays are to be conducted.
we have not explicitly considered additional effects of single-strand secondary structure  formation among full-length cdna or probes
nor, for that matter, have we considered the possibility of hybridization between cdna molecules in solution
the hybridization between probes, or folding thereof, is customarily not a substantial problem in oligonucleotide microarrays.
probe design software packages such as oligoarray  <dig>  routinely screen probes for secondary structure, and inter-probe hybrids are precluded due to spatial separation.
the remaining interactions among cdnas  serve only to sequester single-stranded cdna from the probes  <cit> .
inasmuch as this can only further degrade the sensitivity of probe-target interactions, we have restricted our algorithm to consider only probe-cdna hybridization dynamics.
hence, the method we have presented for probe and array characterization  gives the performance of a microarray under a "best case" scenario, even as it fully accounts for the interactions between oligonucleotide probes and full length cdna targets at the system-scale.
for the end user we note that simulation-based characterization of putative microarray probe sets requires no more information than that contained in mage-tab  or mage-ml  formatted array information  <cit> , which are required for publication of microarray data to ensure miame compliance.
the array design format  component of mage-tab contains all of the probe sequence information and its target, whereas the investigation description format  contains the experimental protocols, including the hybridization temperature and salt conditions.
additionally, the adf contains information regarding the relationship between reporter sequences and features: there may be multiple features with the same reporter sequence, or alternately, signals from several reporters may be combined to produce a signal for a single gene .
as the total cdna concentration is ostensibly included in the idf, and the individual cdna concentrations are randomly generated, the system-scale selectivity of any array can be computed via simulation from a proposed experimental protocol.
in this work, we have not explicitly considered the prediction of the time evolution of the hybridization process, focusing instead on the equilibrium state.
simulations employing our algorithm are most accurate when they utilize experimentally-determined rate constants  <cit> .
however, care must be taken to ensure simulations are conducted for the same surface densities and ionic strengths employed in the experiments employed to estimate the rate constants, for the reasons previously discussed.
given accurate estimates of the rate constants  and , the time required for microarray hybridizations may be estimated via the method illustrated in fig.
<dig>  as the time variable would represent the actual hybridization time.
CONCLUSIONS
in this work, we have developed an algorithm for the stochastic simulation of exceptionally large and complex probe-cdna hybridization reaction networks that underlie microarray assays.
using the method of partial sums in conjunction with the data structure we denote the "hybridization table", our algorithm requires o operations per reaction event.
this is substantially fewer than gillespie's direct method  operations per event) and the next reaction method of gibson and bruck ) operations per event).
moreover, our algorithm requires less the data storage than others, obviating the need for pointers that track of the dependencies of reactions.
for instance, the next reaction method requires  pointers for each of its 2np nt reactions, which can consume a vast amount of memory for genome-scale simulations.
as a result of these innovations, our algorithm permits system-level simulation of the complete reaction network composed of all potential probe-target hybridizations  without the need for high-performance computing.
furthermore, such simulations are now possible within a reasonable amount of time.
thus, given robust thermodynamic predictions of the free energies of dna hybridization, one may obtain a conservative estimate of the reliability of a candidate probe set in silico.

 presently, there are several high throughput methods of quantifying changes in gene expression including oligonucleotide microarrays, quantitative realtime pcr and "next generation sequencing". of these, high density oligonucleotide microarrays are arguably the most important tools for genomic investigation. although next generation sequencing is a promising alternative to microarrays for genome-scale expression profiling, and exhibit more sensitivity in the low-expression limit, microarray technology is substantially less expensive and the resulting data sets require much less information processing. moreover, microarrays have substantially higher throughput than qpcr. microarrays consist of dna probes which are orderly arranged on a glass slide. probes may be attached to the slide surface via mask-dependent or maskless photolithographic dna synthesis technology, or robotic printing of pcr products or synthetic oligomers. the first two of these methods yield oligonucleotide arrays, which are more reliable than pcr product-based arrays; they comprise the majority of commercial arrays by market share. thus, we restrict our considerations to microarrays of these types. each probe is designed to hybridize with a specific cell-derived and fluorolabeled dna species such as cdna or crna. if these targets originate from mrna, then the fluorescence associated with each microarray feature is assumed to be proportional to the amount of each transcript. however, in recent years, this assumption has been called into question. several studies have shown that different microarray assays yield different results when used to quantify differences in expression. moreover, significant differences have also been reported among the results of microarray and qpcr assays. although the microarray quality control consortium confirmed a "high level of interplatform concordance in terms of genes identified as differentially expressed" for several commercial microarrays targeting the human transcriptome in, the reliability of many other array designs and experimental protocols have not been characterized systematically. theoretical analyses of smaller systems have suggested that problems with probe reliability are thermodynamic in origin. sequence and gc content heuristics are commonly employed in the design of probes. however, these heuristics do not preclude the possibility of finite lengths of complementary sequence between probe candidates and non-target ssdna. consequently, all probe-ssdna interactions will exhibit favorable energetics of interaction to some extent, and there will always be a finite amount of cross-hybridization between probes and non-target cdna. however, experimental identification and quantification of the relative amounts of correctly-hybridized and cross-hybridized probes is impractical, since the only measurement of an array reader is the fluorescence associated with each microarray feature. chemical dynamics of dna hybridization however, the relative amounts of these species may be quantified using population-balances combined with an appropriate thermodynamic model of hybridization. consider a microarray with np oligonucleotide probes that target many if not all of the nt solution-phase ssdnas the hybridization network may be written as where ℓ ∈. most arrays have multiple reporters per target. the deterministic time evolution of the process is completely specified by the following chemical population balance equations in eq. v is the volume of cdna solution added to the array, is the population of unhybridized probes of type m, is the population of unhybridized transcripts of type ℓ, and is the population of hybrids composed of probe m and transcript ℓ. before the addition of cdna to the array, there are n unhybridized probes per feature and molecules of each cdna ℓ. thus, the initial conditions of this system of np + nt + np nt ordinary differential equations are = n, = and = numerical solution of eq. until t ~∞ yields the equilibrium populations of desired hybrids as well as cross-hybrids. unfortunately, this deterministic approach exhibits certain practical pitfalls. since the reaction rate constants and cdna populations vary over many orders of magnitude, eqs. are "stiff". the size of the hybridization reaction network compounds the problem; typical genomic assays are designed to measure the expression levels of thousands of transcripts. for example, baker's yeast possesses approximately genes and humans possess approximately. since most microarrays feature one or more distinct reporters for each target, the size of the hybridization network will be enormous. for these genomes, eq. represents millions to billions of stiff ordinary differential equations. stochastic simulation alternately, one may utilize the stochastic approach to chemical kinetics, which underlies the aforementioned rate equations. to begin, let us consider the general case where a volume of solution containing nt cdnas at populations is added to an array with np surface-bound probes at populations = n, where again, the superscript "0" denotes the initial state. upon addition of the solution to the slide, unhybridized probes and targets will randomly hybridize in accordance with eq. assuming perfect mixing and isothermal hybridization - both ostensibly achieved with most assays - the state of this system may be defined in terms of the populations of unhybridized probes, unhybridized target dna, hybrids, and the hybridization volume, v. the probability of a transition from one state to another is defined by the stoichiometry of eq. since ssdna molecules directly interact to form a dsdna hybrid, is the probability that probe m will hybridize with cdna ℓ within the imminent time interval δt, and is the probability that any of the hybrids composed of probe m and cdna ℓ will dehybridize within the imminent time interval δt. we recognize and as the rates of the forward and reverse reactions of eq. respectively. eqs. and are microphysically valid and are in fact the basis of the validity of the aforementioned rate equations. we refer interested readers to gillespie's paper on this subject. eqs. and are the bases of the stochastic simulation algorithms. ssas simulate the time evolution of a chemical process by repetitively selecting a reaction μ among a set m of potential reactions, selecting the time τ until that reaction occurs, updating the state of the system to reflect the occurrence of the selected reaction, and updating the time. exact ssas differ only in how the first two steps are implemented. each methodology features its own memory and speed enhancements and restrictions. direct method the direct method samples two exact density functions to obtain the quiescence time and imminent reaction event. the selection rule for the quiescence time is where a is the sum of all reaction rates and r ~u. the selection rule for the imminent event μ is in a microarary hybridization network, these rates are defined by eqs. and as discussed, such that μ ∈. although the direct method is the faster and more memory efficient of the two ssas first developed by gillespie, it is numerically intensive when applied to large chemical networks. since o operations are required for eq. simulations with m reaction types require o calculations per time step. since each probe may potentially bind every target, m = 2np nt for the reaction system described by eq. thus, there are o operations per time step with the direct method. considering typical values for np and nt, there will be hundreds of millions of operations per time step in dm simulations of conventional microarrays. next reaction method in gibson & bruck proposed a new exact ssa called next reaction method. this approach purportedly reduces the number of required calculations per time step from the o of the direct method to o, where k is tantamount to the number of chemical species with which the average chemical species will react. the speed enhancement of this algorithm is most prominent when the reaction network is sparse, however, it should be much faster than the direct method for reaction networks as coupled as eq. where k ~np and m = 2npnt. the next reaction method is significantly different than the direct method in both its data handling and mc selection rules. first, the absolute times at which all reactions might occur ~ exponential, ν ∈ ) are selected by mc, by contrast to the mc selection of just one quiescence time in the direct method. one may show that the smallest of these times is the time at which the next reaction occurs, and that this selection rule is an exact mc selection from the reaction probability density function, like eqs. and. another noteworthy difference between the direct method and the next reaction method is the employment of a "dependency graph" by the latter. this data structure reduces the number of calculations per reaction event, both for event selection and updating the reaction times. however, the dependency graph requires o objects to store the dependencies of reaction rates upon the populations of reactants and products shared with other reactions. for microarray hybridization, this translates to a storage requirement of o() reaction objects, which is prohibitively large for genome-sized microarray simulations; even current supercomputers with terabytes of memory cannot meet these requirements.  algorithm in simulations of microarray hybridization, the extreme computational burden of the direct method and the memory burden of the next reaction method may be alleviated by judicious storage and summation of the terms in eq. our algorithm employs a data structure called a hybridization table that stores partial sums of the terms in eq. for this reason, we call our approach the "method of partial sums".    initialize    calculate the auxiliary variables, ϕj, j =... np, α and ϕ    repeat       calculate the total reaction rate α and quiescence interval τ       if the next reaction is a hybridization then          select the hybridizing probe m and cdna ℓ       else          select the dehybridizing probe m and cdna ℓ       end if       update populations       update hybridization rates: αm and αj ≠ m, j =... np       update dehybridization rates: ϕm, ϕ       t ← t + τ    until p > for = algorithm 1: method of partial sums for the simulation of the coupled reaction network composed of the hybridizations of np oligonucleotide probes and nt cdna species. hybridization table the structure of the hybridization table is outlined in fig. in addition to storing the populations of targets {}, probes {}, hybrids {} and rate constants and, the table also contains the rates at which probes of type j will hybridize with any target, and the rates at which all hybrids composed of probes j will dissociate, two additional quantities are stored in the table: the total rate of hybridization and the total rate of dehybridization the total rate of reaction may then be calculated from eqs. and effectively subdivide the running sum in eq. into independently manageable groups of information. this in turn reduces the number of operations required for event selection reaction rate updates. reaction event selection the selection of the imminent hybridization or dehybridization event is accomplished by summing the terms of eq. such that the total number of operations in eq. is minimized. we begin by ordering the reactions according to column and row in the ht such that. eq. is then implemented by summing the reaction rates corresponding to each entry in the ht, by column and then by row: first, by hybridization events, and then by dehybridization events. one begins by identifying whether or not the imminent event will be a hybridization. since the imminent event is defined by the reaction whose rate causes the sum of αν to exceed the quantity r2α and all hybridization reactions precede the dehybridization reactions in the ht, it follows that if eq. results in the selection of a hybridization event, α >r2a and one needn't consider the dehybridization rates in the selection of the event to come. the next step is the selection of the probe that will hybridize in the imminent event. using the notation of the direct method and the order of summation, the index of the event to come must be less than or equal to np nt, where is the last of the hybridization rates in the table. eq. becomes in essence, the event is selected by summing the rates corresponding to each column of the "hybridization" section of the ht, column by column, followed by row, until the quantity r2a is exceeded. noting that the quantity is equal to αj, eq. may be simplified to which defines the probe that will hybridize in the imminent event. the equation defining the selection of the target may be derived similarly first, we express the sums on the right and left sides of eq. as and simplifying these expressions, we obtain thus, the target in the imminent event may be obtained by summing the quantities until the quantity is exceeded. if the event to come is a dehybridization, then α <r2a hence, the sum of all hybridization rates may be subtracted from each term in eq. leaving where {aν } are the rates of the dehybridization reactions and we have explicitly noted the fact that all np nt hybridization rates are contained within α. expressing this in terms of the hybrid indices, we obtain again, the selection of the the event may be simplified using the auxiliary variables in the ht. if the probe m is released in the imminent event, the value of m may be defined by sequential addition of the values of {Φj} until the quantity is exceeded note that this also defines the identity of the probe that will dissociate from the target. the cdna to be released in the dehybridization event may be calculated by subtraction of eq. from eq. 19: the selection rules for hybridization and dehybridization substantially reduce the number of operations required by eq. whereas the direct method may require 2np nt operations to select a reaction, our selection rules require at most operations. for a genome-sized microarray simulations, with np ~nt ~ our reaction selection approach will be several orders of magnitude faster than those of other algorithms. system state accounting upon selection of the imminent event, the populations of the species involved must be updated in accordance with the stoichiometry of the reaction. the structure of the ht facilitates this procedure: the row and column of the selected event designate the identities of both the reactants and products as well as the stoichiometric changes in population. however, other quantities must be updated, most notably the partial sums of reaction rates that facilitate event selection. the first quantity to be updated is the hybridization rate of the probe m where i = - for hybridization and + for dehybridization. since each partial sum αj is a function of the population of the ℓth cdna, these quantities must also be updated: subsequently, a must be recalculated using eq. collectively, 2np operations are required for the update of the hybridization section of the ht. after a hybridization or dehybridization event, only one of the dehybridization rates must be updated. thus, only one partial sum requires modification: moreover, inasmuch as Φm is the only affected partial sum, Φ may also be updated in one operation, thus, only two operations are required to update the dehybridization section of the ht. in summary, after selection of a reaction, the ht may be updated in o operations, a substantial improvement over the direct method. this is largely a result of the fact that reaction rates are not stored in the ht, precluding the need to update np rates of events featuring probe m and nt rates featuring target ℓ. the same argument applies to the dehybrization reaction rates {Φj}. ultimately, the method of partial sums is efficient so long as the partial sums can be updated easily. determination of equilibrium microarray analysis is predicated upon the assumption that the probes and solution-phase cdna have equilibrated prior to scanning the slide and measuring the fluorescence associated with each feature. we follow this experimental convention in silico. the hybridization process is at equilibrium when the rates of change of all chemical populations in eq. are zero, implying as straightforward as this criterion appears, it is difficult to employ in practice. eq. represents millions of comparisons for genome-scale microarrays, requiring o operations per time step. this many operations would also be required if one determined the steady states of all molecular species, which would additionally require storage of the current state as well as previous states. clearly, this is memory-prohibitive. furthermore, eq. is exact only on average. hence, it will never be exactly satisfied at any point in time within a single simulation. to circumvent these issues, we propose an alternate approach that employs the average total rates of hybridization and dehybridization. considering the definitions of these quantities summation of eq. over all ℓ and m yields the result that α = Φ at equilibrium. this criterion may be established using student's t test for two populations with unknown means and standard deviations. strictly speaking, this is necessary but insufficient for the specification of thermodynamic equilibrium. however, it is remarkably effective as a heuristic. we implement it as follows:    repeat       save α and Φ to disk every o reaction steps       maintain the last ten saved values of α and Φ in memory as vectors α and Φ       keep running averages of the numbers in these vectors, and.       if <then          if p < for h0: = then             reinitialize α and Φ          end if       end if    until p > for h0: = testing heuristic for the determination of equilibrium to evaluate the heuristic approach to the determination of the equilibrium state, we performed simulations of the hybridization of modified agilent probes with yeast cdna. a typical example of the time evolution of a and f is presented in fig. where the time at which equilibrium is attained is marked by the red circle. at and after this point, the populations of all hybrid species were at steady state, fluctuating in accordance with the predictions of equilibrium statistical mechanics. subsequent comparison of hybrid populations with eq. confirmed the findings of our "equilibrium heuristic" for all simulations. in standard hybridizations, imperfect mixing causes the transport of cdna or crna to be diffusion limited. this in turn presents an obstacle to hybridization, as the time required for a target to diffuse to it's probe is large. this, in turn affects the kinetics. the most common solution to this problem is to increase the concentration of target cdna or crna, which results in an abundance of these molecules at equilibrium. as we have discussed, our simulations feature perfect mixing. thus, we may use substantially less solution-phase cdna. as a result, free cdna is sparse at equilibrium, which introduces fluctuation into a. by contrast, f exhibits little fluctuation because hybrid populations are fairly large compared to the change in population accompanying a reaction. effect of reaction rate constants upon the steady state as discussed, the values of the kinetic rate constants should not affect the equilibrium state of our simulations provided that their ratios are constant. this assumption may be formulated as a testable hypothesis as follows: if the values of and affect the equilibrium state of simulation, then they will affect the fractional occupancy of each probe m ∈ defined by in this expression is the total population of probes m hybridized with any type of cdna at equilibrium. we test these hypotheses at the system level by performing simulations using rate constants that is, we perturb the rate constants generated via eqs. and by multiplying their results by a uniform random number rℓm. since a unique random number is generated for each probe m and cdna ℓ, one may quantify the effect of kinetic perturbations on the equilibrium state of a simulation via the hypothesis if kinetic rate constants significantly affect any of the fractional occupancies at the equilibrium state, then this hypothesis will fail. simulations of the hybridization of all modified agilent probes with full-length yeast cdna molecules were conducted for two types of perturbations. in the first study, r ~u, where u is a uniform random number on the interval. this allowed us to evaluate the effect of the widest variation of the rate constants. in the second study, r ~u. rate constants for all, hybridization events were independently perturbed. the results of our hypothesis tests are illustrated in fig. in both cases, our results clearly indicate that the equilibrium state is unaffected by the values of the rate constants, as expected. interestingly, the time required to reach equilibrium correlates with the heterogeneity of the dehybridization rate constants: deviations of the results for the two cases at hours of cpu time indicate that the hybridization of many probes with solution-phase cdna was incomplete in cases where their rate constants were perturbed by factors less than 10- analyses of the timeseries of the overall rate of reaction revealed that the progression to equilibrium is considerably slower if r ~u than if r ~u. this conforms with the experimental observations of dai and coworkers, which demonstrated differences between the kinetics of specific and nonspecific hybridizations. comparison of stochastic simulation algorithms results of all ssas applied to the process described by eq. should yield statistically indistinguishable results since they share a common stochastic process. in his seminal work, gillespie showed that the "first reaction method" and "direct method" were equivalent. subsequently, gibson and bruck demonstrated that their "next reaction method" is equivalent to gillespie's algorithms. since the method of partial sums shares the mathematical underpinnings of the direct method, its results should be indistinguishable from those generated by the next reaction method or gillespie's algorithms, as well as the law of mass action. we initially tested this hypothesis by conducting simulations of the hybridization of ten full length cdna species from yeast to ten probes for those species. five simulations were conducted via both the method of partial sums and the next reaction method. additionally, we solved the corresponding population balance equations for this illustrative hybridization process. our results clearly show that all methods yield equivalent results. average populations for all hybrid species could not be distinguished by simulation or calculation method via the t-test. interestingly, several common differential equation solvers could not integrate eqs. for this ten by ten system from t = until steady state, including the matlab packages ode and ode the two hundred equations necessary to model this small array required the use of the ode15s, which is designed for stiff sets of odes. by contrast, the stochastic simulation algorithms were unimpeded in their numerical progress. although the stochastic simulation algorithms give indistinguishable results for the population states both in time and at equilibrium, their computational performances are significantly different. the differences in the computational speeds of the ssas were evaluated by conducting simulations of the hybridization of full length yeast cdnas to microarrays featuring and probes targeting those cdnas. additionally, a simulation with probes and targets was conducted with the method of partial sums. the same initial populations of targets and the same set of probes were used in both the next reaction and mps simulations. estimation of the time at which equilibrium is attained was determined as discussed, and then used in simulations conducted with the next reaction method. our results are illustrated in fig. and the contrast is stark. the method of partial sums outperforms the next reaction method for simulations of microarrays of all sizes. it requires one hour to complete a simulation for an array large enough to be used for genomic characterization. by contrast, the next reaction method is not capable of performing such simulations due to its memory requirements. if np ~nt ~ the pointers required by the next reaction method will consume approximately three terabytes of bit computer memory by themselves; the ht requires no such pointers. next reaction simulations with as few as probes required hrs of cpu time to reach equilibrium, and for a computer with extraordinary memory, we forecast that simulations using the next reaction method would require a month to simulate the hybridization to the agilent yeast array. in addition, the scaling of the cpu time with respect to the number of probes differs among the two stochastic simulation algorithms. both algorithms feature power-law scaling, however, the method of partial sums scales as, whereas the next reaction method scales as. these differences arise due to the differences in the data structures underlying the two methods, which in turn affect the number of calculations required per time step. our algorithm permits two concurrent simulations of an million-reaction system to reach equilibrium within two hours using ghz dual g apple servers with gb memory. gb of memory are sufficient for simulations featuring probes, e.g. simulations of yeast arrays with one perfect match probe and one mismatch probe for each transcribed gene. simulations of an array designed for the human genome with one probe per gene - at approximately genes - possess hybridization reaction spaces approximately times larger, and an estimated run time of - hours based on cpu time scaling of t ~. such simulations can be achieved on most university shared-resource machines, which commonly feature hundreds of gigabytes of ram. illustrative example: characterization of cross hybridization the in silico characterization of cross hybridization is only as sensitive as the number of probe molecules per feature. agilent arrays have × probe molecules per feature, facilitating the hybridization and measurement of as many target cdnas and giving, in principle, as much resolution. however, the time required for stochastic simulations is proportional to the number of molecules therein. balancing these resolution and population considerations, we selected a hybridization volume of nl corresponding to initial probe populations of molecules/feature and concomitant scaling of the feature diameter to maintain the surface concentration. we also employed a total concentration of ng per μl hybridization volume. at this concentration, only a few if any probes become saturated. this concentration is a tenfold dilution of that recommended by agilent's protocol, however, it is is within the range of commercial oligonucleotide microarray protocols. moreover, the agilent array effectively has a probe population four times higher than the one we used, as it has four redundant features per probe. in fig. we illustrate a small subset of the hybrid populations predicted from our simulations. the average populations of hybrids were calculated from five replicate simulations, each of which required hours of cpu time as discussed. the complete sets of results are provided in additional files and many of the cross-hybridizing cdna species are exceptionally homologous. for instance, a_06_p - the probe designed to target cdna for ypl281c also hybridizes with cdna for yor393w; in this case, the two orfs are identical. the third hybrid, ymr323w, shares all but twelve bases of the other two. other cross hybridizing cdna species have less overlap in sequence, but share the probe sequences completely or in part. all cross-hybrids listed in fig. share at least 90% of the target sequence. a probe's proclivity to cross hybridize in the presence of thousands of potentially competitive cdnas may be expressed in terms of its selectivity, s, defined as the fraction of fluorescence intensity associated with a microarray feature that originates from the corresponding target. mathematically, where is the population of hybrids composed of probe m and its target ), and is calculated via eq. for example, the selectivity of the aforementioned agilent probe a_06_p is %. in fig. we present a summary of the selectivities for all probes in the agilent set. as these results show, the vast majority of the probes for this commercial microarray are selective, and do not exhibit cross hybridization when they bind yeast cdna at the concentrations specified by the expression state. for other microarrays and probe sets designed by various publicly available software packages, the cross hybridization may be more extreme. since the populations of all hybrids depend upon the populations of potentially cross-hybridizing cdna molecules, the value of s for each probe depends upon the expression state of a cell. interestingly, the distribution of s among the probes is independent of the expression state. simulations with different initial conditions corresponding to four additional mc-generated expression states yielded selectivity distributions that were indistinguishable from the aforementioned distribution. this result suggests a method of characterizing the overall reliability of a proposed probe set. since the distribution of selectivities is invariant with respect to the initial target populations, a statistic that characterizes that distribution should be robust. we propose that the average selectivity can serve as such a metric. the use of such metrics should be used with care, however, inasmuch as they lend themselves to "ecological fallacy". implementation our software is implemented in c++ and executed on apple g processors running macos x. unafold is available for download at the dinamelt server http://dinamelt.bioinfo.rpi.edu/ and compiles on a variety of operating systems. additionally, we have implemented the next reaction method as described by gibson and bruck in c++ for the purposes of comparing its performance with that of our algorithm. discussion in recent years, advances in microarray manufacturing have opened the doors to custom microarray design. individual researchers can now upload their own microarray probe sequences to one of many sites and have a custom array manufactured. in response, many probe design algorithms have been arisen to fill the needs of researchers intent on performing global investigations of gene expression. however, none of the probe sets generated by these algorithms have been evaluated in terms of their selectivity or proclivity to give a linear relationship between feature fluorescence and target concentration. studies of this type carried out by the microarray quality control consortium were costly, involving dozens of participants. such laborious quality control is not feasible for each and every probe set designed by a novel probe tool. however, robust population-based simulations of the hybridization process may be employed to evaluate candidate probe sets, given robust estimates of the thermodynamic free energies of hybridization. in the tests of our simulation algorithm, we have utilized equilibrium constants that were calculated via the nn model of allawi and santalucia. this model predicts free-energies to within three significant digits of experimentally measured values for solution-phase hybridization and has been validated for oligonucleotide arrays in which the oligonucleotides are connected to the surface by linkers, making them more "solution like". however, many studies have demonstrated that nn models do not accurately predict the thermodynamic properties of hybridization between solution-phase and surface-bound oligonucleotides. electrostatic interference in conventional microarrays - which do not feature 3d linkers of the type employed by weckx and coworkers, is a major reason for this discrepancy. solution-phase na+ may shield the phosphate groups of both hybridized and single-stranded probes and targets. cations will also shield the negatively-charged glass surface to which the probes are bound by organizing the formation of a layer of counter ions. the surface density the oligonucleotide probes also plays a role, partly via steric interference and partly via changes in the charge distribution at the glass-water interface. ultimately, surface electrostatic effects cause a location-dependent effect of mismatches. therefore, extra care must be taken when applying our algorithm to characterize cross hybridization in real microarray assays. if equilibrium constants are calculated using nn models, corrections for the effect of the debye layer should be introduced at a bare minimum. theoretical predictions of the effect of the resulting "debye layer" upon the melting temperature tm have been confirmed via experimental measurements. this correction may be applied to nn models for dna-crna duplexes as well, if simulations of crna-dna assays are to be conducted. we have not explicitly considered additional effects of single-strand secondary structure formation among full-length cdna or probes nor, for that matter, have we considered the possibility of hybridization between cdna molecules in solution the hybridization between probes, or folding thereof, is customarily not a substantial problem in oligonucleotide microarrays. probe design software packages such as oligoarray routinely screen probes for secondary structure, and inter-probe hybrids are precluded due to spatial separation. the remaining interactions among cdnas serve only to sequester single-stranded cdna from the probes. inasmuch as this can only further degrade the sensitivity of probe-target interactions, we have restricted our algorithm to consider only probe-cdna hybridization dynamics. hence, the method we have presented for probe and array characterization gives the performance of a microarray under a "best case" scenario, even as it fully accounts for the interactions between oligonucleotide probes and full length cdna targets at the system-scale. for the end user we note that simulation-based characterization of putative microarray probe sets requires no more information than that contained in mage-tab or mage-ml formatted array information, which are required for publication of microarray data to ensure miame compliance. the array design format component of mage-tab contains all of the probe sequence information and its target, whereas the investigation description format contains the experimental protocols, including the hybridization temperature and salt conditions. additionally, the adf contains information regarding the relationship between reporter sequences and features: there may be multiple features with the same reporter sequence, or alternately, signals from several reporters may be combined to produce a signal for a single gene. as the total cdna concentration is ostensibly included in the idf, and the individual cdna concentrations are randomly generated, the system-scale selectivity of any array can be computed via simulation from a proposed experimental protocol. in this work, we have not explicitly considered the prediction of the time evolution of the hybridization process, focusing instead on the equilibrium state. simulations employing our algorithm are most accurate when they utilize experimentally-determined rate constants. however, care must be taken to ensure simulations are conducted for the same surface densities and ionic strengths employed in the experiments employed to estimate the rate constants, for the reasons previously discussed. given accurate estimates of the rate constants and, the time required for microarray hybridizations may be estimated via the method illustrated in fig. as the time variable would represent the actual hybridization time.  in this work, we have developed an algorithm for the stochastic simulation of exceptionally large and complex probe-cdna hybridization reaction networks that underlie microarray assays. using the method of partial sums in conjunction with the data structure we denote the "hybridization table", our algorithm requires o operations per reaction event. this is substantially fewer than gillespie's direct method operations per event) and the next reaction method of gibson and bruck ) operations per event). moreover, our algorithm requires less the data storage than others, obviating the need for pointers that track of the dependencies of reactions. for instance, the next reaction method requires pointers for each of its 2np nt reactions, which can consume a vast amount of memory for genome-scale simulations. as a result of these innovations, our algorithm permits system-level simulation of the complete reaction network composed of all potential probe-target hybridizations without the need for high-performance computing. furthermore, such simulations are now possible within a reasonable amount of time. thus, given robust thermodynamic predictions of the free energies of dna hybridization, one may obtain a conservative estimate of the reliability of a candidate probe set in silico. 
by integrating stochastic simulation with thermodynamic prediction tools for dna hybridization, one may robustly and rapidly characterize of the selectivity of a proposed microarray design at the probe and "system" levels. since hybridization is a chemical process, it may be modeled at a population-level using a combination of material balance equations and thermodynamics. cross-hybridization between microarray probes and non-target ssdna has been implicated as a primary factor in sensitivity and selectivity loss. however, the hybridization reaction network may be exceptionally large for commercial arrays, which often possess at least one reporter per transcript.
 although oligonucleotide microarray technology is ubiquitous in genomic research, reproducibility and standardization of expression measurements still concern many researchers. cross-hybridization between microarray probes and non-target ssdna has been implicated as a primary factor in sensitivity and selectivity loss. since hybridization is a chemical process, it may be modeled at a population-level using a combination of material balance equations and thermodynamics. however, the hybridization reaction network may be exceptionally large for commercial arrays, which often possess at least one reporter per transcript. quantification of the kinetics and equilibrium of exceptionally large chemical systems of this type is numerically infeasible with customary approaches.  in this paper, we present a robust and computationally efficient algorithm for the simulation of hybridization processes underlying microarray assays. our method may be utilized to identify the extent to which nucleic acid targets will cross-hybridize with probes, and by extension, characterize probe robustnessusing the information specified by mage-tab. using this algorithm, we characterize cross-hybridization in a modified commercial microarray assay.  by integrating stochastic simulation with thermodynamic prediction tools for dna hybridization, one may robustly and rapidly characterize of the selectivity of a proposed microarray design at the probe and "system" levels. our code is available at http://www.laurenzi.net. 