Theoretical Study of Structural Changes in DNA under High External

Article pubs.acs.org/JPCB

Theoretical Study of Structural Changes in DNA under High External Hydrostatic Pressure P. Sudheer Kumar, Arnab Mukherjee,* and Anirban Hazra* Department of Chemistry, Indian Institute of Science Education and Research, Dr. Homi Bhabha Road, Pune, Maharashtra 411 008, India S Supporting Information *

ABSTRACT: The study of DNA under high hydrostatic pressure provides fundamental insights into the nature of interactions responsible for its structure and its remarkable stability in extreme conditions. We have investigated the structural changes in DNA under 2000 bar external pressure using electronic structure calculations and molecular dynamics simulations. Both these methods predict very small distortions in the structure; notably, the change in hydrogen bond lengths is an order of magnitude smaller than previously reported experimental values using NMR. The large discrepancy suggests further investigation into the analysis of the experimental data obtained from NMR.

■

INTRODUCTION DNA can undergo structural changes, sometimes leading to mutations, due to a variety of environmental factors which can either be chemical or physical in nature. Chemical factors include interaction with small molecules or enzymes leading to uncoiling of DNA or breaking of hydrogen bonds,1 and mutagens like reactive oxygen species which change the chemical composition of the DNA.2 Physical factors can be heat that can break the hydrogen bonds between the base pairs leading to strand separation in DNA.3 It can also be UV light which can cause chemical degradation in the bases4 or high external pressure which can exert a large enough force on the DNA molecule to distort its structure.5 Among the external factors affecting the structure of DNA, the effect of high pressures on DNA is relatively unexplored. In fact, compared to other biological macromolecules like proteins, there are fewer studies on the effect of pressure on DNA. The recent progress in the area of high-pressure studies of biomolecules has been described comprehensively in a justpublished review.6 The relevance of high-pressure studies on DNA is 2-fold: First, they give insight into the fundamental nature of forces within the DNA and the surrounding water hydration layers. High pressures provide a handle of distorting the molecular structure in a controlled manner (without chemical change or other external perturbations) to obtain information about different interatomic interactions, for example, those stabilizing a specific molecular conformation. This tool has been very effectively used to understand the mechanism of conformational change between the B- and Zforms of DNA by Barciszewski et al.7 These authors found that the effect of pressure on water structure was an important driving force for the conformational change. They related the © XXXX American Chemical Society

decreased compressibility and increased density of the hydration water around DNA, to the idea that similar changes in water properties can be effected by an increase in pressure.7,8 Second, high-pressure studies are important to understand the biology of organisms which live in the deep oceans under high hydrostatic pressures. Such organisms are called barophiles or piezophiles and have been found to survive even at the deepest part of the ocean with depths in excess of 10 000 m.9 Interestingly, certain barophilic bacteria have been isolated that can only grow under hydrostatic pressure greater than 500 bar. Understanding the survival and adaptation mechanisms of these organisms are very exciting questions.9,10 The first high-resolution NMR study of the structural changes in B-DNA under pressure was carried out by Wilton et al. in 2008.5 These authors reported that under 200 MPa (2000 bar) pressure the DNA compresses to a small extent, with a reduction in hydrogen-bond (H-bond) length of approximately 0.11 Å for the G−C base pairs and 0.29 Å for the A−T base pairs. This study is in disagreement with the only theoretical predictions made earlier by Norberg and Nilsson in 1996, who found the compression of the hydrogen bonds in the nonterminal base pairs to be around 0.02 Å at 6000 atm pressure with the help of molecular dynamics (MD) simulations.11 The large deviation between the experimental and theoretical values is striking. While the theoretical prediction based on short simulations (∼200 ps) may be responsible for the discrepancy, the analysis tools used in the experimental study may also contribute to the error. In fact, Received: October 25, 2014 Revised: January 14, 2015

A

DOI: 10.1021/jp5107185 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

ing and impractical for a large DNA macromolecule. Therefore, instead of the full DNA, we considered single DNA base pairs to understand the effect of pressure on the H-bond. We calculated the force on a single base pair by estimating an effective area that may be exposed to the external pressure. The distortion in the hydrogen bonds is due to the pressure acting perpendicular to the axis of the DNA. The effective area per base pair is approximated to be the area of a rectangle, the height of which is taken to be 3.4 Å (a typical rise value),12 while the width is taken to be 5.7 Å for the A−T base pair and 6.5 Å for the G−C base pair. The width values were estimated based on the largest separation between atoms in the direction perpendicular to the direction of the hydrogen bonds in the standard structures of the base pairs.13 Figure 1 schematically explains the above choices. Once the area is defined, the force per DNA base pair can be calculated for a given pressure.

obtaining distortions in structure at the subatomic length scale poses a significant challenge to experimental methods and thus theoretical investigation is vital and appropriate as long as the common drawbacks of theoretical calculations (for MD simulations these would be inaccuracy in the force-field and inadequate configurational sampling) are addressed. In the present study, we have used two completely different theoretical investigations to address this discrepancy and to help resolve the effect of pressure on DNA structure. Of the two different approaches, the first involves performing quantum chemical electronic structure calculations, which are highly accurate but have to be used on a truncated system to maintain computational feasibility. The second involves performing classical MD simulations, which although not as accurate as quantum calculations are used to model a realistic solvated DNA system. In comparison with quantum calculations, MD simulations have the added advantage of sampling the phase space and providing accurate thermodynamic observables at a finite temperature. Interestingly, our results show that both the approaches agree in predicting negligible distortion of DNA at a very high hydrostatic pressure.

■

THEORETICAL METHODS In this section, we describe in detail the two theoretical approaches we have used to study the effect of pressure on the DNA structure. Electronic Structure Calculations. Quantum mechanical electronic structure calculations involve solving the Schrödinger equation for the electrons in a molecule for a given nuclear geometry. This provides us with accurate energy for a given structure. It is not straightforward to use electronic structure calculations to study pressure effects on geometry because there is no way to specify a thermodynamic quantity like pressure in the definition of the molecular Hamiltonian. We present here a way of estimating the effect of pressure on structural distortion using quantum chemical calculations. Suppose that q is a geometrical parameter in the molecule (like a distance or angle) which describes the molecule’s stretching or compression from its quantum mechanical optimized structure. As convention, consider positive q to correspond to stretching and negative q to compression, with q = 0 corresponding to the optimized structure. The quantum mechanical energy E is calculated for different q values and by interpolation between the q values, E(q) can be obtained as a function of q. Clearly E(q) has a minimum at q = 0 corresponding to the optimized structure. By taking the negative derivative of E with respect to q, we obtain the restoring force = −((∂E)/(∂q)) that the molecule exerts for a distortion of q from q = 0. We then have a one-to-one mapping between the restoring force and q. When subject to external pressure, the equilibrium geometry of the molecule is calculated to be that where the external force corresponding to the pressure equals the restoring force. Thus, to find the geometry corresponding to the desired external pressure, we first calculate the external force and subsequently use the mapping of restoring force and q to find the distorted geometry. We have used this general idea to calculate the distortion in the hydrogen bonds in a base pair of the DNA when subject to 2000 bar pressure. The specifics of our method are described below. To obtain E(q) for different q, quantum mechanical calculations have to be performed for a number of molecular geometries. Each such calculation is computationally demand-

Figure 1. Effective area per base pair in B-DNA approximated as a rectangle with height equal to a typical rise value and width estimated as the largest separation between atoms in the direction perpendicular to the direction of the hydrogen bonds in A−T and G−C.

For obtaining energies of the base pairs at different geometries, electronic structure calculations were performed at the Møller−Plesset second order perturbation (MP2) level of theory and the 6-31++G (d, p) basis set. The GAMESS program package was used.14 First, the equilibrium structures of the A−T and G−C dimers were calculated. Then to calculate the energy profile, the bases were translated with respect to each other from the equilibrium position along the line connecting the C6−C4 atoms in the case of A−T and the line connecting the N1−N3 atoms in the case of G−C as shown in Figure 2. Moving the bases closer to each other corresponds to

Figure 2. Double headed arrows indicate the direction along which the base pairs are translated while calculating the energy profile E(q) as a function of expansion, while the opposite direction corresponds to compression.

compression and moving them apart corresponds to expansion. Both of these changes lead to an increase in energy. However, only the compression is relevant for the present high pressure study. In calculating the energy profile, we implicitly assume that the hydrogen bonds are linear, i.e., the H-bond donor, hydrogen atom and H-bond acceptor lie along a straight line (see justification in the Results). Furthermore, we assume that the major change in structure of the dimer is the change in the length of the interbase hydrogen bonds, while the covalent bonds and other structural parameters like angles are not B

DOI: 10.1021/jp5107185 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Table 1. H-Bond Lengths of the A−T and G−C Base Pairs at Their Optimized Geometries Using Electronic Structure Calculations (Labeled 0 bar) and the Bond Lengths with 1 and 2000 bar Pressuresa Base pairs

H-bonds

0 bar (Å)

1 bar (Å)

2000 bar (Å)

Compression (Å)

A−T

N6H···O4 N1···HN3 O6···HN4 N1H···N3 N2H···O2

1.951 15 1.803 18 1.779 25 1.906 69 1.917 57

1.951 12 1.803 15 1.779 18 1.906 61 1.917 49

1.935 24 1.787 21 1.767 36 1.894 69 1.905 59

0.016 0.016 0.012 0.012 0.012

G−C

a

The change in bond lengths between 1 and 2000 bar pressures is shown in the last column.

respectively. All covalent bonds were restrained using the LINCS algorithm,25 allowing for stable simulations with 2 fs integration time steps. The DNA atom positions were saved every 20 fs for analysis of the H-bond lengths and other structural property calculations. All trajectory information (position and velocities of every atom) was saved every 10 ps, and this was used for the analysis of the H-bond lengths in water. As discussed later, a smaller time sampling for the water is good enough because of the large number of water molecules as compared to DNA.

affected. The energy change from the equilibrium structure is a combined effect of the distortion of all the hydrogen bonds in the dimer. These assumptions, which are quite reasonable, give an unambiguous one-to-one mapping between the energy and structure of the dimer which allows us to predict the structure for a particular external force. From the energy profile, a force profile is calculated as the negative derivative of the energy with respect to displacement by doing a finite difference calculation. The force obtained is due to simultaneous distortions of all the hydrogen bonds in a base pair, rather than individual ones. Then by using this correspondence profile between the force and geometry, the displacement corresponding to the force estimated for the 1 and 2000 bar pressures is obtained. Molecular Dynamics Simulations. Unlike quantum calculations, MD simulations allow direct calculations of the effect of external pressure on the H-bond lengths of the base pairs. The DNA molecule used in this study is the B-DNA dodecamer with sequence d(CGCGAATTCGCG)2 . Its structure was taken from the RCSB protein data bank corresponding to the PDB ID, 1BNA.13 Here we performed MD simulations of the B-DNA dodecamer in aqueous solution at 1 and 2000 bar pressure to obtain the average hydrogen bond lengths at these two pressures. Simulations were carried out with periodic boundary conditions within a cubic box using the GROMACS15 program package with AMBER99sb force field16 with parmbsc0 modifications.17 Water was modeled using the TIP3P parameters.18 The simulations were carried out with 758 DNA atoms solvated by a 1.2 nm water layer (6,798 water molecules) to avoid interaction of DNA with its periodic image. Na+ and Cl− ions were included to ensure zero net charge and physiological ion concentration (150 mM). The total number of atoms in the system thus was 21 212. Electrostatic interactions were treated using the particle mesh Ewald method19 with a 1 nm real-space cutoff and cubic spline interpolation onto the charge grid with spacing of 0.16 nm. The cutoff used for the Lennard-Jones interactions was 1 nm. Initial equilibration was performed using a standard protocol20 as briefly mentioned below. Energy minimization was carried out first. This was followed by heating with position restraints on the heavy atoms of the DNA to 300 K at constant pressure of 1 bar using the Berendsen barostat21 and the velocity-rescale thermostat.22 Subsequently, a series of minimization and position restraint simulations at constant temperature of 300 K and pressure 1 bar was performed, again using the velocity-rescale thermostat22 and Berendsen barostat21 with a decrease in the restraint force at each step. Finally, the system was simulated for 1 ns without any restraint. This equilibrated configuration was taken for further simulations. Three independent production simulations of 99 ns, each at 1 and 2000 bar pressures were carried out using an NPT ensemble with the Nose−Hoover thermostat23 and Parrinello− Rahman barostat24 with coupling constants of 0.4 and 1.0 ps,

■

RESULTS Electronic Structure Calculations. There are five different types of hydrogen bonds in the base pairs: An A−T base pair has two types of hydrogen bonds while a G−C base pair has three types. These are denoted as AT(N6H···O4), AT(N1··· HN3), GC(O6···HN4), GC(N1H···N3), and GC(N2H···O2) using the standard atom numbering as shown in Figure 2. The H-bond lengths of the A−T and G−C stable structures for different pressures are presented in Table 1, while the distance between the corresponding H-bond donor and acceptor are presented in Table 2. The calculated H-bond lengths and Table 2. Distances between H-Bond Donor and Acceptor of A−T and G−C Base Pairs at Their Optimized Geometries Using Electronic Structure Calculations (Labeled 0 bar) and the Distances at 1 and 2000 bar Pressuresa Base pairs A−T G−C

D−A

0 bar (Å)

1 bar (Å)

2000 bar (Å)

Compression (Å)

N6−O4 N3−N1 N4−O6 N1−N3 N2−O2

2.966 30 2.847 39 2.809 37 2.937 88 2.934 28

2.966 27 2.847 35 2.809 30 2.937 81 2.934 20

2.950 35 2.831 42 2.797 49 2.925 89 2.922 34

0.016 0.016 0.012 0.012 0.012

a

The change in bond lengths between 1 and 2000 bar pressures is shown in the last column.

donor−acceptor distances in the absence of pressure are in agreement within 0.1 Å with previous theoretical studies on the DNA base pairs.26 They are also in qualitative agreement with X-ray crystallographic studies.27 The small discrepancies between theoretical prediction and the X-ray structures have been analyzed in detail by Guerra et al.26a The angle between the H-bond donor, hydrogen bond and H-bond acceptor of the five hydrogen bonds is between 173.2 and 179.8° (Table S2 in Supporting Information), thus justifying the linearity assumption discussed in the methods section. The calculated energy profiles as a function of the distortion are shown in Figure 3a where the reference energy is equal to the energy of the undistorted structure. The corresponding restoring forces as a function of distortion are shown in Figure C

DOI: 10.1021/jp5107185 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

(three for the GC and two for the AT) in all the base pairs of the dodecamer. While averaging, the terminal base pairs, which happen to be G−C at both ends, are not considered because these bases exhibit large fluctuations as compared to the nonterminal bases (Figure S1 in Supporting Information). Collecting data over a longer time interval than 20 fs does not change the result significantly (Figure S2 in Supporting Information). With longer simulation time, the sampling of the phase-space increases and the average bond length can be expected to converge. We have performed three independent long (100 ns) simulations to calculate the error in the average H-bond lengths and other DNA parameters. The average hydrogen bond distances in DNA for 1 and 2000 bar pressure for the three independent simulations are shown in Figures 4 and 5 respectively. Although the simulations

Figure 3. (a) Energy calculated using the MP2/6-31++G (d, p) method as a function of the distortion coordinate q. (b) Corresponding force as a function of the distortion coordinate q.

3b. The potential energy curve is approximately harmonic for small distortions. Assuming each of the bases to be rigid point masses attached by a spring, the frequencies calculated for the A−T and G−C vibrations are 132 and 168 cm−1 respectively. Interestingly, our calculated normal-mode frequencies corresponding to the intermolecular stretching vibration of the A−T and G−C pairs are 107 and 127 cm−1. The similarity in the values of the frequency corresponding to the one-dimensional energy profile and the normal-mode frequency indicates a resemblance between these two coordinates and validates the choice of the distortion coordinate for calculating the energy profile. The slightly greater frequencies, and thus force constants, in the one-dimensional case are due to the nuclear degrees of freedom other than the distance between the base pairs not being allowed to relax during the distortion. For large distortions, as expected, there is anharmonicity corresponding to a greater increase in energy for a compressive distortion as compared to an equivalent expansive distortion. To find the distortion in structure under high pressure, the external force was calculated using the procedure of estimating effective area described earlier. For 1 bar (2000 bar values in parentheses) pressure the calculated force was 3.88 × 10−14 N (7.75 × 10−11 N) for the A−T base pair while it was 4.42 × 10−14 N (8.84 × 10−11 N) for the G−C base pair. Using Figure 3b, the change in H-bond length was obtained for the different forces and is presented in Tables 1 and 2. For 1 bar pressure, the change in H-bond length with respect to the structure with zero force was of the order of 10−5 Å and is negligible. For 2000 bar pressure, the change in H-bond length was approximately 0.016 Å for the A−T base pair whereas it was approximately 0.012 Å for the G−C base pair. The change in the distance between the H-bond donor and acceptors is similar to the change in the H-bond lengths. Molecular Dynamics Simulations. MD simulations corroborate the results obtained above. We calculated average H-bond lengths and other DNA parameters under normal and high pressure to observe its effect on DNA. The simulated DNA dodecamer contains 8 G−C base pairs and 4 A−T base pairs. Apart from time averaging for each of the H-bonds at every 20 fs, we also averaged over the same type of H-bonds

Figure 4. Cumulative running average of the hydrogen bond lengths of DNA base pairs at 1 bar pressure as a function of time for simulation 1 (black line), simulation 2 (red dotted line), and simulation 3 (green dashed line). Note that the range on the hydrogen bond length axis is 0.02 Å for all the panels except AT (N6H···O4) for which the range is 0.03 Å.

Figure 5. Cumulative running average of the hydrogen bond lengths of DNA base pairs at 2000 bar pressure as a function of time for simulation 1 (black line), simulation 2 (red dotted line), and simulation 3 (green dashed line). Note that the range on the hydrogen bond length axis is 0.02 Å for all the panels except the AT (N6H···O4) for which the range is 0.07 Å. D

DOI: 10.1021/jp5107185 J. Phys. Chem. B XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry B

Table 3. Average H-Bond Lengths at 1 and 2000 bar Pressure and Compressions Averaged over All Three Independent Simulationsa Base pairs A−T G−C

H-bonds N6H···O4 N1···HN3 O6···HN4 N1H···N3 N2H···O2

water (all) water (0.8) has lower value at high pressure (blue line) compared to that at normal pressure (red line). For hydration water, distribution of Q has lower value in the low Q range (