ARTICLE pubs.acs.org/JPCB
DNA Melting in Slit Pores: A Reaction Density Functional Theory Yu Liu,†,‡ Yazhuo Shang,† Honglai Liu,*,† Ying Hu,† and Jianwen Jiang*,‡ †
State Key Laboratory of Chemical Engineering and Department of Chemistry, East China University of Science and Technology, Shanghai 200237, China ‡ Department of Chemical and Biomolecular Engineering, National University of Singapore, 117576, Singapore
bS Supporting Information ABSTRACT: A reaction density functional theory (R-DFT) is developed for chemical reactions in confined space by integrating reaction thermodynamics and DFT for chain fluids. The theory is applied to investigate DNA melting in slit pores, with nucleotides represented by coarse-grained charged Lennard-Jones particles. Three types of slit pores are considered for DNA melting: repulsive pore, attractive pore, and under electric field. In repulsive pores, the melting temperature increases slightly with reducing pore width, and the increase magnitude is nearly the same for DNA of different chain lengths. The double-strand DNA (dsDNA) and singlestrand DNA (ssDNA) are located in the slit center, particularly for long DNA due to the effect of configuration entropy. In attractive pores, the melting temperature increases with increasing wall-fluid interaction. The DNA chains are preferentially adsorbed near the slit walls with a strong wall-fluid interaction. Under electric field, the melting temperature increases slightly and is more distinct for shorter DNA.
1. INTRODUCTION DNA melting is an important biological process, in which double-stranded DNA (dsDNA) is dissociated into singlestranded DNA (ssDNA). The hydrogen bonds of dsDNA that maintain the double-helix structure are broken during melting. The melting process of DNA is usually characterized by melting curve, indicating the degree of melting versus temperature. Melting temperature is defined when half of dsDNA is melted. In the past, a large number of experimental and theoretical studies have been reported on DNA melting. For example, semiempirical formulations were proposed on the basis of experimental data to predict melting temperature.1-4 Molecular simulations were performed with a variety of coarse-grained models such as the two-site,5 bead-pin,6 and three-sites models.7 In our recent study, a molecular thermodynamic model was developed to predict DNA melting in ionic and crowded solutions.8 In many practical applications, DNA exists in confined space and exhibits significantly different properties from bulk solution. There have been numerous studies focusing on the conformation and diffusion of DNA under confinement. Odijk proposed a scaling theory for long DNA in nanochannels.9 Bonthuis et al. studied the shape of DNA in slit channels and estimated the radius of gyration.10 Cifra et al. examined dsDNA elongation in nanochannels using Monte Carlo simulation as a function of DNA chain length and stiffness as well as channel dimension.11 Chen et al. studied DNA diffusion in a slit pore using both Brownian dynamics simulation and single-molecule experiment. r 2011 American Chemical Society
It was revealed that the normalized chain stretch and diffusivity of DNA largely depended on confined conditions.12 Kato et al. examined the conformational transition of a long DNA from coiled to folded state in a cell-sized space surrounded by phospholipid membrane and found the transition was affected by ionic concentration and droplet size.13 Despite these experimental and theoretical studies for DNA in confined space, DNA melting under confinement is scarcely investigated. The objective of this study is to develop a statistical-mechanics theory for DNA melting in slit pores. The theory is based on the latest advance in density functional theory (DFT), which has been recognized as a robust method to describe the thermodynamic and structural properties of nonuniform fluids. DFT was initially used for simple hard-sphere fluids14,15 and later for more complicated fluids such as hard-core attractive fluids,16-19 LJ fluids,20-22 electrolyte solutions,23-25 flexible chain fluids,18,19,26-28 semiflexible chain fluids,29,30 and rod-like fluids.31-33 DFT was also applied for lipid bilayers,34 meso- and microporous materials,21,35,36 biological helical polymers,37 DNA in a bacteriophage,38 and pgRNA in Hepatitis B virus.39 Here we present a reaction DFT (R-DFT) for DNA melting by integrating a molecular thermodynamic model for chemical reaction40 and DFT for chain fluids under confinement. In Section 2, the theoretical framework is descriReceived: September 3, 2010 Revised: January 7, 2011 Published: February 7, 2011 1848
dx.doi.org/10.1021/jp108415x | J. Phys. Chem. B 2011, 115, 1848–1855
The Journal of Physical Chemistry B
ARTICLE
particles-ions are modeled by LJ and electrostatic potentials ! σ 12 σ6ij Zi Zj e2 ij ð1Þ uij ðrÞ ¼ 4εij 12 - 6 þ r r 4πε0 εr r where e = 1.6022 10-19 C, ε0 = 8.8542 10-12 C2 N-1 m-2, Zi is the charge magnitude of segment i, and εij and σij are calculated from the combining rules εij = (εiεj)1/2 and σij = (σi þ σj)/2. The values of εi, σi, and Zi for dsDNA and ssDNA particles and Naþ ion are listed in Table 1. These molecular-based parameters are also used in our previous study.8 A general chemical reaction is considered in a system with K components, X1, X2, ..., XK K X ni Xi ¼ 0 ð2Þ i ¼1
Figure 1. DNA melting in a slit pore. Each pair of nucleotides in the dsDNA and each nucleotide in the ssDNA are, respectively, represented by two different types of coarse-grained charged particles.
Table 1. Lennard-Jones Potential Parameters and Charges of DNA Sites and Ions diameter σ (nm) ssDNA
0.43
dsDNA
0.542
Naþ
0.235
well depth ε/kB (K)
charge Z (e)
558.71
-1
1117.4 65.422
-2 þ1
bed. DNA melting is examined in Section 3 under three different confined conditions, namely, repulsive pore, attractive pore, and under electric field. The concluding remarks are summarized in Section 4.
2. THEORETICAL FRAMEWORK DNA melting dsDNA h 2 ssDNA is examined in a slit pore with width H, as illustrated in Figure 1. The solvent, water, is treated as a background continuum with dielectric constant εr = 78. The dsDNA and ssDNA are modeled as freely jointed tangent polyion chains. Each nucleotide in the ssDNA is represented by a charged LJ sphere. The sphere has a charge -e and a diameter σss = 0.43 nm, which is the distance between two nucleotides in the ssDNA.41 Drukker et al. proposed a simplified backbone-base model to study the dynamics of DNA denaturation.5 They represented each nucleotide by two sites, one was the backbone sugar plus phosphate group and the other was the base. The LJ parameters of the backbone and base sites were εbackbone/kB = 150 K, σbackbone = 0.24 nm; εbase/kB = 225 K, σA = σG = 0.3 nm, and σT = σC = 0.24 nm.5 The subscripts A, G, T, and C refer to adenine, guanine, thymine, and cytosine. In our model, a nucleotide in the ssDNA is considered as a single site. On the basis of geometrical analysis, the LJ well depth of the single site εss is approximated by εbackbone þ εbase þ (εbackbone εbase)1/2. A pair of nucleotides in the dsDNA is also modeled as a charged LJ sphere but with different potential parameters from those in the ssDNA. By assuming that a pair of nucleotides has the same packing fraction as the two constituent nucleotides, the diameter of the dsDNA sphere is σds3 = 2σss3. The well depth εds is the summation of the ssDNA, εds = 2εss. The charges on DNA are neutralized by Naþ counterions. The interactions of particles and
where ni is the stoichiometric coefficient. At chemical equilibrium, we have K X ni μi ¼ 0 ð3Þ i ¼1
where μi is the chemical potential of component i. In addition, thermodynamic equilibrium requires δΩ½FðRÞ ¼0 δFi ðR i Þ
ð4Þ
where Ω[G(R)] is the grand potential functional, G(R) is the density vector defined as G(R) = [F1(R1), F2(R2), ..., FK(RK)], and Fi(Ri) is density distribution of component i. For chain-like molecules, Ri = (r1, r2, ..., rmi) where ri is the position of the ith segment, and mi is the chain length. The grand potential can be expressed as K X Ω½FðRÞ ¼ F½FðRÞμi Ni i ¼1
¼ F½FðRÞ-
K Z X i ¼1
μi Fi ðR i Þ dR i
ð5Þ
where Ni is the molecule number of component i, dRi = dr1dr2...drmi, and F[G(R)] is the free energy functional and can be decomposed into F ¼ F 9o þ F id þ F ext þ F ex
ð6Þ
where F, Fid, Fext, and Fex are the standard, idea gas, external, and excess free energies, respectively. The standard free energy F 9o is given by K K Z X X Fi ðR i ÞFi9o dR i F 9o ¼ Ni Fi9o ¼ ð7Þ i ¼1
i ¼1
where Fi9o is the standard free energy of component i. The idea gas free energy Fid can be derived exactly βFid ½FðRÞ ¼
K Z X i ¼1
Fi ðR i Þfln½Fi ðR i ÞΛ3i -1-Vintra ðR i Þg dR i ð8Þ
where Λi is the de Broglie wavelength, β = 1/kBT, kB is the Boltzmann constant, T is temperature, and Vintra(Ri) denotes the bonding potential. Because DNA is modeled as a freely jointed tangent chain, Vintra(Ri) can be expressed as Vintra ðR i Þ ¼
m i -1 X
vj, j þ 1 ðjrj þ 1 -rj jÞ
ð9Þ
j¼ 1
1849
dx.doi.org/10.1021/jp108415x |J. Phys. Chem. B 2011, 115, 1848–1855
The Journal of Physical Chemistry B
ARTICLE
where vj,jþ1(r) is the bonding potential between segment j and j þ 1 and given by 1 exp½-βvj, j þ 1 ðrÞ ¼ δðr-σj, j þ 1 Þ 4πσ 2j, j þ 1 The external free energy F F
¼
ext
ext
Ks Z X i ¼1
Combining eqs 3 and 17, we have K X
ð10Þ
i ¼1
¼-
can be simply expressed
K X i¼ 1
Fs, i ðrÞViext ðrÞ dr
F ex ¼ F hc þ F attr þ F ele
ð12Þ
where Fhc, Fattr, and Fele represent the hard-core, attractive, and electrostatic contributions. To evaluate these contributions, the weighted density approximation (WDA) is used but with different weight functions Ks Z X ðiÞ Fs, i ðrÞfhs ½Fhc ðrÞ dr ð13Þ F hc ¼
¼-
i ¼1
F
ele
¼e
Z X Ks i ¼1
ðiÞ
Fs, i ðrÞfattr ½Fattr ðrÞ dr
ð14Þ
Z Zi Fi ðrÞΨðrÞ dr þ
Φele ½Fele ðrÞ dr
ð15Þ
K X
ð16Þ μi9o
where is the standard chemical potential. In common DFT without reaction involved, Fi(Ri) canRbe solved directly at a given μi or with normalization condition Fi(Ri) dRi = Ni at a given number of molecules Ni. In a reaction system, if μi is known, for example, on the basis of the bulk phase information, then eq 16 can also be solved as in a nonreaction system.44 For reaction occurring in a closed system, however, neither chemical potential μi nor the number of molecules Ni can be readily derived, and eq 16 cannot be directly solved. Nevertheless, eq 3 holds true in chemical equilibrium condition, and the integration of eq 16 yields R βμi ¼ βμi9o þ ln
R
"
Fi ðR i Þ dR i
exp -βViext ðR i Þ-βVintra ðR i Þ-
δβF δFi
ex
# dR i
R
"
Fi ðR i Þ dR i
# ex δβF exp -βViext ðR i Þ-βVintra ðR i ÞdR i δFi R
ðiÞ
N0 þ ni R
"
# ex δβF exp -βViext ðR i Þ-βVintra ðR i ÞdR i δFi ð18Þ
is the initial number of molecule i and R is the degree where of reaction. Equation 18 can also be converted into K Y K 9o ¼ ai-ni ð19Þ N(i) 0
i¼ 1
with
R
ai ¼
R
" exp
Fi ðR i Þ dR i
δβF ex -βViext ðR i Þ-βVintra ðR i ÞδFi
#
ð20Þ dR i
as the activity of molecule i. K9o = exp(Σki=1βniμi9o) is the equilibrium constant and can be evaluated on the basis of experimentally determined enthalpy and entropy changes during DNA melting.8 The equilibrium constant can be expanded as ln K 9o ¼ -
42
For hard-core contribution, Tarazona weight function is used. For attractive and electrostatic contributions, we use the functions similar to those of Peng and Yu21 and Pizio et. al.24,25 These functions have been shown to be reliable and efficient,21,24,25,36,42,43 and the details are provided in the Supporting Information. From eqs 4-11, the density distribution can be expressed as " # ex δβF Fi ðR i Þ ¼ exp βðμi -μi9o Þ-βViext ðR i Þ-βVintra ðR i ÞδFi
ni ln
i¼ 1
i ¼1
F attr ¼
R
ð11Þ
where Ks represents the types of segments and Fs,i(r) is the density of segment i. In this study, either ssDNA or dsDNA is modeled as a homopolymer, and thus we have Ks = K. The excess free energy Fex cannot be derived exactly. As usual, Fex is decomposed into different contributions
Ks Z X
ni ln
βni μi9o
ΔH 9o -TΔS 9o RT
ð21Þ
where ΔH9o and ΔS9o are assumed to depend on the number of AT and CG base pairs rather than the sequence of bases ( ΔH9o ¼ nAT ΔHATo9 þ nCGo9 ΔHCG ð22Þ ΔS9o ¼ nAT ΔSATo9 þ nCGo9 ΔSCG with ΔH9oAT = 34.5 kJ/mol, ΔH9oCG = 36.8 kJ/mol, and ΔS9oCG = ΔS9oCG = 102.7 J/mol/K. nAT and nCG represent the numbers of AT and CG base pairs in a dsDNA chain, respectively. The details o o o o 9 9 9 9 , ΔHCG , ΔSAT , and ΔSCG can be found of how to estimate ΔHAT 8 in our previous study. Using eqs 16-20, the density distribution and reaction degree can be solved iteratively. In numerical calculation, the segment density Rin addition to i δ(r - r molecular density can be calculated from Fs,i(r) = Σm j=1 j)Fi(Ri) dRi. In a slit pore, eq 16 can be further expressed as45,46 mi X i, j i, j Fs, i ðzÞ ¼ exp½βðμi -μi9o Þ exp½-βλi, j ðzÞGL ðzÞGR ðzÞ j¼ 1
with
¼
8 > < 1 > : 2σ
i, j GL ðzÞ
Z
1 minðH , z þ σ i Þ
i maxð0,z-σ i Þ
ð23Þ j¼1
i, j-1 exp½-βλi, j-1 ðz0 ÞGL dz0 2 e j e mi
ð24Þ
ð17Þ 1850
dx.doi.org/10.1021/jp108415x |J. Phys. Chem. B 2011, 115, 1848–1855
The Journal of Physical Chemistry B
¼
ARTICLE
i, j GR ðzÞ
8 > < 1 > : 2σ i
Z
j ¼ mi
1 minðH , z þ σi Þ
maxð0,z-σ i Þ
i, j þ 1 exp½-βλi, j þ 1 ðz0 ÞGR
dz0 1 e j e mi -1
δF ex δF ext λi, j ¼ þ δFs,ði, jÞ δFs,ði, jÞ
ð25Þ ð26Þ
where Fs,(i,j) denotes the segment density of the jth segment in molecule i. Considering Ks = K for a homopolymer, eqs 17 and 18 can be simplified as R Fs, i ðzÞ dz o βμi ¼ βμi9 þ ln R P mi ð27Þ i, j i, j exp½-βλi, j ðzÞGL ðzÞGR ðzÞ dz j ¼1
and K X i ¼1
βni μi9o ¼
-
K X i ¼1
R ni ln R P mi j ¼1
¼-
K X i ¼1
Fs, i ðzÞ dz
i, j i, j exp½-βλi, j ðzÞGL ðzÞGR ðzÞ dz ðiÞ
ni ln R P mi j ¼1
mi ðN0 þ ni RÞ ð28Þ i, j i, j exp½-βλi, j ðzÞGL ðzÞGR ðzÞ dz
The activity ai is expressed as ðiÞ
a i ¼ R mi P j¼ 1
mi ðN0 þ ni RÞ i, j i, j exp½-βλi, j ðzÞGL ðzÞGR ðzÞ dz
ð29Þ
The numerical calculations can be carried out by two separate methods. Method 1: Calculate the degree of reaction R at given a temperature T. (1) Choose a reasonable initial degree of reaction R(0), for example, from idea-gas approximation. (2) Solve eq 23Zfor F(k) s,i (z) with the normalization condition. ðkÞ
ðiÞ
Fs, i ðzÞ dz ¼ mi ðN0 þ ni RðkÞ Þ
ð30Þ
(3) Solve eq 28 for a new degree of reaction R(kþ1). (4) Stop if |R(kþ1) - R(k)| < given precision; otherwise, go to step (2). Method 2: Calculate temperature T at a given degree of reaction R. (1) Assume upper and lower limits of temperature Tmax and Tmin. (2) Stop if |Tmax - Tmin| < given precision and T = (Tmax þ Tmin)/2. (Z), F(k),min (Z), and F(k),mid (Z) in (3) Solve eq 23 for F(k),max s,i s,i s,i the normalization condition eq 30 at T = Tmax, Tmin, and (Tmax þ Tmin)/2, respectively. (4) Calculate function f(T) f ðTÞ ¼ ln K 9o ðiÞ K X mi ðN0 þ ni RÞ þ ni ln R P mi i, j i, j i ¼1 exp½-βλi, j ðzÞGL ðzÞGR ðzÞ dz j ¼1
ð31Þ
Figure 2. Melting curves for DNA (m = 10) in repulsive slit pores with different widths.
at T = Tmax, Tmin, and (Tmax þ Tmin)/2, respectively. If Tmax and Tmin are reasonably assumed, then f(Tmax) 3 f(Tmin) should be
