(515h) In silico Studies of the solvation Effects As the Main Stabilizer of the Double Helix Structure of DNA

Feng et al. (2019) recently suggested that the DNA stability was due to the hydrophobic cohesion of the base pairs’ coin-pile stacking, requiring abundant water. The hydrophobic stacking indirectly makes the DNA interior dry, making the hydrogen bonding exerts its full potential. Their findings show that introducing PEG-400 and diglyme will significantly affect the DNA’s stability, as demonstrated by their experimental results. This suggests that hydrophobic base stacking is the main contributor to DNA stability, not hydrogen bonding. However, there are difficulties in the theoretical description of stacking associations in the condensed phase and in understanding the stacking role in biomolecules. The solvation thermodynamics analyses were not conducted; thus, the behavior of DNA molecules at different solvent conditions was not determined.

This study aims to elucidate the fundamental hydrophobic interactions between DNA and water with the influence of cosolvents like EG using molecular dynamics (MD) simulations and the 3-dimensional reference interaction site model (3D-RISM) theory.

All system preparations were done using GROMACS. The 40-bp DNA was prepared in a triclinic PBC with 90° angle sides. The PBC was set at a 1.0 nm distance. The forcefield used for the DNA structure is parmbsc1 (Ivani et al., 2016). The DNA structure was then solvated in TIP3P water (Price & Brooks, 2004). After solvation, sodium (Na⁺) and chloride (Cl^-) ions were introduced for system neutralization.

The Avogadro software (Hanwell et al., 2012) was used to construct the EG PDB file. The molecular structure was optimized using the General AMBER Force Field (J. Wang et al., 2004; Zhang et al., 2018). Structural and topology parameters were generated using the Antechamber and acpype tools (Sousa da Silva & Vranken, 2012). GAFF was used during Antechamber, and the AMBER Force field format was used in acpype to generate the structural and topology files. The partial charge calculation used was Gasteiger (Mortier et al., 1985). The generated files were then added to the topology files of DNA structure for water-cosolvent system runs. EG molecules were added to the system at different concentrations. The number of EG in the system was based on the data density data obtained (Egorov et al., 2010; Lago et al., 2009; Trivedi et al., 2010).

A Verlet cut-off scheme was used for the particle neighbor searching for potential interactions. For the electrostatic forces, Particle Mesh Ewald (PME) was used (Essmann et al., 1995). Both electrostatic and van der Waal (vdW) cut-offs were 1.0 nm. For the minimization algorithm, the steepest descent was used. All NVT and NPT equilibrations were done at 310.15 K. LINCS (Hess et al., 1997) was used for the bond constraints, while Particle Mesh Ewald (PME) (Essmann et al., 1995) with an order of four (4) was used for the long-range electrostatic interactions. The DNA was positionally restrained during the NVT ensemble for 100 ps simulation. A velocity-rescale thermostat (Bussi et al., 2007) was used during the NVT equilibration, while a Berendsen barostat was used for the NPT equilibration. The NPT equilibration was done in a 100 ps simulation with DNA positionally restrained. Molecular dynamics (MD) simulations were done for 100 ns with a 2-fs time step. All energy and trajectory data were recorded every 100 ps. All simulations were done in triplicates.

Clustered files were obtained in GROMACS using the method of GROMOS (Daura et al., 1999) with a cut-off value of 0.2 nm. Only the last ten (10) ns simulations were used for clustering. Only the representative of cluster 1 as the clustered structure was used for energetics analysis. The representative structures and corresponding trajectories were converted to AMBER-readable files using the tleap and cpptraj command tools.

For the 1D-RISM, dielectric constant values of the EG-water mixtures were calculated using the method of Wang and Anderko (2001). The physical properties of both water and EG were obtained (Egorov et al., 2010; Malmberg & Maryott, 1956; Perry & Green, 2008; Sengwa, 2003). For the grid properties, the number of grids was set to 16,384 with 0.025 Å grid spacing, while the buffer was set to 10 Å (Case et al., 2021, p. 119). The MDIIS (Kovalenko et al., 1999) number of vectors was set to 20, step size of 0.3 and tolerance of 1.0×10^−12. For the solvent site of EG, the molecule was represented as a coarse-grained (CG) model presented by Go et al. (2019) instead of an all-atomistic (AA) model to have a simulation convergence.

The 3D-RISM was performed using the rism3d.snglpnt of AmberTools21 (Case et al., 2021). For the MDIIS automatic setup, the step size was set to 0.7, the number of vectors was 5, the number of restart vectors was 10, and the tolerance was 1.0×10^−5. The closure function used was KH. The buffer was set to 10 Å.

Figure 1A to E shows the conformer of the 40-bp DNA at different EG mass concentrations at the last ten (10) ns of MD simulation. There is no significant difference between the DNA structures at 0%, 10%, and 20% EG concentrations. Significant distortions exist at the DNA structures when solvated at 30% and 40% EG concentrations, especially at the middle base pairs. The EG layer was present in the middle part of the DNA structure. The molecular crowding effect induced this EG layer.

Figure 1F shows the average Helmholtz free energy difference as a solute and excess chemical potential energy contribution. The figure shows the downward trend of the Helmholtz free energy as the EG mass concentration increases (-403.4853 to -2127.1421). The decrease in the free energy means that the solute molecule in the water-ions-EG solvent system was spontaneous. The decrease in the Helmholtz free energy difference was mainly due to the contribution of the excess chemical potential (-415.7741 to -2368.5599 kcal/mol). There was little contribution of the solute potential energy to the Helmholtz free energy difference. There was a slight increase in the solute potential energy difference at 30 and 40% EG mass concentration (-47.4962 to 241.4177 kcal/mol) due to the base-pair flipping at these concentration values at the EG-crowded layer. Nevertheless, the solute potential energy increase was not enough to overcome the excess chemical potential, thus, making the DNA molecule stay solvated.

Figure 1G shows the average excess chemical potential contribution as the sum of solvation energy and entropy. The figure shows entropy as the main contributor to EG mass concentrations of 10 and 20% (-475.8239 to -928.9530 kcal/mol). The decrease in entropy means that the entropy of the solution increases with increasing EG molecules. The excess chemical potential difference decreased at 30 and 40% EG mass concentration. This layer formed the solvent in an ordered state, reducing the entropic penalty. However, the solvation energy was the main contributor (-749.1495 to -1313.5432 kcal/mol) to the decrease of excess chemical potential rather than the entropic penalty. In order to overcome the randomness of the arrangement of the solvent system, the solvation energy must compensate to maintain the decreasing excess chemical potential difference values.

Figure 1H shows that the solvation energy difference term can be broken down as a sum of solute-solvent interaction potential energy and solvent reorganization energy. At 10 to 20% EG mass concentrations, the solute-solvent interaction energy difference increased (738.9582 and 1687.6058 kcal/mol for 10% and 20%, respectively) while the solvent reorganization energy decreased (-678.9804 to -1530.3733 kcal/mol). As the water molecule decreased, the interaction between the water and DNA also decreased, thus, increasing the interaction potential energy value. The decrease in the solvent reorganization energy means that the water-ions-EG solvent system requires less energy to reorganize itself for the accommodation of DNA as the EG mass concentration increases up to 20%. The solute-solvent interaction potential energy difference suddenly decreased at 30% EG mass concentration and had a downward trend until 40% EG mass concentration (-734.0113 to -1052.5248 kcal/mol). The solvent reorganization energy difference had decreased at 40% EG mass concentration (-261.0184 kcal/mol for 40%). Consequently, there was a sudden increase in the solvent reorganization energy difference as the solvent system required more energy to create a cavity and to form an EG-crowded layer (-15.1382 kcal/mol for 30%).

These results show that the stability of the DNA molecule is based on the favorable excess chemical potential, i.e. solvation effects. Breaking down the excess chemical potential showed that polar contribution of water was responsible for the significantly favorable solvation entropy and solvation energy with the solute-solvent interaction energy as the main contributor. A decrease in the number of available water molecules resulted in the denaturation if the DNA, further showing the importance of hydrophobicity in the stability of the DNA structure.