Atomistic and coarse-grained simulations of DNAelectrostaticsbyShahzad Ghanbarian AlavijehM.Sc., Physics, Isfahan University of Technology, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Physics)The University of British Columbia(Vancouver)April 2015c© Shahzad Ghanbarian Alavijeh, 2015AbstractThis thesis aims to understand the nanoscale effects of water as a structured sol-vent on the phenomenon of counterion-induced DNA condensation. We presentresults of molecular dynamics simulations of the electrostatic interaction betweentwo DNA molecules in the presence of divalent counterions in different solvationmodels. We also develop a coarse-grained implicit solvent model for investigatingthe dynamics on long time and length scales.In the first project, we investigate the role of the solvation effects on the interac-tion between like-charged cylindrical rods as simplified model for DNA molecules.We obtain the average force between two parallel charged rods in simulations thatdiffer only by their representation of water as a implicit or explicit solvent, buthave otherwise identical parameters. We find that the presence of water moleculeschanges the structure of the counterions and results in both qualitative and quanti-tative changes of the force between highly charged polyelectrolytes.In the second project, we explore the importance of the DNA geometry on theelectrostatic forces by considering two rigid helical models. The simulation resultsindicate that the DNA shape is an essential contributor to the interaction. A regimeof attractive interaction, which disappeared in the cylindrical model, is recovered inthe explicit solvent model in both types of helical models. The results also confirmthat the behaviour of the interactions between two DNA molecules in the explicitsolvent model are different from the implicit solvent models.In the third project, we develop a coarse-grained (CG) representation of thesesolvation effects. This CG model is constructed from explicit simulations and sig-nificantly reduces the computational expense. Short-ranged corrections are addedto the pair-wise interaction potentials in the implicit solvent model such that theiistructure of counterions in the system is consistent with the results from the ex-plicit solvent simulations. This CG model succeeds in reproducing the like-chargeattraction effect between DNA molecules in explicit simulations.In a final project, we apply the CG model developed previously to study threeDNA strands in the presence of divalent counterions as a starting point for investi-gating many-body effects in the mechanism of DNA bundling.iiiPrefaceThe research described in Chapter 3 forms the basis of a co-authored, peer-reviewedjournal article, published as S. Ghanbarian and J. Rottler, Solvation effects on like-charge attraction, J. Chem. Phys. 138,084901(2013). The results of Chapter 4-6are currently being prepared for publication in manuscript form. Parts of the re-sults of chapter 4 and 5 have been presented in APS March meeting 2015. Mycontributions to the work described in this dissertation are as follows:• Designed the research project with my supervisor.• Performed all the simulations.• Performed all the data analysis.• Prepared all of the figures in the publications.• Wrote the first draft of the publication.• Writing of the text for the publication was shared with my supervisor.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixLists of Abbreviations . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 The DNA molecule . . . . . . . . . . . . . . . . . . . . . . . . . 21.2.1 DNA models with different resolution . . . . . . . . . . . 31.3 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . 51.3.1 Analytical approaches . . . . . . . . . . . . . . . . . . . 51.3.2 Implicit solvent simulations of cylindrical models . . . . . 61.3.3 Implicit solvent simulations of helical models . . . . . . . 71.3.4 Fully atomistic models . . . . . . . . . . . . . . . . . . . 81.3.5 Experimental results . . . . . . . . . . . . . . . . . . . . 91.4 Challenges for molecular simulations of charged biomolecules . . 91.5 Thesis outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11v2 Polyelectrolytes in solution: theory and simulations . . . . . . . . . 142.1 Poisson-Boltzmann theory . . . . . . . . . . . . . . . . . . . . . 142.2 Molecular dynamics simulations . . . . . . . . . . . . . . . . . . 172.2.1 Interaction potentials . . . . . . . . . . . . . . . . . . . . 172.3 Solvation models . . . . . . . . . . . . . . . . . . . . . . . . . . 182.3.1 Implicit solvent model . . . . . . . . . . . . . . . . . . . 182.3.2 Explicit solvent model . . . . . . . . . . . . . . . . . . . 192.4 Velocity Verlet . . . . . . . . . . . . . . . . . . . . . . . . . . . 223 Cylindrical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.2 Models and simulations . . . . . . . . . . . . . . . . . . . . . . . 253.2.1 Implicit solvent model . . . . . . . . . . . . . . . . . . . 263.2.2 Explicit solvent model . . . . . . . . . . . . . . . . . . . 303.2.3 Correlation functions . . . . . . . . . . . . . . . . . . . . 353.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384 Helical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 394.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 394.2 DNA models . . . . . . . . . . . . . . . . . . . . . . . . . . . . 404.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 424.3.1 Implicit solvent model . . . . . . . . . . . . . . . . . . . 424.3.2 Explicit solvent model . . . . . . . . . . . . . . . . . . . 444.3.3 Correlation functions . . . . . . . . . . . . . . . . . . . . 464.3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . 525 Coarse-grained model . . . . . . . . . . . . . . . . . . . . . . . . . . 555.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 555.2 Iterative Boltzmann inversion . . . . . . . . . . . . . . . . . . . . 565.2.1 Coarse-grained potential for ionic species . . . . . . . . . 585.3 Simulation details . . . . . . . . . . . . . . . . . . . . . . . . . . 615.3.1 Coarse-grained simulations . . . . . . . . . . . . . . . . . 615.3.2 Fully detailed atomistic simulations . . . . . . . . . . . . 625.4 Results and discussion . . . . . . . . . . . . . . . . . . . . . . . 65vi5.4.1 Local structure from initial guess . . . . . . . . . . . . . 655.4.2 Model refinement . . . . . . . . . . . . . . . . . . . . . . 705.4.3 Interaction between double stranded DNA molecules . . . 705.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 726 Three-DNA model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 736.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 736.2 Pairwise additivity of interactions in the three-DNA system . . . . 746.3 Coarse-grained three-DNA model . . . . . . . . . . . . . . . . . 776.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 787 Conclusions and future work . . . . . . . . . . . . . . . . . . . . . . 797.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 797.2 Future work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94A Interaction potentials . . . . . . . . . . . . . . . . . . . . . . . . . . 94B Blocking method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97viiList of TablesTable 2.1 Parameters of SPC/E model . . . . . . . . . . . . . . . . . . . 20Table 3.1 Parameters for the Lennard-Jones interaction of different atoms( rod ion(RI), counter ion (CI) and oxygen (O)). There is noLennard-Jones interaction between Hydrogen and other parti-cles in the system. The potential is truncated at rc. . . . . . . . 27Table 3.2 Properties of the systems. . . . . . . . . . . . . . . . . . . . . 27Table 4.1 Parameters for the Lennard-Jones interaction of different atoms(central atoms (CA), phosphate (P), counter ion (CI) and oxy-gen (O)). There is no Lennard-Jones interaction between Hy-drogen and other particles in the system. The potential is trun-cated at rc. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43Table 5.1 Parameters for the interaction of different atoms (central atoms(CA), phosphate (P), counter ion (CI)) of initial guess whichcorrespond to Fig. 5.1-5.5. . . . . . . . . . . . . . . . . . . . . 56Table 5.2 Parameters for the hydration contribution in the potential of dif-ferent atoms of initial guess which correspond to Fig. 5.1-5.5. 57Table 5.3 Parameters for the interaction of different atoms in final step,corresponds to Figs. 5.6- 5.9. . . . . . . . . . . . . . . . . . . 68Table 5.4 Parameters for the hydration contribution in the potential of dif-ferent atoms of final step which corresponds to Figs. 5.6- 5.9. . 68viiiList of FiguresFigure 1.1 Schematic structure of a double strand of DNA in the B-form(sequence: CGCGAATTCGCG). Hydrogen atoms are not shownexplicitly. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2Figure 1.2 Different DNA models, (a) charged cylinder, (b) cylindricalmodel (CM), (c) extended cylindrical model (ECM), (d) Montoro-Abascal model (MAM). Phosphate charges are shown as bluespheres and the neutral cylindrical core is coloured in yellow.For the MAM, another neutral sphere is added between thecylindrical core and the charged phosphate sphere which iscoloured in green. . . . . . . . . . . . . . . . . . . . . . . . 4Figure 2.1 Schematic of the SPC/E model. . . . . . . . . . . . . . . . . 21Figure 3.1 Schematic depiction of the simulation geometry. (a) Two chargedpolyelectrolytes are modeled as rigid rods, which consist of rodions separated by b = 1.7A˚. Each rod ion is modeled as a LJsphere with diameter σ1 and charge q1 = −e at the center ofsphere. Each counterion is modeled as a LJ sphere with diame-ter σ2 and charge q2 = 2e at the center of sphere. (b) Thin rodssystem, where rod ions and counterions have the same size,σ1 = σ2 = 3.17A˚. (c) Thick rods system (DNA like), wherethe size of rod ions is almost six times greater than the size ofcounterions, σ1 = 20A˚ and σ2 = 3.166A˚. . . . . . . . . . . . 26ixFigure 3.2 Average force between (a) thin and (b) thick rods system as afunction of rod separation (surface to surface) in implicit sol-vent simulations. The crosses and heavy dots represent forcein implicit system with and without counterions, respectively.(c) Ratio of forces with and without counterions shown in pan-els (a) and (b). The crosses and open circles show this ratio forthick and thin rods system, resp. . . . . . . . . . . . . . . . . 28Figure 3.3 (a) Force in three different models of thin rods as a function ofrod separation (surface to surface). Squares show the force inan explicit model with counterions, crosses show the force inan explicit model without counterions and dots show the forcein an implicit model with counterions. (b) Ratio of the forcesin explicit solvent simulations with and without counterions. . 29Figure 3.4 (a) Force in three different models of thick rods system as afunction of rod separation (surface to surface). Squares showthe force in an explicit model with counterions, crosses showthe force in an explicit model without counterions and dotsshow the force in an implicit model with counterions. (b) Ratioof the forces in explicit solvent simulations with and withoutcounterions. . . . . . . . . . . . . . . . . . . . . . . . . . . . 30Figure 3.5 Cross-sectional picture of distribution of counterions aroundthe rod ions (purple spheres) for explicit model of thin rodssystem with rod separation, (a) d = 0.9A˚ and (b) d = 1.8A˚and (c) d = 4A˚ . . . . . . . . . . . . . . . . . . . . . . . . . . 31Figure 3.6 (a) Counterion-rod ion pair correlation functions in thin rodssystem with implicit solvent for three rod separations (surfaceto surface), d = 0.5,2,4A˚, represented by squares, heavy dotsand crosses. (b) Counterion-rod ion pair correlation functionsin thin rods systems with explicit solvent for three rod sep-arations (surface to surface), d = 0.9,1.8,4A˚ represented bysquares, heavy dots and crosses. . . . . . . . . . . . . . . . . 32xFigure 3.7 (a) Counterion-rod ion pair correlation functions in thick rodsystems with implicit solvent for three rod separations (sur-face to surface), d = 2,4,6A˚ , represented by squares, heavydots and crosses. (b) Counterion-rod ion pair correlation func-tions in thick rod systems with explicit solvent for three rodseparations (surface to surface), d = 1.3,1.8,4A˚ , representedby squares, heavy dots and crosses. . . . . . . . . . . . . . . 33Figure 3.8 Water molecule-rod ion pair correlation function in explicitthin rods system (a) with and (b) without counterions. The re-sults for three rod separations (surface to surface), d= 0.9,1.8,4A˚are represented by squares, heavy dots and crosses, respec-tively, for both systems. . . . . . . . . . . . . . . . . . . . . . 34Figure 3.9 Water molecule-rod ion pair correlation function in explicitthick rods system (a) with and (b) without counterions for threerod separations (surface to surface), d = 0.2,1.8,4A˚ are repre-sented by squares, heavy dots and crosses, respectively. . . . . 36Figure 3.10 Number density plot of water molecules in explicit solventmodel for two different thin rod separations (surface to sur-face), d = 0A˚ for systems (a) with and (b) without counterionsand d = 4A˚, for systems (c) with and (d) without counterions. 37Figure 4.1 Extended cylinder model. Phosphate charges are shown asblue spheres and the cylindrical core is coloured in yellow. . . 40Figure 4.2 Grooved or Montoro-Abascal model. Phosphate charges areshown as blue spheres. The cylindrical core is coloured in yel-low and the neutral particles as green spheres. . . . . . . . . . 41Figure 4.3 A schematic picture explaining the positions of DNA moleculesand the definition of the different azimuthal angles φs, φ0 andφ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42Figure 4.4 Reduced DNA-DNA interaction force 〈F〉/F0 vs separationdistance R with an implicit solvent in the presence of diva-lent counter ions. The circles and squares represent force forthe ECM and MAM, respectively. . . . . . . . . . . . . . . . 44xiFigure 4.5 Reduced DNA-DNA interaction force 〈F〉/F0 as a function ofDNA separation (centre to centre) with an explicit solvent. Thecircles and squares represent force for the ECM and MAM,respectively. . . . . . . . . . . . . . . . . . . . . . . . . . . . 45Figure 4.6 Reduced DNA-DNA interaction force 〈F〉/F0 as a function ofDNA separation (centre to centre), considering extended cylin-der model (ECM) for the DNA structure. The crosses and cir-cles represent force in explicit and implicit solvent model, re-spectively. The inset shows a zoom for the reduced force withimplicit model. . . . . . . . . . . . . . . . . . . . . . . . . . 47Figure 4.7 Reduced DNA-DNA interaction force 〈F〉/F0 as a functionof DNA separation (centre to centre), considering Montoro-Abascal model (MAM) for the DNA structure. The crosses andcircles represent force in implicit and explicit solvent model,respectively. The inset shows a zoom for the reduced forcewith implicit model. . . . . . . . . . . . . . . . . . . . . . . 48Figure 4.8 Phosphate-counterion radial pair correlation function of the ex-tended cylinder model in (a) explicit and (b) implicit systemfor three DNA separations (centre to centre), R = 23,25,28[A˚], represented by heavy dots, squares and triangles. . . . . . . . 49Figure 4.9 Phosphate-counterion radial pair correlation function of theMontoro-Abascal model in (a) explicit and (b) implicit systemfor three DNA separations (centre to centre), R= 22.5,25.25,28[A˚], represented by heavy dots, squares and triangles. . . . . . . . 50Figure 4.10 Central atom-counterion pair correlation functions of the ex-tended cylinder model in (a) explicit and (b) implicit systemfor three DNA separations (centre to centre), R = 23,25,28[A˚], represented by heavy dots, squares and triangles. . . . . . . 51Figure 4.11 Central atom-counterion radial pair correlation function of theMontoro-Abascal model for three DNA separations (centre tocentre), R = 22.5,25.25,28[A˚] , represented by heavy dots,squares and triangles in the explicit solvent model and circles,pluses and crosses in the implicit solvent model. . . . . . . . 52xiiFigure 4.12 Counterion-counterion radial pair correlation function of theextended cylinder model model in (a) explicit and (b) implicitsystem for three DNA separations (centre to centre), R= 23,25,28[A˚], represented by heavy dots, squares and triangles. . . . . . . . 53Figure 4.13 Counterion-counterion radial pair correlation function of Montoro-Abascal model in (a) explicit and (b) implicit system for threeDNA separations (centre to centre), R = 22.5,25.25,28[A˚] ,represented by heavy dots, squares and triangles. . . . . . . . 54Figure 5.1 Development of the central atom-counterion potential. Solidlines show the radial distribution function g(r) (panel (b)) andits Boltzmann inverse (panel (a)) obtained from a fully atom-istic simulation of one DNA strand. Dashed lines show theparameterized coarse grained potential (best fit to the potentialof mean force) and the rdf that is obtained in a coarse grainedsimulation of two DNA strands separated by 28 A˚. . . . . . . 59Figure 5.2 Development of the grooved atom-counterion potential. Solidlines show the radial distribution function g(r) (panel (b)) andits Boltzmann inverse (panel (a)) obtained from a fully atom-istic simulation of one DNA strand. Dashed lines show theparameterized coarse grained potential (best fit to the potentialof mean force) and the rdf that is obtained in a coarse grainedsimulation of two DNA strands separated by 28 A˚. . . . . . . 60Figure 5.3 Development of the phosphate-counterion potential. Solid linesshow the rdfs in full atomistic simulations (see text, panel (b))and the Boltzmann inverse (panel (a)). Dashed lines show theinitial guess of the coarse grained potential obtained from theBoltzmann inverse (panel (a)) and the rdf that is observed whenthe CG simulation of 2 DNA strands is carried out (panel (b)). 62xiiiFigure 5.4 Development of the counterion-counterion potential. Solid linesshow the rdfs in full atomistic simulations (see text, panel (b))and the Boltzmann inverse (panel (a)). Dashed lines show theinitial guess of the coarse grained potential obtained from theBoltzmann inverse (panel (a)) and the rdf that is observed whenthe CG simulation of 2 DNA strands is carried out (panel (b)). 63Figure 5.5 Development of the phosphate-phosphate potential. Solid linesshow the rdfs in full atomistic simulations (see text, panel (b))and the Boltzmann inverse (panel (a)), while dashed line showthe coarse grained potential obtained from the Boltzmann in-verse (panel (a)). . . . . . . . . . . . . . . . . . . . . . . . . 64Figure 5.6 a) Phosphate-counterion potential and b) Radial pair correla-tion function. Solid lines, dashed lines, dotted lines, dot-dashlines and circles represent the function in the explicit model,initial guess and iteration 1, 2, 3 of the CG model, respectively. 65Figure 5.7 Comparison of central atom-counterion radial pair correlationfunction in different solvation models. Solid line, dashed lineand circles represent this function in the explicit model, initialguess and final iteration of CG model, respectively. . . . . . . 66Figure 5.8 Comparison of groove atom-counterion radial pair correlationfunction in different solvation models. Solid line, dashed lineand circles represent this function in the explicit model, initialguess and final iteration of CG model, respectively. . . . . . . 67Figure 5.9 Comparison of counterion-counterion radial pair correlationfunction in different solvation models. Solid line, dashed lineand circles represent this function in the explicit model, initialguess and final iteration of CG model, respectively. . . . . . . 69Figure 5.10 Reduced force between double stranded DNA molecules indifferent solvation models. Solid line, circles, squares andtriangles represent this interaction in the old implicit solventmodel, explicit solvent model, and initial guess and final itera-tion of the CG model, respectively. . . . . . . . . . . . . . . . 71xivFigure 6.1 Three DNA configuration (shown here for the Montoro-Abascalmodel). Phosphate charges are shown as blue spheres. Thecylindrical core is coloured in yellow and the neutral particlesas green spheres. The DNA molecules are located at the ver-tices of an equilateral triangle. Each two DNA molecules areseparated by R. . . . . . . . . . . . . . . . . . . . . . . . . . 74Figure 6.2 Reduced DNA-DNA interaction force 〈F〉/F0 as a function ofDNA separation (centre to centre) for the extended cylindricalmodel (ECM) in primitive electrolyte approximation. The unitof the force is F0 = kBT/P, where P = 34A˚ is the DNA pitchlength. The circles represent force in the three-DNA-system.The triangles represent the force on DNA in the two-DNA-system times√3. . . . . . . . . . . . . . . . . . . . . . . . . 75Figure 6.3 Reduced DNA-DNA interaction force 〈F〉/F0 as a function ofDNA separation (centre to centre) for the Montoro-Abascalmodel (MAM) in primitive electrolyte approximation. Thecircles represent force in the three-DNA-system. The squaresrepresent the force on DNA in the two-DNA-system times√3. 76Figure 6.4 Reduced DNA-DNA interaction force 〈F〉/F0 in three-DNAmodel in the presence of divalent counterions, considering dif-ferent solvation models. Circles and triangles represent thisforce in the coarse-grained model in two- DNA-system times√3 and in three-DNA model, respectively. Solid line repre-sents the force on DNA for the primitive solvent model. . . . 77Figure B.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99xvLists of AbbreviationsBD Brownian dynamicsCA Central atomCG Coarse-grainedCI CounterionCM Cylindrical modelDH Debye-Hu¨ckelDNA Deoxyribonucleic acidds-DNA Double stranded DNAECM Extended cylindrical modelIBI Iterative Boltzmann inversionLJ Lennard-JonesLAMMPS Large-scale Atomic/Molecular Massively Parallel SimulatorMD Molecular dynamicsMAM Montoro-Abascal modelO OxygenP PhosphatePB Poisson-BoltzmannPPPM Particle-particle-particle-meshRI Rod ionSPC Simple point chargeSPC/E Simple point charge/ExtendedxviAcknowledgmentsFirst and foremost, I would like to offer my sincerest gratitude to my supervi-sor, Dr. Jo¨rg Rottler for invaluable advice, insightful guidance, and continuouspatience. He has supported me and encouraged me throughout my PhD study.Honestly without him, this research would not be possible. I would like to thankhim for his consideration and understanding of the conditions and problems that Iwas involved with as a new international student.I would also like to thank the members of my doctoral committee, Dr. GrenfellPatey, Dr. Mona Berciu and Dr. Joshua Folk for their guidance and support duringmy research project. Special thanks to Dr. Grenfell Patey for his constructive sug-gestions and technical comments throughout my research.I would like to express my gratitude to my former supervisor, Dr. KeivanAghababaei Samani who helped me with his brightness and wide range of knowl-edge and let me follow my interests.I would like to thank my great friends around the world who inspired and mo-tivated me throughout my research. Special thanks to Hamid Omid and ArmanAkbarian for invaluable discussions and technical comments.Lastly, and most importantly, I feel indebted to my great parents, my lovelysisters Farnaz and Behnaz and my dear twin brother Behzad for their everlastinglove, support and kindness throughout my entire life.xviiChapter 1Introduction1.1 MotivationMany properties of biomolecules can be explained by considering electrostatic in-teractions. One of the best known examples is DNA (deoxyribonucleic acid), thecarrier of genetic information, which is a long flexible polymer that becomes neg-atively charged when it is solvated in an aqueous environment. Other importantexamples include actin filaments, which are much stiffer polymers that contributeto the mechanical behaviour of cells. In both cases, long ranged electrostatic inter-actions play an essential role in the complex structure and dynamics [1, 2]. Over thepast few decades, many efforts have been made to better understand electrostaticeffects of macromolecules in solution. One of the most interesting and counterin-tuitive phenomena is the emergence of attractive interactions between like-chargedpolyelectrolytes in the presence of mobile multivalent counterions, which has firstbeen noted by Oosawa [3]. These attractions are important for the condensationand self-organisation of DNA as well as other highly charged macroions and poly-electrolytes [2, 4–6]. In the phenomenon of DNA condensation, DNA forms a con-figuration in which helical strands are approximately parallel to each other [7, 8].Despite many studies, the detailed mechanism of this process is still not well un-derstood. There is a strong interest in studying the mechanism of DNA bundlingbecause of its key role in packing DNA into the cell nucleus [9–14]. This attractiveforce is also important for gene delivery and gene therapy and therefore relevant1Figure 1.1: Schematic structure of a double strand of DNA in the B-form(sequence: CGCGAATTCGCG). Hydrogen atoms are not shown ex-plicitly.for biochemists and drug researchers [15, 16].1.2 The DNA moleculeDNA is perhaps one of the best studied biopolymers on earth. It carries all biologi-cal information of all living systems and was first observed by Frederich Miescherin the late 1800s [4]. A double stranded DNA molecule consists of two singlestrands which coil around the same axis and form a double helix with a pitch of34A˚ and a radius of 10A˚ [17]. Each strand is made of units called nucleotides.Each nucleotide consist of a phosphate group, a sugar group and one of the fournucleobases (guanine (G), adenine (A), thymine (T), or cytosine (C)), which arelinked to the sugar group. Repeating these units forms the backbone of the DNAwith alternating phosphate and sugar groups. Hydrogen bonds between nucleotides2bind two complementary strands to each other. These hydrogen bonds form the in-ner core of a double strand with a diameter of approx 20A˚. The sequence of thenucleobases along each strand encodes all the genetic information. The strands arecloser to each other on one side of DNA (minor groove) and farther apart on theother side (major groove). The formation of the grooves is a simple consequence ofthe geometry of the base pairs. The phosphate group is the only charged elementof the DNA molecule and has a charge q = −e. As a result, DNA carries a linecharge density of −e/1.7A˚.The helices can take on different forms depending on the DNA’s environment.B-DNA is one of the most common forms of DNA at neutral pH and physio-logical salt concentrations. Some studies show that condensed DNA is in the B-conformation [18, 19]. Fig. 1.1 shows an atomistic representation of three helicalturns of B-DNA with atomic positions taken from the protein data bank (PDB).1.2.1 DNA models with different resolutionAs can be seen in Fig. 1.1, the detailed atomistic structure of DNA is quite compli-cated. In order to study the structure, properties and dynamics of DNA moleculesefficiently, different models have been developed for double stranded DNA (ds-DNA) at different levels of resolution. All of these models treat DNA as a rigidpolyelectrolyte and neglect bending fluctuations. The most simplified model isthe charged cylinder, in which the line charge density of the DNA molecule isdistributed along the cylinder uniformly. This model is shown in column (a) ofFig. 1.2 and may be taken as the simplest possible starting point for a study fo-cusing on electrostatic interactions. In the cylindrical model (CM), the DNA ismodelled as a chain of discrete spheres. The charge of each phosphate in the realDNA molecule is considered as a point charge and is located at the centre of thespheres, see column (b) of Fig. 1.2. The spheres are separated by b = 1.7A˚ whichcorresponds to the separation of phosphate groups on a DNA strand. Progress-ing further in complexity, helical models have been suggested in which the chargepattern is modelled explicitly as a double-helix structure. One of these models iscalled the extended cylindrical model (ECM) and was designed by Lyubartsev etal. [20]. In this model, the phosphate charge distribution is described by discrete3Figure 1.2: Different DNA models, (a) charged cylinder, (b) cylindricalmodel (CM), (c) extended cylindrical model (ECM), (d) Montoro-Abascal model (MAM). Phosphate charges are shown as blue spheresand the neutral cylindrical core is coloured in yellow. For the MAM,another neutral sphere is added between the cylindrical core and thecharged phosphate sphere which is coloured in green.charges in a double helical array corresponding to the geometry of B-DNA, seecolumn (c) of Fig. 1.2. Another helical structure, which resembles the real DNAappearance even more closely, is called Montoro-Abascal model (MAM) [21]. Inthis model, another neutral sphere is added between the cylindrical core and thecharged phosphate sphere to incorporate the accessibility of the DNA grooves tothe counterions. This model gives a fairly realistic description of the molecularstructure of B-DNA and takes to account the molecular shape of DNA by mod-elling major and minor grooves along the strands, see column (d) of Fig. 1.2.41.3 Literature review1.3.1 Analytical approachesThe Poisson-Boltzmann (PB) formalism [22, 23] is one of the most widely usedapproaches to describe electrostatic systems at finite temperature. It is able to re-produce, explain and predict many experimental observations regarding the struc-ture of the counterion cloud around charged objects. However, the experimentaland numerical results for attractive interactions between like-charged polyelec-trolytes can not be explained with PB as it is a mean-field theory [24, 25]. In PBtheory, two parallel polyelectrolytes with the same line charge always repel eachother [26, 27]. A brief summary of this theory is presented in Chapter 2.1. Othermethods are therefore needed to investigate like-charge attraction such as integralequations [28], field theoretic methods [29], and computer simulations [30–33].These studies have revealed that strong electrostatic interactions between multiva-lent counterions and polyions result in a condensed layer of counterions around themacromolecules, and the dynamic correlation between condensed counterions isthe origin of attractive interactions between polyelectrolytes [34–41].There are several theoretical models which have been proposed to explain theorigin of DNA attraction [42]. One of them [37] models DNA as a uniformlycharged cylinder as shown in Fig. 1.2(a) and neglects structural features. Counte-rions are assumed to form a Wigner crystal on a background of a bundle of rods.This model predicts that DNA condensation occurs in an electrolyte for ions ofvalence ≥ 3. The electrostatic zipper model is another theory proposed by Korny-shev et al. [43]. It assumes that counterions bind to the DNA grooves and thereforethe positively charged grooves of one DNA molecule interlock with the negativelycharged helices of the other DNA strand, which results in an effective attractionbetween two DNA molecules. The assumption of this theory is not supported byexperiments which show that counterions (valence< 3) are more attracted towardthe DNA backbone and not to the DNA grooves [44]. Ha et al. present another the-oretical approach for the counterion-mediated interactions between like chargedrigid polymers [35]. They considered DNA as a uniformly charged cylinder anddivided the counterions into two classes, condensed and free. They found that the5charge fluctuation of condensed ions along the length of the rods results in attrac-tive interaction between DNA molecules.1.3.2 Implicit solvent simulations of cylindrical modelsGrønbech-Jensen et al. examined the nature of electrostatic interaction betweentwo rigid polyelectrolytes by performing Brownian dynamics (BD) simulation withthe cylindrical model (CM) of Fig. 1.2(b) [30]. Water molecules are not representedexplicitly, instead their effect is included in a dielectric continuum approximation(sometimes referred to as the primitive model of electrolytes). Since the DNAstrands are assumed to be rigid, only the equations of motion for the mobile coun-terions need to be integrated. They computed the thermally averaged force as afunction of DNA-DNA separation for different counterion radii and DNA chargedensities. They found that an attractive interaction appears for higher rod ion valen-cies (z≥ 2) and the onset of the attractive regime highly depends on the counterionradius. This study reports a maximum attractive force per helical turn of 0.008nNat a DNA center-to-center separation of 30A˚.Deserno et al. performed molecular dynamics (MD) simulation to investigatethe origin of attraction between charged stiff polyelectrolytes in the presence oftrivalent counterions and in the absence of additional salt. Their study also em-ploys a dielectric continuum approximation for the solvent [33] and used a cylin-drical model of the type shown in Fig. 1.2(b), but in a bulk configuration of parallelcharged rods forming a triangular lattice. A negative osmotic pressure indicatespresence of an attractive regime within a certain density range of counterions. Theyfound that the electrostatic contribution (always negative) shows no particular fea-tures in that specific range of counterion density. However, there is a drop in therepulsive excluded-volume contribution such that the electrostatic contribution canovercome the short-range repulsion. Therefore, the origin of the negative osmoticcoefficient is identified as the reduction of short-range repulsions. They also inves-tigated the interaction between a pair of rods as a function of Bjerrum length, whichis the distance at which the electrostatic energy between two elementary chargesequals the thermal energy kBT . They found that at high Bjerrum length, the to-tal force is attractive although the electrostatic contribution can become repulsive.6They also confirmed the presence of a mutually interlocking pattern of counterionsaround the rods.Stevens used MD simulations to investigate the conditions under which bundlingof polyelectrolytes occurs in solution and in the presence of counterions for a sim-ple bead-spring model of semiflexible polyelectrolytes [31, 32]. His simulationsshow that there is no aggregation in the presence of monovalent ions. Bundlingappears in the presence of divalent counterions; however, the chains must be suf-ficiently long and stiff. He showed that the condensate is stable only for triva-lent and tetravalent counterions. The strong electrostatic interactions dominate en-tropic contributions. The semiflexible polyelectrolytes condense into toroidal orrod structures depending on the stiffness of the polymer. Condensation was foundto be independent of DNA chemistry and sequence [32].1.3.3 Implicit solvent simulations of helical modelsIn a series of papers, Allahyarov et al. investigated the effect of charge pattern ofDNA molecules on their mutual electrostatic interactions, all within the primitivemodel of electrolytes for computational efficiency [38, 45, 46]. Using the extendcylindrical (ECM) model Fig. 1.2(c), they calculated force and torque on each DNAmolecule for different salt concentrations. They found that the effective force de-pends on the relative orientation of the DNA molecules and can be both repulsiveand attractive for center-to-center DNA separations of less than R = 26A˚. Thiswork also shows that the pattern of condensed counterions depends on the geome-try of the DNA molecule by considering different models for the DNA molecules[45]. Besides CM and ECM, they also considered the helical model shown inFig. 1.2(d) which has been termed Montoro-Abascal model (MAM). While thecounterions bind mostly on the phosphate strands in the ECM, they move towardsthe minor grooves in the MAM. Increasing the concentration of counterions alsochanges the distribution of counterions around the DNA molecules. The results ofanother set of simulations show that the interaction between DNA molecules in alldifferent geometries is repulsive in the presence of the monovalent counterions andsalt [46]. For the divalent and multivalent counterions there is an attractive regimewhose strength depends on the valency of the counterions. Typical maximum val-7ues of the attractive force per helical turn are 0.012nN at a DNA center-to-centerseparation of 24A˚. For the multivalent counterions and monovalent salt ions, theattractive regime disappears when increasing the salt concentration. In the systemwith tetravalent counterions and monovalent salt, there is a double-minimum struc-ture in the force distance curve. The maximum attraction per helical turn of 0.03nNappears at R = 30A˚ [47].1.3.4 Fully atomistic modelsMore recently, two groups have presented all-atom simulations, in which two DNAstrands are represented in full molecular detail (see Fig. 1.2) and solvated withexplicit water molecules (see Chapter 2 for details). Such simulations are far morecomputationally demanding as the equations of motion for DNA, counterions andwater molecules must be integrated. The study by Dai et al. was able to explainthe attractive force in terms of the formation of transient ion bridges [48]. Theyinvestigated the interaction between two double stranded DNA in the presence ofpolyamines (Putrescine (Pu), spermidine (Sd), and spermine (Sm)), using umbrellasampling. Their results show a weak attraction regime in the presence of Pu whichis a linear polyamine with two cationic nitrogen charges located at the terminalends. The depth of potential is less than 2kBT and it occurs at R= 24A˚. This resultis consistent with the experimental observation that divalent counterions cannotinduce condensation [4] and the experimentally observed weak DNA attraction inthe presence of divalent magnesium ions [49].Luan et al. also report direct observation of DNA attraction from all-atom MDsimulations in the presence of monovalent and divalent electrolytes [50]. Theirsimulations show the existence of a very weak attractive force in the presence ofNaCl at the DNA separation at R = 25A˚. However the attractive regime is muchstronger in the MgCl2 electrolyte which is in agreement with X-ray diffraction ex-periments [49]. They reported maximum attractive force of ∼ 0.042nN per helicalturn at DNA-DNA separation of R = 26A˚ in the presence of MgCl2 electrolyte.Their simulations suggest that two DNA molecules can attract in high concentra-tion of NaCl and MgCl2 electrolytes.81.3.5 Experimental resultsEarly experimental studies indicate that DNA condensation in vitro occurs in thepresence of counterions with the valence, z≥ 3 (e.g. polyamines spermidine3+ andspermine4+ and the inorganic cation Co(NH3)3+6 and basic histones) [4, 5, 51].However, synchrotron X-ray diffraction experiments by Koltover et al. show thatλ -phage DNA molecules confined to a two-dimensional cationic surface can at-tract each other in the presence of divalent counterions. Thus, the presence of di-valent ions facilitates the parallel alignment of DNA fragments [52]. Qiu et al. alsopresent experimental evidence for attraction between multiple DNA molecules inthe presence of Mg2+ counterions [49]. By increasing the counterion concentrationthe inter-DNA interaction changed from repulsive to attractive.1.4 Challenges for molecular simulations of chargedbiomoleculesIn all simulation studies described above, approximations have been made but it isnot obvious how these simplifications influence the results. With the exception ofthe all-atom simulations, the primitive model of electrolytes has been used, wherethe solvent is described as a continuum dielectric with uniform relative dielectricconstant εr ≈ 78. In the primitive model, the nonlocal and spatially inhomogeneousdielectric behaviour of water is completely ignored. In reality, however, water isa structured polar fluid and one expects this discreteness to become important onthe nanoscale. Indeed, the potential of mean force between ions is no longer de-scribed by the continuum approximation of Coulomb’s law between point chargeson length scales less than∼ 1nm, but shows multiple oscillations and even changessign for some ion separations [53]. One possible interpretation of this behaviouris in form of a wavevector-dependent dielectric response function of water, whichcauses an oscillatory behaviour in the interaction at short distances [54, 55]. Theorigin of this behaviour lies in the discrete nature of the solvent, which restructuresaround the ions in the form of solvation shells [56]. Atomistic simulations witha chemically detailed representation of water molecules automatically take theseeffects into account, but remain a formidable computational challenge for all butthe smallest systems.9For this reason, there is currently a great effort underway to develop coarse-graining strategies that preserve essential nanoscale aspects of solvation of bio-molecules while offering greater computational efficiency. Coarse-graining re-quires a tradeoff between averaging out details at small length scales while pre-serving correlations on larger scales. Most coarse-graining approaches suggestedso far remove a certain number of degrees of freedom by grouping sets of atomsinto superatoms with renormalized pairwise interactions [57, 58]. These effectivemesoscale potentials can be derived systematically from atomistic simulations byperforming a Boltzmann inversion of the pair distribution function, yielding a po-tential of mean force (PMF) [57]. Reith et al. then proposed to refine these setsof potentials iteratively and applied the strategy successfully to polymer systemssuch as polyisoprene [58]. In this so-called structure based coarse-graining, pairpotentials are parameterized such that the radial distribution functions (rdfs) of theCG model converge to the target functions of the full atomistic simulation or ex-perimental results.Lyubartsev et al. [59] derived effective potentials for ion-DNA interactionsfrom molecular dynamics simulations data by applying an inverse Monte Carlomethod [60]. They applied the derived effective potentials to simulate a consid-erably large DNA fragment and they could recalculate the ion distributions andthe relative binding affinities of different ions to DNA with good agreement withexperimental data.In the context of ds-DNA solvated by water, the nanoscale effects of ion sol-vation need to be described in an implicit way. DeMille et al. developed a coarse-grained model of NaCl in water by considering only short-ranged potentials [61].To avoid computationally challenging electrostatics, they took the somewhat rad-ical step and ignored long range electrostatic interactions entirely. Their modelwas successful in reproducing solvation structure, orientational and radial distribu-tion functions despite disregarding Coulombic interactions. They also developed acoarse grained model of DNA with explicit solvation by the ion-water model previ-ously introduced [62]. They apply a parametrization with the aim of replicating thestructure of ions and water around the DNA molecules in the full atomistic simula-tions. The results show that the CG model reproduces radial distribution functionswhich are very similar to the reference model.10In our study of counterion-mediated attraction between DNA molecules, it isnot reasonable to discard electrostatics entirely. A more meaningful approach is topreserve long range Coulombic interactions, but to include short range correctionsthat capture the effect of solvation shells on ion-ion interactions. Lenart et al. de-veloped a coarse-grained model for the interionic interaction potentials for LiCl,NaCl, and KCl in solution [63]. They added two analytic terms to the Coulombpotential between ions. The first correction term incorporates a distance depen-dent dielectric permittivity. In fact the distribution of water molecules aroundions is not homogeneous in the system and the uniform dielectric model is nota good approach for the phenomena which appear at short length scales. The sec-ond correction term incorporates the presence of one or more repulsive hydrationcorrection terms, which arise from water molecules rearranging in shells aroundthe charged particles [64–67]. By parametrizing these corrections, the model wasshown to capture the structure and thermodynamics properties of the correspondingfull atomistic model. An alternative approach would be to directly use a numeri-cally tabulated potential of mean force extracted from atomistic simulations of ionpairs in the dilute limit [68].The method described above was applied by Freeman et al. [69] to a coarse-grained DNA model, in which each nucleotide is represented by 3 beads (3SPN)developed by Knotts et al. [70]. Their work improved previous coarse grainedmodels for DNA [62, 71, 72] such that there is a good agreement with experiment,for melting behaviour. In addition, their model captures the local ion structure inthe vicinity of DNA observed in atomistic simulations.1.5 Thesis outlineIn this thesis, we use molecular dynamics (MD) simulations in order to system-atically investigate solvation effects on the interaction between DNA molecules.Our first goal will be to critically evaluate the application of the primitive elec-trolyte model in simulations of counterion induced attraction in models of the typeshown in Fig. 1.2(b-d). To this end, we shall carefully compare results from im-plicit solvent, continuum electrostatic simulations, with those where the solvent isrepresented by explicit water molecules. Technical details of the MD simulations11and the water models are described in Chapter 2.In a first series of calculations, we obtain the average force between two chargedDNA strands represented by the cylindrical model (CM) in simulations that differonly by their representation of water as an implicit or explicit solvent, but haveotherwise identical parameters. In this way, the effect of the presence of watermolecules on the correlation of counterions and the behaviour of the interactionbetween the rods can be revealed. A previous study of the force between like-charged walls that treated solvation effects with integral equation theory alreadyindicated significant deviations from the dielectric continuum approximation [73].We consider two cases: thin rods, where the diameter of the rods is comparableto that of the ions, and thick rods, where the diameter is six times larger. The ge-ometry of the thick rods is representative of an assembly of two parallel ds-DNAstrands. We show in Chapter 3 that in both cases, an explicit representation of thewater alters the results of implicit simulations both qualitatively and quantitatively.In the next step, we want to obtain a more detailed understanding of the physi-cal mechanism of DNA condensation in different DNA structures. Hence we con-sider in Chapter 4 two different helical models (ECM and MAM) for the DNAstructure and study the effects of explicit solvent on the ion-mediated DNA inter-actions. As in the case of the CM, the behaviour of the interactions between twoDNA molecules and the structure of the surrounding counterions in the implicitsolvent models for both DNA structures are significantly different from a formallyequivalent explicit solvent simulation.Having established the importance of the atomistic structure of the solvent onDNA-DNA interactions, we proceed in Chapter 5 to develop a coarse-grained rep-resentation for ds-DNA solvated by water molecules that can replace the primitiveelectrolyte approximation with minimal computational overhead. It is necessary todevelop such models for investigating longer time and length scales, where explicitsolvent simulations are infeasible with typical computational resources. Since pre-vious work on coarse graining strategies successfully captures aspects of the DNA-counterion structure, we adopt a similar approach as a starting point. We use atwo-step approach, where we obtain an initial guess for short range corrections tothe Coulomb interaction from atomistic simulations, and then apply the iterativeBoltzmann inversion method for refinement. The inter-DNA forces and counterion12distributions generated by the newly defined potentials are a very good match tothose observed in the explicit solvent model. This result proves that it is possibleto capture complex multibody interactions between polyelectrolyte strands withrelatively simple pairwise potentials.As a first demonstration of the potential of the coarse-grained model, we con-sider in Chapter 6 a three-DNA-system in a triangular configuration. Simulationsshow the emergence of attractive interactions for multiple rods, but the force is notsimply additive.Chapter 7 presents a concise summary of the thesis’ findings.13Chapter 2Polyelectrolytes in solution:theory and simulations2.1 Poisson-Boltzmann theoryThe Poisson-Boltzmann model is one of the theories describing the electrostaticpotential and equilibrium distribution of ions around molecules in solution. Asa mean-field theory, the PB theory relies on the following assumptions: i) theCoulomb interaction is the only interaction to be considered between charged par-ticles. By considering particles as point charges, it neglects all short-range non-electrostatic interactions. ii) The solution is represented as a continuous dielectricmedium characterized with a dielectric constant. The relative permittivity is takento be εr = constant for water. iii) It neglects all the interactions between groups ofthree, four, five, etc. ions and considers only the interaction between pairs. Thisapproximation is only valid at very dilute concentration in which the number ofclusters of more than two units is negligible. For homogeneous systems with con-stant dielectric in the presence of a charge density ρ(r); the electrostatic potentialφ(r) is described by the Poisson equation:∇2φ(r) =−ρ(r)/ε, (2.1)14where ε is the dielectric coefficient of the medium. The charge density profile ofall ions ρ(r) is a mean-field continuous function of the position r and is assumedto follow the Boltzmann distribution,ρ(r) = ρ0e−zeφ(r)/kBT , (2.2)where ρ0 is the reference density of ions in the absence of interactions betweenions in the system, φ → 0. z is the valency of the ion and z > 0 for cations andz < 0 for anions. Defining n as the number density (per unit volume) of the eachspecies, then ρ(r) = ezn(r) , Substituting Eq. 2.2 into Eq. 2.1, we get the Poisson-Boltzmann (PB) equation for the potential φ(r),∇2φ(r) =−∑i4pien0i ziε e−eziφ(r)/kBT . (2.3)The PB equation constitutes one of the most useful analytical approaches tostudy electrostatic effects in solution. The equation is non-linear and it has an-alytical solutions for limited number of problems with simple charged boundaryconditions. However the equation can be solved numerically or within some fur-ther approximations. Approximations impose limitations on the validity of thesolution of the PB equation. The PB equation can describe the distribution of ionsaround charged surfaces as long as surfaces are not highly charged. For a systemin contact with an infinite reservoir of symmetric monovalent electrolyte, in whichn0+ = n0− = n0, the PB equation reads,∇2φ(r) = 8pien0ε sinh[eφ(r)kBT]. (2.4)This equation can be linearized under the assumption that the potential energy issmall which results in the famous Debye-Hu¨ckel (DH) theory [22]. The lineariza-tion of the equation by considering βφ(r) 1 is equivalent to assuming that thesolution is dilute enough to produce a low potential by increasing the mean distancebetween the particles.∇2φ(r)' 8pie2n0εkBTφ(r) = λ 2Dφ(r) = φ(r)/l2D. (2.5)15The new parameter, λD is the so called Debye-Hu¨ckel screening constant whichdescribes the exponential decay of the potential in the solvent. lD is called theDebye length which varies from about 3A˚ in a strong ionic strength of 1M ofNaCl to about 1µm in pure water [74]. The DH theory states that the interactionbetween pairs of ions is screened by all other cations and anions in the system. Theapproximation says that for r≤ lD, the Coulombic interaction (1/r) is only slightlyscreened but for r > lD, it is exponentially screened.φ ∼ r−1exp(−λDr). (2.6)The PB equation has been solved for various simple geometries. For the caseof two charged rods as a simple model for the two double stranded DNA, the PBequation has been solved analytically [26] and numerically [27]. A detailed ana-lytical study by Brenner et al. [26] presents the repulsive electrostatic interactionbetween two parallel cylinders immersed in an ionic medium. The equations forthe force and energy of interaction has been derived using the linear superposi-tion method for calculating the electrostatic potential outside of the cylinders. Theresults for the repulsive force were found to beF(r) =ξe−λDr√λDr, (2.7)where r is the interaxis distance and ξ is a parameter which only depends on theionic strength of the medium and the pH of the reservoir. ξ can be determinedby a calculation of the force at a single value of r because it is independent of theinteraxis distance. Values of ξ have been calculated for a variety of pH and ionicstrength values and the results show that a repulsion exists between the cylinderfor all ionic strength and pH values considered. The approximation of linearizedPB theory is only valid for low surface charge densities. This approximation is notappropriate for DNA with high charge density (∼−e/1.7A˚) and a radius (∼ 10A˚)comparable to the Debye length (9.61A˚≤ lD ≤ 30.4A˚ for the ionic strengths in therange, 0.01−0.1 mole/liter ). The numerical solution of the nonlinear PB equationfor the interaction between two charged cylinders also shows that the force betweenequally charged cylinders of equal radii will always be repulsive[27].162.2 Molecular dynamics simulationsIn order to study polyelectrolytes beyond the mean-field approximation of Poisson-Boltzmann theory, we use MD within a canonical ensemble. MD simulations cre-ate long trajectories of a system in phase space. By considering ergodicity, one candefine an ensemble average as a time average. In the most common version, thetrajectories of molecules and atoms are determined by integrating Newton’s equa-tions of motion for a system of interacting particles which yields a microcanonicalaverage for the system.mir¨i =−∇U({ri}), (2.8)where mi are the masses of particles with coordinates ri and U is the energy func-tion. Unfortunately the microcanonical ensemble does not correspond to the condi-tions under which most experiments are carried out. In order to obtain a canonicalaverage, the system has to be coupled to a heat bath or new degrees of freedommust be introduced. The type of thermostat depends on the solvation model repre-senting water which is discussed in the next sections.2.2.1 Interaction potentialsTwo kinds of pairwise interactions are considered in our simulations. First, theLennard-Jones (LJ) potential is the most widely used choice in molecular dynam-ics. It represents a van der Waals interaction between two particles separated by adistance r and includes a repulsive part at short distances to prevent particle overlapVLJ(r) = 4εi j[(σi jr)12− (σi jr)6], (2.9)where εi j and σi j are the energy and distance parameters for the different LJ site-site interactions. εi j is the depth of the potential well and σi j is the distance atwhich the inter-particle interaction is zero. For the cross interactions, the Lorentz-Berthelot mixing rules, εi j =√εiiε j j and σi j = (σii+σ j j)/2 are used.Second, the Coulomb potential VC acts between charged particles which gives17rise to long-ranged interactions,VC(r) =∑i6= jqiq j4piε0εrr, (2.10)where qi is the charge of ion i, εo = 8.85×10−12F/m is the dielectric constant ofvacuum and εr is the relative dielectric permittivity of the medium. Electrostaticinteractions strongly influence the dynamical and structural properties of soft ma-terials. This influence is important if the magnitude of the Coulomb interaction iscomparable to kBT which is the typical energy associated with deformations andstructural degrees of freedom. Another way of saying the same is to introducea length for which the Coulombic energy between two unit charges is equal tothe thermal energy. This is so-called the Bjerrum length, lB = βe2/4piεoεr whereβ ≡ 1/kBT . The Bjerrum length, lB ≈ 7A˚ for a dielectric constant of water, εr = 80and room temperature, T = 300K. The Coulomb interaction between two chargesz1e and z2e can be written in term of this new length asVC(r) =z1z2lBβ r . (2.11)2.3 Solvation models2.3.1 Implicit solvent modelThe implicit solvent model is a method of treating water as a continuum dielectricwith uniform dielectric constant, εr = 80. In reality, however, water is a structuredpolar fluid. One might expect this approximation to work well on length scalesmuch larger than the diameter of a water molecule, 3 A˚.In order to obtain canonical averages, a Langevin thermostat is applied in theimplicit solvent model [75]. It is actually imposed by integrating a set of Langevinequationsmir¨i =−∇U({ri})−Γr˙i +ξi(t), (2.12)where mi are the masses of particles with coordinates ri, U is the energy function,18Γ is a friction coefficient and ξi is a stochastic white noise. In the Langevin equa-tion (2.12), the friction term Γr˙i is proportional to the velocity of a particle and isresponsible for a drag force. The noise term ξi(t) is responsible for random colli-sions. It is a force due to solvent atoms at a temperature T randomly bumping intothe particle. ξi(t) is a Gaussian noise source and its first and second moments aregiven by〈ξi(t)〉= 0 and 〈ξi(t).ξj(t)〉= 6kBTΓδi jδ (t− t ′). (2.13)Although the continuum representation of solvent significantly improves thecomputational speed, this is an approximate method with certain limitations andproblems. It is shown in chapter 3 and 4, how these simplifications affect theresults of the simulations.2.3.2 Explicit solvent modelThe interaction of biomolecules with the solvent can be investigated precisely byfully representing the water molecules in the system. Explicit solvent models relyon using thousands of discrete solvent molecules which results in slow equilibra-tion of the system. In computational biophysics, several different models havebeen proposed for water molecules in order to be represented in the explicit sol-vent model. In this thesis, we use the three-site SPC/E, the extended simple pointcharge model [76, 77]. It is characterized by three point charges, representingthree atoms of the water molecule. Oxygen also gets the Lennard-Jones parame-ters. More information about the geometry of this model is given in table 2.1 andin Fig. 2.1. The SPC/E model adds a polarization correction to the effective pairpotential which results in a better diffusion constant and density than SPC model.This model is widely used because of its simplicity and computational efficiency.In the explicit solvent model, the dynamic quantities and hydrodynamic in-teractions are important, therefore the Langevin thermostat is not an appropriatechoice and another thermostat has to be used. Here we use a Nose´-Hoover thermo-stat to keep the temperature around an average value. In the approach of Nose´, theheat bath is considered as an integral part of the system by adding an extra degreeof freedom [78]. For a system of N particles, with coordinates q′i, momentums p′i,19Property ParametersMass of O 15.9994 [g/mol]Mass of H 1.008[g/mol]Charge of O -0.8476 eCharge of H 0.4238 eσLJ,OO 3.166 [A˚ ]εLJ,OO 0.1553 [kcal/mol]σLJ,HH 0εLJ,HH 0rOH 1 [A˚ ]θHOH 109.47◦Table 2.1: Parameters of SPC/E modelmasses mi, potential energy U(q′) and s as an additional degree of freedom, theHamiltonian is given by,H =∑ip2i2ms2+12 ∑i j,i6= jU(q)+p2s2Q+gkT ln(s). (2.14)where g is equal to the total number of degrees of freedom of the system, Q isan effective mass associated to s and ps is the conjugate momentum of s. qi, piand t are the virtual coordinate, momentum and time respectively. These virtualvariables are related to the real variables q′i, p′i, s′ and t ′ by,q′i = qi, (2.15)p′i = pi/s, (2.16)s′ = s, (2.17)t ′ =∫ t dts. (2.18)Assuming that the Hamiltonian formalism can be applied to equation 2.14 withvirtual variables, the equations of motion are.dqidt=∂H∂pi= pi/mis2, (2.19)20Figure 2.1: Schematic of the SPC/E model.dpidt=−∂H∂qi=−∂U∂qi, (2.20)dsdt=∂H∂ ps= ps/Q, (2.21)dpsdt=−∂H∂ s =(∑ip2imis2−gkT)/s. (2.22)In terms of the real variables, these equations can be written as,dq′idt ′= sdqidt= p′i/mi, (2.23)dp′idt ′i=−∂U∂q′i− (s′p′s/Q)p′i, (2.24)211sds′dt ′=dsdt= s′p′s/Q, (2.25)d(s′p′s/Q)dt ′=sQdpsdt=(∑ip′2imi−gkT)/Q. (2.26)To simplify these equations, we can introduce the thermodynamics friction coeffi-cient ξ = s′p′s/Q. The equations of motion then become,q˙′i = p′i/mi, (2.27)p˙′i =−∂U∂q′i−ξp′i, (2.28)ξ˙ =(∑ip′2imi−gkT)/Q, (2.29)s˙/s =d lnsdt= ξ . (2.30)The equations of motion of the Nose´-Hoover scheme lead to a canonical distri-bution.2.4 Velocity VerletThe Velocity Verlet algorithm is a numerical method to integrate the equation ofmotion to determine the trajectory in the molecular dynamics simulations [79].Given the initial position x(t) and velocity v(t), the algorithm calculates positionof each particle at time t+∆t,x(t+∆t) = x(t)+ v(t)∆t+f (t)2m∆t2 (2.31)where ∆t is a time increment which is taken equal to 2 f s. f (t) is the total forceacting on the particle at time t and m is the particle mass. Given the initial posi-tion x(t) and velocity v(t), the algorithm calculates x(t+∆t) by applying equation(2.31). In the next step, it recalculates the force value on the particle by the updatedvalue of the particle position. Finally, it updates the velocity of particle with thenew force value, in the following way,22v(t+∆t) = v(t)+(f (t)+ f (t+∆t)2m)∆t (2.32)All simulations in this thesis have been carried out using the LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) from Sandia National Lab-oratories [80]. Ewald and particle-particle particle-mesh (PPPM) methods are ap-plied for the evaluation of the Coulomb sum in the implicit and explicit solventmodel respectively, and details are given in the appendix A.23Chapter 3Cylindrical model3.1 IntroductionIn this Chapter, we use molecular dynamics (MD) simulations in order to system-atically investigate the role of solvation effects on the interaction between like-charged cylindrical rods as simplified model for DNA molecules in the presenceof divalent counterions. We obtain the average force between two charged rodsin simulations that differ only by their representation of water as an implicit orexplicit solvent, but have otherwise identical parameters. In this way, the effectof the presence of water molecules on the correlation of counterions and the be-haviour of the interaction between the rods can be revealed. A previous study ofthe force between like-charged walls that treated solvation effects with integralequation theory already indicated significant deviations from the dielectric contin-uum approximation[73]. We consider two cases: thin rods, where the diameter ofthe rods is comparable to that of the ions, and thick rods, where the diameter issix times larger. The geometry of the thick rods is representative of an assemblyof two parallel ds-DNA strands. We show that in both cases, an explicit represen-tation of the water alters the results of implicit simulations both qualitatively andquantitatively. 11A version of this Chapter has been published as ” Solvation effects on like-charge attraction”, J.Chem. Phys. 138, 084901 (2013) [81].243.2 Models and simulationsAll simulations consist of two immobile rods composed of 20 ions per rod, whichare placed parallel to the z-axis and separated along the x-axis of a rectangularsimulation box. Rigid rods are modeled as chains of discrete negatively chargedspheres with two different LJ diameters, σ = 3.17A˚ and σ = 20A˚ for thin andthick rod (ds-DNA-like) systems, respectively. The simulation setup is shown inFig. 3.1. In both systems, the rod ions are separated along the rod axis by b =1.7A˚, which corresponds to the vertical separation of phosphate charges on a ds-DNA strand. Periodic boundary conditions make a set of infinite rods along the z-direction, and the x-dimension is chosen so that it is at least three times larger thanthe largest rod separation considered. The total charge of the rods is compensatedby adding divalent counterions in some of the systems as depicted in Fig. 3.1.In the implicit solvent model, water is treated as a continuum dielectric while inthe explicit solvent model, water molecules are fully represented using the SPC/Emodel [76, 77]. The size of the counterions is the same as that of the diameter ofthe water molecules, σ = 3.17A˚.As discussed in Chapter 2, particles interact via two basic interactions, theLennard-Jones potential (LJ), Eq. 2.9 and the Coulomb potential, Eq. 2.10. TheLJ energy scale ε = 0.155 kcal/mole is employed between all pairs of particlesto ensure self-avoidance. The LJ length scale σ is considered as the diameterof the particles, and the potential is truncated and shifted at rc = 10A˚ for all ofthe interactions (rc = 23.16A˚ for rod ion-counterion interactions in simulationswith thick rods). The parameters for the LJ interactions are given in Table 3.1.Both sets of simulations are conducted at the same temperature or, equivalently,value of the Bjerrum length lB = βe2/4piεoεr, which is the distance at which theCoulomb energy of two unit charges is equal to 1/β = kBT . We ignore Lennard-Jones interactions between rod ions to focus on the electrostatic contribution. Wemake use of Ewald and PPPM methods for the evaluation of the Coulomb sum inthe implicit and explicit solvent models respectively. In order to obtain a canonicalaverage, we apply a Langevin thermostat in the implicit model and a Nose´-Hooverthermostat in the explicit model to keep the temperature at T = 300K, so that lB ≈7A˚. With the above values the Manning parameter ξ = 2lB/b ≈ 8, which places25Figure 3.1: Schematic depiction of the simulation geometry. (a) Two chargedpolyelectrolytes are modeled as rigid rods, which consist of rod ionsseparated by b = 1.7A˚. Each rod ion is modeled as a LJ sphere withdiameter σ1 and charge q1 =−e at the center of sphere. Each counterionis modeled as a LJ sphere with diameter σ2 and charge q2 = 2e at thecenter of sphere. (b) Thin rods system, where rod ions and counterionshave the same size, σ1 = σ2 = 3.17A˚. (c) Thick rods system (DNAlike), where the size of rod ions is almost six times greater than the sizeof counterions, σ1 = 20A˚ and σ2 = 3.166A˚.the systems in a regime where the counterions strongly condense onto the rods[82].This parameter can also be interpreted as the ratio of the electrostatic energy to thethermal energy, kBT . All the properties of the systems are shown in Table 3.2.3.2.1 Implicit solvent modelWe start with a system in which the surface-to-surface distance of two rods is zeroand run the simulation for 10 ns in order to reach an equilibrium state in which thecounterions have small fluctuations around the rods. We consider this state as theinitial configuration, increase the separation between the rods and perform simu-26Pair ε[kcal/mole] σ [A˚] rc[A˚]RI-CI 0.1553 3.17 (11.58) 10 (23.16)RI-O 0.1553 3.17 (11.58) 10 (23.16)O-O 0.1553 3.17 10O-CI 0.1553 3.17 10Table 3.1: Parameters for the Lennard-Jones interaction of different atoms( rod ion(RI), counter ion (CI) and oxygen (O)). There is no Lennard-Jones interaction between Hydrogen and other particles in the system.The potential is truncated at rc.PropertyBjerrum length [A˚] 7Charge separation [A˚] 1.7Manning parameter 8Counterion valency 2Number of counterions 20Rod diameter [A˚] 3.17( 20 )Temperature [K] 300Table 3.2: Properties of the systems.lations for 10 ns at each rod separation while measuring the thermally averagedforce 〈F〉 in the x-direction that one rod exerts on the other. A blocking method(see appendix B) was used to ensure that the error bars are smaller than the symbolsizes [83]. We also ensure that the average force on each of the two rods is obeyingNewton’s 3rd law with a percentage error of less than 1%. Moreover, animations ofthe MD trajectories show that all counterions have condensed onto or near the rodsfor the entire simulation time. Panel (a) and (b) of Fig. 3.2 show our simulationresults for the mean force per rod ion in the thin and thick rods systems with animplicit solvent. By comparing the results of the simulations with the force on therod in the absence of counterions, we can first see that the electrostatic interactionbetween rod ions has been screened by the counterions. In both systems, an at-270 1 2 3 4 5 6 7 8−0.100.1/Ion [nN] (a)d [Å]0 1 2 3 4 5 6 7 8−0.0200.02/Ion [nN] (b)d [Å]0 1 2 3 4 5 6 7 8−0.8−0.400.40.8 (c)d [Å]/Figure 3.2: Average force between (a) thin and (b) thick rods system as afunction of rod separation (surface to surface) in implicit solvent simu-lations. The crosses and heavy dots represent force in implicit systemwith and without counterions, respectively. (c) Ratio of forces with andwithout counterions shown in panels (a) and (b). The crosses and opencircles show this ratio for thick and thin rods system, resp.tractive regime (negative force) appears for specific rod separations which agreeswith previous implicit solvent studies by Grønbech-Jensen et al. [30] and Desernoet al. [33]. We note, however, that the magnitude of the force in the thick rodsystem is about 6 times weaker since the absolute center to center rod separationhas increased by the same amount. In order to compare the results of two systems280 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6−50510d [Å]/Ion [nN] (a)0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6−1−0.500.5 (b)d [Å]/Figure 3.3: (a) Force in three different models of thin rods as a function of rodseparation (surface to surface). Squares show the force in an explicitmodel with counterions, crosses show the force in an explicit modelwithout counterions and dots show the force in an implicit model withcounterions. (b) Ratio of the forces in explicit solvent simulations withand without counterions.on the same scale, we show in panel (c) the ratio of the force with counterions andforce without counterions. It shows that for both systems, the attractive regimestarts when the rod separation is roughly equal to the diameter of counterions. Themaximum attraction occurs at rod separation, d = 3.7A˚ and d = 5.75A˚ in thin andthick rod systems, respectively and the attractive regime is wider in the thick rods290 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 600.20.4d [Å]/Ion [nN] (a)0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 600.511.5d [Å]/(b)Figure 3.4: (a) Force in three different models of thick rods system as a func-tion of rod separation (surface to surface). Squares show the force in anexplicit model with counterions, crosses show the force in an explicitmodel without counterions and dots show the force in an implicit modelwith counterions. (b) Ratio of the forces in explicit solvent simulationswith and without counterions.system.3.2.2 Explicit solvent modelWe now study the exact same systems in an explicit solvent model by filling thesimulation box with 840 and 5160 water molecules for the thin and thick rods30(a)(b)(c)Figure 3.5: Cross-sectional picture of distribution of counterions around therod ions (purple spheres) for explicit model of thin rods system with rodseparation, (a) d = 0.9A˚ and (b) d = 1.8A˚ and (c) d = 4A˚ .systems, respectively. The size of the box is determined by considering the densityof water to be approximately 1g/cm3. First, we study the interaction between twosolvated charged rods without counterions (green ×) in panels (a) of Figs. 3.3 and3.4. The force between the rods is always repulsive, but now shows some structuredue to the formation of solvation shells (see also Fig. 3.10). When compared to theforce experienced in the implicit solvent model (red •), we observe that the forcesare more than an order of magnitude greater, which suggests that the implicit modelstrongly overscreens at short distances.312 2.5 3 3.5 4 4.5 5 5.5 60204060r [Å]g(r)(a)2 2.5 3 3.5 4 4.5 5 5.5 60100200300 (b)r [Å]g(r)Figure 3.6: (a) Counterion-rod ion pair correlation functions in thin rodssystem with implicit solvent for three rod separations (surface to sur-face), d = 0.5,2,4A˚, represented by squares, heavy dots and crosses.(b) Counterion-rod ion pair correlation functions in thin rods systemswith explicit solvent for three rod separations (surface to surface),d = 0.9,1.8,4A˚ represented by squares, heavy dots and crosses.By adding divalent counterions to the thin rods system, an attractive interactionregime arises as in the implicit model (blue  in Fig. 3.3(a)). However, the force-separation curve now features three minima at d = 0.9A˚, d = 1.8A˚, and d = 4.0A˚,respectively. This confirms the role of counterions in the like-charged attractionphenomenon. At the first minimum, the magnitude of the attractive force is more3210 11 12 13 14 15 16 17 180102030r [Å]g(r)(a)10 11 12 13 14 15 16 17 1801020304050 (b)r [Å]g(r)Figure 3.7: (a) Counterion-rod ion pair correlation functions in thick rod sys-tems with implicit solvent for three rod separations (surface to sur-face), d = 2,4,6A˚ , represented by squares, heavy dots and crosses.(b) Counterion-rod ion pair correlation functions in thick rod systemswith explicit solvent for three rod separations (surface to surface),d = 1.3,1.8,4A˚ , represented by squares, heavy dots and crosses.than 80 times greater than in the implicit solvent model, since there is no spacefor water molecules to be placed directly between rod ions and counterions. Inrelation to the force without counterions, however, we find a comparable variationin amplitude (see Fig. 3.3(b)). Each of the three minima correspond to distinctcounterion configurations: at d = 0.9A˚ ions arrange in a quadrupolar pattern, at332 3 4 5 6 7 805101520r [Å]g(r)(a)2 3 4 5 6 7 80102030r [Å]g(r)(b)Figure 3.8: Water molecule-rod ion pair correlation function in explicit thinrods system (a) with and (b) without counterions. The results for threerod separations (surface to surface), d = 0.9,1.8,4A˚ are represented bysquares, heavy dots and crosses, respectively, for both systems.d = 1.8A˚ one layer of ions fits between the two rods, and at d = 4.0A˚ more thanone layer of counterions condenses between the rods, see Fig. 3.5.For the explicit solvent model of the thick rods system with counterions, theinteraction between two rod ions is represented by blue squares in Fig. 3.4(a). Wecan see that the force is repulsive for all rod separations, and the attractive regimewhich exists in the implicit solvent model disappears. In fact, the presence ofcounterions barely seems to affect the interaction between the solvated charged34rods. As a result, the ratio of the forces with and without counterions is always oforder one, see Fig. 3.4(b). Considering the present geometries for the rigid rods, inthe explicit solvent model, attraction is observed only for the thin rods system andfor the thick rods system, the force is always repulsive, see Figs. 3.3(a) and 3.4(a).3.2.3 Correlation functionsIn order to obtain more direct insight into the origin of the changes in behaviour,we investigate the counterion-rod ion and water molecule-rod ion pair correlationfunctions normal to the rods. The radial pair correlation function g(r) is the proba-bility density of finding the center of a particle (counterion or water molecule) at adistance r away from the center of one of the rods relative to the probability densityof finding a particle at the same distance in a random system with the same averagedensity. By inspecting the counterion-rod ion pair correlation functions, Figs. 3.6and 3.7, we find that almost all of the counterions condense in one layer very closeto the rods. The position of this layer depends on the model and geometry of thesystem. For all rod separations in the implicit model of thin rods system, counteri-ons condense at r≈ 3.2A˚, see Fig. 3.6(a). This corresponds to the surface to surfacecontact of counterions and rod ions. Increasing the separation between two rodsdecreases the probability density of condensation of counterions at this particularposition. In the explicit model of the thin rods system, counterions condense closerto the center of rod at r ≈ 2.3A˚ which is smaller than the LJ diameter of the ions,see Fig. 3.6(b). This shows that the unscreened electrostatic interactions betweenions overcome part of the Lennard-Jones interactions and allow the counterions toget close to the rod ions. For the rod separation d = 1.8A˚, where the maximum ofattraction occurs, there is another favorable position for counterions at r ≈ 2.5A˚.In the implicit model simulations of thick rods, the electric fields at the surfaceof the rods are much weaker than for the thin rods since the small free ions cannotenter the much larger excluded rod volume. Fig. 3.7(a) shows that the counteri-ons condense at a distance r ≈ 12A˚ away from the center and the counterion cloudis less localized. In the explicit solvent model, the condensed layer of counterionsmoves further away to r≈ 12.6A˚, see Fig. 3.7(b). This behaviour explains the qual-itative difference between thin and thick rod systems: the counterion correlations3510 11 12 13 14 15 16 170123 (a)r[Å]g(r)10 11 12 13 14 15 16 170123 (b)r[Å]g(r)Figure 3.9: Water molecule-rod ion pair correlation function in explicit thickrods system (a) with and (b) without counterions for three rod separa-tions (surface to surface), d = 0.2,1.8,4A˚ are represented by squares,heavy dots and crosses, respectively.are no longer strong enough and the attractive regime disappears. The ultimatecause of this behaviour is actually not the discrete nature of the solvent, but in factthe dielectric contrast between excluded rod volume and the regions accessible bythe water molecules. In the explicit model, the counterions see the rods as a di-electric with εr = 1. It is therefore energetically more favorable for the ions tobe solvated in water than getting very close to the rod surface. There is no suchdielectric contrast in the implicit solvent model.36Figure 3.10: Number density plot of water molecules in explicit solventmodel for two different thin rod separations (surface to surface), d =0A˚ for systems (a) with and (b) without counterions and d = 4A˚, forsystems (c) with and (d) without counterions.We also investigate the formation of solvent-shell configurations by calculatingwater-rod ion pair correlation functions. The position of the oxygen atom is takenas the center of the water molecule. In the thin rods system, r≈ 2.5A˚ is the distanceat which the first layer of water molecules appears, see Fig. 3.8(b). We can see inpanel (a) that adding counterions to the system does not change the position ofthe solvation shells but reduces the probability density of finding water moleculesin the shells. In the thick rods system, the first shell of water molecules occursat r ≈ 10.9A˚, see Fig. 3.9. The distribution of water molecules in this system isweakly affected by the presence of counterions compared to the thin rods system.37Number density plots of water molecules in the plane perpendicular to the rodaxes shown in Fig. 3.10 further illustrate the distributions of water molecules inthe thin rods system. For rod separation d = 0.2A˚, the screening effect of thecounterions removes the third solvation shell, see panels (a) and (b) of Fig. 3.10.For rod separation d = 4A˚, where water molecules and counterions can be placedbetween two rods, adding counterions smears out the localized positions of watermolecules, see panel (c) and (d) of Fig. 3.10.3.3 ConclusionsWe have investigated the role and importance of representing water as explicit sol-vent in molecular dynamics simulations of the like charge attraction effect betweentwo rigid charged rods. We can see from the results of both geometries that at shortdistances, the behaviour of interactions between rod ions in the explicit solventare not only qualitatively but also quantitatively different from implicit solventmodels[81]. Solvating thin rods, where the rod and ion LJ diameters are equal,with SPC/E water predicts multiple attractive minima with much larger forces andat much smaller rod separations than an equivalent implicit solvent simulation.These changes clearly arise from the structured nature of the solvation shells andthe lack of screening at short distances. In explicit solvent simulations of thick rodswith a diameter comparable to that of ds-DNA, we observed the disappearance ofthe attractive regime present in implicit solvent models. This effect is due to thedielectric contrast between the interior of the rods and the regions filled with water,which is naturally accounted for in atomistic simulations but not in implicit solventmodels with a uniform dielectric constant.The above results indicate that implicit solvent models with a continuum di-electric description are not appropriate for studying electrostatic phenomena whichoccur on the nanoscale. The possibilities to improve on the primitive solvent modelwill be explored in Chapter 5.38Chapter 4Helical models4.1 IntroductionThe goal of this chapter is twofold: first, we want to obtain a detailed understand-ing of the physical mechanism of DNA condensation in different DNA structures.In this chapter two different helical models for the DNA structure are considered.These helical models are the refinement of the cylindrical model (CM) which westudied in the chapter 3. In one of these models, the phosphate charge distribu-tion is described by discrete charges in a double helical array corresponding toB-DNA. This model is called extended cylindrical model (ECM) which is shownin Fig. 4.1 and was first designed by Lyubartsev et al [20]. The other helical struc-ture, which resembles the real DNA appearance, is called Montoro-Abascal model(MAM) [21] and is shown in Fig. 4.2. In this model, another neutral sphere isadded between the cylindrical core and the charged phosphate sphere to incorpo-rate the accessibility of the DNA grooves to the counterions. This model gives afairly realistic description of the molecular structure of B-DNA. Second, we aimto study nanoscale solvation effects of water molecules on the ion-mediated DNAinteractions.39Figure 4.1: Extended cylinder model. Phosphate charges are shown as bluespheres and the cylindrical core is coloured in yellow.4.2 DNA modelsThe system that we want to model consists of two DNA molecules which are sepa-rated from each other by a distance R, with mobile ions which are randomly placedaround them. In the ECM, each DNA is treated as a central cylinder and phos-phate groups situated at the sites corresponding to the B-form of DNA. The centralcylinder composed of 10 spheres of radius rCA = 9A˚ which are separated alongthe z-axis by b = 3.4A˚. There are two phosphate groups per central atom (basepair) which makes the number of phosphate charges per one helical turn, Np= 20.Each phosphate group has a charge q = −e and a radius rp = 2.1A˚. In the MAM,there is another neutral sphere between the cylindrical core and the charged phos-phate sphere. The cylindrical core in the MAM has a radius rCA = 3.9A˚, the stringof neutral spheres is located at a radial distance r = 5.9A˚, and the outer string ofphosphates is centred at a radial distance r = 8.9A˚. The spheres in both stringshave the same φ and z coordinates and radius rp = 2.1A˚.The width of the helices in one DNA molecule in the xy plane is given by φs,40Figure 4.2: Grooved or Montoro-Abascal model. Phosphate charges areshown as blue spheres. The cylindrical core is coloured in yellow andthe neutral particles as green spheres.see Fig. 4.3. φs is 144◦ and 154◦ for the ECM and MAM respectively. The mutualconfiguration of the two DNA molecules is given by two azimuthal angles φ , φ0which are shown in Fig. 4.3. φ determines the the rotation of the second DNAmolecule around its z-axis with regard to the first DNA molecule. φ0 is the anglebetween the phosphate charge and the DNA-DNA separation (centre to centre)vector R. In all figures hereafter we show the interaction forces starting from thedistances R= 22A˚ and for the azimuthal angles; φ0 = 0◦ and φ = 180 for the ECMand φ0 = 18◦ and φ = 0 for the MAM. The torque is zero for both DNA moleculesat these particular azimuthal angles.The dimensions of the box are Lx = 90A˚, Ly = 57A˚ and Lz = 34A˚which cor-responds to one full turn of B-DNA. Periodic boundary conditions make a set ofinfinite DNA molecules along the z-direction. In the implicit solvent model, wateris treated as a continuum dielectric (ε = 72) while in the explicit solvent model,water molecules are fully represented using the SPC/E model [76, 77]. The di-41Figure 4.3: A schematic picture explaining the positions of DNA moleculesand the definition of the different azimuthal angles φs, φ0 and φ .ameter of the divalent counterions is the same as that of the diameter of the watermolecules, σ = 3.17A˚. All particles interact via LJ and Coulomb interaction.Theparameters for the LJ interactions are given in the Table 4.1 for the both DNAmodels. For the cross interactions, the Lorentz-Berthelot rules, εi j =√εiiε j j andσi j = (σii+σ j j)/2 are used.4.3 Results4.3.1 Implicit solvent modelWe start with a system in which the centre-to-centre distance of two double strandedDNAs is R = 22A˚ and run the simulation for 10ns in order to reach an equilibriumstate. We consider this state as the initial configuration, increase the separationbetween the DNA strands. At each separation, simulations are performed for 10nsand 60ns for the ECM and the MAM respectively. We record the force on eachphosphate every 10−5ns and measure the thermally averaged force 〈F〉 in the x-42ECM MAMPair ε[kcal/mol] σ [A˚] rc[A˚] ε[kcal/mol] σ [A˚] rc[A˚]CA-CA 0 18 ... 0 7.8 ...CA-P 0 9 ... 0 8.9 ...CA-CI 0.1553 10.583 21.166 0.1553 5.483 10.966CA-O 0.1553 10.583 21.166 0.1553 5.483 10.966P-P 0 4.2 ... 0 4.2 ...P-CI 0.1553 3.683 10 0.1553 3.683 10P-O 0.1553 3.683 10 0.1553 3.683 10CI-CI 0.1553 3.166 10 0.1553 3.166 10O-O 0.1553 3.166 10 0.1553 3.166 10O-CI 0.1553 3.166 10 0.1553 3.166 10Table 4.1: Parameters for the Lennard-Jones interaction of different atoms(central atoms (CA), phosphate (P), counter ion (CI) and oxygen (O)).There is no Lennard-Jones interaction between Hydrogen and other par-ticles in the system. The potential is truncated at rc.direction in the presence of divalent counter ions. We ensure that the average forceon each of the two rods is obeying Newton’s 3rd law.Fig. 4.4 shows our simulation results for the reduced DNA-DNA interactionforce per helical turn, 〈F〉/F0 vs separation distance R in the ECM and MAM withan implicit solvent. The unit of the force is F0 = kBT/P, where P= 34A˚ is the DNApitch length. In the ECM, an attractive regime (negative force) appears for specificDNA strands separations which agrees with previous implicit solvent studies byAllahyarov et al. [38, 45–47]. It shows that for the ECM, the attractive regimestarts when the separation is equal to R= 22.5A˚ and the maximum attraction occursat separation, R = 24.25A˚. For the MAM with an implicit solvent, the interactionbetween two DNA strands is represented by filled squares in Fig. 4.4. We can seethat the force is repulsive for all separations. The interaction force depends onthe shape of DNA which is modeled differently in the ECM, and MAM. Differentbehaviour of the interaction between DNA strands could be due to the differentcounterion condensation patterns on the DNA surface in the ECM and MAM whichwill be investigated in the next sections by calculating pair correlation functions.4322 24 26 28 30−20−1001020−20R [Å]/F0 Figure 4.4: Reduced DNA-DNA interaction force 〈F〉/F0 vs separation dis-tance R with an implicit solvent in the presence of divalent counter ions.The circles and squares represent force for the ECM and MAM, respec-tively.4.3.2 Explicit solvent modelWe now study the exact same systems in an explicit solvent model by filling thesimulation box with 5093 and 5117 water molecules for the ECM and MAM re-spectively. The size of the box is determined by considering the density of water tobe approximately 1g/cm3. We start with the system in which the centre-to-centredistance of two DNAs equals to R= 22A˚ and let the system run for t = 10ns in or-der to reach an equilibrium state at room temperature T = 300K. We consider thisstate as the initial configuration and increase the separation between DNA strands.We use two different methods for the ECM and MAM in order to be sure that weare sampling the configurations correctly.In the ECM, we choose three different DNA strands separation, R = 23,25,28and for each separation, we heat up the system to T = 2000 K, let the systemrun for 10ns, and save the configuration every 2ns. Then all the configurations4422 23 24 25 26 27 28 29 30−1000100200R [Å]/F0 Figure 4.5: Reduced DNA-DNA interaction force 〈F〉/F0 as a function ofDNA separation (centre to centre) with an explicit solvent. The circlesand squares represent force for the ECM and MAM, respectively.are cooled down from T = 2000K to T = 300K in t = 0.01ns and run for t =5ns. For each DNA separation, we check the force between two DNA strands,the radial pair correlation function, the averaged number of counterions betweentwo DNA strands, the mean square displacement of counterions, the energy andthe pressure of the system for all different initial configurations. The results showthat by considering the statistical error, all the measured parameters are the same,which ensures that the sampling has been done correctly.In the MAM, a similar protocol is employed: the system is heated up to T =1000K for each separation, equilibrated for 1ns so counterion have the chance toexplore the entire box, then cooled down from T = 1000K to T = 300K in t =0.01ns. For each value of the inter-DNA distance, a 30ns simulation is performedto average out thermal fluctuations of the force. The force on each phosphate isrecorded every 10−5ns and we measure the thermally averaged force 〈F〉 in thex-direction. For both DNA models, we also ensure that the average force on each45of the two DNA strands is obeying Newton’s 3rd law with a percentage error ofless than 5%.Fig. 4.5 shows our simulation results for the reduced DNA-DNA interactionforce 〈F〉/F0 vs separation distance R in the ECM and MAM with an explicit sol-vent. In the ECM, the force is repulsive for all DNA separations less than 26A˚. Atthe separation, R= 25A˚, there is a weak attractive force between two DNA strands.In the MAM, the force-separation curve features a multiple minima structure withmaximum attractive forces at R ≈ 22.25A˚ and R ≈ 25.25A˚. The appearance of anattractive regime in the MAM at shorter distances compared to the ECM confirmsthe role of the accessibility of the DNA grooves to counterions through addingneutral atoms in the MAM.For the explicit solvent model of the ECM, the interaction between two DNAstrands is represented by blue heavy dots in Fig. 4.6 and shows some structure dueto the formation of solvation shells. When compared to the force experienced inthe implicit solvent model (red crosses), we observe that the forces are more thanan order of magnitude greater, which suggests that the implicit model stronglyoverscreens at short distances. In the explicit solvent model, the attractive regimeappears at larger distances compared to the implicit solvent model.Fig. 4.7 shows that the repulsive behaviour of the interaction between DNAstrands in the MAM with implicit model changes qualitatively and quantitativelyby explicitly considering the molecular structure of water. In the MAM with ex-plicit water, an attractive interaction regime arises at short distances.4.3.3 Correlation functionsWe study the origin of the different behaviour of the interactions by investigat-ing the counterion condensation pattern around the DNA strands. We calculate thecounterion-phosphates, counterion-central atom and counterion-counterion pair cor-relation functions. The pair correlation function g(r) is the probability density offinding the centre of a counterion at a distance r away from the center of a parti-cle (phosphate or central atom or counterion) relative to the probability density offinding a particle at the same distance in a random system with the same averagedensity.46Figure 4.6: Reduced DNA-DNA interaction force 〈F〉/F0 as a function ofDNA separation (centre to centre), considering extended cylinder model(ECM) for the DNA structure. The crosses and circles represent forcein explicit and implicit solvent model, respectively. The inset shows azoom for the reduced force with implicit model.Counterion-phosphate pair correlation functionsBy studying the counterion-phosphate pair correlation functions in the ECM withexplicit model, Fig. 4.8(a), we find that the counterions condense in two layers.The first condensed layer occurs at r = 3A˚ which is smaller than the LJ diameterof the ions, see Table 4.1. This shows that the bare Coulomb interactions betweenions overcomes part of the Lennard-Jones interactions and allows the counterionsto get close to the phosphate ions. This is also called the contact ion pair configu-ration. The second layer at r = 5A˚ shows the distance of the counterion from theadjacent phosphate along the strands. In the implicit model of the same system,almost all the counterions condense in one layer at r = 3.8A˚ for all separations,see Fig. 4.8(b). This distance is a bit greater than the surface to surface contact ofcounterion and phosphate ion. In fact, this distance is the optimum position for the47Figure 4.7: Reduced DNA-DNA interaction force 〈F〉/F0 as a function ofDNA separation (centre to centre), considering Montoro-Abascal model(MAM) for the DNA structure. The crosses and circles represent forcein implicit and explicit solvent model, respectively. The inset shows azoom for the reduced force with implicit model.counterion to be placed, considering the adjacent phosphates. Increasing the sepa-ration between DNA strands decreases the probability density of the condensationof counterions at this particular position.Fig. 4.9 shows the behaviour of the counterion-phosphate pair correlation func-tions for the MAM with implicit and explicit model. Considering the same LJ pa-rameters for the phosphate-counterion interaction results in similar behaviour forthe counterion-phosphate pair correlation functions in both DNA structures.Counterion-central atom pair correlation functionsFigs. 4.10 and 4.11 show the counterion-central atom pair correlation functions inthe ECM and MAM, respectively. In the explicit model simulations of both DNAstructures, counterions condense in more than one layer around the DNA. The firstpeak of the correlation function occurs at r = 12.2A˚ for the ECM and r = 10.7A˚482 3 4 5 6 7 8 9 10051015r [Å]g(r)(b) 2 3 4 5 6 7 8 9 1002040 (a)r [Å]g(r) Figure 4.8: Phosphate-counterion radial pair correlation function of the ex-tended cylinder model in (a) explicit and (b) implicit system for threeDNA separations (centre to centre), R = 23,25,28[A˚] , represented byheavy dots, squares and triangles.for the MAM, see Figs. 4.10(a) and 4.11(a). However, there is only one condensedlayer of counterions in the implicit model. The counterions condense at a distancer = 12.2A˚ and r = 8.7A˚ away from the central atoms for the ECM and MAM, seeFigs. 4.10(b) and 4.11(b). The position of the condensed layers depends on thesolvent model and geometry of the system. In the explicit solvent model of theMAM, the condensed layer of counterions moves further away to r = 10.7A˚, seeFig. 4.11(a) . The cause of this shift is the dielectric contrast between excludedcentral cylinder volume and the regions accessible by the water molecules. In theexplicit model, the counterions see the central atoms as a dielectric with εr = 1.It is therefore energetically more favourable for the ions to be solvated in waterthan getting very close to the DNA surface. The accessibility of the DNA groovesto the counterions causes the counterions to be condensed closer to the centralatoms in the MAM compared to the ECM for both solvent models. For both DNAstructures, the counterion clouds are more localized in the explicit model compared492.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5 10051015 (b)r [Å]g(r) 2.5 3 3.5 4 4.5 5 5.5 6 6.5 7 7.5 8 8.5 9 9.5 1002040 (a)r [Å]g(r) Figure 4.9: Phosphate-counterion radial pair correlation function of theMontoro-Abascal model in (a) explicit and (b) implicit system for threeDNA separations (centre to centre), R= 22.5,25.25,28[A˚] , representedby heavy dots, squares and triangles.to the implicit model.In the ECM, increasing the separation between DNA strands does not changethe position of the peaks for both solvent models. However, in the explicit modelof the MAM, the peaks move slightly closer to the central atoms as the separationbetween DNA strands increases. The separation between the first and third peakin the explicit model of the ECM and the first and second peak in in the explicitmodel of the MAM is almost equal to the diameter of one water molecule.Counterion-counterion pair correlation functionsThe presence of explicit water molecules changes the counterion-counterion cor-relation function in both DNA structures, see Figs. 4.12 and 4.13. Panel (a) ofFig. 4.12 shows that in the explicit model of the ECM, there are two favoured sep-aration distances for the counterions at r = 6.5A˚ and r = 8.5A˚. The probabilitydensity of finding counterions at these two particular positions decreases by in-509 10 11 12 13 14 15 160510r [Å]g(r)(b) 9 10 11 12 13 14 15 16051015r [Å]g(r) (a)Figure 4.10: Central atom-counterion pair correlation functions of the ex-tended cylinder model in (a) explicit and (b) implicit system for threeDNA separations (centre to centre), R = 23,25,28[A˚] , represented byheavy dots, squares and triangles.creasing the separation between two DNA strands. In the implicit model of theECM, there is a broad peak in the counterion-counterion pair correlation functionfor each DNA separation. By increasing the separation between DNA strands, thebroad peak moves away, see panel (b) of Fig. 4.12.In the explicit model of the MAM, for all the DNA separations, there is a pre-ferred separation distance for the counterions at r = 7A˚. For the specific DNA sep-aration R = 22.5A˚ in which the maximum attraction appears in the MAM (explicitmodel), there is another preferred position at r = 5.5A˚, see panel (a) of Fig. 4.13.The behaviour of the correlation function is totally different in the implicit modelof the MAM compared to the explicit model, panel (b) of Fig. 4.13. There is nolocalized peak in the pair correlation functions of the implicit model.514 6 8 10 12 14 160246r [Å]g(r)(a)4 6 8 10 12 14 160246 (b)r [Å]g(r)Figure 4.11: Central atom-counterion radial pair correlation function of theMontoro-Abascal model for three DNA separations (centre to centre),R = 22.5,25.25,28[A˚] , represented by heavy dots, squares and trian-gles in the explicit solvent model and circles, pluses and crosses in theimplicit solvent model.4.3.4 ConclusionsWe have investigated the interaction forces between a pair of DNA molecules in thepresence of divalent ions using molecular dynamics simulations. The importanceof the DNA geometry on the electrostatic forces has been studied by consideringtwo models for the DNA structure in our simulations. We have also studied theimportance of representing water as explicit solvent in the simulations of the likecharge attraction effect between two DNA strands. The simulation results indicatethat the DNA shape is an essential contributor to the interaction forces. Consider-ing helical models, we have found attractive regimes for the explicit and implicitmodel of ECM and the explicit model of the MAM. The results of different paircorrelation function show different counterion condensation patterns on the DNAsurface in two DNA structures.522 3 4 5 6 7 8 9 10 11012r [Å]g(r)(b)2 3 4 5 6 7 8 9 10 11012 (a)r [Å]g(r)Figure 4.12: Counterion-counterion radial pair correlation function of the ex-tended cylinder model model in (a) explicit and (b) implicit system forthree DNA separations (centre to centre), R = 23,25,28[A˚] , repre-sented by heavy dots, squares and triangles.We can also see from the results that the behaviour of the interactions betweentwo DNA molecules in the explicit solvent models, for both DNA structures aredifferent from the explicit solvent models. For the MAM, multiple attractive min-ima appear with large forces and at small DNA separations. The differences comefrom the structured nature of the solvation shells and the lack of screening at shortdistances.533 4 5 6 7 8 9 10 1100.51r [Å]g(r)(a) 3 4 5 6 7 8 9 10 1100.511.5 (b)r [Å]g(r) Figure 4.13: Counterion-counterion radial pair correlation function ofMontoro-Abascal model in (a) explicit and (b) implicit system forthree DNA separations (centre to centre), R = 22.5,25.25,28[A˚] , rep-resented by heavy dots, squares and triangles.54Chapter 5Coarse-grained model5.1 IntroductionThere is a strong interest in developing molecular models in order to describe struc-tural and dynamical properties of DNA. Fully atomistic models which are used toinvestigate some specific problems are computationally expensive and thereforeare not appropriate approaches to study phenomena at large length and time scales.On the other hand the simple models which are extensively used in investigatingthe problems, ignore important physics at the atomic scale. One way to overcomethis problem is by reducing the degrees of freedom and keeping only those degreeswhich are relevant and essential in that particular problem.In Chapters 3 and 4, we investigated the role of nanoscale effects of water as astructured solvent on the interaction between DNA molecules. The results showedthat the interactions change strongly as water molecules are added to the systemexplicitly. Therefore considering water as a continuum dielectric is not a properapproach to model the phenomena which occur at short length scale. In fact im-plicit solvent models with a primitive electrolyte approximation are too simplifiedto capture the ion mediated attraction between DNA molecules in the presenceof divalent ions. As modeling large system sizes with explicit solvent moleculesbecomes prohibitively expensive, only simulations with a coarse-grained (CG) rep-resentations of solvent remain feasible. In order to design a coarse-grained model,we need to write down a Hamiltonian which can be parametrized to properly match55Pair ε[kcal/mole] σ [A˚] rc[A˚] re[A˚] σe[A˚]CA-P 0 8.9 ... ... ...CA-CI 0.42 9.2 18.4 ... ...g-CI 0.1553 4.35 10 ... ...P-P 0.1553 4.2 ... 4.5 0.4P-CI 0.1553 3.683 10 3.2 0.55CI-CI 0.1553 3.166 10 5.8 1.3Table 5.1: Parameters for the interaction of different atoms (central atoms(CA), phosphate (P), counter ion (CI)) of initial guess which correspondto Fig. 5.1-5.5.a set of empirical information or quantities calculated from an atomistic model.In this Chapter, we present a coarse-grained representation for ds-DNA sol-vated by water molecules. This model makes a bridge between explicit and implicitsolvent models. This CG model replaces the water molecules by implicit solvationin order to significantly reduce the computational expense. However, short-rangedcorrections are added to the pair potentials of the primitive electrolyte such thatthe local molecular structure of the system is consistent with the results from ex-plicit solvent simulations. This CG model succeeds in reproducing the like-chargeattraction effect between DNA molecules in full atomistic simulations.5.2 Iterative Boltzmann inversionThe main goal is to develop a coarse-grained model which is quite easy to simulatebut reproduces the same physical behaviour as a reference atomistic simulation.The first step toward a CG model is to define an effective potential function for thesystem. The effective potential can be determined by the knowledge of all particlecorrelation functions [84]. In this project, we limit ourselves to single coordinatepairwise correlation functions, g(r). This approximation makes the simulationssimple and computationally efficient but on the other hand, it results in deviationfrom the correct underlying potential. The effective potential for each pair is de-56Pair H1 rmh,1 σh,1 H2 rmh,2 σh,2 H3 rmh,3 σh,3[kcal/mole] [A˚] [A˚] [kcal/mole] [A˚] [A˚] [kcal/mole] [A˚] [A˚]CI-CI -0.14 6.9 0.4 4 7.6 1.5 ... ... ...P-P -0.45 5.1 0.45 0.26 6.2 0.4 -0.14 7.6 0.4P-CI -0.18 2.9 0.25 5.68 3.75 0.55 -0.25 5 0.35g-CI -0.78 5.7 0.7 0.4 7 0.55 ... ... ...Table 5.2: Parameters for the hydration contribution in the potential of differ-ent atoms of initial guess which correspond to Fig. 5.1-5.5.rived by taking a simple Boltzmann inverse of the radial distribution functionF(r) =−kBT log(g(r)), (5.1)where F(r) is a free energy and potential energy only for the case of zero density.Such a ”potential of mean force” is temperature dependent and depends specificallyon state variables such as temperature, pressure and density. Determining and un-derstanding the limits of transferability of such potentials is one of the outstandingchallenges in coarse-graining of soft matter.Iterative Boltzmann inversion (IBI) is a procedure to systematically refine coarse-grained potentials. As we mentioned previously, the interactions between the par-ticles are chosen to be analytic functions with a few free parameters. IBI couplesthe information from explicit solvent model with iterative simulation to optimizethe free parameters in the potential to produce a set of radial distribution functionsthat have features similar to those obtained from explicit solvent simulations.The general approach to determine the CG potentials is the following:571. Determine the design of the CG model including particles and interactions be-tween them.2. Perform one or more reference simulations using the full atomistic model. Com-pute gre f (r) of each pair of particles from the full atomistic simulation.3. Calculate the pairwise potentialV0(r) between each pair of particles using Eq. 5.1.4. Optimize the free parameters in the analytical functions of the interaction be-tween each pair of particles. Alternatively, potentials can be tabulated numerically.5. Perform a simulation in the canonical ensemble for the CG system using theoptimized parameters. Measure the pair distributions the CG system.6. Update the potentials according toVi+1 =Vi− log[gi(r)gre f (r)], (5.2)7. Go back to step 4 and repeat a CG simulation. Update until all pair correlationfunctions measured in the CG model match the target pair correlation functions.This needs to be done for multiple rdfs simultaneously, and it is not clear a prioriif all rdfs can be converged with a single set of potentials.Reith et al. implemented similar iterative method for potential inversion fromdistribution functions for polymer systems. For every trial functions, convergencewas reached in few iterations [58]. Lyubartsev et al. introduced an alternativemethod to calculate effective interaction potentials from radial distribution func-tions by considering a reverse Monte Carlo approach [60].5.2.1 Coarse-grained potential for ionic speciesThe first step to define the Hamiltonian for the CG model is to choose the particlesand set the interactions between them. In all the solvent models, we simulatedso far, there are two types of interactions between particles in the system. TheLennard-Jones potential (LJ)V (r) = 4ε[(σr)12− (σr)6], (5.3)58−2−10U/kT(a) 8 9 10 11 12 130246 (b)r [Å]g(r) Figure 5.1: Development of the central atom-counterion potential. Solid linesshow the radial distribution function g(r) (panel (b)) and its Boltzmanninverse (panel (a)) obtained from a fully atomistic simulation of oneDNA strand. Dashed lines show the parameterized coarse grained po-tential (best fit to the potential of mean force) and the rdf that is obtainedin a coarse grained simulation of two DNA strands separated by 28 A˚.with LJ energy scale ε = 0.155 kcal/mole is employed between all pairs of par-ticles. The LJ length scale σ is considered as the diameter of the particles, andthe potential is truncated and shifted at rc. The values of σ and rc for all pairs ofparticles are given in the Table 2.1. Second, the Coulomb potential VC between twocharges qi and q j separated by a distance r isVqq(r) =qiq j4piε0εrr, (5.4)where ε0 is the permittivity of vacuum, and εr is the dielectric constant of thesolvent. For ions in vacuum or explicit water, we set εr = 1. For ions in implicitwater at T = 300K, it is to set εr = 72 at all distances since this is the bulk value ofthe dielectric constant in the SPC/E model [68, 85].594 5 6 7 8−2−11U/kT(a)4 5 6 7 8 90246 (b)r [Å]g(r)Figure 5.2: Development of the grooved atom-counterion potential. Solidlines show the radial distribution function g(r) (panel (b)) and its Boltz-mann inverse (panel (a)) obtained from a fully atomistic simulation ofone DNA strand. Dashed lines show the parameterized coarse grainedpotential (best fit to the potential of mean force) and the rdf that is ob-tained in a coarse grained simulation of two DNA strands separated by28 A˚.In the CG model, the approach we take to reproduce the effects of the solventon the electrolytes is by adding two contributions to the interaction energy betweenparticles. The first is a distance dependent dielectric permittivity representing thedielectric saturation of water in close vicinity to an ion. From a molecular point ofview, the solvent acts inhomogeneously over all the ions in the solution. Therefore,we replace εr in Eq. 5.4 by εD(r) which is the distant dependent dielectric permit-tivity. At very large separation of the ions, the dielectric permittivity reaches a bulkvalue. At small distances or close to the contact separation between the ions, thelower limit of εD(r→ 0) = 5.2 appears due to dielectric saturation [86]. The tran-sition between these two limiting dielectric regimes is modelled by a hyperbolicfunction which results in the Coulomb correction term for small ion separations60[87].εD(r) =(5.2+ εs2)+(εs−5.22)tanh[−(r− rme)22σ2e], (5.5)Here, εs is the dielectric of the bulk solvent,rme is the inflection point of the functionand σe is the width of transition between the bulk value and saturated dielectric astwo ions approach each other.The second short-range correction is an additional term to account for the pres-ence of one or more hydration shells. There are some favourable regions for watermolecules around the charged particles [64, 67].The hydration shells around ionsin the solution create a repulsive barrier which is modelled by a Gaussian function[63, 88].Vhydr(r) =Hσh√2piexp[−(r− rmh)22σ2h], (5.6)where Hσh√2pi is the potential height, rmh is the location of the hydration shell andσh is the width of the shell. The total ion-ion potential is then calculated as the sumof the individual contributions,Vion−ion(r) =VLJ(r)+Vqq(r)+Vhydr(r). (5.7)5.3 Simulation details5.3.1 Coarse-grained simulationsThe system that we want to model consists of two rigid DNA molecules whichare placed parallel to the z-axis and separated along the x-axis of a rectangularsimulation box. The box dimensions are Lx = 90A˚, Ly = 57A˚ and Lz = 34A˚. Wefocus here on the MAM model for the DNA structure as it provides the most re-alistic description of the shape of the DNA molecule. As before, all atoms onthe DNA strands are frozen. The total charge of the DNA molecules is compen-sated by adding divalent counterions in the system. The system is equilibrated for10ns using a canonical (NVT) ensemble, where the temperature was controlled at613 4 5 6 7−20246r [Å]U/kT(a)3 4 5 6 7 80510 (b)r [Å]g(r)Figure 5.3: Development of the phosphate-counterion potential. Solid linesshow the rdfs in full atomistic simulations (see text, panel (b)) and theBoltzmann inverse (panel (a)). Dashed lines show the initial guess ofthe coarse grained potential obtained from the Boltzmann inverse (panel(a)) and the rdf that is observed when the CG simulation of 2 DNAstrands is carried out (panel (b)).T = 300K applying a Langevin thermostat. The parameters for the LJ interactionsare given in Table 4.1.5.3.2 Fully detailed atomistic simulationsInitial guessEffective potentials are derived from the distributions in an atomistic trajectory.The atomistic simulations which are our starting point for the coarse graining in-clude two different sets of MD simulations. First, we perform a fully atomisticsimulation consisting of one rigid MAM-DNA molecule in a periodic box, with10 divalent free counterions which are placed randomly around it. The tetragonalsimulation box with the dimensions, Lx = Ly = 50A˚ and Lz = 34A˚, corresponds626 6.5 7 7.5 8 8.5369U/kT(a)5.5 6 6.5 7 7.5 8 8.5 900.5(b)r [Å]g(r)Figure 5.4: Development of the counterion-counterion potential. Solid linesshow the rdfs in full atomistic simulations (see text, panel (b)) and theBoltzmann inverse (panel (a)). Dashed lines show the initial guess ofthe coarse grained potential obtained from the Boltzmann inverse (panel(a)) and the rdf that is observed when the CG simulation of 2 DNAstrands is carried out (panel (b)).to one full turn of B-DNA. Periodic boundary conditions make a set of infiniteDNA molecules along the z-direction. The box is filled with 2129 SPC/E watermolecules which corresponds to the water density of 1g/cm3. Simulation parame-ters for water molecules, the DNA model and counterions are given in Tables 2.1and 4.1. The goal of running this simulation is to extract the distribution of counte-rions around the central atoms, the grooves and the phosphates of the DNA strand.These pair correlation functions (central atom-counterion, groove-counterion) areconsidered as the initial guess for construction of the coarse-grained model.A second set of atomistic simulation is designed to derive effective potentialsbetween different pairs of ions (phosphate-counterion, counterion-counterion andphosphate-phosphate). The size of all particles are the same as the two strandsatomistic model and are given in Tables 2.1 and 4.1. The periodic simulation box635 6 7 8 90510U/kT(a)5 6 7 8 9 1000.51 (b)r [Å]g(r)Figure 5.5: Development of the phosphate-phosphate potential. Solid linesshow the rdfs in full atomistic simulations (see text, panel (b)) andthe Boltzmann inverse (panel (a)), while dashed line show the coarsegrained potential obtained from the Boltzmann inverse (panel (a)).contains 40 mobile phosphates, 20 mobile counterions and 5412 water molecules.Particles are placed randomly inside the box and are allowed to explore the boxfreely. The dimensions of the box are Lx = 90A˚, Ly = 57A˚ and Lz = 34A˚. Thesystem is simulated in the NVT ensemble at T = 300K for 10ns. In both full atom-istic simulations, we apply a Nose´-Hoover thermostat in the system to keep thetemperature at T = 300K. Radial distribution functions from these two atomisticsimulations are used to optimize parameters of the coarse-grained model to obtainpotentials of mean force that serve as an initial guess.Final targetThe ultimate pair correlation functions to be reproduced correspond to anotheratomistic model which has the same configuration as the CG model. But in thismodel, two DNA molecules are considered at distance d = 28A˚ from each other643 4 5 6 7 8010203040r [Å]g(r)(b) −505 (a)U/kTFigure 5.6: a) Phosphate-counterion potential and b) Radial pair correlationfunction. Solid lines, dashed lines, dotted lines, dot-dash lines and cir-cles represent the function in the explicit model, initial guess and itera-tion 1, 2, 3 of the CG model, respectively.and are solvated by 5093 water molecules. We compared the distribution of coun-terions around each DNA molecules at d = 28A˚ with the distribution of counterionsin 1DNA model and we found that pair correlation function at this distance is in-dependent of the presence of the other DNA molecule. This atomistic simulationis performed for 30ns to get smooth radial distribution functions.5.4 Results and discussion5.4.1 Local structure from initial guessParametrization of the coarse-grained model starts with 5 pair correlation functionswhich are computed from the first two full atomistic simulations, as explained pre-viously. All pair potentials are extracted from the pair correlation functions usingEq. 5.1. The free parameters in the potential, Eq. 5.7 (σ in the Lennard-Jones658 9 10 11 12 13 1402468r [Å]g(r) Figure 5.7: Comparison of central atom-counterion radial pair correlationfunction in different solvation models. Solid line, dashed line and cir-cles represent this function in the explicit model, initial guess and finaliteration of CG model, respectively.repulsion, H, rmh and σh in the hydration repulsions, re and σe in the distancedependent dielectric function) are fitted to match the potential of mean force ob-tained from the explicit solvent simulations, see panels (a) of Figs. 5.1-5.5. Wekeep the LJ energy scale, ε = 0.155 kcal/mole between all pair of particles exceptthe central atom-counterion interaction in which the energy scale, ε = 0.42 kcal/-mole. We consider a new set of parameters which are represented in Tables 5.1and 5.2 as initial guess and the perform a coarse-grained simulation of two DNAstrands separated by r = 28A˚ for 10ns. Panels (b) of Figs. 5.1-5.5 compares theradial distribution functions measured in these simulations against their atomisticcounterparts. This comparison gives an impression of how well the initial guessreproduces the structure of the counterion cloud.In the explicit solvent model, there is a dielectric contrast between the volumeexcluded by the central atoms and the regions accessible by the water molecules.Therefore, the counterions see the central atoms as a dielectric with εr = 1. Hence663.5 4 4.5 5 5.5 6 6.5 7 7.5 801234567r [Å]g(r) Figure 5.8: Comparison of groove atom-counterion radial pair correlationfunction in different solvation models. Solid line, dashed line and cir-cles represent this function in the explicit model, initial guess and finaliteration of CG model, respectively.it is energetically more favourable for the counterions to be solvated in water thangetting very close to the surface of the central atoms. There is no such dielectriccontrast in the implicit solvent model. In order to capture the same distributionfunction for the counterions around the central atoms in the CG model, the coun-terions have to be pushed away from the surface of the central atoms. This is donein an ad hoc but efficient way by increasing central atom-counterion LJ diameterfrom σ = 5.683A˚ in the primitive model to σ = 9.2A˚ in the CG model. Fig. 5.1shows that counterions condense around central atoms at the same position as theexplicit solvent model. The LJ diameter of the grooved atom-counterion interac-tion is also increased from σ = 3.683A˚ in the primitive model to σ = 4.35A˚ in theCG model for the same reason. Two hydration shells are considered around thegrooved atoms at r = 5.7A˚ and r = 7A˚. These two corrections cause counterions tocondense around grooved atoms at the same positions as the explicit solvent model,see Fig. 5.2.67Pair ε[kcal/mole] σ [A˚] rc[A˚] re[A˚] σe[A˚]CA-P 0 8.9 ... ... ...CA-CI 0.3 9.2 18.4 ... ...P-P 0.155 4.2 ... 4.5 0.4P-CI 0.155 3.683 10 3.2 0.55g-CI 0.155 4.35 10 ... ...CI-CI 0.155 3.166 10 5.8 1.3Table 5.3: Parameters for the interaction of different atoms in final step, cor-responds to Figs. 5.6- 5.9.Pair H1 rmh,1 σh,1 H2 rmh,2 σh,2 H3 rmh,3 σh,3[kcal/mole] [A˚] [A˚] [kcal/mole] [A˚] [A˚] [kcal/mole] [A˚] [A˚]P-P -0.45 5.1 0.45 0.26 6.2 0.4 -0.14 7.6 0.4P-CI -1 2.9 0.35 4.87 3.7 0.45 -0.6 4.9 0.35g-CI -0.2 5.9 0.5 0.1 7.1 0.5 ... ... ...CI-CI -0.02 7 0.4 5.5 7.8 2 ... ... ...Table 5.4: Parameters for the hydration contribution in the potential of differ-ent atoms of final step which corresponds to Figs. 5.6- 5.9.In the explicit solvent model, counterions are distributed closer to the phos-phates when compered to the implicit solvent model because of the dielectric sat-uration, see Fig. 4.9. In the CG model, this feature is captured by consideringa distant dependent dielectric which results in the Coulomb correction term forsmall ion separations. For the pairs of phosphates and counterions, four solva-tion shells are considered and information for the first three solvation shells aregiven in Table 5.2. The forth solvation is located at r4 = 6.2A˚ with a height ofH4 = 0.38kcal/mole and a width of σ = 0.3A˚. For the phosphate-counterion paircorrelation function, the positions of the first and second peak are the same as theexplicit solvent model, see Fig. 5.3. However, the probability of finding them at685.5 6 6.5 7 7.5 8 8.5 900.20.40.60.8r [Å]g(r) Figure 5.9: Comparison of counterion-counterion radial pair correlationfunction in different solvation models. Solid line, dashed line and cir-cles represent this function in the explicit model, initial guess and finaliteration of CG model, respectively.those specific positions is lower in the CG model.There is a strong repulsion between pairs of like charged ions (counterionsor phosphates) in the explicit solvent according to the lack of dielectric screeningat short length scales. It can be modelled by considering the distant dependentdielectric function as explained previously. For both pairs (counterions and phos-phates), the magnitude of re in Eq. 5.3 is greater than the LJ diameter of the ions,see Table 5.1. We also consider two and three hydration shells for the counterion-counterion and phosphate-phosphate potential, respectively. The parameters em-ployed for the hydration shells are given in Table 5.2. Panel (b) of Fig. 5.4 showsthat the counterion-counterion radial distribution function in the CG model hasexcellent agreement with the results in the explicit solvent model. Panel (a) ofFig. 5.5 shows the phosphate-phosphate potential. In the CG model, phosphatesare fixed and located on the helices around the central atoms. Therefore there isno relevant correlation function in the CG model to calculate and compare with the69pair correlation function of mobile phosphates in the explicit solvent model whichis shown in panel (b) of Fig. 5.5.5.4.2 Model refinementThe next step is using the Iterative Boltzmann Inversion method to better reproducethe pair correlation functions of the atomistic model which has the same configu-ration as the CG model. Each individual distribution depends on the full set ofpotentials through higher-order correlations. Therefore we try to iterate all paircorrelation functions at every single step. However, the phosphate-phosphate paircorrelation function can not be improved as phosphates are fixed and rigid in thesystem of two DNA molecules. We start from the Boltzmann inverted potential ofthe fully atomistic simulation focus on the position and amplitudes of the minimaand maxima at short length scales in the target rdfs. The minima of the potentialsdetermine where the ion condensations occurs.The phosphate-counterion potential at each step of the iteration is shown inpanel (a) of Fig. 5.6. The depth of the first and second minima decreases at eachiteration and the phosphate-counterion pair correlation function matches the targetfunction well after 3 iterations, see panel (b) of Fig. 5.6. For the central atom-counterion pair correlation function, the energy scale of the LJ interaction is theonly parameter which is changed at each iteration. Fig. 5.7 shows that after threeiteration the first peak from the two models overlaps very well. For the groovedatoms-counterion rdf, the height and area of the first peak approaches to the targetrdf, see Fig. 5.8. Fig. 5.9 shows that for the counterion-counterion interaction, rdfsmatch at short length scale within some fluctuations. The final set of parametersfor the different pairs of particles are given in Tables 5.3, 5.4.5.4.3 Interaction between double stranded DNA moleculesAfter the third iteration, we stop and check the interaction between two DNAmolecules in the CG model that we are constructing. The parameters in Tables 5.3,5.4 are used to build the potentials for the different pairs of atoms. We increase theseparation between the DNA molecules and calculate the thermally averaged forceon each DNA molecule with the same protocol used in Chapters 3 and 4. Force7022 23 24 25 26 27−100−50050100R [Å]/F0 Figure 5.10: Reduced force between double stranded DNA molecules in dif-ferent solvation models. Solid line, circles, squares and triangles repre-sent this interaction in the old implicit solvent model, explicit solventmodel, and initial guess and final iteration of the CG model, respec-tively.averages were taken over t = 10ns for each DNA separation. Fig. 5.10 shows oursimulation results for the averaged reduced force in the two DNA models with dif-ferent solvation models. By comparing the results of the simulations, we see thatthere is good agreement between the behaviour of the force in the CG model withthe results of the explicit solvent model. In the CG model, two attractive regions(negative force) appear for specific DNA separations which agrees with the ex-plicit solvent model. However, the attractive regime disappears for r > 25.5A˚ inthe CG model. We also find the interaction between DNA molecules using the ini-tial guess potentials (using the parameters in Tables 5.1,5.2). The results are shownby squares in the Fig. 5.10. They show a similarly good match with the result ofexplicit solvent model. We see that most of the features of the behaviour of theinteraction in the explicit solvent model are already captured by considering thepotentials of the initial guess.715.5 ConclusionWe have developed a coarse-grained model for double stranded DNA solvated bywater molecules in the presence of divalent counterions. We applied the iterativeBoltzmann inversion method to construct the CG model in order to investigatethe effective interaction between DNA molecules. Starting with an initial guessaccording to the two fully atomistic simulations, we optimized four different po-tentials with 40+ parameters in three iterations. The IBI algorithm produced nu-merical potentials that lead to pair correlation functions which are very similar tothe target pair correlation functions. The IBI successfully generates the local ionicstructure of the system. The forces generated by new defined potentials are verysimilar to those observed in the explicit solvent model. The CG pair potentialssucceed in reproducing the like-charge attraction effect between DNA moleculesin full atomistic simulations at significantly reduced computational expense. TheCG model is around 2000 times faster than the equivalent explicit solvent model.72Chapter 6Three-DNA model6.1 IntroductionIn this Chapter, we investigate the interaction between three DNA molecules inthe presence of divalent counterions as a starting point for investigating many-body effects in the mechanism of DNA bundling. The induced effective attractiveinteraction between like charged biopolymers (DNA, actin filaments (F-actin), mi-crotubules (MT), tobacco mosaic virus (TMV), and the filamentous bacteriophage)leads to the formation of bundles under specific conditions [6]. Henle et al. inves-tigated the equilibrium bundle size of rodlike polyelectrolytes in the presence ofmultivalent ions, assuming that the attraction between polyelectrolytes are pairwiseadditive [89]. They showed how the size of counterions affects the equilibrium sizeof bundle. In another study, Ha et al. showed that in the presence of multivalentcounterions, like-charged polyelectrolytes attract each other but the effective inter-action is not pairwise additive [90].In this Chapter, we try to remove the ambiguity about the pairwise additivityor non-additivity of the interactions between DNA strands. We consider a three-DNA system with two different geometries for the DNA (ECM and MAM) in theprimitive solvent model. The three DNA strands are placed at the vertices of anequilateral triangle as shown in Fig. 6.1. We calculate the interaction between DNAmolecules and compare it with the force on each DNA in the two-DNA system.We consider the same solvation model and system size and we check whether the73Figure 6.1: Three DNA configuration (shown here for the Montoro-Abascalmodel). Phosphate charges are shown as blue spheres. The cylindricalcore is coloured in yellow and the neutral particles as green spheres.The DNA molecules are located at the vertices of an equilateral triangle.Each two DNA molecules are separated by R.effective interstrand interactions are pairwise additive.In Chapter 5, we found that the coarse-grained model reproduces the structureof counterions around DNA molecules in a manner consistent with the results ofthe full atomistic simulations, and successfully predicts the emergence of attractiveforces between two ds-DNA molecules. In this Chapter, we apply the developedcoarse-grained model to investigate the emergence of attractive interaction betweenthree DNA molecules in the presence of divalent counterions.6.2 Pairwise additivity of interactions in the three-DNAsystemWe start with a system in which the three double stranded DNA molecules are lo-cated at the vertices of an equilateral triangle. We consider the centre-to-centre7422 23 24 25 26 27 28−40−20020R [Å]/F 0Figure 6.2: Reduced DNA-DNA interaction force 〈F〉/F0 as a function ofDNA separation (centre to centre) for the extended cylindrical model(ECM) in primitive electrolyte approximation. The unit of the forceis F0 = kBT/P, where P = 34A˚ is the DNA pitch length. The circlesrepresent force in the three-DNA-system. The triangles represent theforce on DNA in the two-DNA-system times√3.distance of each pair of ds-DNA, R = 22A˚ in a rectangular simulation box withdimensions, Lx = Ly = 90A˚ and Lz = 34A˚. The x and y dimensions of the sim-ulation box are chosen to be three times greater than the maximum separation ofeach two DNA molecules. For simplification, symmetric relative orientations areconsidered for the DNA molecules, see Fig. 6.1. For R = 22A˚, the simulation runsfor 10ns in order to reach an equilibrium state. We consider the final configura-tion of this run as the initial state for the separation simulations and increase allseparations between the DNA strands from R= 22A˚ to R= 30A˚ simultaneously sothat the strands remain in an equilateral triangle. At each separation, simulationsare performed for 10ns and 60ns for the ECM and the MAM, respectively. Werecord the force on each phosphate every 10−5ns and measure the thermally av-eraged force components in the presence of divalent counter ions. We ensure thatthe total averaged forces, 〈F〉 on each of three DNA molecules are the same. Theprimitive electrolyte approximation (uniform dielectric constant) is applied and the7522 23 24 25 26 27 28 29 3010152025R [Å]/F 0Figure 6.3: Reduced DNA-DNA interaction force 〈F〉/F0 as a function ofDNA separation (centre to centre) for the Montoro-Abascal model(MAM) in primitive electrolyte approximation. The circles representforce in the three-DNA-system. The squares represent the force on DNAin the two-DNA-system times√3.DNA molecules are surrounded by divalent counterions so that the system is glob-ally charge neutral. We also run simulations for two DNA molecules (both ECMand MAM) with the same relative orientation in the same simulation box for dif-ferent DNA separations. At each separation, simulations are performed for 40nsand 60ns for the ECM and the MAM, respectively, and we calculate the averagedforce on each DNA molecule.Fig. 6.2 shows the results of the simulations for the system with the ECM ge-ometry for the DNA molecules. Comparing the force on each DNA molecule inthree-DNA and two-DNA systems shows that at short length scale the force is notsimply the vector summation of the interaction between two DNA molecules. Ifthe forces were additive, we would expect |F3−DNA| =√3|F2−DNA| for the trian-gular configuration. However, the forces converge upon increasing the separationbetween DNA molecules.For the MAM, the distribution of counterions around three DNA molecules re-sults in a totally different force on each molecules in the three-DNA system com-7622 23 24 25 26 27 28−200−150−100−50050100R [Å]/F 0Figure 6.4: Reduced DNA-DNA interaction force 〈F〉/F0 in three-DNAmodel in the presence of divalent counterions, considering different sol-vation models. Circles and triangles represent this force in the coarse-grained model in two- DNA-system times√3 and in three-DNA model,respectively. Solid line represents the force on DNA for the primitivesolvent model.pared to the two-DNA system. However, the position of the minimum in both twoand three DNA systems is the same and interactions are repulsive for all DNAseparations, see Fig. 6.3.6.3 Coarse-grained three-DNA modelWe continue to study the same three-DNA-system which was described in the pre-vious section, but now apply the coarse-grained potentials developed in Chapter 5.Parameters for the interaction of different atoms and for the hydration contributionsin the potential are given in the Tables 5.3 and 5.4. The system is equilibrated atR = 22A˚ for 20ns using a canonical (NVT) ensemble, where the temperature wascontrolled at T = 300K applying a Langevin thermostat. Then the separation ofthe DNA molecules is increased from R= 22A˚ to R= 28A˚ and data is collected fort = 40ns to ensure the system has reached the equilibrium.Fig. 6.4 shows the reduced force on each DNA molecule for different solvent77models and two and three-DNA model. Using the coarse-grained model for thepotential of the system results in the emergence of attractive regimes. The force-separation curve features a multiple minima structure with maximum attractiveforces at R≈ 23A˚ and R≈ 25A˚ in the three-DNA model. The maximum attractionin the two-DNA coarse-grained model occurs almost at the same position as thefirst maximum position of three-DNA coarse-grained model. However, the forceon each molecule in the three-DNA system is greater compared to the two-DNAsystem at short length scales. Therefore, we can conclude again that DNA-DNAinteractions are not pairwise additive.6.4 ConclusionsThe results of this Chapter show directly that the effective interstrand interactionsare not pairwise additive. For the ECM, the maximum of attraction regime appearsat shorter distances compared to the two strand system. For the MAM, the mini-mum of the interaction occurs almost at the same position for both two and three-DNA-system. However in the three-DNA-system the distribution of counterionsscreens the repulsive force between like-charged DNA molecules and decreasesthe magnitude of the interaction between DNA molecules compared to two-DNAsystem.The results of the simulations also show that applying the coarse-grained modeldeveloped previously to study three DNA molecules in the presence of divalentcounterions results in the emergence of attractive interactions for multiple DNAconfigurations. Again DNA-DNA interactions are not pairwise additive. Thecoarse grained models appear as a promising starting point for investigating thestability of DNA bundles in the presence of divalent counterions.78Chapter 7Conclusions and future work7.1 ConclusionsThe work presented in this thesis is a systematic investigation to understand theeffects of water as a structured solvent on the phenomenon of counterion-inducedDNA condensation. We performed classical molecular dynamics (MD) simulationsof the electrostatic interaction between two rigid DNA molecules in the presenceof divalent counterions with different DNA and solvation models. Since the de-sired time and length scales required to study biophysical situations such as DNAbundling are beyond atomistic simulations, we also constructed a coarse-grainedimplicit solvent model for investigating the dynamics on long time and lengthscales.In Chapter 3, we investigated the role and importance of representing wateras explicit solvent for the like charge attraction effect between two rigid chargedrods as a highly simplified model for a charged DNA molecule. We obtained theaverage force between two parallel charged rods in simulations that differ only bytheir representation of water as an implicit or explicit solvent. We consider twocases: thin rods, where the diameter of the rods is comparable to that of the ions,and thick rods, with a diameter comparable to that of ds-DNA. The results of thesimulation shows that the presence of water molecules changes the structure of thecounterions at short distances. Therefore, the behaviour of interactions betweenrods in the explicit solvent differs qualitatively and quantitatively from implicit79solvent models for both thin and thick rod systems[81]. We can see from the re-sults that multiple attractive minima appears in the force-distance curve for the thinrods system as the SPC/E water molecules are added explicitly to the system. Thelack of screening at short distances and the structured nature of the solvation shellsresults in much larger forces in the explicit solvent model compared to an equiva-lent implicit solvent simulation. The attractive regime occurs at much smaller rodseparations in the system with full representation of water molecules. In the simu-lations with explicit solvent, the dielectric contrast between the interior of the rodsand the regions filled with water is taken into account, however it is not consideredin implicit solvent models with a uniform dielectric constant. This effect results indisappearance of the attractive regime in the explicit solvent models of the thickrods system.In Chapter 4, we continued the investigation on the interaction forces betweena pair of double stranded model DNA molecules in the presence of divalent ionsusing molecular dynamics simulations. Since the attractive regime disappeared inthe explicit solvation simulation of the cylindrical DNA, we decided to explore theimportance of the DNA charge pattern on the mutual DNA-DNA interaction. Weconsidered two different helical models (ECM and MAM) which represent a real-istic description of the charge distribution in the DNA double strand. The resultsof the molecular dynamics simulations revealed the role of the DNA geometry onthe electrostatic interaction between DNA molecules. Counterion condensationpatterns on the DNA surface are different in the two DNA structures. Consider-ing helical models, the regime of attractive interaction is recovered for the explicitand implicit model of ECM and the explicit model of the MAM. The results alsoconfirm that the behaviour of the interactions between two DNA molecules in theexplicit solvent model are both qualitatively and quantitatively different from theimplicit solvent models. The structured nature of the solvation shells and the lackof screening at short distances results in different behaviour in different solventmodels. For the MAM, multiple attractive minima occur at smaller DNA separa-tions with larger force magnitudes while the interaction is always repulsive in anequivalent implicit solvent simulation.The results of our simulations are consistent with the available experimentalstudies [49, 52]. As we mentioned in the introduction, Koltover et al. confirmed80the existence of attraction between λ -phage DNA molecules confined to a two-dimensional cationic in the presence of divalent counterions [52]. They appliedsynchrotron X-ray diffraction experiments and showed that the electrostatic forcesbetween DNA chains reverse from repulsive to attractive above a critical divalentcounterion concentration. Qiu et al. also presented the evidences for attractionbetween multiple DNA molecules in the presence of Mg2+ counterions [49].For the MAM model with explicit solvent model, we found maximum attrac-tive forces of 0.088nN and 0.048nN at the DNA separations 22.5A˚ and 25.25A˚respectively, see Fig. 4.5. The results of our simulations are quantitatively in goodagreement with the results of two available fully atomistic models in which twoDNA strands are represented in full molecular detail and solvated with explicitwater molecules. Luan et al. reported direct observation of DNA attraction fromall-atom MD simulations in the presence of divalent electrolytes [50]. Their simu-lations showed the existence of an attractive force in the presence of MgCl2 at theDNA separation of 26A˚. They reported maximum attractive force of ∼ 0.042nNper helical turn. The study by Dai et al. was able to explain the attractive forcein terms of the formation of transient ion bridges [48]. Their results show a weakattraction regime in the presence of Pu(divalent ion). The depth of potential is lessthan 2kBT and it’s minimum occurs at DNA separation of 24A˚.The results of Chapters 3 and 4 demonstrate the importance of nanoscale ef-fects of water as a structured solvent on the interaction between DNA molecules.However, modeling large system sizes with explicit solvent molecules becomes ex-cessively expensive and is therefore not appropriate to study phenomena at largelength and time scales. One way to overcome this problem is simulations with acoarse-grained (CG) representations of the solvent.In Chapter 5, we developed a coarse-grained model for double stranded DNAsolvated by water molecules in the presence of divalent counterions. We appliedthe iterative Boltzmann inversion (IBI) method to construct the CG model fromall-atom simulations. This CG model was developed in order to investigate theeffective interaction between DNA molecules on longer time and length scales.Short-ranged corrections were added to the pairwise interaction potentials such thatthe structure of counterions is consistent with results from corresponding explicitsolvent simulations. We considered two fully atomistic simulations as an initial81guess to parameterize the added short-range corrections. Using the IBI algorithm,we optimized four different potentials with more than 40 parameters in three itera-tions. The optimized parameters produced numerical potentials that lead to radialdistribution functions which are very similar to the target rdfs. The coarse-grainedmodel succeeded in reproducing the local ionic structure of the system. The ef-fective force between two DNA strands generated by these potentials provided anexcellent match to that observed in the explicit solvent model. Importantly, this in-teraction features multiple minima and reproduces the like-charge attraction effectbetween DNA molecules observed in fully atomistic simulations at significantlyreduced computational expense. This result proves that it is possible to capturecomplex multibody interactions between polyelectrolyte strands with two-body po-tentials.In Chapter 6, we applied the coarse-grained model to study three DNA strandsin the presence of divalent counterions as a starting point for investigating many-body effects in the mechanism of DNA bundling. The results of this Chaptershowed directly that the effective interstrand interactions are not pairwise additive.Moreover, the results showed that applying the coarse-grained model developedpreviously to study three DNA molecules in the presence of divalent counterionsresults in the emergence of attractive interactions for multiple DNA configurations.7.2 Future workThis work showed that key features of the effective interactions between chargedpolyelectrolytes can be successfully captured with a coarse-graining strategy thataccounts for solvation effects through short ranged pair potentials. This methodnow opens up the door to investigate and reexamine more complex scenarios. Animmediate extension of Chapter 6 would be to investigate formation and stabilityof larger DNA bundles. It would also be interesting to revisit the osmotic pres-sure calculations in bulk systems, where DNA molecules form a hexagonal lattice.The previous studies by Deserno et al. [33] (cylindrical model) and Lee et al. [91](helical model) were done only with the primitive electrolyte model.In all of our simulations so far, DNA was modelled as a stiff rigid polymer.Such a picture is appropriate when the length of DNA is less than the persistence82length (lp = 500A˚), but for longer strands the finite bending stiffness should betaken into account. For a short DNA molecule, the associated entropy contributesonly little to the free energy and can be disregarded. However, for flexible poly-electrolytes, bending fluctuations can become important. It would be interesting toapply the present strategy to a flexible coarse-grained DNA model in a simulationthat also allows for DNA rotation. Such a study would make a further step towardsa fully dynamical simulation of DNA aggregation, which becomes possible withthe coarse-grained model developed in this thesis. Further interesting extensionscould include the application of the coarse-graining technique to DNA adsorptiononto surfaces, where the structure of water and counterions locally differ from thebulk.83Bibliography[1] Y. Levin. Electrostatic correlations: from plasma to biology. Rep. Prog.Phys., 65(11), 2002.[2] G. C. L. Wong and L. Pollack. Electrostatics of strongly charged biologicalpolymers: ion-mediated interactions and self-organization in nucleic acidsand proteins. Annu Rev Phys Chem, 61, 2010.[3] F. Oosawa. Interaction between parallel rodlike macroions. Biopolymers,6(11), 1968.[4] V. A. Bloomfield. DNA condensation. Current Opinion in Structural Biology,6(3), 1996.[5] V. A. Bloomfield. DNA condensation by multivalent cations. Biopolymers,44(3), 1997.[6] J. X. Tang, S. Wong, P. T. Tran, and P. A. Janmey. Counterion induced bundleformation of rodlike polyelectrolytes. Berichte der Bunsengesellschaft fu¨rphysikalische Chemie, 100(6), 1996.[7] L. Guldbrand, L. G. Nilsson, and L. Nordenskio¨ld. A Monte Carlo simula-tion study of electrostatic forces between hexagonally packed DNA doublehelices. The Journal of Chemical Physics, 85(11), 1986.[8] A. P. Lyubartsev and L. Nordenskio¨ld. Monte carlo simulation study ofion distribution and osmotic pressure in hexagonally oriented dna. J. Phys.Chem., 99(25), 1995.84[9] D. Marenduzzo. Computer Simulations of DNA Packing inside Bacterio-phages: Elasticity, Electrostatics and Entropy. Computational and mathe-matical methods in medicine, 9(3-4), 2008.[10] D. E. Smith, S. J. Tans, S. B. Smith, S. Grimes, D. L. Anderson, and C. Bus-tamante. The bacteriophage 29 portal motor can package DNA against a largeinternal force. Nature, 413(6857), 2001.[11] M. A. Robinson, J. P. A. Wood, S. A. Capaldi, A. J. Baron, C. Gell, D. A.Smith, and N. J. Stonehouse. Affinity of molecular interactions in the bacte-riophage 29 DNA packaging motor. Nucl. Acids Res., 34(9), 2006.[12] I. Ivanovska, G. Wuite, B. Jo¨nsson, and A. Evilevitch. Internal DNA pressuremodifies stability of WT phage. Proc Natl Acad Sci USA, 104(23), 2007.[13] D. N. Fuller, J. P. Rickgauer, P. J. Jardine, S. Grimes, D. L. Anderson, andD. E. Smith. Ionic Effects on Viral DNA Packaging and Portal Motor Func-tion in Bacteriophage 29. Proceedings of the National Academy of Sciencesof the United States of America, 104(27), 2007.[14] J. P. Rickgauer, D. N. Fuller, S. Grimes, P. J. Jardine, D. L. Anderson, andD. E. Smith. Portal Motor Velocity and Internal Force Resisting Viral DNAPackaging in Bacteriophage 29. Biophysical Journal, 94(1), 2008.[15] H. Farhood, X. Gao, K. Son, J. S. Lazo, L. Huang, J. Barsoum, R. Bottega,and R. M. Epand. Cationic Liposomes for Direct Gene Transfer in Therapyof Cancer and Other Diseasesa. Annals of the New York Academy of Sciences,716(1), 1994.[16] Y. Liu, D. Liggitt, W. Zhong, G. Tu, K. Gaensler, and R. Debs. CationicLiposome-mediated Intravenous Gene Delivery. J. Biol. Chem., 270(42),1995.[17] M. Mandelkern, J. G. Elias, D. Eden, and D. M. Crothers. The dimensions ofDNA in solution. Journal of Molecular Biology, 152(1), 1981.85[18] N. V. Hud, F. P. Milanovich, and R. Balhorn. Evidence of Novel SecondaryStructure in DNA-Bound Protamine Is Revealed by Raman Spectroscopy.Biochemistry, 33(24), 1994.[19] J. G. Duguid, V. A. Bloomfield, J. M. Benevides, and G. J. Thomas. Ramanspectroscopy of DNA-metal complexes. II. The thermal denaturation of DNAin the presence of Sr2+, Ba2+, Mg2+, Ca2+, Mn2+, Co2+, Ni2+, and Cd2+.Biophys J, 69(6), 1995.[20] A. P. Lyubartsev and L. Nordenskio¨ld. Monte Carlo Simulation Study ofDNA Polyelectrolyte Properties in the Presence of Multivalent PolyamineIons. J. Phys. Chem. B, 101(21), 1997.[21] J. C. G. Montoro and J. L. F. Abascal. Ionic distribution around simpleDNA models. I. Cylindrically averaged properties. The Journal of Chemi-cal Physics, 103(18), 1995.[22] U. P. Strauss. The Collected Papers of Peter J. W. Debye. Interscience, NewYork-London. Science, 120(3118), 1954.[23] F. Fogolari, A. Brigo, and H. Molinari. The Poisson Boltzmann equation forbiomolecular electrostatics: a tool for structural biology. J. Mol. Recognit.,15(6), 2002.[24] R. R. Netz and H. Orland. Beyond Poisson-Boltzmann: Fluctuation effectsand correlation functions. Eur. Phys. J. E, 1(2-3), 2000.[25] A. Naji, A. Arnold, C. Holm, and R. R. Netz. Attraction and unbinding oflike-charged rods. EPL, 67(1), 2004.[26] S. L. Brenner and D. A. McQuarrie. On the theory of the electrostatic inter-action between parallel cylindrical polyelectrolytes. Journal of Colloid andInterface Science, 44(2), 1973.[27] D. Harries. Solving the Poisson Boltzmann Equation for Two Parallel Cylin-ders. Langmuir, 14(12), 1998.86[28] G. N. Patey. The interaction of two spherical colloidal particles in electrolytesolution. An application of the hypernetted chain approximation. The Journalof Chemical Physics, 72(10), 1980.[29] A. G. Moreira and R. R. Netz. Strong-coupling theory for counter-ion distri-butions. Europhysics Letters (EPL), 52(6), 2000.[30] N. Grønbech-Jensen, R. J. Mashl, R. F. Bruinsma, and W. M. Gelbart.Counterion-Induced Attraction between Rigid Polyelectrolytes. Phys. Rev.Lett., 78(12), 1997.[31] M. J. Stevens. Bundle Binding in Polyelectrolyte Solutions. Phys. Rev. Lett.,82(1), 1999.[32] M. J. Stevens. Simple simulations of DNA condensation. Biophys J, 80(1),2001.[33] M. Deserno, A. Arnold, and C. Holm. Attraction and Ionic Correlations be-tween Charged Stiff Polyelectrolytes. Macromolecules, 36(1), 2003.[34] J. Ray and G. S. Manning. An attractive force between two rodlike polyionsmediated by the sharing of condensed counterions. Langmuir, 10(7), 1994.[35] B. Y. Ha and Andrea J. Liu. Counterion-Mediated Attraction between TwoLike-Charged Rods. Phys. Rev. Lett., 79(7), 1997.[36] B. Y. Ha and Andrea J. Liu. Effect of Non-Pairwise-Additive Interactions onBundles of Rodlike Polyelectrolytes. Phys. Rev. Lett., 81(5), 1998.[37] B. I. Shklovskii. Wigner Crystal Model of Counterion Induced Bundle For-mation of Rodlike Polyelectrolytes. Phys. Rev. Lett., 82(16), 1999.[38] E. Allahyarov and H. Lo¨wen. Effective interaction between helicalbiomolecules. Phys. Rev. E, 62(4), 2000.[39] A. Diehl, H. A. Carmona, and Y. Levin. Counterion correlations and attrac-tion between like-charged macromolecules. Phys. Rev. E, 64(1), 2001.87[40] M. Kanducˇ, J. Dobnikar, and R. Podgornik. Counterion-mediated electro-static interactions between helical molecules. Soft Matter, 5(4), 2009.[41] M. Sayar and C. Holm. Equilibrium polyelectrolyte bundles with differentmultivalent counterion concentrations. Phys. Rev. E, 82(3), 2010.[42] I. Rouzina and V. A. Bloomfield. Macroion Attraction Due to ElectrostaticCorrelation between Screening Counterions. 1. Mobile Surface-AdsorbedIons and Diffuse Ion Cloud. J. Phys. Chem., 100(23), 1996.[43] A. A. Kornyshev and S. Leikin. Electrostatic Zipper Motif for DNA Aggre-gation. Phys. Rev. Lett., 82(20), 1999.[44] V. Tereshko, G. Minasov, and M. Egli. The Dickerson-Drew B-DNA Dode-camer Revisited at Atomic Resolution. J. Am. Chem. Soc., 121(2), 1999.[45] E. Allahyarov, H. Lo¨wen, and G. Gompper. Adsorption of monovalent andmultivalent cations and anions on DNA molecules. Phys. Rev. E, 68(6), 2003.[46] E. Allahyarov, G. Gompper, and H. Lo¨wen. Attraction between DNAmolecules mediated by multivalent ions. Phys. Rev. E, 69(4), 2004.[47] E. Allahyarov, G. Gompper, and H. Lo¨wen. DNA condensation and redisso-lution: interaction between overcharged DNA molecules. J. Phys.: Condens.Matter, 17(20), 2005.[48] L. Dai, Y. Mu, L. Nordenskio¨ld, and J. R. C. van der Maarel. Moleculardynamics Simulation of Multivalent-Ion Mediated Attraction between DNAMolecules. Phys. Rev. Lett., 100(11), 2008.[49] X. Qiu, K. Andresen, L. W. Kwok, J. S. Lamb, H. Y. Park, and L. Pollack.Inter-DNA Attraction Mediated by Divalent Counterions. Phys. Rev. Lett.,99(3), 2007.[50] B. Luan and A. Aksimentiev. DNA Attraction in Monovalent and DivalentElectrolytes. J. Am. Chem. Soc., 130(47), 2008.[51] J. Pelta, F. Livolant, and J. Sikorav. DNA Aggregation Induced by Polyaminesand Cobalthexamine. J. Biol. Chem., 271(10), 1996.88[52] I. Koltover, K. Wagner, and C. R. Safinya. DNA condensation in two dimen-sions. Proc Natl Acad Sci U S A, 97(26), 2000.[53] L. R. Pratt, G. Hummer, and A. E. Garcia´. Ion pair potentials-of-mean-forcein water. Biophysical Chemistry, 51, 1994.[54] P. A. Bopp, A. A. Kornyshev, and G. Sutmann. Frequency and wave-vectordependent dielectric function of water: Collective modes and relaxation spec-tra. The Journal of Chemical Physics, 109(5), 1998.[55] J. Rottler and B. Krayenhoff. Numerical studies of nonlocal electrostaticeffects on the sub-nanoscale. J. Phys.: Condens. Matter, 21(25), 2009.[56] A. C. Belch, M. Berkowitz, and J. A. McCammon. Solvation structure of asodium chloride ion pair in water. J. Am. Chem. Soc., 108(8), 1986.[57] R. Faller, H. Schmitz, O. Biermann, and F. Mu¨ller-Plathe. Automatic param-eterization of force fields for liquids by simplex optimization. J. Comput.Chem., 20(10), 1999.[58] D. Reith, M. Pu¨tz, and F. Mu¨ller-Plathe. Deriving effective mesoscale poten-tials from atomistic simulations. J. Comput. Chem., 24(13), 2003.[59] A. P. Lyubartsev and A. Laaksonen. Effective potentials for ion-DNA inter-actions. Journal of Chemical Physics, 111(24), 1999.[60] A. P. Lyubartsev and A. Laaksonen. Calculation of effective interaction po-tentials from radial distribution functions: A reverse Monte Carlo approach.Phys. Rev. E, 52(4), 1995.[61] R. C. DeMille and V. Molinero. Coarse-grained ions without charges: Repro-ducing the solvation structure of NaCl in water using short-ranged potentials.The Journal of Chemical Physics, 131(3), 2009.[62] R. C. DeMille, T. E. Cheatham, and V. Molinero. A Coarse-Grained Modelof DNA with Explicit Solvation by Water and Ions. J. Phys. Chem. B, 115(1),2011.89[63] P. J. Lenart, A. Jusufi, and A. Z. Panagiotopoulos. Effective potentials for 1:1electrolyte solutions incorporating dielectric saturation and repulsive hydra-tion. The Journal of Chemical Physics, 126(4), 2007.[64] J. E. Enderby and G. W. Neilson. The structure of electrolyte solutions. Rep.Prog. Phys., 44(6), 1981.[65] L. X. Dang, J. E. Rice, J. Caldwell, and P. A. Kollman. Ion solvation in po-larizable water: molecular dynamics simulations. J. Am. Chem. Soc., 113(7),1991.[66] L. Degre`ve and F. L. B. da Silva. Structure of concentrated aqueous NaClsolution: A Monte Carlo study. The Journal of Chemical Physics, 110(6),1999.[67] A. Kovalenko and F. Hirata. Potentials of mean force of simple ions in ambi-ent aqueous solution. II. Solvation structure from the three-dimensional refer-ence interaction site model approach, and comparison with simulations. TheJournal of Chemical Physics, 112(23), 2000.[68] B. Hess, C. Holm, and N. Vegt. Modeling Multibody Effects in Ionic So-lutions with a Concentration Dependent Dielectric Permittivity. Phys. Rev.Lett., 96(14), 2006.[69] G. S. Freeman, D. M. Hinckley, and J. J. de Pablo. A coarse-grain three-site-per-nucleotide model for DNA with explicit ions. The Journal of ChemicalPhysics, 135(16), 2011.[70] T. A. Knotts, N. Rathore, D. C. Schwartz, and J. J. de Pablo. A coarse grainmodel for DNA. Journal of Chemical Physics, 126(8), 2007.[71] A. Savelyev and G. A. Papoian. Chemically accurate coarse graining ofdouble-stranded DNA. Proc Natl Acad Sci USA, 107(47), 2010.[72] T. R. Prytkova, I. Eryazici, B. Stepp, S. Nguyen, and G. C. Schatz. DNA Melt-ing in Small-MoleculeDNA-Hybrid Dimer Structures: Experimental Charac-terization and Coarse-Grained Molecular Dynamics Simulations. The Journalof Physical Chemistry B, 114(8), 2010.90[73] F. Otto and G. N. Patey. Forces between like-charged walls in electrolytesolution: Molecular solvent effects at the McMillan Mayer level. Journal ofChemical Physics, 112(20), 2000.[74] W. C. K. Poon and D. Andelman. Soft condensed matter physics in molecularand cell biology. Scitech Book News, 30(1), 2006.[75] G. S. Grest and K. Kremer. Molecular dynamics simulation for polymers inthe presence of a heat bath. Phys. Rev. A, 33(5), 1986.[76] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, and J. Hermans.Interaction Models for Water in Relation to Protein Hydration. In B. Pull-man, editor, Intermolecular Forces, number 14 in The Jerusalem Symposiaon Quantum Chemistry and Biochemistry. Springer Netherlands, 1981.[77] H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma. The missing term ineffective pair potentials. J. Phys. Chem., 91(24), 1987.[78] S. Nose´. A unified formulation of the constant temperature molecular dynam-ics methods. The Journal of Chemical Physics, 81(1), 1984.[79] L. Verlet. Computer ”Experiments” on Classical Fluids. I. ThermodynamicalProperties of Lennard-Jones Molecules. Phys. Rev., 159(1), 1967.[80] S. Plimpton. Fast Parallel Algorithms for Short-Range Molecular Dynamics.Journal of Computational Physics, 117(1), 1995.[81] S. Ghanbarian and J. Rottler. Solvation effects on like-charge attraction. TheJournal of Chemical Physics, 138(8), 2013.[82] G. S. Manning. Limiting Laws and Counterion Condensation in Polyelec-trolyte Solutions. III. An Analysis Based on the Mayer Ionic Solution Theory.The Journal of Chemical Physics, 51(8), 1969.[83] H. Flyvbjerg and H. G. Petersen. Error estimates on averages of correlateddata. The Journal of Chemical Physics, 91(1), 1989.[84] J. Zwicker and R. Lovett. When does a pair correlation function fix the stateof an equilibrium system? The Journal of Chemical Physics, 93(9), 1990.91[85] M. Rami Reddy and M. Berkowitz. The dielectric constant of SPC/E water.Chemical Physics Letters, 155(2), 1989.[86] P. J. Stiles, J. B. Hubbard, and R. F. Kayser. Dielectric saturation and dielec-tric friction in electrolyte solutions. The Journal of Chemical Physics, 77(12),1982.[87] I. Danielewicz-Ferchmin and A. R. Ferchmin. Review: Water at ions,biomolecules and charged surfaces. Physics and Chemistry of Liquids, 42(1),2004.[88] A. Jusufi, A. Hynninen, and A. Z. Panagiotopoulos. Implicit Solvent Modelsfor Micellization of Ionic Surfactants. J. Phys. Chem. B, 112(44), 2008.[89] Mark L. Henle and P. A. Pincus. Equilibrium bundle size of rodlike polyelec-trolytes with counterion-induced attractive interactions. Phys. Rev. E, 71(6),2005.[90] B. Y. Ha and A. J. Liu. Counterion-mediated, non-pairwise-additive attrac-tions in bundles of like-charged rods. Phys Rev E Stat Phys Plasmas FluidsRelat Interdiscip Topics, 60(1), 1999.[91] S. Lee, T. T. Le, and T. T. Nguyen. Reentrant Behavior of Divalent-Counterion-Mediated DNA-DNA Electrostatic Interaction. Phys. Rev. Lett.,105(24), 2010.[92] R. W. Hockney and J. W. Eastwood. Computer Simulation Using Particles.CRC Press, 1988.[93] T. Darden, D. York, and L. Pedersen. Particle mesh Ewald: An Nlog(N)method for Ewald sums in large systems. The Journal of Chemical Physics,98(12), 1993.[94] U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Ped-ersen. A smooth particle mesh Ewald method. The Journal of ChemicalPhysics, 103(19), 1995.92[95] M. Deserno and C. Holm. How to mesh up Ewald sums. I. A theoretical andnumerical comparison of various particle mesh. Journal of Chemical Physics,109(18), 1998.[96] M. P. Allen. Computer simulation of liquids. Oxford science publications.Clarendon Press ; Oxford University Press, Oxford [England], 1989.[97] C. Whitmer. Over-relaxation methods for Monte Carlo simulations ofquadratic and multiquadratic actions. Phys. Rev. D, 29(2), 1984.[98] H. Flyvbjerg and H. G. Petersen. Error estimates on averages of correlateddata. Journal of Chemical Physics, 91(1), 1989.[99] C. Domb. Phase transitions and critical phenomena. Academic Press, 1972.93Appendix AInteraction potentialsIn all of the simulations particles interact via two kinds of pairwise interactions, theLennard-Jones potential, Eq. 2.9 and the Coulomb potential, Eq. 2.10. The inter-action potentials and definition of terms within them were introduced in chapter 2.The Coulomb potential is a long-range interaction which is inversely proportionalto the distance between charges and varies slowly at large distances. Consideringperiodic boundary condition, charges will interact with all images in the periodicboxes. Ewald summation is an efficient method which can be used to computelong-range interactions in periodic systems [92–95]. For a system of N particleswith charges qi and cubic box of length L, the total Coulomb potential energy in asimulation cell with periodic boundary condition is given byUC(r) =12∑nN∑i, j=1qiq j|nL+ ri j|, (A.1)where N is the total number of charges, ri j = ri− r j is a vector which gives thedistance between the charges, n is an integer vector identifying either the originalsimulation box (n = 0) or periodic image and L is the box length. The summationdoes not include the interaction of a charge with itself. A factor of 1/2 accountsfor the double counting of the interaction between the same pair of charges. In theEwald method the 1/r part of the potential is modified by dividing it into short andlong range parts.941r=f (r)r+1− f (r)r(A.2)The first term, f (r)r is a short range function which should be negligible beyond thetruncation limit (rcut). The second part[1− f (r)]r , is a slowly varying function for allr and is calculated efficiently in the reciprocal space. A traditional selection for thefunction f is the complementary error functioner f c(r) =2√pi∫ ∞rdΩexp(−Ω2), (A.3)α is the Ewald parameter that sets the relative weight of the real and reciprocalspace contribution to ensure both terms converge rapidly. The real and reciprocalspace contribution of the potential are given by,U (r) =18piε0 ∑i, j∑mqiq jer f c(α|ri j +mL|)|ri j +mL|, (A.4)U (k) =12ε0L3 ∑k 6=0N∑i, jqiq jk2exp(−k24α2 )cos(k.ri j). (A.5)where k is reciprocal lattice vector. In addition, a self energy term U s, for the inter-action of the cancelling distribution with its own point charge must be subtractedof the potential of the system. For the molecular system, an other self-energy term,accounting for intra-molecular Coulombic interactions is also subtracted from thepotential energy [96].U (s) =α4pi3/2ε0N∑iq2i . (A.6)U (d) =14piε0 ∑i(∑aqa∑bqber f (α|d|)|d|) (A.7)where a and b are different charges in the same molecule i and d is the vectorbetween different charges qa and qb. Thus the total Coulomb potential energy isgiven by,95UCoul =U(r)+U (k)−U (s)−U (d). (A.8)We use Ewald method to compute long-range interactions in the implicit sol-vent model of a system of ∼ 100 particles. But for the explicit solvent model witharound 15000 particles, we use particle-particle particle-mesh (PPPM) solver. ThePPPM method was presented by Hockney [92] and it maps the particle charge to a3D grid and uses a Fast Fourier Transform (FFT) to solve poisson equation on thegrid. Then, it interpolates the electric fields on the mesh points back to the atoms.The PPPM solver is a faster choice for the system with large number of particles.96Appendix BBlocking methodThe blocking method is a renormalization group method to calculate the statisticalerror on an average of correlated data. It was invented by K. Wilson and describedby Whitmer [97] and Flyvbjerg et al. [98]. This renormalization method is inter-estingly applied on a real, discrete, one dimensional time series of correlated datawhich is produced by Monte Carlo or molecular dynamics simulations.If the time series, x1,x2, ...,xn are the results of a molecular dynamics simula-tion in thermal equilibrium, by considering ergodicity, the averaged value of x andvariance are given by,m≡ x =1nn∑i=1xi (B.1)σ2(m) = 〈m2〉−〈m〉2 (B.2)By introducing γi, j ≡ 〈xix j〉 − 〈xi〉〈x j〉 as the correlation between xi and x j,γt = γi, j and t = |i− j|, we find,σ2(m) = 1n2n∑i, j=1γi, j =1n[γ0 +2n−1∑t=1(1−tn)γt ], (B.3)In order to calculate σ2(m), an estimation for γt is required. 〈ct〉 is usually usedas an estimator to calculate γt ,97ct ≡1n− tn−t∑k=1(xk−x)(xk+t −x), (B.4)The expectation value of ct is related to γt in the following way,〈ct〉= γt −σ2(m)+∆t , (B.5)where∆t = 2(1nn∑i=1−1n− tn−t∑i=1)1nn∑j=1γi, j. (B.6)Considering the largest correlation time as τ and choosing an appropriate cut-offparameter T , Binder [99] estimated σ2(m) in the following way,σ2(m)≈〈c0 +2∑Tt=1 ctn−2T −1〉. (B.7)Calculation of ct and choosing the value of T make the estimation complicated.However, blocking method estimates σ2(m) computationally more efficient by re-peatedly blocking the data which avoids computation of ct and the choice of T . Atthe first step the data set x1,x2, ...,xn is transferred into half x′1,x′2, ...,x′n, wherex′i =12(x2i−1 + x2i), (B.8)n′ =12n, (B.9)m′ = m, (B.10)σ2(m′) = σ2(m). (B.11)Which shows that m and σ2(m) are invariant under blocking. Eq. B.3 shows thatσ2(m)≥ γ0n . It is obvious that at each step of the methodγ0n increases and one canshow that there is a fixed point for this transformation in which σ2(m) = γ0n and980 5 10 15 20051015x 10−5Number of block transformation appliedσ [nN]Figure B.1γt = 0 for t > 0. Therefore σ2(m) can be estimated by using Eqs. B.4 and B.5 andconsidering ∆0 = 0,σ2(m) =〈c0n−1〉. (B.12)At each step of blocking transformation, c0n−1 is computed and the process isrepeated until n′ = 2. During the process, the value of c0n−1 increases until thefixed point is reached. At the fixed point, the value of c0n−1 remains constant withinfluctuations and the constant value is the estimation for σ2(m). In the case that itis not reached a constant value, the maximum value which is calculated for c0n−1 isconsidered as a lower bound for σ2(m) and it is a sign that the system is not in theequilibrium.The block transformation was used to estimate standard deviation of the aver-aged force on each DNA molecule in the molecular dynamics simulations. Fig. B.1shows the estimates for σ in the implicit solvent model of cylindrical model withthe surface-surface separation of d = 8A˚, running for t = 30ns (see Chapter 3 for99the details). It presents that after 8 block transformations, the estimates for σ re-main constant within error bars.100