This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.106.090662v1 91/11/4133    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Google Scholar Google Scholar Articles by Horne, M. T. Articles by Benight, A. S. Search for Related Content PubMed PubMed Citation Articles by Horne, M. T. Articles by Benight, A. S.

Statistical Thermodynamics and Kinetics of DNA Multiplex Hybridization Reactions

M. T. Horne ^* , D. J. Fish  and A. S. Benight ^*  

^* Departments of Chemistry, Physics and Mathematics, Portland State University, Portland, Oregon; and Portland Bioscience, Portland, Oregon

Correspondence: Address reprint requests to A. S. Benight, Tel.: 503-725-9513; E-mail: abenight{at}pdx.edu.

	ABSTRACT

TOP ABSTRACT INTRODUCTION THEORETICAL APPROACH EXAMPLE CALCULATIONS COMPARISONS WITH PREVIOUS WORK CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
A general analytical description of the equilibrium and reaction kinetics of DNA multiplex hybridization has been developed. In this approach, multiplex hybridization is considered to be a competitive multichannel reaction process: a system wherein many species can react both specifically and nonspecifically with one another. General equations are presented that can consider equilibrium and kinetic models of multiplex hybridization systems comprised, in principle, of any number of targets and probes. Numerical solutions to these systems for both equilibrium and kinetic behaviors are provided. Practical examples demonstrate clear differences between results obtained from more common simplex methods, in which individual hybridization reactions are considered to occur in isolation; and multiplex hybridization, where desired and competitive cross-hybrid reactions between all possible pairs of strands are considered. In addition, sensitivities of the hybridization process of the perfect match duplex, to temperature, target concentration, and existence of sequence homology with other strands, are examined. This general approach also considers explicit sequence-dependent interactions between targets and probes involved in the reactions. Sequence-dependent stabilities of all perfect match and mismatch duplex complexes are explicitly considered and effects of relative stability of cross-hybrid complexes are also explored. Results reveal several interdependent factors that strongly influence DNA multiplex hybridization behavior. These include: relative concentrations of all probes and targets; relative thermodynamic stability of all perfect match and mismatch complexes; sensitivity to temperature, particularly for mismatches; and amount of sequence homology shared by the probe and target strands in the multiplex mix.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION THEORETICAL APPROACH EXAMPLE CALCULATIONS COMPARISONS WITH PREVIOUS WORK CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Historically, theoretical and experimental studies of the equilibrium melting (or hybridization), and kinetics of short duplex DNAs have primarily focused on "simplex" reactions wherein two short single strands, whose sequences are perfectly complementary, anneal to form a perfectly matched Watson-Crick (w/c) basepaired duplex. Simplex melting experiments of short duplex DNAs, with well-defined sequences, and theoretical descriptions of the melting process, have been an active subject of study for over 25 years (1–11). Systematic studies of the melting stability of short duplex DNAs with different sequences have provided evaluations of nearest-neighbor thermodynamic stability parameters that enable prediction of the thermodynamic stability of a short duplex from its basepair sequence (1,8,9). These parameter sets are commonly used in the design of probes and primers having desired hybridization and stability properties for diagnostic assays. Similarly, kinetic studies of the annealing (or hybridization) of two perfectly matched single strands to form a duplex, have been performed and have provided analytical descriptions of the simplex hybridization process. These have provided evaluations of the kinetic rate constants, and other parameters, describing hybridization kinetics of short DNAs (12–23).

Although parameter values determined from kinetic and equilibrium analysis of simplex reactions have been quite useful in facilitating improved predictions of the sequence-dependent melting and kinetic behaviors of homogeneous solutions of short duplex DNAs and their relative sequence-dependent stability, more accurate parameters alone are insufficient to provide realistic descriptions of multiplex hybridization reactions, i.e., when more than one duplex is present in the same solution. There are several additional, essential components of multiplex hybridization reactions that must be considered, which make them remarkably different from more common simplex reactions. These are important considerations as multiplex hybridization forms the basis of many modern nucleic acid diagnostic reactions.

Nucleic acid diagnostic assays based on multiplex hybridization have the potential to revolutionize genomic research and genetic medicine (24–31). Multiplex assays can be performed on microarrays or in solution via various iterations of the polymerase chain reaction (32–42). These assays offer never-before-imagined capabilities for systematic high throughput screening, discrimination, and analysis of large numbers of DNA (and RNA) sequences. Despite the vast and important applications of multiplex hybridization in nucleic acid diagnostics, there have been relatively few studies aimed at gaining a better understanding of the actual hybridization reactions that can occur in mixtures containing more than two strands.

As stated above, the majority of studies of melting (or hybridization) reactions of short duplex DNAs have been performed primarily on equimolar mixtures of two single strands whose sequences are perfectly complementary. Exceptions are studies of solutions containing a single species of a linear or circular self-complementary single strand for the purpose of evaluating DNA sequence-dependent stability parameters, and examining intramolecular hairpin stability (43–45), and studies performed for the purpose of investigating the thermodynamic parameters of single basepair mismatches (1,9,46), bulged loops (47–49), or multiple strand mixtures, such as triplexes and quadruplex (50). Melting studies of mixtures of more than two strands in solution, which can form both perfect match and mismatch duplex complexes, have suggested that sequence-dependent interactions in tandem mismatches might contribute to the stability of mismatch hybrid duplex complexes (51). For this and other reasons described herein, when multiple strands meant to form perfectly matched duplexes are present together in a reaction mix, it cannot be presumed that each pair of strands will form the appropriate duplex, exclusively, in the same way they would in isolation, in the absence of the influence from other single strands (and duplexes). Several additional essential features of multiplex hybridization must be considered. These considerations are associated with the significant probability of formation of mismatch hybrids (cross-hybridization) brought about by sequence homology with other strands and sequence-dependent stability of mismatch basepairs.

Existing nucleic acid platforms and analysis procedures continue to be improved and new technology platforms promising higher sensitivity and more highly reproducible results are being developed to the point of providing reliable quantitative measurements (16,17,52–54). Along with these developmental improvements, it is essential that an analytical foundation be established that will enable formulation and support of robust and realistic models of multiplex hybridization. In turn, these will facilitate enhanced design, analysis, and optimization of nucleic acid multiplex diagnostic assays.

Several authors have recognized the necessity of considering multiple hybridization reactions simultaneously and demonstrated, in several simple cases, generally dramatic effects associated with the potential of multiple two-strand interactions (55,56). As might be expected, these studies have pointed out that interactions between each probe:target pair cannot be assumed to occur as separate events. These studies suggest the essential importance of considering competitive reactions in a multiplex hybridization reaction mixture in a multichannel-systems manner.

For the case of hybridization on microarrays, a formal approach to modeling the hybridization of targets with probes affixed to a surface involves two distinct phases (16,22,56,57), termed "transport" and "reaction." The transport phase involves diffusion of target strands across the probe surface. The reaction phase involves the reaction (binding and dissociation) kinetics between targets and probes at the surface. For the latter, once targets have diffused to the interaction zone, reaction kinetics dominates the process, and effects of sequence-dependent interactions become quite significant. If the target concentration is sufficiently in excess of the probe concentration, to the point where encounters between all targets and probes are equally likely, then the entire process is dominated by the reaction phase. The development presented here considers only the reaction phase of the multiplex hybridization process.

Model systems studied so far have discovered complications associated with competition and cross-hybridization between two and three strands in the same reaction, and anticipated added complexities associated with competition between multiple-strand interactions in multiplex mixes (55,56). However, general descriptions of the equilibrium and kinetic behaviors of DNA multiplex hybridization reactions, containing any number of probes and targets, have not been presented. To model multiplex hybridization reactions in a completely general way, we have developed a systems-analysis approach involving collective consideration of multiple reactions and the simultaneous solution for individual products of specific reactions. In the development presented here, multiplex hybridization is considered to be a competitive, multichannel, reaction process: a system wherein many species can react both specifically and nonspecifically with one another. General equations are presented supporting both equilibrium and kinetic models of multiplex hybridization systems comprised, in principle, of any number of targets and probes. Practical examples demonstrate clear differences between simplex and multiplex hybridization and sensitivity of the hybridization process of perfect-match duplexes to temperature, target concentration, and existence of sequence homology with other strands.

This article is presented in the following sections. In Theoretical Approach, our general theoretical approaches are presented. Both equilibrium and kinetic models are described, and their analytical solutions are provided. Example Calculations contains the results of a number of comparisons of examples of simplex and multiplex calculations, and their sensitivity to temperature, sequence homology, and stability of mismatch hybrid duplexes. In Comparisons with Previous Work, comparisons to relevant published works are provided. The final section (Conclusions) summarizes our major findings and provides conclusions.

	THEORETICAL APPROACH

TOP ABSTRACT INTRODUCTION THEORETICAL APPROACH EXAMPLE CALCULATIONS COMPARISONS WITH PREVIOUS WORK CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Equilibrium statistical thermodynamics
Consider a population of DNA single strands at constant temperature and pressure. For any molecular configuration adopted by these strands, i.e., single-strand intramolecular hairpin or two-stranded duplex structure, the number of microstates is given by the Gibb's factor for that configuration. For a specific configuration, i, at approximately constant solution volume, the Gibb's factor or statistical weight is given by , where is the standard-state free energy of the i^th configuration. The standard-state free energy is given in terms of the differences of the standard-state chemical potentials of the configuration and the unstructured (melted) single strands from which it formed. For example, for a duplex p_it_j formed from a probe p_i and a target t_j, the standard-state free energy is given by

where , , and are the standard-state chemical potentials for the probe:target complex, and the individual probe and target strands, respectively.

In this development, the reference state of the single strands is defined as the unstructured, unstacked (melted) single strand. The reference state for the duplex can be defined in a number of ways. For example, as the nucleation complex of two single strands in approximately the same orientation and volume as the fully duplex state, in the absence of hydrogen-bonding or stacking interactions between the strands. The system partition function Z is the sum of the statistical weights for all possible configurations, i.e., . The probability of an occurrence of configuration i is then given by P=sw_i/Z.

Consider a multiplex system wherein a set of probes, p_i, is designed to specifically hybridize with particular targets, t_j. Implicit in such a specificity requirement is the potential for error. This is necessarily so because, conceivably, every probe and target strand may form either a perfectly matched w/c duplex or a "degenerate" structure, i.e., a nonperfect match hybrid duplex, or an intramolecular hairpin. The standard-state free energies that define the statistical weights, and therefore the population of error states, depend explicitly on the differences in their respective standard-state chemical potentials.

Chemical system
To model the "system" behavior of multiplex hybridization reactions, a system of equilibrium reactions is assumed (see Table 1) where, for example, p_it_j is the ensemble of duplex states that form from single strands p_i and t_j with equilibrium constant . Similar definitions apply for the other reactions. This formulation also considers, through the index , the possibility of microstates within each configuration, such that x represents a specific microstate of the configuration, e.g., p_it_j could represent an overlap duplex state of a configuration involving single strands p_i and t_j, where the two strands might not be perfectly aligned at their ends. If these microstates are not considered (as in formulation of the kinetic model described later), = 1.

(1a)

(1b)

where and are the total concentrations of the probe and target strands, respectively.

Now, the concentrations of the different reaction products are given in terms of the equilibrium constants for their formation:

(2a)

(2b)

(2c)

(2d)

(2e)

Define the following:

With these expressions, Eqs. 1a and 2b become

(3a)

(3b)

For simplicity, if target:target, probe:probe, and intramolecular states in single-strand targets and probes are ignored, i.e., , Eqs. 3a and 3b can be combined to form the following system of coupled, nonlinear equations:

(4a)

(4b)

If there are N_P probes and N_T targets, and l goes from 1 to N_T, then Eq. 4a or 4b represents a system of N_P+N_T nonlinear coupled equations. Given input values of and , and the equilibrium constants, and , a solution of this system produces values of , which can in turn be used to evaluate the corresponding values of . These are then combined with the particular equilibrium constants, Eq. 2, to determine the duplex concentrations, .

Multiplex hybridization kinetics
Consider a multicomponent system comprised of two species types, i.e., probes p_i, and targets t_j. Each p_i can react with every t_j forming a p_it_j complex or can react with each of its counterparts of the same species, forming p_ip_j and t_it_j complexes (for simplicity, we ignore hairpin formation as its inclusion does not affect the results that follow). Each of these individual reactions occurs with forward rate constants , , and and reverse rate constants, , , and respectively (see Table 1). The forward rates constants are assumed to be the same for all species, i.e., .

It should be noted that this development concerns only the "reaction" kinetics of the hybridization process. That is, all targets are available for interactions with all probes. It is implicitly assumed that there are no diffusional barriers to the reaction process. Under the umbrella of this assumption, in combination with the above assignments, the fundamental rate equations are as follows:

(5a)

(5b)

and

(6a)

(6b)

(6c)

Here, is the concentration of the i^th probe (i N), is the concentration of the j^th target (j M), and , , and are the concentrations of the probe:probe (i, j N), probe:target (i N, j M), and target:target (i, j M) duplexes. The factor (_ab+1) is either 1 or 2, depending on whether a b or a = b, respectively. After ignoring repetition due to symmetries (C_pp and C_tt are symmetric), a system of d = (1/2)(N+M)(N+M+5) equations is formed.

Observe that by combining Eqs. 5 and 6, the rate equations for the single-strand probe and target concentrations can be written as

If, for i N and j M, the functions and are defined as

(7a)

(7b)

then

Thus, the functions and are constants with respect to the time of reaction for the system, . In fact, these functions form a set of N+M independent constants. This allows the reduction of the system to the smaller system obtained by performing the following substitution.

Let a solution to be given with initial conditions , , , , and . Then for all time, , and . According to the definitions in Eqs. 7a and 7b, and suppressing the variable ,

(8a)

(8b)

For convenience, denote the right-hand sides of Eqs. 8a and 8b by and , respectively. Then Eqs. 6a–6c can be replaced by

(9a)

(9b)

(9c)

Together with Eqs. 8a and 8b, these define an equivalent system of N+M fewer equations.

Since the functions and are constant along solutions to , Eqs. 8a and 8b form a system of N+M independent nonlinear equations, which is satisfied by the solution to with initial conditions of

By continuity, any equilibrium point of also satisfies this system of equations. Thus, the search for equilibria of the system should begin with a search for solutions to this system. The following relations are seen to hold at equilibrium points of ,

(10a)

(10b)

(10c)

or, if , where x and y are either p_i or t_j, then

(11a)

(11b)

(11c)

Thus, any equilibrium solution to must be contained in the set of common solutions to both systems of Eqs. 8 and 11. These equations comprise the formal foundation of the so-called equilibrium model for the hybridization of probes and targets. In summary, the basic expressions are

(12a)

(12b)

(12c)

(12d)

(12e)

Examination of these expressions reveals they are identical to those in Eqs. 3a and Eq. 3b derived from equilibrium considerations alone. Solutions of these determine equilibrium concentrations for each species.

Recall that p_it_j represents the duplex formed from single strands p_i and t_j, with similar definitions for the other reactions in Table 1 (neglecting hairpins in single strands). Following are alternate forms of the expressions in Eq. 5:

(13a)

(13b)

Likewise, the following are alternate forms of expressions in Eq. 6:

(14a)

(14b)

(14c)

Initial conditions for the system are C_xy = 0 and C_x = 0 (i.e., all strands are initially in single-strand conformations) and C_z = C⁰ (i.e., initial single-strand concentrations are given), where x, y, z represent p_i or t_j. If misaligned, overlapping states are ignored, then = 1, in the expressions in Eq. 13. Assuming that for all states x, then , in analogy with the expressions in Eqs. 5 and 6.

It should be noted that, with little difficulty, the preceding development can be extended to include hairpin formation in single strands. To do this, the expressions in Eq. 5 are amended to include the effect of hairpin formation, and rate equations for hairpin configurations are added to the set of expressions in Eq. 6. As a result, Eqs. 12a and 12b include additional terms, and the expressions in Eqs. 12c–12e are expanded. This extension has been omitted in the above discussion for the sake of simplicity.

Further, the system in Eq. 5 allows for potential interactions between distinct species of probes, i.e., the formation of probe:probe duplexes is not precluded. Thus, the models developed above are not specifically restricted to microarray hybridization, but can also be applied to multichannel simulations of DNA hybridization in solution. Of course, elimination of probe:probe interactions from the simulation can be achieved by the appropriate assignment of rate constants (i.e., for all i and j).

Numerical models
Equilibrium model
Estimates of the solutions to the system of Eq. 12 can be obtained using any standard numerical optimization program. For example, the nonlinear least-squares optimizer provided by MatLab's Optimization Toolbox (58), as well as a similar algorithm provided by the software OCTAVE (59) have been employed with comparable degrees of success. These algorithms were applied to the above equations, written as a vector-valued function F on R^N+M as follows (hairpin reactions are not considered here, but can be included without difficulty):

To numerically find the roots of this system of equations, an initial guess (or seed) must first be made that approximates the true solutions. If seed values are not chosen appropriately, the algorithm can terminate without returning acceptable values for concentrations of probes and targets. To overcome this difficulty, an iterative process was developed that reseeds the algorithm until an acceptable level of error is reached in the approximate solutions. Using this process, the equilibrium concentrations of hybridization systems of N probes and M targets with N and M as high as 1000 have been successfully computed. On a 3.6 GHz Pentium-IV machine with 1 Gigabyte of RAM, the algorithm required 45 min to complete a 1000 x 1000 simulation.

For more accurate values, equilibrium concentrations calculated from the Kinetic Model (see next section) can be used as seed values for the equilibrium algorithm. Unfortunately, the kinetic algorithm involves a much larger number of equations, and due to hardware limitations this method was not feasible for systems higher than 80 x 80. However, for smaller systems, the combination of the two algorithms produces highly accurate results in a reasonable amount of computing time.

Kinetic model
The set of differential expressions in Eqs. 5 and 6 above (again, hairpin reactions are not considered) comprise a model of the kinetics for the reaction phase of microarray hybridization. An ODE solver written in FORTRAN was implemented to find solutions to this set of equations for kinetic simulations up to N = M = 80. The number of equations (when N = M) grows like N², and the size of the Jacobian matrix grows like N⁴, making simulations of larger systems highly memory-intensive. In Fig. 1, the computation times for the two types of simulations (top, kinetic model; bottom, equilibrium model) are shown. Kinetic simulations up to N = M = 35 and equilibrium simulations up to N = M = 100 were run on a 3.6 GHz Pentium-IV CPU and 1 Gigabyte of RAM, running Windows XP. Similar results were obtained when the simulations were carried out on a 1.8 GHz/512 RAM machine.

View larger version (8K): [in this window] [in a new window]   FIGURE 1  Computation time versus number of probes and targets for the kinetic model (top) and equilibrium model (bottom).

View larger version (20K): [in this window] [in a new window]   FIGURE 2  Flowcharts of the model calculations. (a) Flowchart of the kinetic model calculation. (b) Flowchart of the equilibrium model calculation.

	EXAMPLE CALCULATIONS

TOP ABSTRACT INTRODUCTION THEORETICAL APPROACH EXAMPLE CALCULATIONS COMPARISONS WITH PREVIOUS WORK CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Model system: definition and rationale
Our general models can, in principle, be applied to multiplex hybridization of any number of probes and targets. In the examples that follow, a simple model system was defined to demonstrate features of the calculation and test sensitivity to a variety of experimental and model parameters. Although described methods are generally applicable to systems of arbitrary dimension, for simplicity, details of a system comprised of three probes and three targets (3 x 3 system) are considered. Although relatively simple, this system demonstrates the features of added complexity of multiplex hybridization and provides adequate resolution to reveal typical behaviors. Two sets of sequences, set I and set II, were the subject of model calculations. In this regard, this system is an extension of the 2 x 2 competitive hybridization analysis that has been reported by others (55). An added feature of our system is the inclusion of the specific sequence-dependent stability of perfect-match and cross-hybrid duplexes, and corresponding consideration of the effects of sequence homology on the hybridization process. As seen for set I in Fig. 3 a, there is no sequence homology between the three probes or targets comprising the set. In contrast, for set II, two probes P1 and P2 share 50–54% sequence homology, while the sequence P3 has only 25% homology with P1 and 30% homology with P2. All model probe:target duplexes contain 24 basepairs.

View larger version (34K): [in this window] [in a new window]   FIGURE 3  Sequences used in model calculations. Mismatched sequences are underlined. (a) Set I is comprised of three probe and target sequences. The probe and target sequences are independent and do not share any homology with one another. (b) Set II sequences are comprised of three probe and target sequences. Probes P1 and P2 are 50% and 54% homologous with target strands T2 and T1, respectively. Probes P1 and P2 share 25 and 30% homology, respectively, with target T3.

This duplex contains eight nearest-neighbor interactions, including single-base 5' dangling ends. The nearest-neighbor dependent parameters, , , , , , etc. for the appropriate sequences and interactions as tabulated in Tables 2–4 were utilized. Experimentally derived parameters for w/c perfect matched doublets, single base dangling ends, and single basepair mismatches were available from the published literature (46). Nearest-neighbor tandem mismatches were assigned values as described in Eq. 15 (see below). The actual parameter values employed are summarized in Tables 2–4. In addition, two initiation factors, and were assigned depending on the particular identities of the end basepairs. Values of the initiation thermodynamic parameters that were employed are

Furthermore, recall the formulas for total free energy, G = HTS, and T_m = H/S, G = H(1 – T/T_m).

View this table: [in this window] [in a new window]   TABLE 2  Nearest-neighbor thermodynamic parameters for w/c doublets (61–64)

View this table: [in this window] [in a new window]   TABLE 4  Sequence-dependent thermodynamic parameters for single basepair mismatches (72)

View this table: [in this window] [in a new window]   TABLE 3  Sequence-dependent parameters for dangling ends (60)

(15)

where G_PM, H_PM, and S_PM are the free energy, enthalpy, and entropy, respectively, for melting a hydrogen-bonded w/c basepair doublet. The factor was introduced as a means of scaling values of thermodynamic parameters of mismatch basepairs in tandem mismatches as a relative fraction of the stability of w/c perfect matches. In examples that were performed, tandem mismatches were treated in two ways. They were either assumed to be minimal, = 0, or assigned a value of = 0.5. Although consideration of tandem mismatches in this manner is clearly an oversimplified generalization, it provided a convenient means of universally weighting non-w/c tandem mismatch pair interactions differently than w/c basepairs, and discerning potential effects of tandem mismatch stability on multiplex hybridization.

Theoretical results
Single-channel versus multichannel hybridization
Our approach adds several new features to modeling of the reaction-phase of the multiplex hybridization process. For a multiplex mixture comprised of any numbers of probes and targets, duplexes formed between perfect match probes and targets, as well as all cross hybrids formed from mismatched probes and targets, i.e., cross hybrids, are all considered collectively in a multichannel, system approach. There are multiple reaction channels available to every probe and target. In addition, sequence-specific effects are considered for both perfect match duplexes and mismatched duplexes containing some amount of basepairing. Consideration of the potential fractional stability of tandem mismatches and associated stronger weighting of mismatch complexes is also included.

The multichannel reaction model is inherently more complex than the more common single channel approach, where hybridization reactions between perfect match probe:target pairs are considered in isolation. Practically, when the single channel model is employed, potential cross-reactions between mismatched probe:target pairs are either screened and filtered according to their degree of sequence homology, or ignored entirely (66). For the model system comprised of three probes and three targets that will be examined in detail, the single-channel calculation (as so defined) considers independently the hybridization behaviors of only the three perfect match probe:target complexes. With added complexity the multichannel model considers the interdependent behaviors of nine different duplexes, i.e., three perfectly matched, and six cross-hybrids that contain mismatches. As will be seen, both the kinetic and equilibrium behaviors of the three perfect match probe:target duplexes are influenced significantly by added considerations inherent in the multichannel model.

Plots in Fig. 4, a–d, show the results of the single channel (top panels) and multichannel (bottom panels) calculations for a simple 3 x 3 system. Kinetic curves (duplex concentration versus time) that are displayed were calculated for set I (Fig. 4, a and c) and set II (Fig. 4, b and d) at 37°C for the single-channel and multichannel models. Plots in Fig. 4 also show effects of assigning differential stability to tandem mismatches and essentially giving more weight to cross-hybrids. Plots in Fig. 4, a and b, were obtained by assuming the free energy of tandem mismatches as 0, i.e., = 0, while those in Fig. 4, c and d, were obtained assuming a universal constant factor, = 0.5 (to scale tandem-mismatch free energies relative to normal w/c basepairs).

View larger version (22K): [in this window] [in a new window]   FIGURE 4  Kinetic plots (time versus concentration) for the 3 x 3 systems comprised of sequence sets I and II (Fig. 4) at constant temperature, 37°C. Curve assignments are given in the individual figures. (a) Set I. Single-channel (top) and multichannel calculation (bottom). (b) Set II. Single-channel (top) and multichannel calculation (bottom). Curves in panels a and b were generated assuming a minimal contribution of tandem mismatches to cross-hybrid duplex stability ( = 0). (c) Set I. Single-channel (top) and multichannel (bottom) calculations assuming tandem mismatch stability is comparable to w/c basepair doublet stability ( = 0.5). (d) Set II. Single-channel (top) and multichannel (bottom) calculations. All curves were generated assuming negligible stability for tandem mismatches ( = 0.5).

Calculated kinetic curves for set II are shown in the upper and lower panels in Fig. 4 b. Recall that, in set II, probes 1 and 2 share at least 50–54% sequence homology with targets 2 and 1, respectively, and probe 3 is 25% homologous with targets 1 and 2. Kinetic curves in the upper panel of Fig. 4 b are from the single-channel calculation and exhibit typical exponential behavior.

In the multichannel calculation for set II, behavior of the curves displayed in the lower panel of Fig. 4 b contrasts to what was observed for set I. In the multichannel calculation, when = 0, kinetic curves for the perfect match duplexes of set II differ considerably from those obtained in the single-channel calculation. In this case, concentrations of the perfect match complexes, P1:T1 and P2:T2, do not increase as rapidly as in the single-channel calculation (upper panel, Fig. 4 b) and take longer to level off to equilibrium values. Apparently, due to the reduction of available single strands from competition with cross-hybrids, particularly at early times, the perfect match complexes reach equilibrium more slowly, and achieve lower equilibrium values at comparable times. This is evident from the kinetic curves in Fig. 4 b for the dominant mismatch complexes, P1:T2 and P2:T1 (black and orange lines in the lower panel of Fig. 4 b), which increase at early times attaining concentrations comparable with the perfect matches, then decrease exponentially to their equilibrium concentrations. The extent to which the competition for resources affects formation of the perfect match hybrids depends directly on the amount of sequence homology among the strands. Note that, for set II, the perfect-match complexes, P1:T1 and P2:T2, approach equilibrium much more slowly than the P3:T3 complex. Consistent with this scenario, the P3:T3 complex suffers less dramatic effects of cross-hybridization. The P3 and T3 strands share only 25% homology with the other strands. Consequently, more of them are available to form the perfect-match complex resulting in much faster approach to equilibrium than the other perfect match complexes. The curves for the P1:T2 and P2:T1 cross-hybrids (black and orange lines in plot in lower panel, Fig. 4 b) bracket the values for the perfect matches P1:T1 and P2:T2 at very early times, but then decrease exponentially with time to reach equilibrium. Some of these curves demonstrate a stark departure from the kinetic behavior seen thus far. Apparently, competition from mismatches reduces the amount of corresponding perfect matches containing the same strands. In Fig. 4 b the curves for the P1:T2 and P2:T1 mismatches are the only ones visible on the scale of perfect matches. In summary, when tandem mismatches have marginal stability ( = 0), cross-hybrid formation is largely dominated by the amount of sequence homology shared by the probe and target strands.

When tandem mismatches have a significant fraction of the stability of normal w/c basepairs, = 0.5 (Fig. 4, c and d), results are much more dramatic. For sets I and II, in the single-channel calculation, assigning more stability to tandem mismatches has no effect, since mismatch cross-hybrids are not considered. As a result, for the single-channel calculations, when = 0.5 for sets I and II (upper plots in Fig. 4, c and d), plots are identical to those in Fig. 4, a and b, respectively. However, in the multichannel calculation, when = 0.5, cross-hybrids are given more stability. Because of increased competition from mismatches, all the perfect matches reach lower equilibrium values much more slowly than in the single-channel model calculation. When = 0.5, some mismatches are actually found to reach higher concentrations than the perfect matches (which seems unlikely). Note that all perfect-match complexes reach lower equilibrium values than observed in the single-channel calculation. The scale of the plot is different in order to help visualize the qualitative behavior. Plots in the bottom panel of Fig. 4 c for set I underscore the significant effects of tandem mismatch stability on the multiplex hybridization process, even when no sequence homology exists between the probes and targets. In effect, stable cross-hybrids compete for the perfect-match probe and target strands. As a consequence of this competition and resulting resource depletion, the perfect-match complexes approach equilibrium much more slowly in the multichannel calculation compared to the single-channel calculation.

For set II, when = 0.5 (lower panel of Fig. 4 d), competition from additional cross-hybrids formed with the aid of sufficient sequence homology acts to decrease levels of the perfect matches and the most significant cross-hybrids, compared to what was observed in the bottom panel of Fig. 4 b. In a dynamic sense, sequence homology and heightened stability of tandem mismatches acts to increase concentrations of all mismatch cross-hybrids. Only one hybrid complex (P3:T1) displays contrasting behavior, first attaining levels of the perfect-match duplexes, then decreasing exponentially to equilibrium. Note that there are only two mismatch complexes observable in Fig. 4 b, bottom panel. This enhanced cross-hybridization acts to decrease rates of hybridization for perfect match complexes (which must compete for resources with the cross-hybrids). Under this environment the equilibrium concentrations of the perfect match duplexes are also depressed.

In summary, assigning significant stability to tandem mismatches ( = 0.5) has a much larger effect on the kinetic behaviors and distributions of equilibrium concentrations of the perfect-match duplexes than sequence homology, which also enhances cross-hybridization. Although probably somewhat unrealistic in that the magnitude of tandem-mismatch interactions is surely not so large in all cases (as has been assumed), this example helps us discern the interrelated effects of sequence homology and tandem-mismatch stabilization. It further demonstrates relative sensitivity of the calculations to influences of different potential sources of cross-hybridization. Examples in Fig. 4, b and d, also reveal a different type of kinetic behavior of some mismatch complexes. In contrast to the typical exponential behavior, these curves increase rapidly then decrease exponentially. This behavior was further revealed when the temperature-dependence of the kinetic behavior for the perfect-match and mismatch complexes present in the 3 x 3 system was examined.

Effects of temperature
Effects of temperature on the kinetics of two-strand complex formation are depicted in Fig. 5. Each of the nine plots shown in the panels of these figures are of duplex concentrations versus time for the different probe:target complexes (probes are rows, targets are the columns) at different temperatures. Predictions from the single-channel (top) and multichannel (bottom) model calculations are shown. Plots on the diagonal are for the perfect matches, and off-diagonal plots are cross-hybrids. For example, the plot in the upper-left corner is for the P1:T1 complex (concentration versus time of the perfect-match complex formed from probe 1 and target 1). The plot in the first row, second column corresponds to the P1:T2 complex, and is the concentration as a function of time of the cross-hybrid complex formed from probe 1 and target 2. Likewise, the plot in the second row, third column is for the P2:T3 complex, the concentration versus time of the cross-hybrid formed from probe 2 and target 3. In each case, temperature was varied over the 60° range from 310 K (yellow) to 370 K (red). Note that the scales in the plots of Fig. 5 are not all equal. Different scales were used to emphasize and compare qualitative behaviors.

View larger version (36K): [in this window] [in a new window]   FIGURE 5  Kinetic plots (time versus concentration) for the 3 x 3 systems comprised of sequence sets I and II (Fig. 4) as a function of temperature from 37 (yellow) to 97°C (red). (a) Set I. Single-channel calculation (top panel) and multichannel calculation (bottom panel). Columns are probes, rows are targets. Thus, the upper-left plot in either panel corresponds to the P1:T1 complex, while the third column, second-row plot is for the P2:T3 complex. (b) Set II. Single-channel calculation (top panel) and multichannel calculation (bottom panel). All curves were generated assuming negligible stability for tandem mismatches ( = 0).

As shown in Fig. 5 b, results are different for the set II sequences that share from 25 to 54% homology with one another. For the single-channel calculation, results are essentially identical to what was observed for set I. However, in the multichannel calculation, sequence homology has a significant effect on the cross-hybrids containing that homology. In particular, the curves in Fig. 5 b for the P1:T2 and P2:T1 that share at least 50% homology are considerably different than seen in Fig. 5 a, where the sequences share no sequence homology. In essence, the kinetic behavior and level of cross-hybridization seen in these plots indicates the level of sequence homology present in the complexes.

A notable observation emerges from this analysis. The kinetic plots for the P1:T2 and P2:T1 complexes obtained from the multichannel calculation display two distinctly different behaviors. Similar observations were also made for these mismatches in the multichannel calculation (bottom panels of Fig. 4, b and d). These contrasting behaviors are referred to as "overdamped" and "critically damped" because of their similarity with the behavior of damped harmonic oscillators. For overdamped behavior, the concentration of hybridized complexes first increases at early times, then decreases exponentially with time, and levels off smoothly to some equilibrium value. In contrast, for critically damped behavior (what is typically observed) the hybridized concentration increases with time exponentially and converges to a constant value. In certain cases, depending on the temperature, both types of behavior can be observed for the same complex. The transition from critically damped to overdamped behavior seems to occur as temperature is increased, relative to the melting temperature T_m, for that particular complex. It can be shown that the maximum value of an overdamped duplex (one that increases at first, then decreases to an equilibrium value) is an unstable local equilibrium for the subsystem consisting of only that duplex and its probe:target pair. The exact mechanism underlying these observations, and when each type of kinetic behavior might be expected, is currently under investigation (D. J. Fish and A. S. Benight, unpublished).

In summary, effects of temperature on hybridization kinetics of mismatch complexes are manifest in two distinguishable kinetic regimes, termed "critically damped" and "overdamped". Which regime is occupied depends strongly on the relative stability of the hybrid duplex, i.e., the melting temperature of the complex relative to the assay temperature, and homology (cross-reactivity) with other sequences of the system.

Effects of differential strand concentrations
The influence of strand concentration and sequence homology on both the single channel and multichannel calculations, for sets I and II, is shown in Fig. 6. To explore the interdependent effects of concentration and sequence homology, three different scenarios were considered. In the first case (Fig. 6 a), concentration of target 1 was set to be 100 times higher than that for targets 2 and 3. In the second case (Fig. 6 b), concentration of target 2 was 100 times higher than targets 1 and 3. And in the third case (Fig. 6 c), concentration of target 3 was 100 times higher than targets 1 and 2. In each case, probe concentrations were equal to the concentration of the two lesser concentrated targets. On the plots in Fig. 6, multichannel calculations for set I (left) and set II (right) are displayed. Upper curves, provided for comparison, are the same in Fig. 6, a–c, and correspond to the situation in the multichannel calculation when concentrations of all strands are equal. In all cases, temperature was 315 K.

View larger version (57K): [in this window] [in a new window]   FIGURE 6  Kinetic plots (time versus concentration) for the 3 x 3 systems comprised of sequence sets I and II (Fig. 3) as a function of differential target concentrations at constant temperature, 42°C. Results from the multichannel calculation for set I are on the left, set II on the right. In the panels that are shown there are nine different plots. In these plots, columns are targets and rows are probes. Thus, the plot in the second column, third row corresponds to the P3:T2 complex, etc. Top panels are results when concentrations of the target and probe strands are all equivalent and are replicated in each figure for comparison. Bottom panels correspond to results obtained in the multichannel calculation for set I (left) and set II (right) when: (a) concentration of target 1 is 100 times that of the other target and probe strands, i.e., [T1] = 100 [T2] = 100 [T3]; (b) [T2] = 100 [T1] = 100 [T3]; and (c) [T3] = 100 [T1] = 100 [T2].

The situation is similar when [T2] = 100 [T1] = 100 [T3] (Fig. 6 b). For set I, the remarkable effect is to speed up duplex formation with increased equilibrium concentration of the P2:T2 perfect-match complex and the P3:T2 cross-hybrid complex. For set II, the P2:T2 perfect-match complex and the P2:T1 and P1:T2 cross-hybrids reach the highest equilibrium values. The P3:T1 and P3:T2 cross-hybrids also increase. These results are entirely analogous to those displayed in Fig. 6 a, except the complexes dominated by T2 (P1:T2 and P2:T2) have the highest rates of formation and reach the highest equilibrium values. The P3:T2 cross-hybrid also increases. Meanwhile, the P1:T1 and P2:T1 complexes, although they form faster, reach lower equilibrium values. Here, depletion of P1 and P2, through formation of the P1:T2 cross-hybrid and P2:T2 perfect match, respectively, driven by the relatively higher concentration of T2, results in lower equilibrium values of the P1:T1 and P2:T2 complexes. The above interpretations are validated by the plots in Fig. 6 c. Here, T3=100 [T1]=100 [T2]. Recall, in this case, T3 and P3 share only 25% homology with the other probes and targets in set II. As observed for the other situations, the effect of increasing T3 is to increase the rate of formation and final equilibrium value of the P3:T3 perfect-match complex. However, there is apparently little effect on formation of the cross-hybrids (regardless of sequence homology). Effects of resource depletion due to cross-hybrid formation are not observed for the higher (more nonhomologous) T3 target concentration.

Results of increasing the concentration of one of the target strands can be summarized as follows. Sequence homology alone has the largest effect on the cross-hybrid kinetics, but higher concentration of one strand enhances formation of perfect-match and cross-hybrid complexes. Formation of these complexes, driven by the strand in excess concentration, then acts to deplete the strands it binds from the reaction, thus reducing formation of other cross-hybrids or even perfect-match complexes.

	COMPARISONS WITH PREVIOUS WORK

TOP ABSTRACT INTRODUCTION THEORETICAL APPROACH EXAMPLE CALCULATIONS COMPARISONS WITH PREVIOUS WORK CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
In similar ways from different perspectives, several authors have approached the problem of defining hybridization error (41,61–74). The following provides brief descriptions and comparisons of some of these approaches, and how results generated from them compare with ours. In some cases, our results are identical; for others, some differences are observed.

Quantification of degeneracy and hybridization error
Rose and co-workers introduced the concept of hybridization error (67–72). Their approach is based on definition of the "computational incoherence." However, their definition did not take into account certain effects of concentration. It will be seen that consideration of concentration effects can have significant consequences on the hybridization error. This is clearly evident from the plots shown below in Fig. 7, where the calculated error determined with and without concentration considerations is plotted versus temperature (explained below). For multiplex hybridization reactions, the probability of error was defined (68,72) as

(16)

View larger version (9K): [in this window] [in a new window]   FIGURE 7  Comparison of estimates on hybridization error. Plots of the predicted error obtained from Eq. 17 of the text, under four different assumptions as a function of temperature from 300 to 400 K. (Solid line, full model; dashed line, full model minus hairpins; dot-dashed line, T2, Eq. 18; and dotted line, T3, Eq. 19.)

(17)

where is the equilibrium constant for the reaction involving formation of the duplex i