We study the following systems using all-atom molecular dynamics simulations: (1) EcoRI and DNA containing cognate sequence (EcoRI-DNA), (2) EcoRI-DNA in presence of Mg2+ ions (Mg2+-EcoRI-DNA), (3) EcoRI-DNA in presence of Ca2+ ions (Ca2+-EcoRI-DNA) (4) EcoRI with mutated DNA where A5:T20 base pair in the DNA is mutated by G5:C20, in presence of Mg2+ ions (Mg2+-EcoRI-mtDNA) and (5) Free DNA containing the cognate sequence.
The co-ordinations of Mg2+ in Mg2+-EcoRI-DNA complex, Ca2+ in Ca2+-EcoRI-DNA and Mg2+ in Mg2+-EcoRI-mtDNA obtained from MD trajectories are shown in Figs. 1(b)-(d) respectively. As the co-ordinations of the divalent metal ions in the two active sites in each complex are identical, we only show those of the divalent metal ion in the first active site. If the average distance between any heavy atom of DNA or water or EcoRI and metal ion is less or equal to 3.0 Å, coordination is considered. OD1 of ASP75, OE1 of GLU95 and OP1 of P5 participate in coordination with Mg2+. Similarly, OD1 and OD2 of ASP75, OE1 and OE2 of GLU95 and OP1 of P5 participate in coordination with Ca2+. The co-ordinations of Mg2+ in Mg2+-EcoRI-mtDNA are same as those of Mg2+ in Mg2+-EcoRI-DNA. The rest of the co-ordinations of Mg2+ and Ca2+ are satisfied by the O atoms of water molecules.
Equilibrium hydration pattern
We calculate the water density around a given chemical moiety. First, we consider the metal ions. Figure 2(a) shows the density of water, \({{\rho }}^{\text{M}\text{g}}\)(r) and \({{\rho }}^{\text{C}\text{a}}\)(r) around Mg2+ in Mg2+-EcoRI-DNA complex and Ca2+ in Ca2+-EcoRI-DNA complex respectively, where r is distance between the metal ion and water oxygen. For all cases, the distributions show two well defined peaks. The first peak in \({{\rho }}^{\text{M}\text{g}}\)(r) appears at r = 0.22 nm, while the second peak is at r = 0.29 nm. The density is observed to be maximum of 3.1 gm/cc at 0.22 nm. As the radius of Ca2+ is greater than that of Mg2+, the peaks in \({{\rho }}^{\text{C}\text{a}}\)(r) is further away from Ca2+ than those in \({{\rho }}^{\text{M}\text{g}}\)(r). The first peak arises at 0.25 nm, while the second peak at 0.32 nm. Here, the maximum density is 1.98 gm/cc which occurs at 0.25 nm. Thus, Mg2+ is better hydrated than Ca2+ due to smaller size and consequently, larger charge density.
Next, we consider the hydration of the phosphate groups of the DNA in various cases. We calculate the density of water molecules, \({{\rho }}_{{\alpha }}^{0}\left(r\right)\), \({{\rho }}_{{\alpha }}^{\text{c}}\left(r\right)\), \({{\rho }}_{{\alpha }}^{\text{M}\text{g}}\left(r\right)\), \({{\rho }}_{{\alpha }}^{\text{C}\text{a}}\left(r\right)\) and \({{\rho }}_{{\alpha }}^{Mut}\left(r\right)\) around phosphate group (\(\alpha\)) in free DNA, EcoRI-DNA, Mg2+-EcoRI-DNA, Ca2+-EcoRI-DNA and Mg2+-EcoRI-mtDNA complexes respectively where \(r\) is the distance between phosphorus atom and the oxygen of water. We compare the cases for α = P5 in the cleavage region and α = P7 in the non-cleavage region in Figs. 2(b) and (c) respectively. \({{\rho }}_{\text{P}5}^{0}\left(r\right)\) in the free DNA shows a peak at 0.4 nm with maximum density 0.82 gm/cc, while the other peak at 0.6 nm has similar maximum density. \({{\rho }}_{\text{P}5}^{\text{c}}\left(r\right),\) in EcoRI-DNA complex without any metal ion shows a well-defined maximum at 0.43 nm, while it shows another broad peak around 0.48 nm. Similarly, \({{\rho }}_{\text{P}5}^{\text{M}\text{g}}\left(r\right)\) shows two well defined maxima. The position of the first maximum in \({{\rho }}_{\text{P}5}^{\text{M}\text{g}}\left(r\right)\) is smaller than that in \({{\rho }}_{\text{P}5}^{0}\left(r\right)\). On the other hand, the first peak in \({{\rho }}_{\text{P}5}^{\text{C}\text{a}}\left(r\right)\) is comparable to that in \({{\rho }}_{\text{P}5}^{0}\left(r\right)\). This indicates that P5 is more hydrated in Mg2+-EcoRI-DNA complex than both EcoRI-DNA and Ca2+-EcoRI-DNA due to larger hydration around Mg2+ ion. On the other hand, the density of water around P7 show two different peaks at 0.3 nm and 0.4 nm for all the cases. Thus, the hydration level of the phosphate groups at the cleavage site is sensitive to the type of metal ions but not so at the non-cleavage site.
We calculate the fraction of water molecules in the first coordination shell, defined by the first peak in the density profile, having hydrogen bonds with O atoms of the phosphate groups in the recognition region. We denote the fraction by \({f}_{{\alpha }}^{0}\), \({f}_{{\alpha }}^{\text{c}}\), \({f}_{{\alpha }}^{\text{M}\text{g}}\) and \({f}_{{\alpha }}^{\text{C}\text{a}}\) in free DNA, EcoRI-DNA, Mg2+-EcoRI-DNA and Ca2+-EcoRI-DNA respectively. As the two strands are symmetric, the results shown in Fig. 2(d) in the recognition region are averaged over both strands. For free DNA, the values of \({f}_{{\alpha }}^{0}\) for all the phosphate groups in recognition region are similar. However, for the complexes, the results show sequence specificity. In EcoRI-DNA complex, \({f}_{\text{P}5}^{\text{c}}\) shows a minimum value. The values of \({f}_{\alpha }^{\text{c}}\) for the other phosphate groups are similar to those in free DNA. In Mg2+-EcoRI-DNA complex, the minimum of \({f}_{\text{P}5}^{\text{M}\text{g}}\) is even deeper and a maximum at P7 is observed unlike in the EcoRI-DNA complex. In Ca2+-EcoRI-DNA complex, \({f}_{\text{P}5}^{\text{C}\text{a}}\) has a minimum comparable to that in EcoRI-DNA complex. The maximum of \({f}_{\text{P}7}^{\text{C}\text{a}}\) is similar to that of \({f}_{\text{P}7}^{\text{M}\text{g}}\). Thus, the data suggest that the hydrogen bond formation is specific to the kind of the metal ion. The number of hydrogen bonded water molecules at the cleavage site is the least in case of Mg2+-EcoRI-DNA complex.
We compare the fluctuations of base pair step parameters in EcoRI-DNA, Mg2+-EcoRI-DNA and Ca2+-EcoRI-DNA complexes with those in free DNA. A given Watson-Crick paired bases is taken to define a base pair. The base pair step parameters define the relative displacement (shift, slide and rise) and orientations (tilt, roll and twist) of two successive base pairs.[22] For a given base pair step (i, i + 1), the distribution is denoted by \({\text{H}}_{\text{i},\text{i}+1}^{\text{f}}\)(b) in the free state and \({\text{H}}_{\text{i},\text{i}+1}^{\text{c}}\)(b) in the complex for a given base pair step parameter b where i and i + 1 are two successive bases in a given strand, omitting the WC paired bases. The distributions in free DNA, Mg2+-EcoRI-DNA and Ca2+-EcoRI-DNA complexes are denoted as \({\text{H}}_{\text{i},\text{i}+1}^{0}\)(b), \({\text{H}}_{\text{i},\text{i}+1}^{\text{M}\text{g}}\)(b) and \({\text{H}}_{\text{i},\text{i}+1}^{\text{C}\text{a}}\)(b) respectively. The base pair parameters are taken to be conformational variables for the DNA in our previous analysis [2].
Here, we illustrate the data for base pair step, (A5:T20, A6:T19) in the cleavage region. We show the histograms in Supplementary Information (SI), Figs. S1(a)-(f). τ for all cases are broad single peaked (see SI, Fig. S1(a)). However, the peak heights of \({\text{H}}_{\text{5,6}}^{0}\)(τ) and \({\text{H}}_{5,6}^{\text{C}\text{a}}\)(τ) are smaller than that of \({\text{H}}_{5,6}^{\text{M}\text{g}}\)(τ). The distributions of ρ (see SI, Fig. S1(b)) and ω (see SI, Fig. S1(c)) are all single peaked in all the cases. However, the peak height is smaller and broader in \({\text{H}}_{5,6}^{0}\)(ω) compared to the other complexes. Although, Dx shows single peaked distributions for all cases, the peak height for \({\text{H}}_{5,6}^{\text{M}\text{g}}\)(Dx) is larger than that of other complexes (see SI, Fig. S1(d)). Similarly, the distributions of Dy (see SI, Fig. S1 (e)) and Dz (see SI, Fig. S1(f)) are single peaked in all complexes. However, the peak height of the distributions of Dy is highest for \({\text{H}}_{5,6}^{\text{M}\text{g}}\)(Dy). Thus, we observe that the distributions of rotational parameters of the base pair steps are broader in all cases than those of translational parameters.
We quantify the fluctuations in terms of the changes in conformational free energy and entropy of DNA in complexes with respect to the free DNA.[23] A negative change in the free energy denotes stabilization and that in entropy denotes enhanced order in the complexed DNA with respect to the reference free DNA. For a given step parameter, the changes in free energy and entropy of (\(\text{i}, \text{i}+1\)) step in the Mg2+-EcoRI-DNA complex with respect to the free DNA are denoted as \({{\Delta }\text{G}}_{\text{i},\text{i}+1}^{\text{M}\text{g}}\)(b) and \({\text{T}{\Delta }\text{S}}_{\text{i},\text{i}+1}^{\text{M}\text{g}}\)(b) respectively. Similarly, the analogous quantities for Ca2+-EcoRI-DNA complex with respect to the free DNA are denoted by \({{\Delta }\text{G}}_{\text{i},\text{i}+1}^{\text{C}\text{a}}\)(b) and \({\text{T}{\Delta }\text{S}}_{\text{i},\text{i}+1}^{\text{C}\text{a}}\)(b). The total free energy (\({{\Delta }\text{G}}_{\text{i},\text{i}+1}^{\text{M}\text{g}}\) and \({{\Delta }\text{G}}_{\text{i},\text{i}+1}^{\text{C}\text{a}}\)) and entropy changes (\({\text{T}{\Delta }\text{S}}_{\text{i},\text{i}+1}^{\text{M}\text{g}}\) and \({\text{T}{\Delta }\text{S}}_{\text{i},\text{i}+1}^{\text{C}\text{a}}\)) of a base pair step parameters of DNA in the complexes with respect to those in free DNA, adding the individual contributions, are shown in Figs. 2(e) and (f) respectively. The base pair steps are highly stabilized in Mg2+-EcoRI-DNA. Figure 2(e) shows that the base pair step \({{\Delta }\text{G}}_{\text{4,5}}^{\text{M}\text{g}}\) in cleavage region of Mg2+-EcoRI-DNA shows maximum conformational stabilization. Figure 2(f) shows that \({\text{T}{\Delta }\text{S}}_{\text{4,5}}^{\text{M}\text{g}}\) has maximum order compared to the other base pair steps and that at \({\text{T}{\Delta }\text{S}}_{\text{6,7}}^{\text{M}\text{g}}\) pair has a minimum. In Ca2+-EcoRI-DNA complex (Fig. 2(e)), the stabilization in \({{\Delta }\text{G}}_{\text{4,5}}^{\text{C}\text{a}}\) is maximum but it is lower compared to that in Mg2+-EcoRI-DNA. Here, \({{\Delta }\text{G}}_{\text{6,7}}^{\text{C}\text{a}}\) shows destabilization. The base pair steps in Ca2+-EcoRI-DNA are all ordered with a similar value of entropy change − 3.0 kJ/mol, except \({\text{T}{\Delta }\text{S}}_{\text{6,7}}^{\text{C}\text{a}}\) (Fig. 2(f)).
Next, we consider a mutated DNA sequence, denoted by Mg2+-EcoRI-mtDNA where we mutate the cleavage base pair A-T in the first active site by G-C in Mg2+-EcoRI-DNA complex. First, we show the density profile of water molecules, \({{\rho }}_{\text{P}5}^{\text{M}\text{u}\text{t}}\) and \({{\rho }}_{\text{P}7}^{\text{M}\text{u}\text{t}}\) in the first strand of the mutated complex in Figs. 2(b) and (c) respectively. We observe that the density of water around P5 and P7 in the mutated complex are similar to those in Mg2+-EcoRI-DNA. The fraction of hydrogen bonds, \({\text{f}}_{{\alpha }}^{\text{M}\text{u}\text{t}}\) for different \({\alpha }\) of Mg2+-EcoRI-mtDNA is shown in Fig. 2(d). It is observed that \({\text{f}}_{{\alpha }}^{\text{M}\text{u}\text{t}}\) has a minimum at P5. However, \({\text{f}}_{\text{P}5}^{\text{M}\text{u}\text{t}}\) >\({\text{f}}_{\text{P}5}^{\text{M}\text{g}}\). Thus, the hydrogen bonding network is sequence specific and sensitive to the mutation of the DNA base.
Time dependent hydration pattern
We investigate the changes in hydration pattern, while the DNA gets unbound from the binding pocket of EcoRI. We show the forces \({\text{F}}^{0}\), \({\text{F}}^{\text{M}\text{g}}\) and \({\text{F}}^{\text{C}\text{a}}\) experienced by the systems due to unbinding of DNA from EcoRI for EcoRI-DNA, Mg2+-EcoRI-DNA and Ca2+-EcoRI-DNA complexes respectively in Fig. 3. For all the cases, the forces increase gradually until they reach a maximum, and after that they gradually decrease. The force is maximum when all the nucleotides of the DNA come out from the interaction regime of protein. This is in agreement to an earlier report [21]. The maximum of the force is larger in Mg2+-EcoRI-DNA than all other cases. This indicates that the interactions between EcoRI and DNA is strongest in Mg2+-EcoRI-DNA.
We also follow the unbinding of the nucleotides. As time progresses, the Nβ nucleotide of the DNA fragment unbinds. The distances of the centre of mass of a nucleotide which contains a phosphate group, a sugar group and one of the four bases, with respect to the centre of mass of the binding residues of EcoRI are denoted by \({\text{R}}_{{\beta }}^{0}\), \({\text{R}}_{{\beta }}^{\text{M}\text{g}}\) and \({\text{R}}_{{\beta }}^{\text{C}\text{a}}\) for EcoRI-DNA, Mg2+-EcoRI-DNA and Ca2+-EcoRI-DNA respectively. Here Nβ = N4-N9 for the first strand and Nβ = N16-N21 for the second strand. The data are shown in SI, Fig. S2. The jump in the displacement denote the unbinding of the corresponding nucleotide.
Based on jumps in the distances with time, we show in Table 1 the times at which different nucleotides are liberated from the binding pocket. We first consider EcoRI-DNA, shown in SI, Figs. S2(a) and (b) for the first and second strands respectively. \(\text{N}20\) in the second strand breaks away first from the interaction regime of EcoRI around 1.35 ns. \(\text{N}19\)and \(\text{N}21\) break away at ⁓ 2.05 ns and 2.9 ns respectively. After that around 3.55 ns, N4 in the first strand leaves the binding site. In the range 4-5.5 ns, N5-N7 in the first strand and N18 in the second strand leave the binding site. All the nucleotides in the recognition region of first strand leave around 6.9 ns, while those in the second strand leave around 6.6 ns, yielding the total rapture time 6.9 ns. The sequence of nucleotide rupture is similar to the earlier study.[21] However, the rupture times are somewhat less in our system which may be due to differences in force fields in two cases.
Next, we consider the Mg2+-EcoRI-DNA complex (SI, Figs. S2(c) and (d) for the first and second strands respectively). In EcoRI-DNA, the nucleotide starts breaking away at 1.35 ns, whereas in Mg2+-EcoRI-DNA complex, this time is 3 ns. Both N20 and N21 in the second strand are ruptured at 3 ns. N4 and N19 break away at 3.4 ns and 3.9 ns respectively. N18 breaks away at 4.95 ns. N4-N6 rupture gradually in a gap of ⁓ 1 ns, while both N7 and N8 break away approximately at the same time (⁓ 5.8 ns). N16 and N17 in the second strand leave the binding site around ⁓ 6.5 ns, while the remaining nucleotide N9 in the first strand breaks away around ⁓ 8 ns. Thus, the rupture time in Mg2+-EcoRI-DNA complex is ~ 8 ns. It indicates that the interactions between EcoRI and DNA is stronger in Mg2+-EcoRI-DNA than in EcoRI-DNA.
The sequence of events in Ca2+-EcoRI-DNA complex is different. The data for the first and second strands are shown in SI, Figs. S2 (e) and (f) respectively. In this complex, N4, N20 and N21 leave the binding site approximately at the same time ⁓ 3.7 ns. N19 breaks away at ⁓ 4.05 ns. After that all four nucleotides N5, N6, N7 and N9 in the first strand break away at approximately same time (⁓ 5.1 ns) (SI, Fig. S2(e)). N18 breaks away at ⁓ 5.6 ns, while N8 breaks away around ⁓ 6 ns. Here, N9 is ruptured before N8 (SI, Fig. S2(e)). N16 comes out from the binding pocket around ⁓ 6 ns. N17 breaks away last at ⁓ 6.55 ns, while in Mg2+-EcoRI-DNA the last nucleotide is N9 that breaks away at ⁓ 8 ns, indicating Mg2+-EcoRI-DNA complex is more stable than Ca2+-EcoRI-DNA.
Now we consider the mutated complex (SI, Figs. S2 (g) and (h) for the first and second strands respectively). N21 is ruptured first at around ⁓3.2 ns, while N4 at 3.6 ns in the first strand. In the range 4–6 ns, all the nucleotides except nucleotide 9 in the first strand, and N19 and N20 in the second strand break away. The nucleotides N9, N16-N18 all break away at approximately same time ⁓6.35 ns. The rapture time in Mg2+-EcoRI-mtDNA is ⁓6.35 ns. This indicates higher stability of the Mg2+-EcoRI-DNA complex compared to the Mg2+-EcoRI-mtDNA complex.
We show histogram, H(t) of the breaking times, t of the nucleotides in different cases in Fig. 3(b). The histogram shows three different maxima at 3.5 ns, 5 ns and 6 ns. Hence, we carry out stopped SMD simulations by considering snapshots at these times from the SMD trajectory for all the complexes. In addition to these times, we perform stopped SMD simulations of the snapshots at 4 ns corresponds to the first minimum of the histogram and at 6.5 ns where most of the nucleotides come out from the binding pocket of EcoRI. We calculate the fraction of hydrogen bonds, \({\text{f}}_{\alpha }^{\text{c}}\), \({\text{f}}_{{\alpha }}^{\text{M}\text{g}}\), \({\text{f}}_{{\alpha }}^{\text{C}\text{a}}\) and \({\text{f}}_{{\alpha }}^{\text{M}\text{u}\text{t}}\) at these times, shown in Fig. 4. Let us consider the scissile phosphate P5 in Fig. 4(a). \({\text{f}}_{\text{P}5}^{\text{c}}\) does not show significant changes with time. The values of \({\text{f}}_{\text{P}5}^{\text{C}\text{a}}\) at different times are similar to those of \({\text{f}}_{\text{P}5}^{\text{c}}\). However, \({\text{f}}_{\text{P}5}^{\text{M}\text{g}}\) shows much lower values at lower times compared to those of \({\text{f}}_{\text{P}5}^{\text{c}}\) and approaches \({\text{f}}_{\text{P}5}^{\text{c}}\) at larger times. Similar time dependent growth is observed in \({\text{f}}_{\text{P}5}^{\text{M}\text{u}\text{t}}.\) In both cases the time dependences are linear as shown by the fitted curves with slopes 0.14 and 0.07 for Mg2+-EcoRI-DNA and Mg2+-EcoRI-mtDNA cases respectively. Thus, the growth for the mutated complex is not as rapid as that in the case of the wild type complex. On the other hand, \({\text{f}}_{\text{P}7}^{\text{c}}\), \({\text{f}}_{\text{P}7}^{\text{M}\text{g}}\),\({\text{f}}_{\text{P}7}^{\text{C}\text{a}}\) and \({\text{f}}_{\text{P}7}^{\text{M}\text{u}\text{t}}\) do not show significant changes with time as shown in Fig. 4(b).
It may be noted that the phosphodiester bond cleavage cannot be observed from classical MD simulations. However, our studies based on hydration pattern from classical MD trajectories give significant insight to the role of the metal ion and DNA sequence. The smaller metal ion like Mg2+ is solvated more efficiently than larger ion like Ca2+. This is favourable from the stand point of activation of water molecules for cleavage [6, 24, 25]
The overall fluctuations in the microscopic conformational variables in the cleavage region is less in Mg2+-EcoRI-DNA compared to the Ca2+-EcoRI-DNA. This scenario is quite analogous to that of EcoRV. The enhanced stability and order are pre-requisite for positioning the metal ion in the vicinity of the scissile phosphate and is more efficiently done in Mg2+-EcoRI-DNA than in Ca2+-EcoRI-DNA. It may be noted that the maximum stabilization and ordering in (4,5) step of Mg2+-EcoRI-DNA complex has the presence of lowest number of hydrogen bonded water molecules. On the other hand, the most destabilized and disordered (6,7) step has the maximum number of hydrogen bonded water molecules. Less abundance of H-bonds indicates that the scissile phosphate group in this complex is more accessible to have the nucleophilic attack by a non-hydrogen bonded water molecule. A single base pair mutation in Mg2+-EcoRI-DNA significantly reduces the fraction of hydrogen bonds around the scissile phosphate group.
Time dependent changes in hydrogen bond fraction during the unbinding gives idea about the kinetics of the cleavage process. The increase in the fraction \({\text{f}}_{\text{P}5}^{\text{M}\text{u}\text{t}}\) and \({\text{f}}_{\text{P}5}^{\text{M}\text{g}}\) while the unbinding may also be taken to indicate decreases in the fraction as the DNA enters into the binding pocket. Thus, the slower rate of increase in mutated complex than that of Mg2+-EcoRI-DNA suggests a lower cleavage reaction in the mutated complex.