Copyright © 2006 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 91, Issue 11, 4133-4153, 1 December 2006
doi:10.1529/biophysj.106.090662
Nucleic Acids
Statistical Thermodynamics and Kinetics of DNA Multiplex Hybridization Reactions
M.T. Horne*, ‡, D.J. Fish‡ and A.S. Benight*, †, ‡, 
, 
* Department of Chemistry, Portland State University, Portland, Oregon
† Departments of Physics and Mathematics, Portland State University, Portland, Oregon
‡ Department of Portland Bioscience, Portland, Oregon
Address reprint requests to A. S. Benight, Tel.: 503-725-9513.
Abstract
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
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,2,3,4,5,6,7,8,9,10,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,13,14,15,16,17,18,19,20,21,22,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,25,26,27,28,29,30,31. Multiplex assays can be performed on microarrays or in solution via various iterations of the polymerase chain reaction 32,33,34,35,36,37,38,39,40,41,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,44,45, and studies performed for the purpose of investigating the thermodynamic parameters of single basepair mismatches 1,9,46, bulged loops 47,48,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,53,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
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 ith 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 pitj formed from a probe pi and a target tj, 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=swi/Z.
Consider a multiplex system wherein a set of probes, pi, is designed to specifically hybridize with particular targets, tj. 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, pitjσ is the ensemble of duplex states that form from single strands pi and tj 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., pitjσ could represent an overlap duplex state of a configuration involving single strands pi and tj, 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.
| | |
| Reaction | Equilibrium constant | Rate constants | |
|---|
|  |  |  | |
|  |  |  | |
|  |  |  | |
|  |  |  | |
|  |  |  | |
| | |
For the above equilibrium system, strand conservation for the probe and target strands is given by
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:
Define the following:
With these expressions, Eqs.
(1a) become
For simplicity, if target:target, probe:probe, and intramolecular states in single-strand targets and probes are ignored, i.e.,
, Eqs.
(3a) can be combined to form the following system of coupled, nonlinear equations:
If there are
NP probes and
NT targets, and
l goes from 1 to
NT, then Eq.
(4a) represents a system of
NP+
NT 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.
(2a), to determine the duplex concentrations,
.
Multiplex hybridization kinetics
Consider a multicomponent system comprised of two species types, i.e., probes pi, and targets tj. Each pi can react with every tj forming a pitj complex or can react with each of its counterparts of the same species, forming pipj and titj 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:
and
Here,
is the concentration of the
ith probe (
i≤
N),
is the concentration of the
jth 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 (
Cpp and
Ctt are symmetric), a system of
d=(1/2)(
N+
M)(
N+
M+5) equations is formed.
Observe that by combining Eqs. (5a), 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
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 suppressing the variable τ,
For convenience, denote the right-hand sides of Eqs. (8a) by 
and 
, respectively. Then Eqs. (6a) can be replaced by
Together with Eqs.
(8a), these define an equivalent system of
N+
M fewer equations.
Since the functions 
and 
are constant along solutions to Σ, Eqs. (8a) 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 Σ,
or, if
, where
x and
y are either
pi or
tj, then
Thus, any equilibrium solution to Σ must be contained in the set of common solutions to both systems of Eqs.
(8a). 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
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 pitj represents the duplex formed from single strands pi and tj, with similar definitions for the other reactions in Table 1 (neglecting hairpins in single strands). Following are alternate forms of the expressions in Eq. (5a):
Likewise, the following are alternate forms of expressions in Eq.
(6a):
Initial conditions for the system are
Cxyσ=0 and
Cxσ=0 (i.e., all strands are initially in single-strand conformations) and
Cz=
C0 (i.e., initial single-strand concentrations are given), where
x,
y,
z represent
pi or
tj. If misaligned, overlapping states are ignored, then
σ=1, in the expressions in Eq.
(13a). Assuming that
for all states
x, then
, in analogy with the expressions in Eqs.
(5a).
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. (5a) are amended to include the effect of hairpin formation, and rate equations for hairpin configurations are added to the set of expressions in Eq. (6a). As a result, Eqs. (12a) 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. (5a) 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. (12a) 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 RN+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.6GHz Pentium-IV machine with 1 Gigabyte of RAM, the algorithm required ∼45min to complete a 1000×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×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. (5a) 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 N2, and the size of the Jacobian matrix grows like N4, 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.6GHz Pentium-IV CPU and 1 Gigabyte of RAM, running Windows XP. Similar results were obtained when the simulations were carried out on a 1.8GHz/512 RAM machine.
Flow chart of the calculation scheme
The ultimate aim of developing this analytical foundation is to have the ability to diagnose each and every hybridization reaction that can conceivably occur in a multiplex reaction scheme. The equilibrium and kinetic models ultimately lead to quantitative predictions of the concentrations of probe:target complexes. Flow charts diagramming steps in the calculational schemes for both the kinetic and equilibrium models are shown in Figure 2ab.
Although some aspects of these models are very similar, the output produced by the two models is quite different. The kinetic model produces a time series of all concentrations involved in hybridization reactions, while the equilibrium model produces only the equilibrium values of these concentrations (i.e., τ=∞). Ideally, the kinetic model could be applied to any system, but computational complexity and hardware limitations prohibit this option for very large systems. On the other hand, for large systems, full knowledge of the concentration levels at all times is perhaps unnecessary, and concentrations at different times before equilibrium and the final equilibrium values may be sufficient for many desired applications. In general, combined use of both models is probably most appropriate and expected to produce the best results.
Example calculations
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×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×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 Figure 3a, 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.
Calculation of duplex thermodynamics and stability
For each duplex that can form from any pair of the three probes and three targets in sets I and II, thermodynamic transition parameters, ΔH, ΔS, and ΔG, used in kinetic and equilibrium model calculations, were determined from published sequence-dependent thermodynamic parameters 9,10,46,60,61,62,63,64. For each pair of strands all overlap complexes were considered and their corresponding thermodynamic quantities were calculated for the particular duplex sequence at each overlap. Of these states at different alignments, the one having the lowest calculated free energy was selected. For this state, the total helix-coil transition thermodynamics were calculated from the sum of appropriate nearest-neighbor sequence-dependent parameters, where available 65. For example, consider the following hybrid duplex sequence and its decomposition into nearest-neighbor components of the enthalpy, ΔH (mismatches are underlined):
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
Table 2,
Table 3,
Table 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
Table 2,
Table 3,
Table 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=Δ
HTΔ
S, and
Tm=Δ
H/Δ
S, Δ
G=Δ
H(1−
T/
Tm).
| | |
| W/C doublet | Enthalpy(cal/mol) | Entropy (cal/kmol) | |
|---|
| AA | −7900 | −22.2 | |
| AC | −8400 | −22.4 | |
| AG | −7800 | −21.0 | |
| AT | −7200 | −20.4 | |
| CA | −8500 | −22.7 | |
| CC | −8000 | −19.9 | |
| CG | −10,600 | −27.2 | |
| GA | −8200 | −22.2 | |
| GC | −9800 | −24.4 | |
| TA | −7200 | −21.3 | |
| | |
| | |
| Dangling end | Enthalpy (cal/mol) | Entropy (cal/Kmol) | Dangling end | Enthalpy (cal/mol) | Entropy (cal/Kmol) | |
|---|
| TA/-T | −6900 | −20.0 | CG/G- | −3200 | −10.4 | |
| AC/-G | −6300 | −17.1 | AG/-C | −3700 | −10 | |
| CA/G- | −5900 | −16.5 | AT/-A | −2900 | −7.6 | |
| GT/-A | −4200 | −15.0 | CC/G- | −2600 | −7.4 | |
| CT/G- | −5200 | −15.0 | -C/AG | −2100 | −3.9 | |
| GC/-G | −5100 | −14.0 | TG/A- | −1600 | −3.6 | |
| TG/-C | −4900 | −13.8 | GA/-T | −1100 | −1.6 | |
| AG/T- | −4100 | −13.1 | AA/T- | −500 | −1.1 | |
| -C/TG | −4400 | −13.1 | TA/A- | −700 | −0.8 | |
| CT/-A | −4100 | −13.0 | TT/-A | −200 | −0.5 | |
| CC/-G | −4400 | −12.6 | -C/CG | −200 | −0.1 | |
| AT/T- | −3800 | −12.6 | AA/-T | 200 | 2.3 | |
| CG/-C | −4000 | −11.9 | CA/-T | 600 | 3.3 | |
| -C/GG | −3900 | −11.2 | TT/A- | 2900 | 10.4 | |
| GG/-C | −3900 | −10.92 | AC/T- | 4700 | 14.2 | |
| TC/-G | −4000 | −10.9 | TC/A- | 4400 | 14.9 | |
| | |
| | |
| MM | Enthalpy (cal/mol) | Entropy (cal/Kmol) | MM | Enthalpy (cal/mol) | Entropy (cal/Kmol) | MM | Enthalpy (cal/mol) | Entropy (cal/Kmol) | |
|---|
| GC/GG | −6000 | −15.8 | CC/GT | −800 | −4.5 | GA/GT | 1600 | 3.6 | |
| CT/GT | −5000 | −15.8 | AC/CG | −700 | −3.8 | CA/GC | 1900 | 3.7 | |
| CG/GG | −4900 | −15.3 | AG/TA | −700 | −2.3 | AA/TC | 2300 | 4.6 | |
| GC/TG | −4400 | −12.3 | CA/GG | −700 | −2.3 | GC/CT | 2300 | 5.4 | |
| CG/GT | −4100 | −11.7 | AA/TG | −600 | −2.3 | AA/GT | 3000 | 7.4 | |
| AG/GC | −4000 | −13.2 | GA/CG | −600 | −1.0 | GG/CT | 3300 | 10.4 | |
| AG/TG | −3100 | −9.5 | TA/GT | −100 | −1.7 | CA/AT | 3400 | 8.0 | |
| AC/AG | −2900 | −9.8 | AC/TC | 0 | −4.4 | CC/CG | 3600 | 8.9 | |
| CT/GG | −2800 | −8.0 | TA/TT | 200 | −1.5 | GT/TG | 4100 | 9.5 | |
| AT/TT | −2700 | −10.8 | AC/GG | 500 | 3.2 | AA/AT | 4700 | 12.9 | |
| AT/TG | −2500 | −8.3 | AG/CC | 600 | −0.6 | CC/AG | 5200 | 14.2 | |
| GT/CT | −2200 | −8.4 | AC/TT | 700 | 0.2 | CC/TG | 5200 | 13.5 | |
| CC/GC | −1500 | −7.2 | GA/AT | 700 | 0.7 | GA/CC | 5200 | 14.2 | |
| CG/TC | −1500 | −6.1 | AG/TT | 1000 | 0.9 | AC/TA | 5300 | 14.6 | |
| TT/AG | −1300 | −5.3 | TT/AC | 1000 | 0.7 | GG/TT | 5800 | 16.3 | |
| AT/TC | −1200 | −6.2 | AA/TA | 1200 | 1.7 | CA/CT | 6100 | 16.4 | |
| AG/AC | −900 | −4.2 | TA/CT | 1200 | 0.7 | AA/CT | 7600 | 20.2 | |
| | |
For strands resulting in hybrid complexes containing tandem mismatches, quantitative prediction of their stabilities is a little less certain. Although there are nearest-neighbor parameters for single basepair mismatches for nearly all of the possible nearest-neighbor combinations (Table 4 below) 10, a parameter set for tandem mismatches does not currently exist. In some of the example calculations that were performed, the influence of the relative thermodynamic stability of tandem mismatches was investigated. For the purpose of example, tandem mismatches were assigned thermodynamic parameter values that were a fraction of the corresponding w/c basepair doublet values, i.e., the free-energy of a mismatch basepair doublet in a tandem mismatch complex was assigned according to
where Δ
GPM, Δ
HPM, and Δ
SPM 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 Figure 4ad, show the results of the single channel (top panels) and multichannel (bottom panels) calculations for a simple 3×3 system. Kinetic curves (duplex concentration versus time) that are displayed were calculated for set I (Figure 4ac) and set II (Figure 4bd) 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 Figure 4ab, were obtained by assuming the free energy of tandem mismatches as 0, i.e., κ=0, while those in Figure 4cd, were obtained assuming a universal constant factor, κ=0.5 (to scale tandem-mismatch free energies relative to normal w/c basepairs).
Examination and comparison of the plots in Fig. 4 reveals specific features of the model calculations and how they are influenced by sequence homology among the probes and targets, and tandem mismatch stability. First, compare the curves in Figure 4ab (κ=0). Figure 4a is for set I, where there is no sequence homology between the different probe and target sequences. In the top panel for the single-channel calculation, three curves are shown corresponding to the three perfect match duplexes. These curves exhibit typical exponential behavior. Plots in the bottom panel correspond to the multichannel calculation where colored lines are the perfect matches and gray lines are the mismatched cross-hybrids. Although there are nine curves generated in the multichannel calculation for set I, the curves for mismatches (or cross-hybrids) lie along the horizontal axis and are not visible. As seen by comparison of the curves in the upper and lower panels of Figure 4a, for set I, kinetic curves for the perfect match duplexes are essentially identical in both the single-channel and multichannel calculations, and display the typical exponential behavior. The comparison in Figure 4a shows when there is no sequence homology between the different probe and target strands, and there is no extra weighting for tandem mismatches (κ=0), the single-channel and multichannel results for the perfect match duplexes are indistinguishable.
Calculated kinetic curves for set II are shown in the upper and lower panels in Figure 4b. 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 Figure 4b 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 Figure 4b 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, Figure 4b) 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 Figure 4b for the dominant mismatch complexes, P1:T2 and P2:T1 (black and orange lines in the lower panel of Figure 4b), 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, Figure 4b) 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 Figure 4b 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
