The H$_{i}$ is the elastic energy of the single strand i (=1 or 2), which, in a discrete represen- tation, is given by where r ( i ) n is the three-dimensional position vector of n th bead in strand i ( n = 1 , 2 , · · · , N ) and b ( ∼ = 0 . 34 nm) is an average distance between neighboring beads within the ss. The first term accounts for bending energy with the bending modulus κb 3 for the ss, which is the related to its persistence length ( L$_{p}$ ) via κb 3 = L$_{p}$k$_{B}$T [8]. The second term, the stretching energy of each strand, is introduced to impose chain inextensibility condition. An appropriate value of k is numerically found by matching simulation of either force- extension curves or the mean end-to-end distances of ssDNA resulting from this model, with corresponding theoretical expression of the inextensible WLC without the stretching term. The V$_{12}$ is the pairing energy between complimentary bases. To describe the bp openings due to thermal excitation, namely, thermal bubbles and forks, the interaction is represented by V$_{12}$ = ∑ N n =1 V$_{n}$ ( r$_{n}$ ) where is the Morse potential [9, 10], r$_{n}$ = | r (1) n − r (2) n | is the distance between n th bp, r$_{0}$ and a are the bond distance and range which correspond to mean and fluctuation of DNA diameter respectively, and D is the potential depth. Whenever bps are unbound, the duplex becomes no more than two single strands with net persistence length as short as twice of L$_{ss}$ = 1 ∼ 4 nm [11], while for bound bps it takes L$_{ds}$ = 50 nm, the persistence length of long DNA, due to stacking interaction in the double strand. In order to include the stacking and destacking that cause such variation of the persistence length in the model, we consider κ ( n ) = L$_{p}$ ( n ) k$_{B}$T/ 2 b $^{3}$, and propose the persistence length per single strand L$_{p}$ ( n ) takes the form Here ˜ θ ( r ) = [1 − erf( r/c )] / 2 is a steplike form function which smoothly increases from 0 to 1 over the width c and r$_{1}$$_{/}$$_{2}$ is the distance at which L$_{p}$ is ( L$_{ds}$/ 2 + L$_{ss}$ ) / 2 (see Fig. 1). L$_{p}$ depends on three consecutive bp distances naturally as it represents the cooperativity of the stacking interaction [12]. The variation of L$_{p}$ depending on the bp distance r is 3 L$_{p}$ ( n ) = L$_{ds}$/ 2 − ( L$_{ds}$/ 2 − L$_{ss}$ ) ˜ θ ( r$_{n}$$_{−}$$_{1}$ − r$_{1}$$_{/}$$_{2}$ ) ˜ θ ( r$_{n}$ − r$_{1}$$_{/}$$_{2}$ ) ˜ θ ( r$_{n}$$_{+1}$ − r$_{1}$$_{/}$$_{2}$ ) . (4) V$_{n}$ ( r$_{n}$ ) = De − ( r$_{n}$ − r 0 ) $^{/a}$[ e − ( r$_{n}$ − r 0 ) /a − 2] (3) H$_{i}$ = N − 1 ∑ n =2 κ 2 ( r ( i ) n − 1 − 2 r ( i ) n + r ( i ) n $_{+1}$) 2 +$^{N}$$^{−}$$^{1}$ ∑ n =1 k 2 ( | r ( i ) n +1 − r ( i ) n | − b ) $^{2}$, (2)