Structural Characterization of Bubbles Formed in DNA Melting: A Monte Carlo Simulation Study.

Bubbles in DNA are involved in many important biological processes. In this work, a coarse-grained model is used for characterizing bubbles formed in DNA melting. The model resorts only to electrostatic interactions at the Debye-Hückel level, in combination with a short-ranged attractive interaction within a base pair. In spite of also omitting atomistic details, the model is able to capture experimentally established trends in persistence length and radius of gyration. By applying it on different systems, it is possible to conclude that there is a minimum size for stable bubbles, in the interval between 12 and 20 bp, which agrees well with previously published experimental findings. Simulated scattering data distinguishes between different bubbles and detects conformational changes in the melting process. Therefore, SAXS is regarded as useful for bubble formation studies, while simulations provide a molecular understanding.

Introduction

Accessing the genetic code is fundamental to important biological processes such as transcription and replication. It requires local openings of the DNA helix such that bubbles are formed. In cells, proteins are involved in this process, although DNA properties like the stability of the different conformations and the properties of the bubbles also play a substantial role. This process can, for example, be studied by thermal denaturation of DNA (melting).

DNA melting is an entropy-driven process, where there is a transition from a rigid double-stranded (ds) state to a much more flexible single-stranded (ss) state, i.e., the entropy gain overrules the attractive hydrogen bonds. The adenine–thymine (A–T) pair and the guanine–cytocine (G–C) pair involve two and three hydrogen bonds, respectively. Hence, the melting and the formation of bubbles depend on the sequence, and the process generally starts in A–T-rich regions.

There are several factors that influence the melting process, for example, pH,[1] ionic strength,[2,3] DNA concentration,[4,5] and DNA characteristics such as composition,[2,3] base sequence,[6] and length.[5,7] Experimentally, DNA melting can be studied by UV spectroscopy because disruption of base stacking leads to an increase in absorption. Another commonly used technique is calorimetry because it provides the enthalpy change associated with melting. Both of these methods measure the fraction of open base pairs, f. However, with the mentioned standard methods, no information regarding the fraction of completely separated molecules, p, is obtained. The latter can, for example, be accessed for special sequences using a quenching procedure,[8] whereas small-angle scattering of X-rays or neutrons can be used for the determination of p for any sequence.[9] The combined knowledge of f and p provides information about the intermediate states in melting, such as the average relative bubble length.[10,11] Our results presented in this article will be compared to the experimental findings of Zeng et al.[10,11]

There are many available theoretical models for studying DNA properties and melting, varying in their level of detail; from statistical Poland and Sheraga models[12−14] and thermodynamic nearest-neighbor models[15] to atomistic models.[16−18] Due to the time required for atomistic simulations, melting studies are restricted to short oligonucleotides.[19] For this reason, there is a large number of coarse-grained models, such as a two-site model by Drukker et al.,[20] the three-sites-per-nucleotide model (3SPN),[21−25] and oxDNA, where each rigid nucleotide is composed of four interaction sites.[26−28] The mesoscopic Peyrard–Bishop–Dauxios model[29,30] has reproduced experimental observations of bubble stability.[31]

In this work, a coarse-grained model for characterizing the bubbles formed in DNA melting is presented, with the purpose of understanding bubble formation on the molecular level. Here, electrostatic interactions on the Debye–Hückel level are used in combination with a short-ranged attractive interaction within a base pair. Even though it is a gross model and omits atomistic details, it gives new physical insight into bubble formation in DNA melting and the factors influencing it.

Theoretical Methods

Model

The structural properties of one dsDNA molecule with monovalent counterions have been studied using a simple model, in which electrostatic interactions at the Debye–Hückel level are used in combination with a short-ranged attractive interaction within a base pair. Positively charged hard spheres represent the counterions, whereas the solvent enters the system Hamiltonian only through its permittivity. Additional salt is treated implicitly. The molecule is constructed by two strands, where each strand is described as a chain of hard negatively charged spheres, connected by harmonic bonds. To represent the hydrogen bonds, which keep the strands together, a short-ranged attractive interaction was applied between complementary nucleotides. Each strand contains four bead types: two that correspond to the A–T pair and two that correspond to the G–C pair. The two bead types representing a pair are not distinguishable; hence, no difference is made between A and T nucleotides, and analogously, there is no difference between G and C nucleotides. This implies that within this model a GCGG sequence is identical to a CGGC sequence. Due to modeling purposes, a sequence of only one pair type is constructed by alternating the two bead types; hence, other beads of the same type in close proximity are not expected to contribute to the total short-ranged attractive interaction. For a schematic description of the model, see Figure . The radius of the sphere corresponding to a nucleotide was set to 2 Å to give a realistic representation of the electrostatic interactions. Notice that this value is not expected to correspond to the actual excluded volume effect of a nucleotide. All particles were enclosed in a spherical cell with hard walls in an attempt to mimic the cell environment.

Figure 1

Schematic description of the DNA model. (a) Grouping of the atoms for an adenine nucleotide. (b) Model of an oligonucleotide with the following sequence: four G–C pairs (beads with the two darkest shades), three A–T pairs (the two brightest shades), and four G–C pairs. Each bead has the charge −e. The arrows describe the attractive interactions, which are 3/2 times stronger for G–C pairs than A–T pairs. The attraction is applied only between beads of the same shade. The attractions represented by dashed arrows, presented only for the A–T block, are negligible for all beads due to the steep decay of the potential.

All interactions are assumed to be pairwise additive, and the system Hamiltonian has five contributions, according toThe hard sphere potential, Uhs, is given bywhere the summation extends over all strand beads and ions. Here, r is the center-to-center distance between particles with indices i and j, and uhs represents the hard sphere potential between two particles, according towhere R is the radius of particle i. The electrostatic potential energy, Uel, is calculated aswhere the summation extends over all particles in the system. Z is the valency of particle i, e is the elementary charge, and κ = [2NAIe2/(ε0εrkBT)]1/2, where NA is the Avogadro constant and I is the ionic strength. ε0 refers to the vacuum permittivity, and εr refers to the relative permittivity for water (78.4 at 298 K), which is assumed in this study to be temperature-independent.

The bond energy, Ubond, and the angular energy, Uangle, presented below, apply only to the strandsIn eqs and 6, Nseg corresponds to the number of segments (beads) in the strand, kbond is the force constant = 0.4 N/m, and r is the center-to-center distance between two connected segments, having the equilibrium separation r0 = 5 Å. The angular force constant kang = 8.6 × 10–24 J/deg2, and α is the angle formed by the vectors r – r and r – r from three consecutive segments, with the equilibrium angle α0 = 180°.

The hydrogen bonds are modeled as a short-ranged attraction with a 1/r9 decaywhich is sufficient for our purposes. Here, ε determines the interaction strength. The energy at closest contact (4 Å) is 24.4 and 16.2 kJ/mol for each base pair, which corresponds to 9.8 and 6.5 kT, respectively. The attraction is applied only between beads of the same type.

To mimic and induce bubble formation in dsDNA, different triblock sequences have been investigated, with each sequence having G–C pairs in the ends of the molecule and a central block of A–T pairs (see Figure ).

Systems

The structural and conformational changes of the dsDNA upon melting, as well as the effect of the bubble size, have been studied by using a model consisting of 48 bp, confined to a spherical cell of radius 150 Å, including counterions. The cell radius was approximately 2.5 times larger than the radius of gyration in the associated state, and increasing the cell radius did not affect the radius of gyration. Five different middle block lengths, presented in Table , were studied, and the simulations were performed with 150 mM implicit salt to mimic physiological conditions. The B20 sequence, i.e., a dsDNA molecule consisting of a middle block of 20 base pairs and two end blocks of 14 base pairs each, was regarded as the reference system. Hence, this system was also simulated at 50 and 10 mM salt to investigate the effect of ionic strength.

Table 1

Systems Studieda

	G–C:A–T:G–C
system	block lengths (bp)	NAT/Ntot
B6	21:6:21	0.125
B12	18:12:18	0.250
B20	14:20:14	0.417
B28	10:28:10	0.583
B36	6:36:6	0.750

The strand length of the molecule is 48 bp in all systems. The number of base pairs in each block is given by the G–C:A–T:G–C block length, and NAT/Ntot is the number of A–T pairs divided by the total number of base pairs. The reference system is marked in bold.

Although the sequences presented in Table are not very common experimentally, the systems were chosen to enable a systematic study of the effect of block length. However, they bear resemblance to the sequences in the experimental study by Zeng et al.[10,11]

Simulation Details

The equilibrium properties of the model systems were obtained by employing Monte Carlo simulations in the canonical ensemble according to the Metropolis algorithm,[32] using the simulation package MOLSIM, version 4.0.8.[33] All particles were confined in a spherical cell with hard boundaries, and the simulations were initiated in a state of completely associated DNA strands. Hence, all nucleotides are involved in base pairs. After an equilibration run, typically 1 × 107 trial moves/particle, where the strands reached a typical configuration for the system, a production run was performed. The latter was divided into sub-batches sufficiently long for achieving convergence of all properties of interest, between 2 × 107 and 6.6 × 108 trial moves/particle. More specifically, the criteria of converged base pair properties required the observation of approximately symmetric distributions of base pairing over the whole strand, and the radius of gyration should display a distribution close to Gaussian and a stable average value.

The strands were subjected to four types of trial displacements: (i) translation of a single bead, (ii) slithering, (iii) pivot rotation, and (iv) translation of an entire chain. In the slithering move, a randomly selected end bead is moved to a random position within the bond length of the other end and thereafter the segments are shifted such that the sequence order is not changed. In the pivot rotation, the shorter end of the chain is rotated around a randomly selected bond. The rotational angle was uniformly sampled between −90 and 90°. The selection among the trial moves was made randomly using the weights 0.85, 0.05, 0.05, and 0.05, respectively. The values of the translation parameters were set to keep the acceptance rate around 30%.

Analyses

Melting curves are represented as the average number of base pairs as a function of temperature. A geometrical definition of base pairing has been applied such that two nucleotides of different strands were defined to be in contact and form a base pair if their center-to-center distance was within 6 Å. Different values of the contact distance have been evaluated, and it was evident that the melting curve was not very sensitive to the exact value of the cutoff distance. Additionally, the average contact probability has been analyzed. Furthermore, the average distance within native base pairs has also been analyzed.

The strand stiffness was determined by the persistence length lp. The calculation was based on the projection of angles between bond vectorswhere ⟨Rbb2⟩1/2 is the root-mean-squared bead–bead separation and α = π – θ1, where θ1 is the first directional angle, further described by Akinchina and Linse.[34]

The effect of melting on the structural and the conformational properties was studied by calculating the radius of gyration of each strandHere, Rg is the radius of gyration and ⟨···⟩ refers to an ensemble average, whereas N is the number of beads of the strand and rcom corresponds to the center of mass. The radius of gyration was also determined from simulated scattering data by employing the Guinier approximation.

Reported uncertainties of the simulated quantities are one standard deviation of the mean, estimated from the deviation among the means of the subdivisions of the total simulation, according toHere, ⟨x⟩s is the average of x from one subdivision, ⟨x⟩ is the average of x from the total simulation, and ns is the number of subdivisions.

Scattering profiles and pair distance distribution functions, P(r), were evaluated from the simulations. At temperatures where the strands were associated, the scattering curves were calculated by taking into account both strands. Hence, the P(r) includes all distances within and between the strands.

Results and Discussion

Characterization of Reference System

System B20, with a total length of 48 bp and a 20 bp block of A–T nucleotides in the middle, was used as the reference system. The melting curve in Figure shows that the process occurs in two steps: (i) opening of the weaker A–T pairs, which induces a stable bubble (see snapshots in Figure ), and (ii) opening of G–C pairs, which causes a full separation of the two strands. Such stepwise melting is often difficult to capture with UV absorption measurements,[35] but it has been observed for DNA molecules with a length of several thousands of base pairs.[36] However, such behavior has been found for shorter molecules[10,11] when calculating the average fractional length of the bubble from UV absorption experiments. In the simulations performed in this study, the end separation and the complete strand separation are negligible at temperatures up to the highest plateau temperature, proving that the base pair separations observed in this study originate from the middle block. Hence, up to that temperature, the fraction of open base pairs is directly comparable to the average fractional length of the bubble.

Figure 2

Melting curve for B20, depicted as the average number of base pairs, Nbp, versus temperature, T. Representative snapshots at 273, 343, and 410 K are included, with counterions excluded for clarity.

Figure shows that the persistence length (lp) decreases as a function of temperature. More drastic changes in stiffness are induced at the same temperatures as the melting, indicating that the strand separation indeed results in a decreased lp; hence, it is not solely a temperature effect. This originates in an energy penalty for bending when the strands are associated, making the strands more rigid than in the separated state. This can be ascribed to an increased electrostatic persistence component. However, in reality, the helical structure of dsDNA makes the difference in lp much larger, cf. lp ≈ 22 Å for ssDNA[37] and lp ≈ 500 Å for dsDNA.[38] In spite of not including helicity, the model still captures an almost 40% change in lp due to melting. Melting is also associated with a decrease in radius of gyration (Rg) of approximately 20% for the individual strands (Figure ). Regarding Rg, the decrease is more drastic in the first melting step than in the second.

Figure 3

Left axis: Persistence length, lp, and radius of gyration, Rg, both normalized with respect to the lowest temperature studied, as a function of temperature, T, for the sequence B20. The solid lines are guides for the eye. Right axis: Melting curve expressed as the probability of base pairing, Pbp, versus T. Uncertainty, expressed as one standard deviation of the mean, is also included.

Simulated scattering intensity functions, I(q), at three different temperatures for the B20 system are presented in Figure a. The lowest temperature corresponds to completely associated strands, whereas at the intermediate temperature, a bubble is present. At the highest temperature studied, the strands are completely separated; hence, the scattering curve corresponds to one ssDNA. The featureless scattering pattern exhibited at all temperatures is typical for flexible polymers that can be in many different conformations. The I(q) of the lowest temperature decays more rapidly at low q compared to the I(q) at higher temperatures, in agreement with a larger Rg. The Rg determined from the Guinier approximation is in good agreement with values calculated directly from simulations (eq ), although it is slightly smaller (6%).

Figure 4

(a) Simulated scattering intensities and (b) dimensionless Kratky representation for sequence B20 at temperatures corresponding to fully associated strands, a bubble state, and complete separation. (c) Pair distance distribution function at 273 and 343 K.

The unitless Kratky representation (Figure b) provides information on the overall conformational state. For a rigid rod, (qRg)2I(q) increases linearly with q, whereas for a Gaussian chain, it initially increases and then reaches a plateau value. dsDNA behaves almost like a rigid rod at low temperatures, whereas with increased temperature, the flexibility increases and the behavior becomes more like a Gaussian chain. Hence, the analysis captures the decrease in persistence length. The extra bend on the curve at the bubble temperature is due to the induced regions with different flexibility, as the connected end blocks are more rigid than the middle part. The bubble also gives rise to a peak in P(r), which is visible when comparing the curves corresponding to fully associated strands and a bubble in Figure c. The high, narrow peaks at short distances correspond to highly repeated distances between beads, such as the bonds between the nearest neighbors. The contraction of the strands corresponding to a smaller Rg at higher temperatures is also detectable in P(r), as longer pair distances become less probable with increased temperature.

It is generally established that the melting temperature increases with increased ionic strength,[2,3] mainly due to screening reducing the repulsive electrostatic interactions between the backbones of the strands. The melting curves for the reference system with 10, 50, and 150 mM monovalent salt follow this trend (see Figure ). Notice also that the temperature interval for melting is smaller when the ionic strength is lower. Since the model captures these experimentally established trends despite its simplicity, it shows the importance of the interplay between electrostatic repulsion and hydrogen bonding.

Figure 5

Melting curves for the B20 system with 10, 50, and 150 mM salt, depicted as the average number of base pairs, Nbp, versus temperature, T.

Effect of Middle Block Length on Bubble Properties

Five systems with different lengths of the A–T middle block have been studied to capture the melting behavior and the effect of bubble stability. For the two shortest middle blocks, i.e., 6 and 12 bp, no plateaus are detected in the melting curve, which indicates that the fraction of separated base pairs increases very slowly with temperature until complete separation starts (see Figure ). The plateau appears for the sequences with a middle block length ≥ 20 bp, and with increasing middle block length, the plateau decreases in size, which is in agreement with experimental studies.[11,39] This implies that when the A–T block fraction increases the melting curve approaches the behavior of homogeneous sequences, which melt in a single step. As expected,[2,3] the melting temperature shifts toward lower temperatures as the A–T content increases.

Figure 6

Melting curves for the five systems with different middle block lengths, shown as the fraction of open base pairs, fsep, versus temperature, T. The total length is 48 bp, and the length of the bubble forming A–T middle block is given by B.

Figure a displays the average probability of base pair separation as a function of the segment number. Here, it is shown that bubbles are always present at the plateau temperatures if the middle block length is between 20 and 28 bp, i.e., Psep = 1 for A–T pairs and Psep = 0 for G–C pairs. For the sequence with only six A–T pairs, the bubble is more often closed than opened because the probability of separation is less than 50%. Simultaneously, Psep in the ends has increased to approximately 20%, indicating that fraying might occur even before the bubble is completely established. The experimental procedure described in references (8), (10), and (11) is not able to distinguish between fraying and bubbles, which makes it evident that simulations can provide additional information. The Psep analysis suggests that bubbles of the size ≤ 12 bp are not stable because Psep < 1, which is in agreement with the conclusion in reference (11), i.e., that there is a minimum length for a stable bubble on the order of 20 bp. For the largest A–T block, 36 bp, Psep > 0.1 in the end blocks. This shows that complete separation occurs and hence the bubble is not stable at the plateau temperature. The G–C end clamping regions appears too short to prevent complete separation; thus, the bubble is not stable. It is therefore possible that the stable bubble size is dependent on the length of the G–C clamping regions. The average separation within native base pairs, presented in Figure b, shows that as the length of the middle block increases, the bubble size increases, since a larger A–T fraction allows for larger separation within base pairs. Data for system B36 is not shown due to the complete separation of the strands.

Figure 7

(a) Probability of base pair separation, Psep, versus segment number, N, for the systems B6 (bubble forming A–T block length B = 6), B12, B20, B28, and B36 at their plateau temperatures (398, 383, 343, 323, and 298 K, respectively). (b) Average base pair separation, rbp, versus segment number, N, for systems B6, B12, B20, and B28 at the same temperatures.

The I(q) for the four smallest bubbles (data not shown) as well as their radii of gyration determined from the Guinier approximation appear very similar. The unitless Kratky representation (Figure a) is in agreement with the earlier conclusion, i.e., that a larger fraction of open base pairs corresponds to a higher flexibility because the flexibility increases with increased bubble size. The small unstable bubble systems, i.e., B6 and B12, show the same behavior as fully associated dsDNA (cf. Figure ). The other bubbles exhibit a bend and less steep slope, corresponding to higher flexibility, as larger parts of the strands are not base paired.

Figure 8

(a) Dimensionless Kratky representation of simulated scattering intensities for the systems with different lengths of the bubble forming A–T middle block: B6, B12, B20, and B28 at their plateau temperatures (398, 383, 343, and 323 K, respectively). (b) Pair distance distribution function for the stable bubble systems B20 and B28 at their plateau temperatures.

The P(r) for the stable bubbles shows a shift in peak position toward larger pair distances when the bubble size increases (Figure b); hence, P(r) analysis can distinguish between bubbles of different sizes. An experimental study showed that SAXS can be used for obtaining the fraction of completely separated molecules.[9] Hence, in combination with another technique for obtaining the fraction of open base pairs, such as regular UV spectroscopy, the average bubble length can be extracted for sequences known to form bubbles, according to references (8), (10), and (11). This opens up for experimental studies of bubbles in sequences that do not form hairpins, which was a restriction of the procedure described in references (8), (10), and (11). Additionally, our simulated scattering data suggests that SAXS can also provide information on bubble size and flexibility changes.

Conclusions

We have presented a coarse-grained model for DNA melting, in which only a few relevant variables are imposed. We have shown that the model is able to qualitatively capture experimentally established trends in lp and Rg, as well as effects of ionic strength and composition. Applying the model on systems with different sequences, it provided information about the bubble formation in melting. It was concluded that there is a minimum size for stable bubbles somewhere between 12 and 20 bp, which agrees well with experimental observations.[11]

Furthermore, it is shown that simulated scattering data distinguishes between different bubble sizes and that it can detect overall conformational changes. As has been shown experimentally,[9]I(0) can provide information about the fraction of completely separated molecules, which also opens up for studies of bubble length. For these reasons, SAXS can be useful in studying bubbles in DNA. However, for obtaining more information about the bubble stability and its dimensions, simulations are important because they provide a molecular understanding.

Future studies of interest include, for example, separation at one end of the molecule as well as simulations of longer molecules to be able to compare with experimental SAXS data.