Thermal conductivity of B-DNA
aa r X i v : . [ phy s i c s . c o m p - ph ] S e p D R A F T Thermal conductivity of B-DNA
Vignesh Mahalingam a,1 and Dineshkumar Harursampath a,2 a Department of Aerospace Engineering, Indian Institute of Science, Bangalore, 560012, IndiaThis manuscript was compiled on October 1, 2020
The thermal conductivity of B-form double-stranded DNA (dsDNA) ofthe Drew-Dickerson sequence d(CGCGAATTCGCG) is computed us-ing classical Molecular Dynamics (MD) simulations. In contrast toprevious studies, which focus on a simplified 1D model or a coarse-grained model of DNA to improve simulation times, full atomistic sim-ulations are employed to understand the thermal conduction in B-DNA. Thermal conductivity at different temperatures from 100 to 400K are investigated using the Einstein Green-Kubo equilibrium andMüller-Plathe non-equilibrium formalisms. The thermal conductivityof B-DNA at room temperature is found to be 1.5 W/m · K in equilib-rium and 1.225 W/m · K in non-equilibrium approach. In addition, thedenaturation regime of B-DNA is obtained from the variation of ther-mal conductivity with temperature. It is in agreement with previousworks using Peyrard-Bishop Dauxois (PBD) model at a temperatureof around 350 K. The quantum heat capacity ( C vq ) has given the ad-ditional clues regarding the Debye and denaturation temperature of12-bp B-DNA. B-DNA | thermal conductivity | heat capacity | molecular dynamics | Green-Kubo | Müller-Plathe P robing thermal conduction in shorter length scales us-ing computations gives us a fundamental understandingof the link between structure and phenomena. DNA satis-fies the low thermal conductivity requirements for buildingmolecular thermoelectric devices (1), This was a motivationfor this study. In an experiment to find the thermal conduc-tion of DNA-gold composite where λ -DNA coated with goldnanoparticles (2), the thermal conductivity ( κ ) was found tobe 150 W/m · K. The later works have cited that gold primar-ily contribute to this high value of thermal conductivity andhave mentioned that the ultra-low thermal conductivity ofDNA (3–5). An experiment using a new transient electro-thermal technique with crystalline DNA composite fibers inNaCl solution, has suggested that this might indeed be thecase as the thermal conductivity found for the fibers are lowaround 0.25-0.85 W/m · K (3, 4).The earlier computational works use the Peyrad-Bishop-Dauxois (PBD) model which is a simplified 1-D non-linearbead spring model (6). This model has been used to under-stand DNA dynamics and to find the denaturation point ofDNA (7). The denaturation regime of DNA was found tobe around 350 K and thermal conductivity of the DNA in(8) at 300 K is 1.8 mW/m · K. This value is much lower thanthat of poly(G) DNA ( ≈ · K) obtained using a 12-coarse grained (12-CG) model (5). The disagreement raisesthe question of the validity of thermal conductivity computedusing PBD model. A few infinite chains formed by permutingAdenine (A) and Guanine (G) sequences-poly(A), poly(G),poly(AG), poly(A G ), poly(A G ) are investigated us-ing PBD model and it was found that the denaturation pointcan be shifted depending on the sequence (9). Additionally,the thermal conductance ratio R = κ H /κ L , which is definedas the high thermal conductivity to low thermal conductivity Fig. 1.
The simulated B-DNA structure is in a water box with Na + ions. The B-DNAhas a total length of 40.8 Å and a radius of 10 Å and is colored according to differentatom types. Also shown is the non-equilibrium Müller-Plathe scheme where the redregion with higher temperature swaps kinetic energies with blue regions with lowertemperatures. of a various sequences is quantified to analyze thermal switch-ing in DNA. The sequence poly(A G ) seemed to have higherthe thermal conductance ratio than poly(A) and poly(G) se-quences. Again, the thermal conductance values were ex-tremely lower, probably owing to the PBD model (9). ThePBD model has also been used to understand the improvingheat conduction through external force (10). As an exten-sion of this, DNA switching was studied by considering oneend of DNA helical turn as drain, another end as source andthe central region as gate(11). An all-atom picture resolvesthe discrepancies in understanding the fundamental nature ofthermal conduction in a DNA. Hence, the water and Na + ionsare considered as the stabilizing media and only the thermalconductivity of B-DNA is calculated... Significance Statement
The thermal conductivity of B-form double-stranded DNA (ds-DNA) of the Drew-Dickerson sequence d(CGCGAATTCGCG)is computed using classical Molecular Dynamics (MD) simula-tions. In contrast to previous studies, which focus on a simpli-fied 1D model or a coarse-grained model of DNA to improvesimulation times, full atomistic simulations are employed to un-derstand the thermal conduction in B-DNA. Thermal conductiv-ity at different temperatures from 100 to 400 K are investigatedusing the Einstein Green-Kubo equilibrium and Müller-Plathenon-equilibrium formalisms.
Vignesh Mahalingam performed research and Vignesh Mahalingam, Dineshkumar Harursampathwrote the paper togetherThe authors declare no conflict of interest. To whom correspondence should be addressed. E-mail: vigneshmiisc.ac.in
October 1, 2020 | vol. XXX | no. XX | R A F T k z ( W / m K ) Temperature(K)
Fig. 2.
Thermal conductivity of B-DNA along the length of DNA as a function oftemperature (4).
Results and discussion
Thermal conductivity using Green-Kubo equilibrium method.
The thermal conductivity of B-DNA is calculated using theequilibrium GK method as in equation (2). This requiresthe calculation of heat-heat auto-correlation function (HCAF)from heat current using equation (3) and this settles to anequilibrium value. The thermal conductivity value obtainedat 300 K along the length of DNA after it has settled toan equilibrium value for 10 ns. All equilibrium GK thermalconductivity values at different temperature are obtained sim-ilarly. The thermal conductivity of the dsDNA sequence alongits length as a function of temperature is shown in Figure 2.The thermal conductivity increases as a function of tempera-ture and eventually saturates around 350 K.The equilibrium method also allows us to calculate theheat flux along the other directions such as those betweenthe strands as shown in Figure 3. These thermal conductiv-ities are inaccessible to experiments at such shorter strandlengths. Interestingly, the heat transfer along base pairs be-tween backbone is higher than that along the strand. Noparallels exist in reported literature about the heat conduc-tion between strands, although (5) has mentioned that theheat conduction along the length is primarily due to sugar-phosphate backbone. The classical treatment here has notaccounted for the transfer of heat between base pairs by tun-neling. Nevertheless, it is evident from Figure 3 that moreheat can be transferred along base pairs than along the lengthof phosphate backbone.
Power density spectrum.
To understand the molecular originof the temperature dependence of the thermal conductivity ofthe DNA, the Density of States (DoS) of the 12-bp B-DNAhave been calculated using a 2-point (2pt) code (12–14) andhave been plotted in Figure 4. Only continuous low frequencymodes can be seen till 800 cm − (Debye frequency, ω D ). Inan earlier work on poly-G DNA (5), the DoS spectrum had agap with no modes between 200 cm − and 300 cm − . More-over, few modes existed beyond this gap till 400 cm − . Nosuch gap is seen in the spectrum between optical and acousticmodes. No high frequency modes has been observed in both k y ( W / m K ) Temperature(K)
Fig. 3.
Thermal conductivity of B-DNA along base-pairs κ y as a function of temper-ature. Fig. 4.
Power spectrum density of B-DNA at 300 K, 100 K and 400 K. this work and in earlier work (5), meaning that the phononmodes have large wavelengths and hence are scattered at thedsDNA boundaries. Around the denaturation temperature,the DNA strands separate, and the thermal conductivity sat-urates. It is possible to compute the Debye temperature,T D from this spectrum by using T D = ~ ω D / k B , where ~ and k b are reduced Planck’s constant and Boltzmann constant, re-spectively. Substituting ω D to be 723 cm − as it the lastavailable frequency mode, T D ≈
165 K.Using the DoS shown in Figure 4, also computed is thequantum molar specific heat capacity ( C vq ) of the dsDNA(Figure 5). The resulting heat capacity is plotted as a func-tion of temperature. The peaks in heat capacity at 200 K and273 K in Figure 5 are probably due to water as similar featurescan be seen in water heat capacity at the same temperaturesin Figure 6. The peak at 150 K is close to the calculatedDebye temperature (165 K) of DNA. The peak around 350K ought to correspond with the denaturation regime, wherethe transition from double strand to two single strands hap-pens (15). It is till this point that the thermal conductivityincreases and beyond which thermal conductivity saturates. et al. R A F T c v q ( d i m e n s i on l ess ) Temperature (K)(a) Bezier fitDebye temperature (T D )Denaturation temperature (T m ) Fig. 5.
Quantum heat capacity ( C vq ) of B-DNA as a function of temperature. Fig. 6.
Quantum heat capacity ( C vq ) of TIP3P water as a function of temperature. A similar study (8) has described the same phenomenon andit is mentioned that the increase in thermal conductivity isstrongly correlated with the anharmonicity in the bond be-tween the complementary base pairs till there is effectively nocontact between the complementary base pairs.
Thermal conductivity using Reverse Non-Equilibrium Molec-ular Dynamics (RNEMD) method.
A non-equilibrium MP ap-proach (16) is also used to understand the low thermal conduc-tivity obtained earlier using equilibrium formulation. Here, atemperature gradient can be set along the length of B-DNAand surrounding water box. The temperature profile acrossthe surrounding water clearly has a gradient between the cen-ter hot region and cold regions on either side. Figures 7 and 8shows the water and DNA temperature profiles, respectively.The thermal conductivity κ due to a linear temperaturegradient between the DNA ends (between 22.5 and 67.5 Å) isgiven as κ z = (cid:0) QAt (cid:1) ( ∂T /∂z ) , [1]where Q is the heat exchange between hot and cold regions, Ais the cross-sectional area of the water box and t is the time for T e m p e r a t u r e ( K ) DNA−length (Å) 100K150K200K250K300K350K400K D T ¶ T/ ¶ x Fig. 7.
Temperature profile across the water box which has B-DNA along z-direction T e m p e r a t u r e ( K ) Water box length (Å) 100K150K200K250K300K350K400K
Fig. 8.
Temperature profile across the 12-bp B-DNA in RNEMD simulation alongz-directionMahalingam et al.
PNAS |
October 1, 2020 | vol. XXX | no. XX | R A F T k z ( W / m K ) Temperature(K)
Fig. 9.
Thermal conductivity versus temperature for 12-bp B-DNA in RNEMD simu-lation. Inset shows the Green-Kubo values plotted in Figure 2. heat exchange. The temperature gradient is computed acrossB-DNA from the temperature profile in Figure 8. Cautionwas exercised in calculating the temperature profile as con-strained SHAKE atoms were excluded and a mild Berendsenthermostat was used (17). Figure 9 shows the temperature de-pendence of the thermal conductivity using this method andthe profile similar to GK thermal conductivity (Figure 2).Till denaturation there is an increase in thermal conductiv-ity κ and heat capacity (C) of B-DNA with respect to temper-ature. This is attributed to increase of phonon density withrespect to temperature (4) as κ = Cv τ r in the single timerelaxation approximation. Moreover, it can be seen from Fig-ure 4 that only low frequency (long wavelength) soft modes areavailable in DNA allowing this classical approximation to bevalid. The heat capacity (C) of B-DNA also increases as De-bye law states: C vq ∝ T (Refer Figure 5). It is being assumedthat the phonon velocity (v) remains almost constant withtemperature and the relaxation time τ r ∝ T − (4). Hence,the thermal conductivity initially increases almost quadrati-cally with respect to temperature till the DNA strands sepa-rate. The thermal conductivity of the DNA-gold compositeis found to be 150 W/m · K (2), whereas recent estimates ofthe thermal conductivity of DNA-water composite mixturevia equilibrium and non-equilibrium MD were 0.381 W/m · Kand 0.373 W/m · K, respectively (18). Our results suggest thatthere is a definite contribution of gold and water to the ther-mal conductivity of the composite mixture in these works andthe thermal conductivity of DNA is somewhere close to thereported values in this work.
1. Discussion
We have examined the thermal conductivity of 12-bp B-DNA,the most common form of DNA from GK calculation fromRNEMD calculation. Both calculations show an increase inthermal conductivity till denaturation temperature. The fullatom description as opposed to coarse grained or 1D non-linear chain models indicate the regions where the modelssucceed and fail. The thermal conductivity obtained usingPBD models needs to be refined as they seem to be very low.Nevertheless, all models predict the denaturation regime close to 350 K, where the thermal conductivity saturates with in-crease in temperature. The 2pt-calculations show that the De-bye temperature is consistent with the earlier works (5, 8, 9).The engineering of thermal conductivity, based on the base-pair sequences along the long lengths, might play a role in itsusage as a molecular thermoelectric device operating at roomtemperature. This work lays the foundation for an all-atomstudy of DNA thermal conductivity. Building further usingthe methods here would give us insight into the dependenceof thermal conductivity on base-pair sequences from all-atomperspective. The thermal conductivity computed here mightbe a necessary validity check for coarse-grained DNA models.
Materials and Methods
All the calculations are performed on the 12-base pair (bp) B-DNAof Drew-Dickerson sequence d(CGCGAATTCGCG) (19). NucleicAcid Builder (NAB) module of AMBERTOOLS18 (20) is used tobuild the initial structures of the double stranded (ds) DNA. ThedsDNA is then placed in a bath of TIP3P water box (21) usingxleap module of AMBERTOOLS18 software package. A water boxwith dimensions of 65 Å ×
65 Å ×
90 Å ( x × y × z ) is chosen to ensure15 Å solvation shell around the B-DNA. 22 Na + ions are added atthe lowest electrostatic potential locations to the solvated dsDNAsystem. DNA OL15 force-field (22) is used. This has parmbsc0 (23)and OL15 (22) corrections to the ff99 force-field (24) to considerthe bonded and non-bonded interactions of the dsDNA. Joung-Cheathem parameter (25) set is used to consider the interaction ofmonovalent Na + ions with TIP3P water and dsDNA. After prepar-ing the initial system using AMBERTOOLS18 (20), LAMMPS (26)software module is used for all further simulations. The whole sol-vated dsDNA system is energy minimized using first 5000 stepsof steepest descent and 5000 steps of conjugate gradient keepingthe B-DNA restrained with a force of 500 kcal/molÅ. The DNAis then slowly released into water by reducing the restraint from20 kcal/molÅ to 0 kcal/molÅ in 5 cycles of the steepest descentand conjugate gradient minimization steps. All the atoms are thenassigned velocities according to Maxwell-Boltzmann distribution.Throughout the MD simulation, the DNA has a small restraint of 1kcal/molÅ to prevent the same from changing its orientation whilstmeasuring thermal conductivity. SHAKE constraints (27) are ap-plied to the hydrogen atoms, bond and angles of DNA and waterwith a tolerance of 10 − (24). The system is equilibrated for 10ns with Nosé-Hoover thermostat and barostats with coupling con-stants 0.1 ps and 1.0 ps, respectively (28–31). Finally, a productionrun of 20 ns for the calculation of thermal conductivity ensures thatthe thermal conductivity values converge. The thermal conductiv-ity is computed using the equilibrium Green-Kubo (GK) method,where the heat-heat auto-correlation function is used to computethe thermal conductivity as (32–34): κ x,y,z = Vk B T Z ∞ h J x,y,z (0) · J x,y,z ( t ) i dt, [2]where thermal conductivity κ x,y,z at a temperature T in a direc-tion x, y or z is obtained from heat current J x,y,z in that direction. k B is the Boltzmann constant. The heat current is obtained as J = 1 V "X i e i v i − X i σ i v i [3]where e i is the total energy of an atom, v i is the velocity of anatom, σ i is the virial stress per atom and V is the volume of thetotal group of atoms. The exact volume (V) of B-DNA is computedfrom the atomic volumes of adenine (136.1 Å ), guanine ((143.8Å ), cystosine (113.2 Å )and thymine (132.6 Å ) base-pair groupsand sugar-phosphate (174.8 Å ) groups (35, 36). Each strand (leftor right of B-DNA symmetrical axis) contains 4 cytosine, 4 guanine,2 adenine and 2 thymine and 12 sugar-phosphate groups giving usthe total volume of the 12-bp B-DNA to be 7326 Å . The powerspectrum density or Density of States (DoS) of the 12-bp B-DNA et al. R A F T is obtained from a fast Fourier transform of velocity-velocity auto-correlation, C(t) as DoS ( ν ) = lim t →∞ k B T τ Z − τ C ( t ) e − πνt dt, [4]where t is correlation time window of 200 ps and ν is the fre-quency. Only the solid component of DoS is considered as liquidand gaseous states are not relevant for B-DNA. The canonical par-tion function (Q) can be constructed from DoS, with a harmonicoscillator assumption:ln Q = ∞ Z DoS ( ν ) q HO ( ν ) dν, [5]where q HO = e − βhν − e − βhν is the harmonic oscillator partition func-tion, β = k B T and h is the Planck’s constant. The entropy S andthe heat capacity C vq are then found using the partition functionand DoS as S = k ln Q + β − (cid:16) ∂ ln Q∂T (cid:17)
N,V = k ∞ Z DoS ( ν ) W s ( ν ) dν [6] C vq = (cid:18) ∂S ∂T (cid:19) = k (cid:16) ∂ ln Q∂T (cid:17)
N,V + β − (cid:18) ∂ ln Q∂T (cid:19) N,V = k β − ∞ Z DoS ( ν ) W C ( ν ) dν [7]with weighting functions W s ( ν ) = βhνe βhν − (cid:2) − (cid:0) e − βhν (cid:1)(cid:3) , [8] W C ( ν ) = e − βhν (cid:2) − (cid:0) e − βhν (cid:1)(cid:3) . [9] ACKNOWLEDGMENTS.
The authors thank Dr. Navaneetha Kr-ishnan, Dr. Prabal Maiti and Abhishek Aggarwal for their fruitfuldiscussions and suggestions.
1. LL Nian, W Liu, L Bai, XF Wang, Spin caloritronics in a chiral double-strand-dna-based hybridjunction.
Phys. Rev. B , 195430 (2019).2. T Kodama, A Jain, KE Goodson, Heat conduction through a dna-gold composite. NANOLETTERS , 2005–2009 (2009).3. Z Xu, S Xu, X Tang, X Wang, Energy transport in crystalline dna composites. AIP ADVANCES (2014).4. Z Xu, X Wang, H Xie, Promoted electron transport and sustained phonon transport by dnadown to 10 k. Polymer , 6373 – 6380 (2014).5. AV Savin, MA Mazo, IP Kikot, LI Manevitch, AV Onufriev, Heat conductivity of the dna doublehelix. Phys. Rev. B , 245406 (2011).6. M Peyrard, Nonlinear dynamics and statistical physics of DNA. Nonlinearity , R1–R40(2004).7. M Peyrard, AR Bishop, Statistical mechanics of a nonlinear model for dna denaturation. Phys.Rev. Lett. , 2755–2758 (1989).8. KA Velizhanin, CC Chien, Y Dubi, M Zwolak, Driving denaturation: Nanoscale thermal trans-port as a probe of dna melting. Phys. Rev. E , 050906 (2011).9. CC Chien, KA Velizhanin, Y Dubi, M Zwolak, Tunable thermal switching via DNA-based nano-devices. Nanotechnology , 095704 (2013).10. S Behnia, R Panahinia, Ballistic induced pumping of hypersonic heat current in DNA nanowire. EUROPEAN PHYSICAL JOURNAL B (2016).11. S Behnia, R Panahinia, Molecular thermal transistor: Dimension analysis and mechanism. CHEMICAL PHYSICS , 40–46 (2018).12. ST Lin, M Blanco, WA Goddard, The two-phase model for calculating thermodynamic prop-erties of liquids from molecular dynamics: Validation for the phase diagram of lennard-jonesfluids.
The J. Chem. Phys . , 11792–11805 (2003).13. ST Lin, PK Maiti, WA Goddard, Two-phase thermodynamic model for efficient and accurateabsolute entropy of water from molecular dynamics simulations. The J. Phys. Chem. B ,8191–8198 (2010). 14. TA Pascal, ST Lin, WA Goddard III, Thermodynamics of liquids: standard molar entropiesand heat capacities of common solvents from 2pt molecular dynamics.
Phys. Chem. Chem.Phys. , 169–181 (2011).15. A Wildes, et al., Structural correlations and melting of b-dna fibers. Phys. Rev. E , 061923(2011).16. F Müller-Plathe, A simple nonequilibrium molecular dynamics method for calculating the ther-mal conductivity. The J. Chem. Phys . , 6082–6085 (1997).17. M Zhang, E Lussetti, LES de Souza, F Müller-Plathe, Thermal conductivities of molecularliquids by reverse nonequilibrium molecular dynamics. The J. Phys. Chem. B , 15060–15067 (2005).18. NA Jolfaei, et al., Investigation of thermal properties of dna structure with precise atomicarrangement via equilibrium and non-equilibrium molecular dynamics approaches.
Comput.Methods Programs Biomed . , 105169 (2020).19. R Dickerson, et al., The anatomy of a-, b-, and z-dna. Science , 475–485 (1982).20. D Case, et al., Amber 2018 (2018).21. WL Jorgensen, J Chandrasekhar, JD Madura, RW Impey, ML Klein, Comparison of simplepotential functions for simulating liquid water.
The J. Chem. Phys . , 926–935 (1983).22. M Zgarbova, et al., Refinement of the Sugar-Phosphate Backbone Torsion Beta for AMBERForce Fields Improves the Description of Z- and B-DNA. JOURNAL OF CHEMICAL THEORYAND COMPUTATION , 5723–5736 (2015).23. A Pérez, et al., Refinement of the amber force field for nucleic acids: Improving the descriptionof α / γ conformers. Biophys. J . , 3817 – 3829 (2007).24. J Wang, P Cieplak, P Kollman, How well does a restrained electrostatic potential (RESP)model perform in calculating conformational energies of organic and biological molecules? JOURNAL OF COMPUTATIONAL CHEMISTRY , 1049–1074 (2000).25. IS Joung, TE Cheatham, Determination of alkali and halide monovalent ion parameters foruse in explicitly solvated biomolecular simulations. The J. Phys. Chem. B , 9020–9041(2008).26. S Plimpton, Fast parallel algorithms for short-range molecular dynamics.
J. Comput. Phys . , 1 – 19 (1995).27. JP Ryckaert, G Ciccotti, HJ Berendsen, Numerical integration of the cartesian equations ofmotion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys . ,327 – 341 (1977).28. W Shinoda, M Shiga, M Mikami, Rapid estimation of elastic constants by molecular dynamicssimulation under constant stress. Phys. Rev. B , 134103 (2004).29. GJ Martyna, DJ Tobias, ML Klein, Constant pressure molecular dynamics algorithms. The J.chemical physics , 4177–4189 (1994).30. M Parrinello, A Rahman, Polymorphic transitions in single crystals: A new molecular dynam-ics method.
J. Appl. physics , 7182–7190 (1981).31. ME Tuckerman, J Alejandre, R López-Rendón, AL Jochim, GJ Martyna, A liouville-operatorderived measure-preserving integrator for molecular dynamics simulations in the isothermal–isobaric ensemble. J. Phys. A: Math. Gen . , 5629 (2006).32. M GREEN, MARKOFF RANDOM PROCESSES AND THE STATISTICAL MECHANICS OFTIME-DEPENDENT PHENOMENA. JOURNAL OF CHEMICAL PHYSICS , 1281–1295(1952).33. M GREEN, MARKOFF RANDOM PROCESSES AND THE STATISTICAL MECHANICS OFTIME-DEPENDENT PHENOMENA .2. IRREVERSIBLE PROCESSES IN FLUIDS. JOUR-NAL OF CHEMICAL PHYSICS , 398–413 (1954).34. R Kubo, Statistical-mechanical theory of irreversible processes. i. general theory and simpleapplications to magnetic and conduction problems. J. Phys. Soc. Jpn . , 570–586 (1957).35. K Nadassy, I Tomás-Oliveira, I Alberts, J Janin, S Wodak, Standard atomic volumes in double-stranded dna and packing in protein-dna interfaces. Nucleic Acids Res . , 3362–3376(2001).36. N Voss, M Gerstein, Calculation of standard atomic volumes for rna and comparison withproteins: Rna is packed more tightly. J. Mol. Biol . , 477 – 492 (2005). Mahalingam et al.
PNAS |
October 1, 2020 | vol. XXX | no. XX || vol. XXX | no. XX |