- Email: [email protected]
Application of a coarse-grained model for DNA to homo- and heterogeneous melting equilibria
Chemical Physics Letters 485 (2010) 354–359
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Application of a coarse-grained model for DNA to homo- and heterogeneous melting equilibria Nicholas B. Tito, John M. Stubbs * Department of Chemistry and Physics, University of New England, 11 Hills Beach Road, Biddeford, ME 04005, USA
a r t i c l e
i n f o
Article history: Received 22 August 2009 In final form 29 December 2009 Available online 4 January 2010
a b s t r a c t Configurational-bias Monte Carlo simulations were carried out on deoxyribonucleic acid (DNA) decamers using a coarse-grained molecular model. The effects of single mutations on the melting transition were investigated as were heterogeneous systems with immobilization of one strand on a surface, both with and without a spacer. The destabilizing effect of an internal mutation is attributed to a lack of cooperativity, which acts through a hydrogen bonding nucleotide’s restriction of the conformational freedom of neighboring bases. A surface-oligomer spacer is necessary for duplex stability with the destabilizing effect of the surface coinciding with the volume it excludes. Ó 2009 Elsevier B.V. All rights reserved.
1. Introduction Hybridization and denaturation or ‘melting’ transitions between single- and double-stranded forms of deoxyribonucleic acid (DNA) are an important property of DNA behavior. The recent development of DNA sensor microarrays has the potential to significantly benefit important applications, such as mutation detection, genetic sequencing, species discrimination and RNA expression [1–5]. These applications rely on knowledge of the thermodynamics and kinetics of DNA melting and hybridization transitions and are sensitive to many experimental factors such as concentration, solvent environment, ionic strength, sequence and temperature. Although these systems have been studied extensively, several questions remain unanswered including how changes in physical environment for DNA oligomers between solution and surface-bound transitions influence hybridization, helix–coil equilibria or competitive adsorption of closely-related sequences [6]. DNA sensors are composed of single-stranded DNA oligomers, or probes, bound to a surface which hybridize with a complementary oligomer, or target, in solution. DNA melting is a phenomenon that occurs when hybridized double-stranded DNA helices dissociate into single-stranded DNA coils. The thermodynamic melting temperature, Tm, is defined as the temperature at which the double strand is exactly half-dissociated and can be measured quite accurately experimentally [7] via hyperchromicity in ultraviolet absorption or differential scanning calorimetry. However, the mechanisms and structures governing the melting process are not as well understood and while there are extensive measurements of melting temperatures of duplexes in solution, there is * Corresponding author. Fax: +1 207 602 5993. E-mail address: [email protected] (J.M. Stubbs). 0009-2614/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2009.12.079
much less data available about heterogeneous melting transitions on a surface, which are complicated by additional factors such as length of spacer between the surface and probe oligomer [8,9] and the density of probe sites on the surface [10,11]. Although in solution it is in general possible to accurately predict melting temperatures for a given sequence [12–14], this approach cannot be extended to surface-bound equilibria without the input of significant, additional experimental data. Recently Zhang et al. [15,16] have applied weighting factors to a nucleotide’s hybridization interaction depending upon its position in the sequence of a surface-bound strand relative to the surface itself. They were able to improve predictions of binding thermodynamics and hence amount of target binding, though certain poorly characterized outlier sequences existed, attributed to predominantly G-rich sequences and multiplex formation [16]. Another problem for binding prediction is the effect of base pair mismatches, or mutations, on Tm. For example, a 25-base oligomer with a single mutation is expected to have very similar hybridization thermodynamics compared to the perfectly complementary strand due to the rather small energetic differences. However, it has been observed that for these same sequences – as targets in a surface hybridization experiment – the single mutant strand does not achieve the same efficiency (100%) as the perfectly complementary strand does even after long hybridization times, potentially indicating the presence of a kinetic barrier to hybridization [17]. It has also been shown that the location of a mismatch along the strand plays a role in the efficiency with mismatches in the center having the greatest destabilizing effect on hybridization [18]. Levicky and Horgan [6] have recently outlined the challenges for microarray technology from a physical chemistry perspective and conclude that a more thorough understanding is needed for many aspects including the origin of differences between solution and surface-bound DNA hybridization equilibrium constants as
355
N.B. Tito, J.M. Stubbs / Chemical Physics Letters 485 (2010) 354–359
well as how mutation discrimination is affected by solution conditions and surface probe construction. Molecular simulation can provide valuable microscopic insight into these systems which can lead to better understanding of the differences between solution- and surface-based DNA hybridization, which could in turn ultimately lead to more robust predictive models. Although a few surface simulation studies have been carried out investigating melting of a single strand [19], optimum surface strand density [20], and molecular recognition [21], a systematic simulation study of hybridization and melting transitions at equilibrium has not yet been carried out, which is the primary focus of this work. 2. Methodology 2.1. Molecular models The coarse-grained molecular model is largely identical to that of Drukker et al. [22] which is a representation of a B-DNA helix. The principal modification is that a reduced hydrogen bonding interaction was necessary to reproduce estimated Tm values based on initial results, which showed that reduction of the hydrogen bonding interaction by 20% was necessary for agreement [23]. Additionally, non-complementary nucleotides are taken to not interact via hydrogen bonding. For surface simulations, a hard surface model was employed (i.e. any Monte Carlo move that caused a DNA interaction site to overlap with the surface was rejected). The surface-bound strand had its movement restricted by not allowing the tethered backbone or spacer bead to move from its position 5 Å above the surface but otherwise allowing all other moves, thus fixing it in place on the surface. For simulations with spacers between the surface and strand, each spacer unit had the same interaction parameters as a backbone unit and thus represents a unit consisting of phosphate and deoxyribose subunits, approximately the same size as the ð—CH2 CH2 CH2 —PO 4 —Þn units described in the work of Shchepinov et al. [8].
available for hydrogen bonding, additional surface-bound simulations were carried out without the second, target, strand i.e. only the probe strand is present. The same temperature range was investigated for all three immobilization scenarios of no spacer, a 5-unit spacer and a 10-unit spacer. Each simulation was equilibrated for at least 108 Monte Carlo (MC) cycles after which at least 3.2 108 MC cycles were carried out for production where one MC cycle consists of Nmolecule MC steps. The MC simulations used three different types of trial moves: translation of a molecule, rotation of a molecule, and conformational moves via coupled-decoupled configurational-bias Monte Carlo (CBMC) [26–28] and self-adapting fixed endpoint CBMC [29] to ensure adequate sampling of interior degrees of freedom. The type of MC move was selected at random using a fixed set of probabilities equally divided between translational, rotational and conformational moves. The maximum translational and rotational displacements were adjusted during the equilibration part in order to yield 50% acceptance rates for each type of molecule in the simulations. All error bars depict the standard error of the mean of block averages from a single simulation trajectory with blocks 1.6 107 MC cycles in size. 3. Results and discussion 3.1. Hydrogen bonding For this model a pair of complementary bases was determined to be hydrogen bonded if the energy of their interaction was less than 200 K 1.7 kJ mol1. This value was chosen to exclude those pairs that have some interaction but are significantly further apart then the more energetically stable hydrogen-bonded configurations. Note that the interactions have been parameterized to reproduce the experimental melting temperature; for comparison, changing the criterion to more strict values of 700 K for AT and 1300 K for GC pairs had no effect on the resulting Tm values but did not allow some helical structures to be counted as completely hydrogen bonded even at low temperatures.
2.2. Simulation details Three duplex systems were studied in detail as listed below: a fully complementary (perfect match or PM) decamer with moderate GC content (40%) for which an experimental Tm is known [24] as well as two possible single mutations from G to A of the PM sequence, with either a terminal mismatch (M1) or interior mismatch (M2). The mismatch is given in bold below:
PM : 50 CCAACTTCTT30 30 GGTTGAAGAA50
3.2. Solution results Fig. 1 shows the average number of hydrogen bonding base pairs, NHB, as a function of temperature for the PM, M1 and M2 sequences. Values of Tm are determined by estimating the temperature at which half the base pairs are hydrogen bonding, i.e. when
10
M1 : 50 CCAACTTCTT30 8
30 AGTTGAAGAA50 0
0
M2 : 5 CCAACTTCTT3
30 GGTTAAAGAA50
N HB
Simulations were carried out in the canonical (NVT) ensemble using periodic boundary conditions [25] at various temperatures between 285 K and 355 K and a concentration of 1.66 mM for each strand, corresponding to a cubic simulation cell 100 Å on each side. For surface-bound studies, the upper sequence (as given above) was immobilized at its 50 end, and the surface was taken to have a thickness of 14 Å. The z-dimension was extended by 14 Å to maintain a constant volume for the strands thus corresponding to a simulation cell with the dimensions 100 Å 100 Å 114 Å. Long range interactions were truncated at 14 Å so that no interactions occurred from one side of the surface to the other. To investigate the effect of the surface on the conformations of the bound strand that are
6
4
2
0 280
300
320
340
360
T (K) Fig. 1. Average number of hydrogen bonding base pairs as a function of temperature for solution simulations. Circles, diamonds and triangles correspond to PM, M1 and M2 sequences, respectively. Lines are drawn as guides for the eye.
356
N.B. Tito, J.M. Stubbs / Chemical Physics Letters 485 (2010) 354–359
the curves cross NHB = 5. As the range of temperature values was covered in 5 or 10 K increments, the uncertainty in the simulation Tm values is likely approximately 5 K. Estimated Tm values were determined using the unified parameters of the empirical nearest-neighbor model [13] for both the simulated and experimental (single strand) concentrations of 1.66 mM and 1 lM, respectively, and Table 1 summarizes this data along with the experimental value. As shown in Table 1, the empirical model value matches the experimental value to within 2 K, and the simulation values match the empirical model values to within 4 K for all 3 sequences. Thus, these solution-based results are validated directly against a robust predictive model and indirectly against experiment. Fig. 2 shows the fraction of configurations where a nucleotide is hydrogen bonded, fHB, as a function of its location in the strand where the nucleotides are numbered starting from the 30 end of the lower sequence as given previously, and as a function of temperature. For clarity, error bars are omitted. Note that the Tm can be seen in this figure as the sharp drop-off between 325 and 330 K. As can be seen, the two ends have on average noticeably fewer hydrogen bonds even at low temperatures, which corresponds to end fraying of the duplex. It is more pronounced for AT pairs which is expected since their interaction is not as strong as that of GC pairs. That this behavior is present in the simulation results despite having the same interaction parameters as interior AT or GC pairs indicates that an important factor in maintaining nucleotide hydrogen bonding is for a nucleotide to have its available conformations restricted, i.e. be held in place by hydrogen bonding neighbors. Because the ends lack a neighbor on one side, they are somewhat destabilized. Another aspect of this result is that this model implicitly captures cooperative effects due to hydrogen bonding neighbors, a key idea of the empirical nearest neighbor model [12,13]. This same effect can also be seen for the neighbors of the mismatched bases, where the effect of being adjacent to a mismatch lowers fHB by between approximately 5–10% (Fig. 3a and b) at 285 K as compared to the PM sequence (Fig. 2). The interior mutation (M2) is shown to be more destabilizing than the terminal mutation (M1), as is expected. These results show that the destabilizing effect of a mutation on neighbors can be reproduced purely through a mutant’s inability to hold its neighbors in place, even when the neighbors have the same interaction potential, a situation analogous to terminal bases with no neighbor on one side. 3.3. Surface-bound results Graphs of NHB as a function of temperature for all three sequences where the 50 end of the upper strand (as given previously)
Table 1 Summary of Tm where [expt] and [sim] denote the experimental and simulation concentrations, respectively, and Expt, NN and Sim denote experimental, empirical nearest-neighbor model and simulation results. Sequence
Method
System
Tm (K)
PM PM PM M1 M2 PM M1 M2 PM M1 M2 PM M1 M2
Expt NN NN NN NN Sim Sim Sim Sim Sim Sim Sim Sim Sim
Solution [Expt] [Sim] [Sim] [Sim] Solution Solution Solution 5-Spacer 5-Spacer 5-Spacer 10-Spacer 10-Spacer 10-Spacer
305.3 303.5 323.3 318.9 300.3 327 319 299 331 326 305 328 326 311
fHB
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
280 290 300 310 T (K) 320 330 340 350 360
1
2
3
4
5
6
7
8
9
10
Base Index
Fig. 2. Fraction of configurations in which a particular nucleotide is hydrogen bonded, with the base index starting at the 50 end of the upper strand, as a function of temperature for the upper strand in the PM solution results.
is immobilized directly above a hard surface are given in Fig. 4, and an analogous plot corresponding to the 5-unit spacer between the immobilized strand and the surface is given in Fig. 5. Results for the 10-unit spacer are similar to those of the 5-unit spacer and are not shown. Fig. 4 indicates that without a spacer, none of the three sequences forms more than 4–6 hydrogen-bonded pairs even at low temperatures. Plots of fHB (not shown) show that without a spacer the distal half of the strand hydrogen bonds normally (as in the solution case) but that the half adjacent to the surface has negligible hydrogen bonding. On the other hand, Fig. 5 indicates that with a spacer of at least 5 units long the corresponding Tm values are almost identical to the non-immobilized solution-based results (Table 1). To investigate the effect of the surface on the probe strand’s ability to form a complete duplex with the target, an analysis was performed on the probe-only simulations to find the percentage of probe conformations that disallowed complete hybridization to a target strand due to the surface. The analysis was carried out by taking each conformation of the probe and determining the hypothetical location of the target strand that would have a position consistent with optimum hydrogen bonding to the probe. If this hypothetical target position overlapped with the surface, then this probe conformation was counted as non-viable with respect to duplex formation. Averaging over all conformations and temperatures, those directly adjacent to the surface (no spacer) were non-viable in 27.6% of all configurations, whereas the 5- and 10-unit spacers were only non-viable in 8.5% and 4.8% of all configurations, respectively. These values show no temperature dependence over the range investigated here, 285–340 K. This analysis indicates that the destabilizing effect of the surface is due primarily to the decrease in number of probe conformations that allow the target strand to completely hybridize with it, i.e. the presence of the surface excludes a set of probe conformations from forming a duplex. In contrast, distancing the probe from the surface by a spacer reduces the number of forbidden conformations considerably, with little difference between the two spacers investigated here. Another interesting aspect is that the spacer simulations show that the PM and M1 sequences have slightly more similar melting curves than their solution results, Figs. 1 and 5. This could be due to the surface having a destabilizing effect on an immobilized strand’s nucleotides closest to the surface, and is more clearly shown in the fHB versus nucleotide plots given in Fig. 6a and b compared to Figs. 2 and 3a. Thus, the destabilizing influence of a
357
N.B. Tito, J.M. Stubbs / Chemical Physics Letters 485 (2010) 354–359
A
B
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
fHB
fHB
280 290 300 310 T (K) 320 330 340 350 360
1
2
3
4
5
6
7
8
9
10
Base Index
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
280 290 300 310 T (K) 320 330 340 350 360
1
2
3
4
5
6
7
8
9
10
Base Index
Fig. 3. Fraction of configurations in which a particular nucleotide is hydrogen bonded, with the base index starting at the 50 end of the upper strand, as a function of temperature for the upper strand in the M1 (left, A) and M2 (right, B) solution results.
For simulations of surface-bound probes with either 5- or 10unit spacers, the amount of end fraying present (from fHB plots) is approximately the same as in the solution simulations. Consistent with experimental data [8], the overall effect of the spacers is that while complete lack of one prevents full duplex formation even at low temperatures, there is relatively little difference between 5-units and 10-units, and they both give similar melting temperatures to the solution results.
10
8
N HB
6
4
3.4. Thermodynamic quantities 2
0 280
For the association between two single strands A and B into a double-stranded duplex AB the equilibrium constant can be written as 300
320
340
360
T (K)
Kc ¼
Fig. 4. Average number of hydrogen bonding base pairs as a function of temperature for surface-bound probes. Styles as in Fig. 1.
10
ðcAB =c0 Þ ðcA =c0 ÞðcB =c0 Þ
where cAB, cA and cB are the concentrations of the duplex AB, strand A and strand B, respectively, and c0 is the standard state concentration and is taken to be 1 mol L1. This can then be combined with the relationships
DA0 ¼ RT ln K c 8
ð1Þ
ð2Þ
and
DA0 ¼ DU 0 T DS0
N HB
6
ð3Þ
to yield 4
ln K c ¼ 2
0 280
300
320
340
360
T (K) Fig. 5. Average number of hydrogen bonding base pairs as a function of temperature for surface-bound probes with a 5-unit spacer. Styles as in Fig. 1.
mismatch is decreased if that nucleotide is less able to hydrogen bond in the first place due to surface proximity. However, the effect is not large and could also be due to the uncertainties in the melting curve data.
DU 0 1 DS0 þ R T R
ð4Þ
where the constant volume conditions imply the use of the Helmholtz (A) and internal (U) energies. Thus, under the assumption that DU0 and DS0 are temperature independent, ln Kc will vary linearly with 1/T and will have a coefficient of DU0/R and constant of DS0/R. Carrying out this analysis for PM duplexes in both solution and surface simulations results in the values listed in Table 2 with error estimates based on the standard error of the linear regression analysis. As can be seen the interaction energy is in general substantially decreased for all surface-bound systems with the no spacer case the most affected at approximately half the solution interaction value. There is little difference between the 5- and 10-unit spacer results which recover much of the solution interaction. The destabilizing effect of the surface is thus most severe when there
358
N.B. Tito, J.M. Stubbs / Chemical Physics Letters 485 (2010) 354–359
A fHB
B
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
280 290 300 310 T (K) 320 330 340 350 360
fHB
1
2
3
4
5
6
7
8
10
9
Base Index
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
280 290 300 310 T (K) 320 330 340 350 360
1
2
3
4
5
6
7
8
9
10
Base Index
Fig. 6. Fraction of configurations in which a particular nucleotide is hydrogen bonded, with the base index starting at the 50 end of the upper strand, as a function of temperature for the upper strand in the PM (left, A) and M1 (right, B) 5-unit spacer results.
Table 2 Simulation results for thermodynamic quantities of DNA hybridization for the PM sequence. System
DU0 (kJ mol1)
DS0 (J K1 mol1)
Solution No Spacer 5-Spacer 10-Spacer
217 ± 19 100 ± 11 160 ± 11 165 ± 11
559 ± 61 264 ± 34 359 ± 33 372 ± 36
is no spacer and partially ameliorated with the presence of at least a 5-unit spacer. The effect on the entropy of association, DS0, is similar with the largest value for the solution simulations, the smallest value for the no spacer surface-bound system and the 5- and 10-unit spacer systems at intermediate values. This is consistent with expectations based on the amount of constraint placed upon the probe strand where direct immobilization on a surface results in the largest change in conformational, rotational and translational entropy relative to the unbound system and the presence of a spacer partially mitigating these changes. As DS0 is less negative for the surface-bound systems, the entropic penalty to form a duplex if one strand is bound to a surface is less than the penalty when both single strands are in solution. Since for the association of A and B
DS0 ¼ S0AB ðS0A þ S0B Þ
4. Summary The application of a coarse-grained DNA model to three DNA decamer systems is shown to give a consistent Tm value with respect to experiment. Additionally, the model rationalizes experimental phenomena in terms of molecular-level events. Simulation results allow for the microscopic interpretation of phenomena such as nucleotide neighbor destabilization due to mutations through increased conformational freedom. The destabilizing nature of the surface in this context corresponds well to the decrease in conformations of the surface-tethered probe strand that allow for complete duplex formation without overlap with the surface. Finally, both the energy and entropy of duplex formation were found to decrease in magnitude when the probe strand is surfacebound. Acknowledgments Financial support from a UNE Presidential Faculty Scholarship Mini-Grant and Grinnell College Faculty Scholarship Grants is gratefully acknowledged. We thank J.I. Siepmann for helpful discussions and suggestions and the Minnesota Supercomputing Institute where the calculations were partially carried out.
ð5Þ
where strands A and B are the probe and target, respectively, if S0B is assumed to be constant regardless of whether strand A is surfacebound or not, comparing two associations gives
DDS0 ¼ DS0AB DS0A
for the duplex when surface-bound and unable to completely associate along its entire length due to the presence of the surface.
ð6Þ
where
DDS0 ¼ DS0surface-bound DS0solution
ð7Þ
DS0AB ¼ S0AB;surface-bound S0AB;solution
ð8Þ
DS0A ¼ S0A;surface-bound S0A;solution
ð9Þ
Since DDS0 > 0 for all surface-bound associations (Table 2), based on these results the effect of surface binding on the probe strand’s entropy, S0A , is much greater than the effect of surface binding (via one strand) on the entropy of the duplex, S0AB . This is no doubt in part due to the increased number of conformational states possible
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
E. Southern, K. Mir, M. Shchepinov, Nat. Genet. (Suppl. 2) (1999) 5. P.O. Brown, D. Botstein, Nat. Genet. 21 (1999) 33. R.J. Lipshutz, S.P.A. Fodor, T.R. Gingeras, D.J. Lockhart, Nat. Genet. 21 (1999) 20. D.J. Lockhart, E.A. Winzeler, Nature 405 (2000) 827. A. Pozhitkov et al., Nucl. Acids Res. 34 (2006) e66. R. Levicky, A. Horgan, Trends Biotechnol. 23 (2005) 143. W. Saenger, Principles of Nucleic Acid Structure, Springer-Verlag, New York, 1984. p. 143. M.S. Shchepinov, S.C. Case-Green, E.M. Southern, Nucl. Acids Res. 25 (1997) 1155. E.L.S. Wong, E. Chow, J.J. Gooding, Langmuir 21 (2005) 6957. A.W. Peterson, R.J. Heaton, R.M. Georgiadis, Nucl. Acids Res. 29 (2001) 5163. M.L. Del Giallo, F. Lucarelli, E. Cosulich, E. Pistarino, B. Santamaria, G. Marrazza, M. Mascini, Anal. Chem. 77 (2005) 6324. J. SantaLucia Jr., Proc. Natl. Acad. Sci. USA 95 (1998) 1460. J. SantaLucia Jr., D. Hicks, Annu. Rev. Biophys. Biomol. Struct. 33 (2004) 415. R. Owczarzy, P.M. Vallone, F.J. Gallo, T.M. Paner, M.J. Lane, A.S. Benight, Biopolymers (Nucl. Acid Sci.) 44 (1997) 217. L. Zhang, M.F. Miles, K.D. Aldape, Nat. Biotechnol. 21 (2003) 818. L. Zhang, C. Wu, R. Carta, H. Zhao, Nucl. Acids Res. 35 (2007) e18. A.W. Peterson, L.K. Wolf, R.M. Georgiadis, J. Am. Chem. Soc. 124 (2002) 14601.
N.B. Tito, J.M. Stubbs / Chemical Physics Letters 485 (2010) 354–359 [18] L.M. Wick, J.M. Rouillard, T.S. Whittam, E. Gulari, J.M. Tiedje, S.A. Hashsham, Nucl. Acids Res. 34 (2006) e26. [19] K. Qamhieh, K.-Y. Wong, G.C. Lynch, B.M. Pettitt, Int. J. Numer. Anal. Model. 6 (2009) 474. [20] L. Yao, J. Sullivan, J. Hower, Y. He, S. Jiang, J. Chem. Phys. 127 (2007) 195101. [21] A. Jayaraman, C.K. Hall, J. Grenzer, Biophys. J. 91 (2006) 2227. [22] K. Drukker, G. Wu, G.C. Schatz, J. Chem. Phys. 114 (2001) 579. [23] M. Buchanan, Summer Chemistry Research Program Final Report, Grinnell College, Grinnell, IA, 2006.
359
[24] R. Owczarzy, Y. You, G. Moreira, J.A. Manthey, L. Huand, M.A. Behlke, J.A. Walder, Biochem 43 (2004) 3537. [25] M.P. Allen, D.J. Tildesley, Computer Simulations of Liquids, Oxford University Press, Oxford, 1987. [26] J.I. Siepmann, Mol. Phys. 70 (1990) 1145. [27] J.I. Siepmann, D. Frenkel, Mol. Phys. 75 (1992) 59. [28] M.G. Martin, J.I. Siepmann, J. Phys. Chem. B 103 (1999) 4508. [29] C.D. Wick, J.I. Siepmann, Macromolecules 33 (2000) 7207.
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Application of a coarse-grained model for DNA to homo- and heterogeneous melting equilibria Nicholas B. Tito, John M. Stubbs * Department of Chemistry and Physics, University of New England, 11 Hills Beach Road, Biddeford, ME 04005, USA
a r t i c l e
i n f o
Article history: Received 22 August 2009 In final form 29 December 2009 Available online 4 January 2010
a b s t r a c t Configurational-bias Monte Carlo simulations were carried out on deoxyribonucleic acid (DNA) decamers using a coarse-grained molecular model. The effects of single mutations on the melting transition were investigated as were heterogeneous systems with immobilization of one strand on a surface, both with and without a spacer. The destabilizing effect of an internal mutation is attributed to a lack of cooperativity, which acts through a hydrogen bonding nucleotide’s restriction of the conformational freedom of neighboring bases. A surface-oligomer spacer is necessary for duplex stability with the destabilizing effect of the surface coinciding with the volume it excludes. Ó 2009 Elsevier B.V. All rights reserved.
1. Introduction Hybridization and denaturation or ‘melting’ transitions between single- and double-stranded forms of deoxyribonucleic acid (DNA) are an important property of DNA behavior. The recent development of DNA sensor microarrays has the potential to significantly benefit important applications, such as mutation detection, genetic sequencing, species discrimination and RNA expression [1–5]. These applications rely on knowledge of the thermodynamics and kinetics of DNA melting and hybridization transitions and are sensitive to many experimental factors such as concentration, solvent environment, ionic strength, sequence and temperature. Although these systems have been studied extensively, several questions remain unanswered including how changes in physical environment for DNA oligomers between solution and surface-bound transitions influence hybridization, helix–coil equilibria or competitive adsorption of closely-related sequences [6]. DNA sensors are composed of single-stranded DNA oligomers, or probes, bound to a surface which hybridize with a complementary oligomer, or target, in solution. DNA melting is a phenomenon that occurs when hybridized double-stranded DNA helices dissociate into single-stranded DNA coils. The thermodynamic melting temperature, Tm, is defined as the temperature at which the double strand is exactly half-dissociated and can be measured quite accurately experimentally [7] via hyperchromicity in ultraviolet absorption or differential scanning calorimetry. However, the mechanisms and structures governing the melting process are not as well understood and while there are extensive measurements of melting temperatures of duplexes in solution, there is * Corresponding author. Fax: +1 207 602 5993. E-mail address: [email protected] (J.M. Stubbs). 0009-2614/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2009.12.079
much less data available about heterogeneous melting transitions on a surface, which are complicated by additional factors such as length of spacer between the surface and probe oligomer [8,9] and the density of probe sites on the surface [10,11]. Although in solution it is in general possible to accurately predict melting temperatures for a given sequence [12–14], this approach cannot be extended to surface-bound equilibria without the input of significant, additional experimental data. Recently Zhang et al. [15,16] have applied weighting factors to a nucleotide’s hybridization interaction depending upon its position in the sequence of a surface-bound strand relative to the surface itself. They were able to improve predictions of binding thermodynamics and hence amount of target binding, though certain poorly characterized outlier sequences existed, attributed to predominantly G-rich sequences and multiplex formation [16]. Another problem for binding prediction is the effect of base pair mismatches, or mutations, on Tm. For example, a 25-base oligomer with a single mutation is expected to have very similar hybridization thermodynamics compared to the perfectly complementary strand due to the rather small energetic differences. However, it has been observed that for these same sequences – as targets in a surface hybridization experiment – the single mutant strand does not achieve the same efficiency (100%) as the perfectly complementary strand does even after long hybridization times, potentially indicating the presence of a kinetic barrier to hybridization [17]. It has also been shown that the location of a mismatch along the strand plays a role in the efficiency with mismatches in the center having the greatest destabilizing effect on hybridization [18]. Levicky and Horgan [6] have recently outlined the challenges for microarray technology from a physical chemistry perspective and conclude that a more thorough understanding is needed for many aspects including the origin of differences between solution and surface-bound DNA hybridization equilibrium constants as
355
N.B. Tito, J.M. Stubbs / Chemical Physics Letters 485 (2010) 354–359
well as how mutation discrimination is affected by solution conditions and surface probe construction. Molecular simulation can provide valuable microscopic insight into these systems which can lead to better understanding of the differences between solution- and surface-based DNA hybridization, which could in turn ultimately lead to more robust predictive models. Although a few surface simulation studies have been carried out investigating melting of a single strand [19], optimum surface strand density [20], and molecular recognition [21], a systematic simulation study of hybridization and melting transitions at equilibrium has not yet been carried out, which is the primary focus of this work. 2. Methodology 2.1. Molecular models The coarse-grained molecular model is largely identical to that of Drukker et al. [22] which is a representation of a B-DNA helix. The principal modification is that a reduced hydrogen bonding interaction was necessary to reproduce estimated Tm values based on initial results, which showed that reduction of the hydrogen bonding interaction by 20% was necessary for agreement [23]. Additionally, non-complementary nucleotides are taken to not interact via hydrogen bonding. For surface simulations, a hard surface model was employed (i.e. any Monte Carlo move that caused a DNA interaction site to overlap with the surface was rejected). The surface-bound strand had its movement restricted by not allowing the tethered backbone or spacer bead to move from its position 5 Å above the surface but otherwise allowing all other moves, thus fixing it in place on the surface. For simulations with spacers between the surface and strand, each spacer unit had the same interaction parameters as a backbone unit and thus represents a unit consisting of phosphate and deoxyribose subunits, approximately the same size as the ð—CH2 CH2 CH2 —PO 4 —Þn units described in the work of Shchepinov et al. [8].
available for hydrogen bonding, additional surface-bound simulations were carried out without the second, target, strand i.e. only the probe strand is present. The same temperature range was investigated for all three immobilization scenarios of no spacer, a 5-unit spacer and a 10-unit spacer. Each simulation was equilibrated for at least 108 Monte Carlo (MC) cycles after which at least 3.2 108 MC cycles were carried out for production where one MC cycle consists of Nmolecule MC steps. The MC simulations used three different types of trial moves: translation of a molecule, rotation of a molecule, and conformational moves via coupled-decoupled configurational-bias Monte Carlo (CBMC) [26–28] and self-adapting fixed endpoint CBMC [29] to ensure adequate sampling of interior degrees of freedom. The type of MC move was selected at random using a fixed set of probabilities equally divided between translational, rotational and conformational moves. The maximum translational and rotational displacements were adjusted during the equilibration part in order to yield 50% acceptance rates for each type of molecule in the simulations. All error bars depict the standard error of the mean of block averages from a single simulation trajectory with blocks 1.6 107 MC cycles in size. 3. Results and discussion 3.1. Hydrogen bonding For this model a pair of complementary bases was determined to be hydrogen bonded if the energy of their interaction was less than 200 K 1.7 kJ mol1. This value was chosen to exclude those pairs that have some interaction but are significantly further apart then the more energetically stable hydrogen-bonded configurations. Note that the interactions have been parameterized to reproduce the experimental melting temperature; for comparison, changing the criterion to more strict values of 700 K for AT and 1300 K for GC pairs had no effect on the resulting Tm values but did not allow some helical structures to be counted as completely hydrogen bonded even at low temperatures.
2.2. Simulation details Three duplex systems were studied in detail as listed below: a fully complementary (perfect match or PM) decamer with moderate GC content (40%) for which an experimental Tm is known [24] as well as two possible single mutations from G to A of the PM sequence, with either a terminal mismatch (M1) or interior mismatch (M2). The mismatch is given in bold below:
PM : 50 CCAACTTCTT30 30 GGTTGAAGAA50
3.2. Solution results Fig. 1 shows the average number of hydrogen bonding base pairs, NHB, as a function of temperature for the PM, M1 and M2 sequences. Values of Tm are determined by estimating the temperature at which half the base pairs are hydrogen bonding, i.e. when
10
M1 : 50 CCAACTTCTT30 8
30 AGTTGAAGAA50 0
0
M2 : 5 CCAACTTCTT3
30 GGTTAAAGAA50
N HB
Simulations were carried out in the canonical (NVT) ensemble using periodic boundary conditions [25] at various temperatures between 285 K and 355 K and a concentration of 1.66 mM for each strand, corresponding to a cubic simulation cell 100 Å on each side. For surface-bound studies, the upper sequence (as given above) was immobilized at its 50 end, and the surface was taken to have a thickness of 14 Å. The z-dimension was extended by 14 Å to maintain a constant volume for the strands thus corresponding to a simulation cell with the dimensions 100 Å 100 Å 114 Å. Long range interactions were truncated at 14 Å so that no interactions occurred from one side of the surface to the other. To investigate the effect of the surface on the conformations of the bound strand that are
6
4
2
0 280
300
320
340
360
T (K) Fig. 1. Average number of hydrogen bonding base pairs as a function of temperature for solution simulations. Circles, diamonds and triangles correspond to PM, M1 and M2 sequences, respectively. Lines are drawn as guides for the eye.
356
N.B. Tito, J.M. Stubbs / Chemical Physics Letters 485 (2010) 354–359
the curves cross NHB = 5. As the range of temperature values was covered in 5 or 10 K increments, the uncertainty in the simulation Tm values is likely approximately 5 K. Estimated Tm values were determined using the unified parameters of the empirical nearest-neighbor model [13] for both the simulated and experimental (single strand) concentrations of 1.66 mM and 1 lM, respectively, and Table 1 summarizes this data along with the experimental value. As shown in Table 1, the empirical model value matches the experimental value to within 2 K, and the simulation values match the empirical model values to within 4 K for all 3 sequences. Thus, these solution-based results are validated directly against a robust predictive model and indirectly against experiment. Fig. 2 shows the fraction of configurations where a nucleotide is hydrogen bonded, fHB, as a function of its location in the strand where the nucleotides are numbered starting from the 30 end of the lower sequence as given previously, and as a function of temperature. For clarity, error bars are omitted. Note that the Tm can be seen in this figure as the sharp drop-off between 325 and 330 K. As can be seen, the two ends have on average noticeably fewer hydrogen bonds even at low temperatures, which corresponds to end fraying of the duplex. It is more pronounced for AT pairs which is expected since their interaction is not as strong as that of GC pairs. That this behavior is present in the simulation results despite having the same interaction parameters as interior AT or GC pairs indicates that an important factor in maintaining nucleotide hydrogen bonding is for a nucleotide to have its available conformations restricted, i.e. be held in place by hydrogen bonding neighbors. Because the ends lack a neighbor on one side, they are somewhat destabilized. Another aspect of this result is that this model implicitly captures cooperative effects due to hydrogen bonding neighbors, a key idea of the empirical nearest neighbor model [12,13]. This same effect can also be seen for the neighbors of the mismatched bases, where the effect of being adjacent to a mismatch lowers fHB by between approximately 5–10% (Fig. 3a and b) at 285 K as compared to the PM sequence (Fig. 2). The interior mutation (M2) is shown to be more destabilizing than the terminal mutation (M1), as is expected. These results show that the destabilizing effect of a mutation on neighbors can be reproduced purely through a mutant’s inability to hold its neighbors in place, even when the neighbors have the same interaction potential, a situation analogous to terminal bases with no neighbor on one side. 3.3. Surface-bound results Graphs of NHB as a function of temperature for all three sequences where the 50 end of the upper strand (as given previously)
Table 1 Summary of Tm where [expt] and [sim] denote the experimental and simulation concentrations, respectively, and Expt, NN and Sim denote experimental, empirical nearest-neighbor model and simulation results. Sequence
Method
System
Tm (K)
PM PM PM M1 M2 PM M1 M2 PM M1 M2 PM M1 M2
Expt NN NN NN NN Sim Sim Sim Sim Sim Sim Sim Sim Sim
Solution [Expt] [Sim] [Sim] [Sim] Solution Solution Solution 5-Spacer 5-Spacer 5-Spacer 10-Spacer 10-Spacer 10-Spacer
305.3 303.5 323.3 318.9 300.3 327 319 299 331 326 305 328 326 311
fHB
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
280 290 300 310 T (K) 320 330 340 350 360
1
2
3
4
5
6
7
8
9
10
Base Index
Fig. 2. Fraction of configurations in which a particular nucleotide is hydrogen bonded, with the base index starting at the 50 end of the upper strand, as a function of temperature for the upper strand in the PM solution results.
is immobilized directly above a hard surface are given in Fig. 4, and an analogous plot corresponding to the 5-unit spacer between the immobilized strand and the surface is given in Fig. 5. Results for the 10-unit spacer are similar to those of the 5-unit spacer and are not shown. Fig. 4 indicates that without a spacer, none of the three sequences forms more than 4–6 hydrogen-bonded pairs even at low temperatures. Plots of fHB (not shown) show that without a spacer the distal half of the strand hydrogen bonds normally (as in the solution case) but that the half adjacent to the surface has negligible hydrogen bonding. On the other hand, Fig. 5 indicates that with a spacer of at least 5 units long the corresponding Tm values are almost identical to the non-immobilized solution-based results (Table 1). To investigate the effect of the surface on the probe strand’s ability to form a complete duplex with the target, an analysis was performed on the probe-only simulations to find the percentage of probe conformations that disallowed complete hybridization to a target strand due to the surface. The analysis was carried out by taking each conformation of the probe and determining the hypothetical location of the target strand that would have a position consistent with optimum hydrogen bonding to the probe. If this hypothetical target position overlapped with the surface, then this probe conformation was counted as non-viable with respect to duplex formation. Averaging over all conformations and temperatures, those directly adjacent to the surface (no spacer) were non-viable in 27.6% of all configurations, whereas the 5- and 10-unit spacers were only non-viable in 8.5% and 4.8% of all configurations, respectively. These values show no temperature dependence over the range investigated here, 285–340 K. This analysis indicates that the destabilizing effect of the surface is due primarily to the decrease in number of probe conformations that allow the target strand to completely hybridize with it, i.e. the presence of the surface excludes a set of probe conformations from forming a duplex. In contrast, distancing the probe from the surface by a spacer reduces the number of forbidden conformations considerably, with little difference between the two spacers investigated here. Another interesting aspect is that the spacer simulations show that the PM and M1 sequences have slightly more similar melting curves than their solution results, Figs. 1 and 5. This could be due to the surface having a destabilizing effect on an immobilized strand’s nucleotides closest to the surface, and is more clearly shown in the fHB versus nucleotide plots given in Fig. 6a and b compared to Figs. 2 and 3a. Thus, the destabilizing influence of a
357
N.B. Tito, J.M. Stubbs / Chemical Physics Letters 485 (2010) 354–359
A
B
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
fHB
fHB
280 290 300 310 T (K) 320 330 340 350 360
1
2
3
4
5
6
7
8
9
10
Base Index
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
280 290 300 310 T (K) 320 330 340 350 360
1
2
3
4
5
6
7
8
9
10
Base Index
Fig. 3. Fraction of configurations in which a particular nucleotide is hydrogen bonded, with the base index starting at the 50 end of the upper strand, as a function of temperature for the upper strand in the M1 (left, A) and M2 (right, B) solution results.
For simulations of surface-bound probes with either 5- or 10unit spacers, the amount of end fraying present (from fHB plots) is approximately the same as in the solution simulations. Consistent with experimental data [8], the overall effect of the spacers is that while complete lack of one prevents full duplex formation even at low temperatures, there is relatively little difference between 5-units and 10-units, and they both give similar melting temperatures to the solution results.
10
8
N HB
6
4
3.4. Thermodynamic quantities 2
0 280
For the association between two single strands A and B into a double-stranded duplex AB the equilibrium constant can be written as 300
320
340
360
T (K)
Kc ¼
Fig. 4. Average number of hydrogen bonding base pairs as a function of temperature for surface-bound probes. Styles as in Fig. 1.
10
ðcAB =c0 Þ ðcA =c0 ÞðcB =c0 Þ
where cAB, cA and cB are the concentrations of the duplex AB, strand A and strand B, respectively, and c0 is the standard state concentration and is taken to be 1 mol L1. This can then be combined with the relationships
DA0 ¼ RT ln K c 8
ð1Þ
ð2Þ
and
DA0 ¼ DU 0 T DS0
N HB
6
ð3Þ
to yield 4
ln K c ¼ 2
0 280
300
320
340
360
T (K) Fig. 5. Average number of hydrogen bonding base pairs as a function of temperature for surface-bound probes with a 5-unit spacer. Styles as in Fig. 1.
mismatch is decreased if that nucleotide is less able to hydrogen bond in the first place due to surface proximity. However, the effect is not large and could also be due to the uncertainties in the melting curve data.
DU 0 1 DS0 þ R T R
ð4Þ
where the constant volume conditions imply the use of the Helmholtz (A) and internal (U) energies. Thus, under the assumption that DU0 and DS0 are temperature independent, ln Kc will vary linearly with 1/T and will have a coefficient of DU0/R and constant of DS0/R. Carrying out this analysis for PM duplexes in both solution and surface simulations results in the values listed in Table 2 with error estimates based on the standard error of the linear regression analysis. As can be seen the interaction energy is in general substantially decreased for all surface-bound systems with the no spacer case the most affected at approximately half the solution interaction value. There is little difference between the 5- and 10-unit spacer results which recover much of the solution interaction. The destabilizing effect of the surface is thus most severe when there
358
N.B. Tito, J.M. Stubbs / Chemical Physics Letters 485 (2010) 354–359
A fHB
B
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
280 290 300 310 T (K) 320 330 340 350 360
fHB
1
2
3
4
5
6
7
8
10
9
Base Index
1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0
280 290 300 310 T (K) 320 330 340 350 360
1
2
3
4
5
6
7
8
9
10
Base Index
Fig. 6. Fraction of configurations in which a particular nucleotide is hydrogen bonded, with the base index starting at the 50 end of the upper strand, as a function of temperature for the upper strand in the PM (left, A) and M1 (right, B) 5-unit spacer results.
Table 2 Simulation results for thermodynamic quantities of DNA hybridization for the PM sequence. System
DU0 (kJ mol1)
DS0 (J K1 mol1)
Solution No Spacer 5-Spacer 10-Spacer
217 ± 19 100 ± 11 160 ± 11 165 ± 11
559 ± 61 264 ± 34 359 ± 33 372 ± 36
is no spacer and partially ameliorated with the presence of at least a 5-unit spacer. The effect on the entropy of association, DS0, is similar with the largest value for the solution simulations, the smallest value for the no spacer surface-bound system and the 5- and 10-unit spacer systems at intermediate values. This is consistent with expectations based on the amount of constraint placed upon the probe strand where direct immobilization on a surface results in the largest change in conformational, rotational and translational entropy relative to the unbound system and the presence of a spacer partially mitigating these changes. As DS0 is less negative for the surface-bound systems, the entropic penalty to form a duplex if one strand is bound to a surface is less than the penalty when both single strands are in solution. Since for the association of A and B
DS0 ¼ S0AB ðS0A þ S0B Þ
4. Summary The application of a coarse-grained DNA model to three DNA decamer systems is shown to give a consistent Tm value with respect to experiment. Additionally, the model rationalizes experimental phenomena in terms of molecular-level events. Simulation results allow for the microscopic interpretation of phenomena such as nucleotide neighbor destabilization due to mutations through increased conformational freedom. The destabilizing nature of the surface in this context corresponds well to the decrease in conformations of the surface-tethered probe strand that allow for complete duplex formation without overlap with the surface. Finally, both the energy and entropy of duplex formation were found to decrease in magnitude when the probe strand is surfacebound. Acknowledgments Financial support from a UNE Presidential Faculty Scholarship Mini-Grant and Grinnell College Faculty Scholarship Grants is gratefully acknowledged. We thank J.I. Siepmann for helpful discussions and suggestions and the Minnesota Supercomputing Institute where the calculations were partially carried out.
ð5Þ
where strands A and B are the probe and target, respectively, if S0B is assumed to be constant regardless of whether strand A is surfacebound or not, comparing two associations gives
DDS0 ¼ DS0AB DS0A
for the duplex when surface-bound and unable to completely associate along its entire length due to the presence of the surface.
ð6Þ
where
DDS0 ¼ DS0surface-bound DS0solution
ð7Þ
DS0AB ¼ S0AB;surface-bound S0AB;solution
ð8Þ
DS0A ¼ S0A;surface-bound S0A;solution
ð9Þ
Since DDS0 > 0 for all surface-bound associations (Table 2), based on these results the effect of surface binding on the probe strand’s entropy, S0A , is much greater than the effect of surface binding (via one strand) on the entropy of the duplex, S0AB . This is no doubt in part due to the increased number of conformational states possible
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
E. Southern, K. Mir, M. Shchepinov, Nat. Genet. (Suppl. 2) (1999) 5. P.O. Brown, D. Botstein, Nat. Genet. 21 (1999) 33. R.J. Lipshutz, S.P.A. Fodor, T.R. Gingeras, D.J. Lockhart, Nat. Genet. 21 (1999) 20. D.J. Lockhart, E.A. Winzeler, Nature 405 (2000) 827. A. Pozhitkov et al., Nucl. Acids Res. 34 (2006) e66. R. Levicky, A. Horgan, Trends Biotechnol. 23 (2005) 143. W. Saenger, Principles of Nucleic Acid Structure, Springer-Verlag, New York, 1984. p. 143. M.S. Shchepinov, S.C. Case-Green, E.M. Southern, Nucl. Acids Res. 25 (1997) 1155. E.L.S. Wong, E. Chow, J.J. Gooding, Langmuir 21 (2005) 6957. A.W. Peterson, R.J. Heaton, R.M. Georgiadis, Nucl. Acids Res. 29 (2001) 5163. M.L. Del Giallo, F. Lucarelli, E. Cosulich, E. Pistarino, B. Santamaria, G. Marrazza, M. Mascini, Anal. Chem. 77 (2005) 6324. J. SantaLucia Jr., Proc. Natl. Acad. Sci. USA 95 (1998) 1460. J. SantaLucia Jr., D. Hicks, Annu. Rev. Biophys. Biomol. Struct. 33 (2004) 415. R. Owczarzy, P.M. Vallone, F.J. Gallo, T.M. Paner, M.J. Lane, A.S. Benight, Biopolymers (Nucl. Acid Sci.) 44 (1997) 217. L. Zhang, M.F. Miles, K.D. Aldape, Nat. Biotechnol. 21 (2003) 818. L. Zhang, C. Wu, R. Carta, H. Zhao, Nucl. Acids Res. 35 (2007) e18. A.W. Peterson, L.K. Wolf, R.M. Georgiadis, J. Am. Chem. Soc. 124 (2002) 14601.
N.B. Tito, J.M. Stubbs / Chemical Physics Letters 485 (2010) 354–359 [18] L.M. Wick, J.M. Rouillard, T.S. Whittam, E. Gulari, J.M. Tiedje, S.A. Hashsham, Nucl. Acids Res. 34 (2006) e26. [19] K. Qamhieh, K.-Y. Wong, G.C. Lynch, B.M. Pettitt, Int. J. Numer. Anal. Model. 6 (2009) 474. [20] L. Yao, J. Sullivan, J. Hower, Y. He, S. Jiang, J. Chem. Phys. 127 (2007) 195101. [21] A. Jayaraman, C.K. Hall, J. Grenzer, Biophys. J. 91 (2006) 2227. [22] K. Drukker, G. Wu, G.C. Schatz, J. Chem. Phys. 114 (2001) 579. [23] M. Buchanan, Summer Chemistry Research Program Final Report, Grinnell College, Grinnell, IA, 2006.
359
[24] R. Owczarzy, Y. You, G. Moreira, J.A. Manthey, L. Huand, M.A. Behlke, J.A. Walder, Biochem 43 (2004) 3537. [25] M.P. Allen, D.J. Tildesley, Computer Simulations of Liquids, Oxford University Press, Oxford, 1987. [26] J.I. Siepmann, Mol. Phys. 70 (1990) 1145. [27] J.I. Siepmann, D. Frenkel, Mol. Phys. 75 (1992) 59. [28] M.G. Martin, J.I. Siepmann, J. Phys. Chem. B 103 (1999) 4508. [29] C.D. Wick, J.I. Siepmann, Macromolecules 33 (2000) 7207.