Monte Carlo simulation of melting transition on DNA nanocompartment

DNA nanocompartment is a typical DNA-based machine whose function is dependent of molecular collective effect. Fundamental properties of the device have been addressed via electrochemical analysis, fluorescent microscopy, and atomic force microscopy. Interesting and novel phenomena emerged during the switching of the device. We have found that DNAs in this system exhibit a much steep melting transition compared to ones in bulk solution or conventional DNA array. To achieve an understanding to this discrepancy, we introduced DNA-DNA interaction potential to the conventional Ising-like Zimm-Bragg theory and Peyrard-Bishop model of DNA melting. To avoid unrealistic numerical calculation caused by modification of the Peyrard-Bishop nonlinear Hamiltonian with the DNA-DNA interaction, we established coarse-gained Monte Carlo recursion relations by elucidation of five components of energy change during melting transition. The result suggests that DNA-DNA interaction potential accounts for the observed steep transition.


Introduction
Studies on the physical chemistry of DNA denaturation have been lasted for almost forty years [1][2][3]. In 1964, Lifson proposed that a phase transition exists in one-dimensional polymer structure. He introduced several pivotal concepts, like sequence partition function, sequence generating function, etc., and established a systematic method to calculate the partition function [1]. These allow us to derive important thermodynamic quantities of the system. In 1966, Poland and Scherage applied Lifson's method to conduct research on amino acid and nucleic acid chains. They built Poland-Scherage (PS) model for calculating the sequence partition function and discussing the behavior of polymers in melting transitions.
Another excellent progress would be the building of Peyrard-Bishop (PB) model [4,5] for DNA chains. In PB model, the Hamiltonian of a single DNA chain, which is constructed by phonon calculations, is given so that we can obtain the system properties through statistical physics method. The PB model has introduced mathematical formula of stacking energy, as well as the kinetic energy and potential energy of each base pair. By theoretical calculation, one can show the entropy-driven transition that leads DNA to shift from ordered state to disorder one [6,7].
However, all these works have not involved the DNA-DNA interactions because the subject investigated is DNAs in bulk solution, and the interaction between them has ever been neglected. The main idea of this paper is to inspect the influence of collective effect on the DNA melting process, primarily motivated by the experiment results of DNA nanocompartment [8,9]. Under the enlightenment of Poland-Scherage model and Zimm-Bragg model [10], we simplify Peyrard-Bishop model to meet a reasonable Monte Carlo simulation by the elucidation of five components of energy changes during melting transition. The result shows that the melting temperature and transition duration depend on whether we take into account the DNA-DNA interactions among columnar assemblies of DNA. Figure 1. (a) Schematic drawing of a dual-state ADNC. By adding or removing 'fuel' strands (shorter segment), the ADNC can be switched between on (right) and off (left) state. Inset, a typical sequence used to fabricate ADNC. The 'fuel' strands is a segment of human p53 gene containing one site of most frequent mutation. (b) Isotherms of Γ for fluorescein encaged in a closed ADNC (h eff = 15 bp). The isotherm fits well to the Langmuir model: x/Γ nc = (1/Γ nc,max ))x + (1/KΓ nc,max ), where x is the concentration of the reporter and K is the association constant per site, Γ nc the surface concentration of encaged molecules. (c) Melting curves using the encaged fluorescein molecules as indicators. Filled circles or squares are corresponding to perfect complementary strands, and hollow circles or squares to single-base mismatched strands. The unit of relative intensity of fluorescein is defined as the light intensity of 5 µl 10 nM fluorescein on a spot size with 5 mm diameter. Inset shows the length of the ssDNA skeleton of nanocompartment and the complementary type. The background noise is within 1 unit.

Experiment
Recently, we found that specially designed DNA array can form a molecular cage on surfaces [8,9]. This molecular cage is switchable due to allosteric transformation driven by the collective hybridization of DNA. We named it "active DNA nanocompartment (ADNC)". Typical DNA motif designed to fabricate ADNC comprises two contiguous elements (inset to figure 1A): a double-stranded DNA (dsDNA) whose array is responsible for a compact membrane (figure 1A, right), and a single-stranded DNA (ssDNA) serving as skeleton supporting the dsDNA membrane, which is terminated on its 5′ end by a surface linker such as an alkanethiol group that can be tethered to gold surface with a sulphur-gold bond [9] or an amino group that can be tethered to SiO 2 substrate with specific surface attachment chemistry [11]. Because the diameter of ssDNA is much smaller than that of dsDNA, a compartment with designable effective height (h eff , 5~50 nm, commensurate with the length of ssDNA skeleton) can form between the dsDNA membrane and substrate surface.
Since ADNC is reversibly switchable, it is able to encage molecules with suitable size. We name this phenomenon molecular encaging effect. Both electrochemical methods [12] and fluorescent microscopy are used to substantiate the molecular encaging effect and the reversibility of switching. Once the closed ADNC entraps some chemical reporters, the surface concentration (Γ nc ) of the encaged reporters can be determined by cyclic voltammetry or fluorescent microscopy. Figure 1B shows the isotherms of the molecular encaging effect for fluorescein (C 20 H 10 Na 2 O 5 ). Figure 1C presents the melting curves of ADNC. Using the encaged molecules as indicator greatly sharpens the melting profiles for the perfectly complementary targets, and flattens denaturation profiles for the strands with a wobble mismatch. The observation shows that single-base mismatched strands are incapable of closing ADNC on surfaces. The result is highly consistent to our observation by electrochemical analysis [12]. These observations bring up an intriguing question: why the melting curves exhibit so steep transition compared to the case of DNA in bulk solutions or on a loosely packed microarray? We try to address this question in this paper.
Worthy of mention is that the steepness of melting transition is useful when the ANDC is applied to DNA detection [8,9]. First, it greatly enhances the discrepancy of perfect targets and single mismatches. This provides much enhanced specificity in DNA recognition, 100:1~10 5 :1 of our system versus 2.7:1 of conventional system. Second, more sensitivity is obtained with optimally decreased ambiguity. Therefore, the clarification of the origin of the steep shape should help us to further extend the experience to related fields or generate new techniques.

Modeling
Taking into account the directional specificity of the hydrogen bonds, the Hamiltonian of a single DNA chain is obtained as following form according to PB model [4][5][6], where the y n is the component of the relative displacement of bases along the direction of hydrogen bond. The stacking energy The Morse potential ( ) n V y describes the potential for the hydrogen bonds (3) However, in this study, the Hamiltonian in equation (1) is not sufficient; it neglects the structure of close-packing of DNA in ADNC. In our system, one should take into account the interactions between the nearest neighboring molecules [13,14]. To model the interaction, one envisions the molecules as rigid cylinders, carrying helical and continuous line charges on their surfaces. Each DNA duplex carries the negative charge of phosphates plus a compensating positive charge from the adsorbed counterions. Let 0 < θ < 1 be the degree of charge compensation, f 1 , f 2 and f 3 the fractions of condensed counterions in the minor and major grooves (f 1 + f 2 + f 3 = 1). The mobile counterions in solution screen the Coulomb interactions between the two molecules, causing at large separations an exponential decay of the latter with the Debye screening length κ -1 . The solvent is accounted for by its 34 H ≈ Α (g = 2π/H), and hard-core radius a = 9 Å. We take the following form for the pair interaction potential [15][16][17][18]: where R (> 2a) is the distance between the two parallel DNA molecules, ∆z a vertical displacement, equivalent to a "spin angle" φ = g∆z. Here, u 0 = 8πσ 2 /εκ 2 (about 2.9k B T/Å at physiological ionic strength), and κ n = with the modified Bessel functions K n (x) and I j (y). The primes denote derivatives. The sum rapidly converges, and it can be truncated after |n| = 2. Since κ n R > 3 and g～ κ, each of the terms in the sum decreases exponentially at increasing R with the decay length κ n -1 ∝ 1/n. Figure 2 present a scheme of interaction between two neighboring columnar DNA molecules charged with counterions on its surface. The distance between two DNA columns in our simulation is about 30 Å and the helical pitch of DNA molecule is about 36 Å. For brevity, we take the mean-field approximation that the pair interactions mainly exist between charges in the same height.

Monte Carlo Simulation
Let  The simulation begins at t = 0, T = 0 ˚C. At each step, t increases by 1, and state of base pair at (x 0 , y 0 , i) is inverted, i.e., ψ (x 0 , y 0 , i, t+1) =1-ψ (x 0 , y 0 , i, t). We assume that by changing the state of the system for M N P Z × × × times, the system will approximate the equilibrium state infinitely. The change will be applied to each base pair for average Z times. Z is determined by experience and should be reasonable. We increase T by ∆T during the simulation. Therefore we have the Whether the state inversion is permitted depends on the energy change (∆E t ) in each step. The possibility of the state change at each step is If the current state change is permitted, we keep up changing the system state at t+1. If the state change is forbidden by the possibility, the system state remains unchanged at t and waits for another change at t+1.
To achieve a relatively precise simulation, the change of the total energy at time t+1 relative to that at the time t is analyzed by five components. The recursion relation of energy change in each step is written as: where E(t+1) and E(t) are the system energy for the instant t+1 and t respectively, and ( 1,2,3,4,5) l E l ∆ = is the variation of the lth component. The energy change depends on both the recursion relation of the base pair at (x 0 , y 0 , i) and the states of its nearest neighbors. The global energy variation is determined by the local states around the base pair (x 0 , y 0 , i). Following analysis presents the recursion relation of energy changes.

The hydrogen-binding energy (∆E 1 )
This component of the energy consists of the Morse potential (equation (3)) and kinetic energy along the orientation of the hydrogen bonds. The binding energy is independent of the states of its neighboring base pairs.
where J (J<0) is the binding energy for each base pair that is in '1' state, while the binding energy for the '0' state is zero to be reference..

The stacking energy (∆E 2 )
To simplify the calculation of the stacking energy shown in equation (2), we take into account the states of base pairs at (x 0 , y 0 , i-1, t) and (x 0 , y 0 , i+1, t). Their states remain unchanged during the interval from t to t+1. We employ the periodic boundary condition (PBC) listed below.
where w is the stacking energy stored in two nearest neighboring base pairs in '1' state.

Morse potential away from equilibrium point (∆E 3 ).
We set ∆E 1 = 0 for the uncoupled hydrogen bond at base pairs. However, for a '0' state is next near to a '1' state in the same DNA strand, the distance between two base pairs is so close that the Morse potential should be taken into account. We assigned energy E to every two nearest neighboring base pairs that are in different states in the same DNA.

The effect of excluded volume (∆E 4 )
The effect of excluded volume in the nature of DNA phase transition is discussed in Fisher's work [20]. The excluded volume effect is connected to the system entropy variation. The effect is prone to separate two complementary strands in a double helix. We use F to represent the energy change corresponding to this effect. One should notice The energy changes discussed above are summarized in the figure 4 below, which does not take into account the DNA-DNA interactions so far.

DNA-DNA interaction potential (∆E 5 )
We  (12) where G is the interaction energy between each pair of ions. Adding ∆E 5 to ∆, we will get the energy variation including the DNA-DNA interaction.

Results and Discussions
Following the equation (5) - (12), we could achieve a coarse-gained simulation of the melting curves of ADNC as well as that of DNA in bulk solutions. To perform the task, we choose suitable scale parameters to carry out the simulation: M = 100, N = 100, P = 20. The values of M and N chosen are much smaller than ones of the actual situation, which is up to 4 10 in the experiment. Since we take the periodic boundary condition, the values of M and N used do not change our result. The starting temperature is 0˚C, and the final temperature is 100˚C, with increment of 0.01˚C for each step.
To guarantee the system reaches equilibrium state, we take 6 10 state changes under a specific temperature. Each base pair has average 5 times to be changed. At each step, we count the number of DNAs that is still hybridized and calculate the percentage for dsDNA in ADNC. The simulation result shown in figure 5 shows a steep melting transition (hollow circles), consistent to the experimental observations. The simulated result without considering DNA-DNA interaction show in filled circles in figure 5 also agrees with the DNA melting curves in bulk solution. Comparison between the two cases suggests that the DNA-DNA interaction greatly increases the melting point of dsDNA chains.
In conclusion, we have established a simple coarse-gained model to simulate the melting transition of DNA in ADNC. The result provides a reasonable explanation for our experimental observations. Although the simulation method discretizes the Morse potential and stacking energy proposed in Peyrard-Bishop model, the result still present a comparable approximation to experimental data due to our fine treatment of energy changes during melting transition. However, this work is only the beginning of insightful theoretical investigation for the rationality of ADNC. In future work, we will establish a more precise model to employ an extensive investigation of phase transition occurring in ADNC as well as its derived DNA machines.