Article Information

• PDF (151 kb)
• Full Text with Large Figures

PubMed

• Articles by Srikanta Sen
• Articles by Lennart Nilsson

Related Articles

• Molecular computing: Does DNA compute?

  Molecular computing: Does DNA compute? Current Biology, Volume 6, Issue 3, 1 March 1996, Pages 254-257 Daniel E. Rozen, Steve McGrew and Andrew D. Ellington Summary The demonstration that DNA molecules can act as parallel processors to solve hard problems has excited interest in the possibility of developing molecular computers based on recombinant DNA techniques. Summary | Full Text | PDF (83 kb)

• Coupling of Fast and Slow Modes in the Reaction Pathway of t...

  Coupling of Fast and Slow Modes in the Reaction Pathway of the Minimal Hammerhead Ribozyme Cleavage Biophysical Journal, Volume 93, Issue 7, 1 October 2007, Pages 2391-2399 Ravi Radhakrishnan Abstract By employing classical molecular dynamics, correlation analysis of coupling between slow and fast dynamical modes, and free energy (umbrella) sampling using classical as well as mixed quantum mechanics molecular mechanics force fields, we uncover a possible pathway for phosphoryl transfer in the self-cleaving reaction of the minimal hammerhead ribozyme. The significance of this pathway is that it initiates from the minimal hammerhead crystal structure and describes the reaction landscape as a conformational rearrangement followed by a covalent transformation. The delineated mechanism is catalyzed by two metal (Mg) ions, proceeds via an in-line-attack by CYT 17 O2′ on the scissile phosphorous (ADE 1.1 P), and is therefore consistent with the experimentally observed inversion configuration. According to the delineated mechanism, the coupling between slow modes involving the hammerhead backbone with fast modes in the cleavage site appears to be crucial for setting up the in-line nucleophilic attack. Abstract | Full Text | PDF (504 kb)

• Studies of Proton Translocations in Biological Systems: Simu...

  Studies of Proton Translocations in Biological Systems: Simulating Proton Transport in Carbonic Anhydrase by EVB-Based Models Biophysical Journal, Volume 87, Issue 4, 1 October 2004, Pages 2221-2239 Sonja Braun-Sand, Marek Strajbl and Arieh Warshel Abstract Proton transport (PTR) processes play a major role in bioenergetics and thus it is important to gain a molecular understanding of these processes. At present the detailed description of PTR in proteins is somewhat unclear and it is important to examine different models by using well-defined experimental systems. One of the best benchmarks is provided by carbonic anhydrase III (CA III), because this is one of the few systems where we have a clear molecular knowledge of the rate constant of the PTR process and its variation upon mutations. Furthermore, this system transfers a proton between several water molecules, thus making it highly relevant to a careful examination of the “proton wire” concept. Obtaining a correlation between the structure of this protein and the rate of the PTR process should help to discriminate between alternative models and to give useful clues about PTR processes in other systems. Obviously, obtaining such a correlation requires a correct representation of the “chemistry” of PTR between different donors and acceptors, as well as the ability to evaluate the free energy barriers of charge transfer in proteins, and to simulate long-time kinetic processes. The microscopic empirical valence bond (Warshel, A., and R. M. Weiss. 1980. . 102:6218–6226; and Åqvist, J., and A. Warshel. 1993. . 93:2523–2544) provides a powerful way for representing the chemistry and evaluating the free energy barriers, but it cannot be used with the currently available computer times in direct simulation of PTR with significant activation barriers. Alternatively, one can reduce the empirical valence bond (EVB) to the modified Marcus’ relationship and use semimacroscopic electrostatic calculations plus a master equation to determine the PTR kinetics (Sham, Y., I. Muegge, and A. Warshel. 1999. . 36:484–500). However, such an approximation does not provide a rigorous multisite kinetic treatment. Here we combine the useful ingredients of both approaches and develop a simplified EVB effective potential that treats explicitly the chain of donors and acceptors while considering implicitly the rest of the protein/solvent system. This approach can be used in Langevin dynamics simulations of long-time PTR processes. The validity of our new simplified approach is demonstrated first by comparing its Langevin dynamics results for a PTR along a chain of water molecules in water to the corresponding molecular dynamics simulations of the fully microscopic EVB model. This study examines dynamics of both models in cases of low activation barriers and the dependence of the rate on the energetics for cases with moderate barriers. The study of the dependence on the activation barrier is next extended to the range of higher barriers, demonstrating a clear correlation between the barrier height and the rate constant. The simplified EVB model is then examined in studies of the PTR in carbonic anhydrase III, where it reproduces the relevant experimental results without the use of any parameter that is specifically adjusted to fit the energetics or dynamics of the reaction in the protein. We also validate the conclusions obtained previously from the EVB-based modified Marcus’ relationship. It is concluded that this approach and the EVB-based model provide a reliable, effective, and general tool for studies of PTR in proteins. Finally in view of the behavior of the simulated result, in both water and the CA III, we conclude that the rate of PTR in proteins is determined by the electrostatic energy of the transferred proton as long as this energy is higher than a few kcal/mol. Abstract | Full Text | PDF (917 kb)

• …more

Copyright © 1999 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 77, Issue 4, 1801-1810, 1 October 1999

doi:10.1016/S0006-3495(99)77025-6

Biophysical Theory and Modeling

Free Energy Calculations and Molecular Dynamics Simulations of Wild-Type and Variants of the DNA-EcoRI Complex

Srikanta Sen and Lennart Nilsson^,  

Center for Structural Biochemistry, Karolinska Institute, Department of Biosciences, Huddinge, Sweden

Address reprint requests to Lennart Nilsson, Center for Structural Biochemistry, Karolinska Institute, Department of Biosciences, Halsovagen 7, Floor 7, S-141 57 Huddinge, Sweden. Tel.: 46-8-608-9228; Fax: 46-8-608-9290.

EcoRI is experimentally the most studied restriction endonuclease and is now used as one of the prototypes in understanding the physical basis of biomolecular recognition (Draper, 1993,Kim et al,Newman et al,Lesser et al,Robinson and Sligar, 1993,Robinson and Sligar, 1994,Venclovas et al,Misra et al). EcoRI binds to DNA specifically to the base sequence d(GAATTC)₂ and, in the presence of Mg⁺² ion as a cofactor, it cuts both the DNA strands by hydrolyzing the sugar-phosphate backbone at the position between G and A in the recognition base sequence. In the companion paper we report the results of a detailed study of the structural, interactional and dynamical aspects of the wild-type DNA-EcoRI complex in aqueous solution by molecular dynamics (MD) simulation (Sen and Nilsson, 1999). However, like the other restriction endonucleases, EcoRI shows extreme selectivity in its interaction with DNA. Alteration in a single basepair in the recognition site can affect its binding affinity and functional activity substantially (Lesser et al,Lesser et al). To characterize the different specific interactions considerable experimental work has been done, and plenty of experimental data on the free energy differences measured in biochemical experiments is available for different mutations made in the recognition site base sequence of the DNA in the DNA-EcoRI complex (Lesser et al,Lesser et al). Experimental data show that a mutation resulting from a chemical modification of a functional group of a base in the recognition site is associated with a free energy difference of 1 to 2kcal/mol, whereas an entire basepair substitution causes a larger free energy difference of about 10kcal/mol (Lesser et al,Lesser et al). On the other hand, theoretical estimation of the free energy differences in biomolecular interactions by molecular simulations (MD or Monte Carlo) is quite common nowadays (Bash et al,Cieplak et al,Härd and Nilsson, 1993,Elofsson et al,Miyamoto and Kollman, 1993,Eriksson and Nilsson, 1995,Essex et al,Soares et al). In the present work, we have performed such calculations by MD simulations for several cases of the mutant variants of the DNA-EcoRI complex. Because calculating the free energy differences due to mutations in biomolecular systems is very time consuming, we have selected only a few specific cases for our study where mutation has been made in a DNA base in the recognition site by altering a functional group that is known to be involved in direct interaction with the protein in the wild-type complex. The objective of the present work is twofold. One aspect is to calculate the free energy differences in the cases of the selected mutants by chemical perturbation and molecular dynamic simulation methods and to compare these values with the corresponding experimental data in order to see how successfully the modeling studies can describe the intermolecular interactions and the stability difference between the wild-type DNA-EcoRI complex and its different mutant variants. The other aspect is to characterize the structural and interactional properties of each of the mutant variants of the complex considered, comparing those to the corresponding wild-type complexes to identify the differences introduced in the properties of the complex due to the individual mutations in each case. For this purpose we have performed additional ordinary MD simulations of each of these mutants. It is particularly interesting to point out that there are experimental data for cases where the same chemical modification made on the same base at two different positions has resulted in opposite effects on the stability of the complex (Lesser et al). Mutation by replacing the NH₂ group at the atomic position 6 of the first adenine base in the recognition site by a hydrogen atom is experimentally found to be unfavorable by a free energy difference of 1.3±0.2kcal/mol, which is consistent with the fact that a H-bond with the protein is deleted by this mutation (Lesser et al,Draper, 1993). On the other hand, the same modification at the next adenine base of the DNA is experimentally found to be energetically favorable, as it is accompanied by a net negative free energy change of −1.0±0.1kcal/mol (Lesser et al), even though one H-bond with the protein is deleted in this case, too. In order to account for this observed preference for the change, it has been suggested that the resulting negative free energy difference in this case may be the consequence of the overall structural relaxation of the kink deformation of the DNA associated with this mutant complex (Lesser et al). In these experiments with different modified bases, it was further assumed that these individual modifications remove only the contribution of the particular functional group of the base which is modified and keep the other interactions between the base and the protein intact. So one of our major objectives in this study was to perform the free energy calculation and ordinary MD simulations of the fully solvated system in this second case and to compare the results with those for the wild-type complex, to obtain better insight into the details of what happens in these two cases and, if possible, to verify the above speculation and assumptions made in this case. It may also be noted that in the free energy calculations for each mutant, we have performed several independent dynamic simulations with different initial velocity conditions to avoid any initial condition-dependent bias in the estimated free energy differences. The estimated values of the free energy differences in the different cases of the mutants (except one case) show good qualitative agreement with the corresponding experimental results. However, the free energy estimates for such complex biomolecular systems by molecular simulations are generally associated with large statistical fluctuations resulting mainly from inadequate sampling of the conformational space. For the analysis and comparison of the structural and interactional aspects of the different mutants, we have looked into different average quantities from the dynamic trajectories and have compared them for all the mutant cases and the wild-type complex. We have also compared the interaction strengths and H-bond lists for all the mutants in order to have further details about the difference in the interaction patterns of the different mutants. The analysis of the local structure and interaction pattern in the mutation locality from ordinary MD simulations indicates that even a small specific alteration in the DNA-protein interaction can affect both the specific (DNA base and protein) and nonspecific (DNA backbone and protein) interactions of the other nearby nucleotides and the altered nonspecific interaction thus can have some role in determining the overall specificity of binding. On the whole, we feel that the results presented here, along with the corresponding experimental data, provide us with much more detailed insight into the DNA-protein interaction at the atomic level in the cases of different mutants of the DNA-EcoRI complex.

As the free energy calculation by chemical perturbation method is now a widely used technique, we will present only an outline of the basic theory (Fleischman and Brooks, 1987,King, 1993; van Gunsteren et al,Straatsma et al). In thermodynamic perturbation method, one starts with a hybrid molecular system consisting of both the reactant and the product parts along with the environmental part. The Hamiltonian is constructed as follows, depending on the molecular coordinates (r) and a coupling parameter λ, as

(1)

(2)

(3)

We have performed calculations for several cases of mutations on the DNA bases of the recognition part in the complex for which we have experimental data. In each of these mutants a functional group or atom of a relevant base has been replaced by another group or atom such that the formation of a crucial hydrogen bond between the DNA and the protein in the complex is eliminated. The different modifications used are described below.

A→⁷A: In this modified base the N7atom of an adenine base is replaced by a CH group (Figure 1a)

Display large version of this figure

Figure 1

Schematic diagrams of the different modifications used in the present study. (a) The N7atom of Adenine base is replaced by a C-H group. (b) The N7atom of Guanine base is replaced by a C-H group. (c) The NH2 group of Adenine base is replaced by a hydrogen atom.

G→⁷G: This represents the replacement of the N7atom of a guanine base by a CH group (Figure 1b).
A→P: This is a modified adenine base where the NH₂ group at the position 6 is replaced by a hydrogen atom (Figure 1c).

Because in these mutations the bases are partially modified, one needs to know the partial atomic charges of the modified base. Calculation of accurate partial atomic charges for a molecular system is a very difficult task and there is no unique, straightforward way of finding them on a rigorous level. In the present work we have adapted the following procedure. In each case, the partial atomic charges of the altered base in the product state were calculated by MOPAC 6.0 package using the AM1 parameter set and the electrostatic potential (ESP) fitting method (Steward, 1990). The atoms nearest to the product atoms in each case were considered the colocated atoms, whose charges were different in the reactant and product states. As a working approximation, we have used the calculated partial atomic charges for the product atoms and distributed the difference in the total charge of the product and colocated atoms between their calculated charges and the corresponding charges in CHARMM topology equally on the calculated charges of the colocated atoms to keep the net charge the same. The charges thus obtained for the product and colocated atoms in the different mutant bases are summarized in Table 1. The partial atomic charges of the other atoms in the base were kept as usual CHARMM charges to maintain an overall consistency with the other parameter set of CHARMM.

We have considered the following cases of mutants with a chemical modification and assigned a name (given in parentheses) to each case: (i) A → P is made at the base position 5 of the first strand (M1); (ii) A → P is made at the base position 6 of the first strand (M2); (iii) G → ⁷G P is made at the base position 4 of the first strand (M3); (iv) A → ⁷A is made at the base position 5 of the first strand (M4); and (v) A → ⁷A is made at the base position 6 of the first strand (M5).

In order to calculate the free energy difference ΔΔA for each mutant we prepared two independent system setups, one for the DNA-protein complex in solution and the other for the free DNA in solution. In the case of the complex, starting from the energy-minimized crystallographic coordinates of the EcoRI-DNA complex, we cut out a part that was within a sphere of 16Å radius with the nitrogen atom (N6 or N7, depending on the case) of the reactant part of the respective base in the mutation site at the center. The respective base in the mutation site was replaced by a hybrid system consisting of the normal base for the reactant state and the modified base for the product state in each mutant case. This system containing the hybrid part was then immersed in a pre-equilibrated TIP3P water sphere of 19Å radius with the NH₂ group of the hybrid system at the center, and the water molecules whose oxygen atoms were within a distance of 2.8Å from any non-hydrogen atom of the complex were removed (Jorgensen et al). To make the system electrically neutral we placed 13 Na⁺ counterions by replacing 13 water molecules whose oxygen atoms had highest electrostatic energies and are >5Å apart from each other. The system was then energy-minimized by 500 steepest descent steps, keeping the reactant and the product atoms of the hybrid part fixed. The energy-minimized system was then used for the MD simulation for free energy perturbation calculations. In all the energy minimization and subsequent dynamic simulations, the free ends of the backbone of the disjointed fragments of the protein created as a result of cutting out a sphere (16Å radius) from the complex as mentioned above were kept harmonically constrained with a force constant of 10kcal/mol/A² to preserve the effect of the natural covalent continuity of the protein backbone.

For the free DNA simulation, in a similar way, from the DNA in standard B-form corresponding to the DNA in the complex and containing the same hybrid part as in the complex, a sphere of 16.0Å radius as mentioned above was cut out. Seventeen Na⁺ counterions were included for electro-neutrality of the system and each counterion was placed at a distance 3.5Å from the phosphorus atom on the line bisecting the line joining the two oxygen atoms of the respective phosphate group of this DNA fragment. The whole system was then solvated in a TIP3P water sphere of 19Å radius and energy-minimized by 500 steepest descent steps, keeping the reactant and the product atoms of the hybrid part fixed, and was subsequently used in the MD simulation.

To investigate the local structural and interactional characteristics of each of the different mutant variants at the modification site, we prepared the solvated complex in a way similar to that described above, except that in these cases we replaced the hybrid base with the modified base. Each such system was then subjected to SDB dynamics simulations.

As the mutant M2 experimentally yields free energy value qualitatively quite different from those for other mutants, we treated the case of this mutant in more detail. For a direct comparison with the properties of the wild-type complex, we prepared the setup with the full complex M2 containing the modified part and fully solvated by a solvent shell 7Å thick as described in the companion paper (Sen and Nilsson, 1999). Ordinary MD simulation was then performed on this system to investigate the details of the interactional and structural aspects of the particular case.

In all the cases of mutants we have performed MD simulation of the solvated DNA-EcoRI complex system by SDB dynamics algorithm for sampling at different intermediate states using the standard CHARMM potential energy function and parameter set (Brooks et al,MacKerell et al,MacKerell et al) and the free energy perturbation implementation (Fleischman and Brooks, 1987). The atoms in the spherical shell (the buffer region) between radii 17Å and 19Å executed dynamics according to Langevin dynamics algorithm, whereas the atoms inside the sphere of 17Å radius were subjected to ordinary MD, following leap-frog algorithm. The solvent molecules were subjected to a deformable boundary force arising from the mean field interaction of water molecules beyond the 19Å boundary (Brooks and Karplus, 1983).

In the cases of free energy calculation for the mutants, the simulations were performed with a time step of 2 fs. We used 13.0Å as the cutoff value for the nonbonded interactions and the nonbonded interaction list was updated every 10 steps. We used force shift option causing the interaction energies and the forces to vanish smoothly at a distance of 12.0Å. We also used the SHAKE algorithm (Ryckaert et al) for constraining all bonds involving hydrogen atoms. We used a dielectric constant of value 1.0. The system was connected to a heat bath at a temperature of 300K. All the non-hydrogen atoms were assigned a friction coefficient of 50ps⁻¹. The system was first equilibrated at 300K with λ=0.5 for 50ps and then simulated in seven consecutive windows at seven λ values of 0.05, 0.125, 0.25, 0.50, 0.75, 0.875, and 0.95, respectively. In each window a 20-ps equilibration followed by a 40-ps production run was performed. We have kept the bond term and the bond angle term in the potential energy function unperturbed to maintain the structure of the perturbed part of the system for λ values close to the limiting values 0 or 1. We also used a double-wide sampling method over the full range of λ to calculate the overall free energy difference in the transformation. Simulation was performed in an identical way for the free DNA.

For each mutant, the same protocol for the dynamics equilibration and the production runs was followed. We performed at least three independent simulations on each of these systems, starting with a different initial velocity assignment, to get an idea of the initial velocity-dependent fluctuation in the computed values.

In the cases of setups with the modified bases in the product state to investigate the local structural and dynamical properties, SDB dynamics was performed following the same protocol as described above, without the perturbation part.

In the case of the mutant M2, for direct comparison with the properties of the wild-type complex we have performed the ordinary molecular dynamics simulations of the fully solvated mutant complex M2 by a solvent shell 7Å thick in the manner described in the companion paper (Sen and Nilsson, 1999).

Table 2 summarizes the results of free energy calculations. It is seen that in all the mutations considered, the calculated free energy difference ΔΔA from MD simulations indicates that the chemical modifications make the resulting DNA-protein complex less favorable compared to the wild-type complex. This is consistent with the fact that in each case, a H-bond between the DNA and the protein is deleted. The computational data in most of the cases (except the case of M2) are thus in good qualitative agreement with the corresponding experimental data, although the quantitative agreements are not that good (Lesser et al). The calculated data also show considerable fluctuations in ΔΔA values obtained from different independent free energy perturbation simulations for the same mutant and it is noticeable that the root mean square fluctuation in the estimated free energy differences (ΔA′) in the case of the complex are, in general, larger than the corresponding ΔA in the free DNA case. However, such statistical fluctuations are generally associated with the free energy calculations of complex macromolecular systems from molecular simulations. Several things seem to be responsible for these observed differences between the calculated average values of the free energy differences, the corresponding experimental data, and the significant fluctuations of the values obtained in different independent simulations. First, the empirical force field and the parameters used to describe such complex molecular systems may not be sufficiently accurate. Second, errors may be introduced due to finite sampling, which may not be adequate in all the cases studied. Finally, in such calculations, because ΔΔA (which is usually a small quantity) is estimated as the difference of two relatively large quantities, ΔA′ and ΔA”, significant errors are generally associated with it. Calculations show that, in general, a large change in free energy (ΔA′ or ΔA”) is associated with such chemical transformations. This large change in free energy actually arises due to the difference in the interaction energies of the perturbed base with the rest of the system in the native and modified cases. Thus, as it basically represents the difference in the average enthalpies of the native and the modified DNA, its actual value is not important for our present purpose and when a difference is taken to obtain ΔΔA, this factor cancels out, as it is present in both the free DNA and the DNA in the complex.

The plot of the ΔA-vs.-λ curve for each mutant as represented in Figure 2 and Figure 3 demonstrate how the cumulative value of ΔA varies with λ in the free DNA (dashed line) and in the complex (solid line) for the different mutants. It may be noted that even in cases of the same chemical modification but at different base positions, there are some differences in the nature of the plots. The reason seems to be the difference in the microscopic nature of the neighborhood. However, these are not physically very relevant, as the states of the hybrid system in the range 0<λ<1 represents unphysical chemical states of the system and, the free energy being a thermodynamic state function, only the end point values are important.

Display large version of this figure

Figure 2

The variation of the accumulated free energy difference (DA) in the different cases of mutations. Solid lines represent the case of the complex and the dashed line corresponds to the free DNA simulation in the cases of mutants M1 (a) and M2 (b).

Display large version of this figure

Figure 3

The variation of the accumulated free energy difference (DA) in the different cases of mutations. Solid lines represent the case of the complex and the dashed line corresponds to the free DNA simulation in the cases of mutants M3 (a), M4 (b), and M5 (c).

Because the DNA in the complex contains a distorted kinked part at the center of the recognition sequence and the modifications are made in the bases in the locality of the kink, the DNA is the most likely molecular component where significant structural relaxation may occur. We have compared the average RMSD of the central five bases (GAATT) of the DNA recognition site (to avoid the end effects) containing the kink for all the different mutants. The results are presented in the Table 3. Comparison shows very similar RMSD values for this part in all the different mutant cases, indicating that there is no exceptional structural relaxation of the kink of the DNA in any of these cases, including the mutant M2, where a significant relaxation was postulated. However, in reference to the wild-type complex the average values of the RMSD in the individual cases of the mutants are found to be larger, as expected.

Fig. 4 compares the average RMSD of the individual nucleotides of the DNA. It is found that generally, for the different mutants, larger structural rearrangements occurred than in the wild-type complex, and the effects of these small modifications extend over a substantial part of the DNA.

Display large version of this figure

Figure 4

Comparison of the average RMSD of the individual bases of the recognition site of the DNA strand DNA1 against the nucleotide number. The symbols ♢, +, □, ×, △, and &cjs3752; represent the cases of the mutant variants M1, M2, M3, M4, M5, and the wild-type complexes, respectively.

Comparison of the specific (DNA base and protein) and nonspecific (DNA backbone and protein) interaction energies in the average structures of the different mutant variants and the wild-type of the complex in aqueous solvent indicates that even these small changes in a specific functional group or atom of a single base can induce significant changes both in the specific (Fig. 5) and nonspecific (Fig. 6) interaction pattern of the individual bases, and the influence is extended over at least a few basepairs around the modification sites. The quantity ΔE=E_mutant−E_wt can be used as a measure of the specificity of a particular modification, where E_mutant and E_wt represent the total energy of interaction between the DNA and the protein in the mutant and the wild-type complex, respectively. If ΔE is large it means that the modification strongly affects the specificity, whereas an ΔE value ∼0 indicates a nonspecific nature of the modification. Again, ΔE can be expressed as ΔE=ΔE (specific)+ΔE (nonspecific), where ΔE (specific) comes from the difference in the direct specific interaction and ΔE (nonspecific) represents the difference resulting from the altered nonspecific interaction part. This clearly indicates that even a difference in nonspecific interaction due to a modification can significantly contribute to the overall specificity of the modification site.

Display large version of this figure

Figure 5

Comparison of the average specific interaction energies of the individual bases of the DNA strand DNA1 with the entire protein, against the nucleotide number. The cases of the mutant variants M1, M2, M3, M4, M5, and the wild-type complexes are represented by the symbols ♢, +, □, ×, △, and respectively.

Display large version of this figure

Figure 6

Comparison of the average nonspecific interaction energies of the individual bases of the DNA strand DNA1 with the entire protein against the nucleotide number. The symbols ♢, +, □, ×, △, and represent the cases of the mutant variants M1, M2, M3, M4, M5, and the wild-type complexes, respectively.

It was also found that both specific and nonspecific interactions of the Gua⁴ and the protein are highly sensitive to the modifications made at other bases in the recognition site. It is, then, quite interesting to identify Gua⁴ as the first base in the recognition site on the DNA strand and one of the two bases closest to the scissile bond where the DNA strand is cut by the protein. It may also be noticed that in some of the mutants, the interaction energies between some DNA nucleotides and the protein are strikingly similar, though for others they are quite different (Figure 5 and Figure 6). However, in all the cases except that of the Ade⁵ → P modification, the interaction strengths at the modified base are changed maximally. It was also found that these modifications not only affected the interactions of the bases of the modified DNA strand with the protein, but also influenced the interactions of the other DNA strand and the protein as well, though to a lesser extent. Thus, even the effects of small modifications are found to be quite complicated as the whole complex behaves as a single integrated system.

Table 4 summarizes the comparison of the specific, nonspecific, and overall intermolecular interaction energies between the GAATT part of the recognition sequence on the first strand (DNA1) of the DNA duplex and the protein. It can be seen clearly that each of these quantities is quite different for the different mutants and the wild-type complex. However, the overall relative stability of a complex is governed by the total free energy difference between the two states of the DNA and the protein where they are free in solution and in bound form as a complex in solution, and thus these values cannot be directly correlated with the relative stabilities of these mutant variants. These results are presented mainly to demonstrate that such small modifications may also cause significant differences in the overall intermolecular interactions between the DNA and the protein.

Comparison of the H-bond list (Table 5A,Table 5B) between the DNA (in the modification site) and the protein for all the mutant variants and the wild-type complex clearly indicates that most of the major H-bonds between the modified base and the protein seen in the wild-type complex are also well maintained in the solution-simulated structures of the different mutant variants. This is also true for the other neighboring bases which preserve their H-bonding with the protein. However, the interaction energies are changed significantly in all these different mutants (Figure 5 and Figure 6). The alteration in the interaction energy strengths in each case is caused by altered interaction geometry of the interacting atoms as a result of small local rearrangements due to the chemical modifications. It may be noted in Table 2 that in the case of the mutant M4, the theoretically estimated values of ΔΔA are small in all the independent calculations and, in the H-bond lists, it is found that the H-bond involving N7atom of Ade⁶ is not strong. Thus, these results are consistent with each other.

The case of this mutant variant is of special interest, as it was mentioned earlier that the experimental value of ΔΔA for this mutant (M2) is negative (−1.0kcal/mol), indicating that the mutated complex is more stable than the wild-type complex. It was then speculated (Lesser et al) that this mutant variant must be associated with a structural relaxation of the kink in the DNA, resulting in an overall favorable change in free energy. However, the theoretically computed value of ΔΔA from MD simulation in this case gives an average value of 1.1±1.6kcal/mol (Table 2), indicating that the mutation in this case is also less stable than the wild-type complex. This result contrasts qualitatively with the experimental data (Table 2), even though the difference (2.1kcal/mol) from the experimental value is similar to the other cases, and it is consistent with the fact that this mutant also lacks a functional group on the DNA base that interacts with the protein through a H-bond in the wild-type complex. The calculated positive values of ΔΔA in most of the independent simulations with this mutant clearly indicate the generality of this result. For better insight into what actually happened in this case, we analyzed the 700-ps trajectory of the whole complex with this mutant. Comparison of the RMSD of the individual nucleotides of the DNA (Fig. 4) does not reveal any indication of any significant difference. The average torsion angles of the nucleotide Ade⁶ (Table 6A) and of the side chain Asn¹⁴¹ (Table 6B), which interacts with the Ade⁶ base in the wild-type complex, also indicates very similar local structures in the wild-type complex and in the mutant variant.

It is possible that the wild-type conformation from which our simulations were started is kinetically stable enough for this mutant that a larger relaxation cannot be achieved in a finite simulation. Thus, we find that the result of the free energy calculation from our simulations in this case does not agree with the corresponding experimental data and in the present dynamic simulation we did not find any evidence of significant relaxation of the DNA kink as suggested by the experimental group. Our theoretical estimate of ΔΔA for this mutant obtained from several independent simulations is, however, completely consistent with the fact that we have not observed any structural relaxation in the DNA.

MD simulations and free energy calculations of the wild-type EcoRI-DNA complex and several mutant variants have been performed in aqueous solvent to gain insight about the interaction and recognition mechanism of the complex. The main results we obtained are summarized here.