Large bubble drives melting in circular DNA

DNA serves as the repository of genetic information, encompassing the instructions essential for the development, survival, and reproduction of all living organisms, ranging from the simplest prokaryotes, such as bacteria [1, 2] and viruses [3, 4], to the more intricate eukaryotes, including humans [5]. Prokaryotic DNA takes the form of a continuous loop without loose ends, characterized by its circular structure. On the other hand, eukaryotes typically exhibit linear DNA with distinct ends. Additionally, recent research findings indicate that extrachromosomal circular DNAs (eccDNAs) are present independently of linear chromosomal DNA in eukaryotes as well [6, 7]. The stability of DNA depends on the a number of physical and chemical parameters, such as the pH value of the environment [8], salt concentration, temperature [9] as well as mechanical forces [10, 11]. The thermal denaturation transition, i.e. melting [12] of DNA, in particular, has been studied for the last 50 years, both in its biological specificity, as well as a statistical and polymer physics problem, due to the greater accessibility of heating experiments in vitro and the importance of DNA melting to the polymerase chain reaction [13]. The melting of linear DNA has been understood theoretically, largely as a sharp first-order phase transition [14]. The linear DNA melting is driven by the replication fork or the Y-fork, since the base-pairs at the ends are entropically favoured to open up. However, the thermal melting of circular DNA is comparatively less understood, although the importance of circular DNA in human pathology [15] is becoming increasingly clear.

In the absence of free ends in circular DNA, the melting behavior is controlled by the presence of bubbles. The majority of theoretical studies have investigated circular DNA within the framework of topological constraints, where the twisting and bending dynamics of DNA are inherently interconnected [16, 17, 18, 19, 20]. This is because the Calugareanu-Fuller-White theorem [21] places a strict constraint on the linking number difference of any circular DNA system. Any flexibility gained by one part of the system by opening up a few base-pairs will have to be compensated by the twist and writhe of other parts of the molecule [22], where the base-pairs are hard to open up. These constraints, relaxed in linear DNA, result in a distinct melting behavior – the transition in topologically constrained DNA, such as closed circular viral DNA, is far less abrupt compared to that in linear DNA [23]. Studies have demonstrated that the competition between melting and local supercoiling induces phase coexistence of denatured and intact phases at the single-molecule level, contributing to the broadening of the transition in circular DNA [24, 18].

The bubble-mediated melting theory shows that the nature of DNA melting is reliant on the non-extensive logarithmic contribution to the entropy, which is of the form $S=Ns_{0}-c\ln{N}$ for a bubble of length $N$ with $s_{0}$ as the bulk entropy per unit length and $c$ is the reunion exponent for bubble formation [25, 26, 27, 28]. In this scheme, a circular DNA under overtwist and supercoiling has been shown to undergo a continuous melting transition with the emergence of a macroscopic loop at $T_{c}$ [29, 30]. The bubble length grows monotonically with temperature $T\to T_{c}^{-}$ until a maximum of half the chain is denatured, as the topological constraints prevent complete melting of the DNA molecule, unlike the linear case. Experiments conducted on supercoiled molecules also indicate that only one denatured bubble occurs per molecule [17]. Moreover, no stable bubbles could be detected when supercoiling was absent due to their short lifetime [31, 17, 32], although these measurements were made away from melting due to experimental aims and constraints. In the absence of supercoiling, one would also expect complete denaturation of the molecule, similar to the prediction of a diverging bubble size at criticality in the original Poland-Scheraga (PS) model.

It is important to acknowledge the existence of circular DNA structures without helical intertwining or topological linkages between strands, such as the singly-nicked polyoma form II DNA [4, 23] and relaxed plasmid DNA [17], and multiple other studies have both proposed and found evidence for topologically unconstrained closed circular DNA [33, 34, 35]. These structures provide an opportunity to unravel the effects of two separate kinds of constraints that could be applied to a closed circular DNA: (a) the helical winding of the two strands and the consequent torsional stress and supercoiling, and (b) the topological constraint of closure of the two covalent bonded strands individually while in the bound state. In the absence of a constraint of the first kind, which may be relaxed either via nicking or topoisomerase, the thermal behaviour and sedimentation analysis of circular DNA was at first seen to be very similar to linear DNA [36, 23]. One would generally expect that singly-nicked circular DNA will undergo melting due to the fraying of its ends originating from the nicked site and results in melting similar to that observed in linear DNA. However, the specific impact of dual covalent bond closure alone in relaxed circular DNA melting, as compared to linear DNA, remains open for investigation.

To address this, we use coarse-grained Brownian dynamics simulations to study the melting transition in non-supercoiled circular and linear DNA. Section II describes our model and the details of our simulations. Section III describes the melting mechanism employed in circular DNA, and the transient involvement of bubbles in the melting process. We summarize our findings in Section IV. Some of the details are provided in the Supplementary Information in the form of figures.

We utilize a simplified, low-resolution model for homopolymeric DNA composed of two complementary strands [37]. The energy function governing the interactions of $N$ base-pair DNA in this model is $E_{\mathrm{DNA}}=E_{\mathrm{strand-I}}+E_{\mathrm{strand-II}}+E_{\mathrm{strand-I,strand-II}}$. Here, $E_{\mathrm{strand-I}}$ and $E_{\mathrm{strand-II}}$ represent the potential functions characterizing each individual DNA strand, given by

	$\displaystyle\left.\begin{array}[]{ll}E_{\mathrm{strand-I}}\\ E_{\mathrm{strand-II}}\end{array}\right\}$	$\displaystyle=$	$\displaystyle\sum_{i=1}^{N-1}K(d_{i,i+1}-d_{0})^{2}+k_{\mathrm{topo}}(d_{1,N}$		(3)
	$\displaystyle-d_{0})^{2}+$		$\displaystyle 4\sum_{\rm N-nat}\bigg{(}\frac{\sigma_{i,j}}{d_{i,j}}\bigg{)}^{12}$		(4)

The distance between beads $d_{i,j}$ is defined as $|\mathbf{r}_{i}-\mathbf{r}_{j}|$, where $\mathbf{r}_{i}$ and $\mathbf{r}_{j}$ are the position vectors of beads $i$ and $j$, respectively with $i,j\in[1,N]$. We use dimensionless distances with $\sigma_{i,j}=1$ and $d_{0}=1.12$. The energy parameters in the Hamiltonian are in units of $k_{B}T$ where $k_{B}$ is the Boltzmann constant (set to $1$) and $T$ is the temperature measured in reduced units. The harmonic spring with dimensionless spring constant $K=100$ connecting adjacent beads along the strand, is given by the first term on the rhs of Eq. 4. The second term in Eq. 4 with spring constant $k_{\mathrm{topo}}$ is responsible for connecting the first and last beads on each strand. The third term introduces a repulsive potential that prevents the overlap of non-native pairs of monomers in strand-I and strand-II [38, 39, 40]. The two complementary strands of DNA interact via the following potential

	$\displaystyle E_{\mathrm{strand-I,strand-II}}$	$\displaystyle=$	$\displaystyle\sum_{\rm N.C.}4{\epsilon_{p}}\bigg{[}\bigg{(}\frac{\sigma_{i,j}}{d_{i,j}}\bigg{)}^{12}-\bigg{(}\frac{\sigma_{i,j}}{d_{i,j}}\bigg{)}^{6}\bigg{]}$		(5)
			$\displaystyle+\sum_{\rm N-nat}4\bigg{(}\frac{\sigma_{i,j}}{d_{i,j}}\bigg{)}^{12}.$

The base pairing between strand-I and strand-II is considered using the first term on the rhs of Eq. 5 with $\epsilon_{p}=1$. The native base-pair contacts (same $i$ of both the chains [38, 37]) give rise to the ladder like structure of DNA. The second repulsive term of the potential energy in Eq. (5) prevents overlapping of non-native pairs of monomers of strand-I and strand-II [37, 38, 39, 40].

The dynamics of this system is governed by the Brownian equations, given by

$$\frac{d\mathbf{r_{i}}}{dt}=\frac{1}{\zeta}(\mathbf{F}_{c}+\Gamma)$$

(6)

Here, $\mathbf{F}_{c}=-\mathbf{\nabla}E_{\mathrm{DNA}}$ is the conservative force and $E_{\mathrm{DNA}}$ is the sum of $E_{\mathrm{strand-I}}$, $E_{\mathrm{strand-II}}$ and $E_{\mathrm{strand-I,strand-II}}$. $\zeta$ is the friction coefficient, here set to 50, and $\Gamma$ is the random force, a white noise with zero mean and correlation $\langle\Gamma_{i}(t)\Gamma_{j}(t^{\prime})\rangle=2\zeta k_{B}T\delta_{i,j}\delta(t-t^{\prime})$. These equations of motion are integrated via the Euler method with time step $\delta t=0.0015$ for $10^{8}$ iterations at every temperature. To investigate how topology affects the DNA melting process, chains were subjected to topological constraints by setting $k_{\mathrm{topo}}=K=100$, effectively creating circular DNA. The melting of linear DNA, which lacks topological constraints ($k_{\mathrm{topo}}=0$) has been studied extensively [41, 42, 43, 44, 45, 46, 47] and provides a platform for investigating the influence of topology on the melting of DNA.

In the absence of a free end, circular DNA can undergo melting through the formation of bubbles. To obtain a microscopic view of the bubbles involved in the melting of circular DNA devoid of supercoiling, we partitioned the entire chain into three equal sections and monitored sector-wise change in average base-pair distance, denoted as $r_{1}$, $r_{2}$ and $r_{3}$, relative to the initial bound state. The sector wise measurements enable the detailed characterization of large or small bubbles within the system. If one or two segments play a substantial role, with $r_{1}$ or $r_{2}$ assuming a significant magnitude, while $r_{3}$ takes on a comparatively smaller value, an anticipated structure involves a large bubble, resembling the configuration in Fig. 1a. Conversely, if each segment contributes roughly equally i.e. $r_{1}\approx r_{2}\approx r_{3}$, the resulting structure is multiple small bubbles distributed uniformly along the chain in Fig 1b. Hence, the melting of circular DNA can be characterized by two types of possible pathways : (i) Type-I, marked by the emergence of a large bubble, and (ii) Type-II, where small bubbles uniformly form across the entire chain. As temperature increases, these bubbles grow, which eventually leads to the melting of the whole circular DNA. It would be interesting to investigate whether the melting of circular DNA (without supercoiling) follows the Type-I pathway (distinguished by a large bubble) as proposed in earlier studies for supercoiled DNA [18, 17, 30] or if it employs the Type-II pathway, especially in the presence of small bubbles.

In this paper, we investigated the thermal melting of non-supercoiled circular and linear DNA of four different lengths using Brownian dynamics simulations across a range of temperatures. Our study revealed that covalently closed circular DNA exhibits a higher melting temperature than linear DNA, particularly at smaller length scales. In circular DNA, we identified a critical minimum length threshold for bubbles, which is crucial for their growth or maintenance. Bubbles shorter than this specific threshold tend to shrink rather than maintain stability or expand. This result has significant implications as shorter DNA chains can only accommodate small bubbles, rendering the latter susceptible to shrinkage. As a result, higher temperatures are required to facilitate the formation of larger bubbles within smaller DNA chains, ultimately leading to the chain’s melting. Unlike circular DNA, there is no such restriction on threshold length scale in linear DNA. Our findings also confirmed that the melting of linear DNA occurs through Y-forks originating from either end. However, there is a weak correlation between the Y-forks at opposite ends of linear DNA, meaning that only one side predominates at any one time. We also ruled out the possibility of melting of non-supercoiled circular DNA via the proliferation of bubbles distributed along the chain. We reaffirmed that the melting of circular DNA predominantly occurs through the formation of one large bubble due to the significantly lower probability of occurrence of multiple threshold-sized bubbles along the chain compared to smaller bubbles.

GM gratefully acknowledges the financial support from SERB India for a start-up grant with file Number SRG/2022/001771. SS acknowledges the financial support from Ashoka University and Ashoka’s High-Performance-Computing Cluster.