INTRODUCTION

Top Abstract Introduction Materials & Methods Results & Discussion Conclusions References

Treatment of long-range electrostatic interactions

Calculation of long-range electrostatic interactions is the most time-consuming part of molecular dynamics (MD) simulations. The usual way to improve the efficiency of the calculation of electrostatic interactions is to truncate Coulombic interactions and neglect those beyond a specified cutoff distance. The truncation discontinuity can be smoothed by applying switching or shifting functions to the energies or forces (Steinbach and Brooks, 1994). The use of truncation methods nevertheless leads to over or underestimation of Coulombic forces acting on the atoms and to significantly reduced accuracy and artificial behavior (Loncharich and Brooks, 1989; Schreiber and Steinhauser, 1992; Guenot and Kollman, 1993). In particular, highly charged molecules like nucleic acids tend to become unstable during MD simulation when a truncation method is applied (Singh et al., 1985; Fritsch et al., 1993; Miaskiewicz et al., 1993).

The particle mesh Ewald (PME) method is one of the approaches to consider all electrostatic interactions over the complete system in a computationally tractable manner (Darden et al., 1993). The algorithm employs interpolation of reciprocal space Ewald sums from a grid and computation of convolutions using fast Fourier transforms. For DNA and RNA molecules of 10-20 basepairs in solution, use of the PME method has enabled stable structures to be obtained during simulations of up to 5 ns (Auffinger and Westhof, 1997; Cheatham et al., 1995; York et al., 1995; Young et al., 1997; Cheatham and Kollman, 1997). These simulations have enabled conformational transitions, the effects of the environment, and the properties of the ordered solvent and counterions to be studied. The robustness of the PME method for DNA simulations was recently shown when the effects of using different equilibration protocols (De Souza and Ornstein, 1997a) and different sizes of periodic box (with minimum solute-box edge distances of 5, 10, and 15 Å) (De Souza and Ornstein, 1997b) for simulation of the same DNA molecule were investigated. The parameters (root mean square deviation (RMSD) values, global helix axes curvature) compared showed no significant differences in these MD simulations.

For MD simulations of proteins, application of the PME method (York et al., 1994; Fox and Kollman, 1996) has produced lower RMSD values than those obtained with truncation of electrostatic interactions. For example, for BPTI, the backbone RMSD values for equivalent simulations with PME and truncation were 0.5 Å and 1.8 Å, respectively (York et al., 1994). Ubiquitin maintained an overall fold very close to the crystallographic structure in a simulation with the PME method: the RMSD for backbone atoms was 1.1 Å after 1.2 ns of simulation. The radius of gyration, the solvent accessible area, and the hydrogen bond network all remained nearly constant (Fox and Kollman, 1996).

These data show the robustness and value of the PME method for simulation of nucleic acids. Its value has also been clearly demonstrated for proteins, but sensitivity to system setup and parametrization has not been studied to so great an extent as for DNA. Here, we investigate sensitivity to the placement and modeling of salt ions in protein simulations.

Counterions

A neutral system is required for use of the PME method (Darden et al., 1993). One way to achieve this is to modify the partial atomic charges of the solute. For example, reducing electrostatic charges to 0.24 e.u. on the phosphate groups to treat the effect of counterions implicitly kept the MD trajectory of DNA stable during a 1-ns simulation (with truncation of long-range interactions), although the conformations derived were ~4.5 Å from that of canonical B-DNA (McConnell et al., 1994). For proteins, many simulations have been carried out with partial atomic charges modified so that the net charge of charged side chains is zero. This is usually done when water molecules are not modeled explicitly so as to implicitly account for the damping of electrostatic interactions by water.

A more realistic way to obtain a neutral system is to add counterions by replacing some of the water molecules solvating the solute that are near to charged groups on the solute. This is the approach adopted for simulations with the PME method. For simulation of nucleic acids, the counter-ions are initially placed close to the phosphate groups (Young et al., 1997). For proteins, charge-balancing counterions are positioned close to charged groups in the regions with the most favorable electrostatic potential for binding on the basis of either geometric or energetic criteria and, sometimes, using quite sophisticated procedures, such as genetic algorithms (Li et al., 1997).

As a neutral system is not necessary for simulations of proteins with truncation of long-range interactions, no counterions are modeled in most such simulations. Neglect of counterions is possible because of the smaller charge density of proteins compared to nucleic acids (for which neutralization of charge or damping of electrostatic interactions is essential for obtaining reasonable stability). The neglect of counterions simplifies the system and prevents potential problems with strong long-range electrostatic interactions and long equilibration times for the counterions. Nevertheless, improved protein stability has been observed when adding counterions (Yelle et al., 1995; York et al., 1993). They screen the charged side chains, preventing unnatural electrostatic interactions between neighboring parts of the macromolecule.

In addition to assigning the positions of counterions for PME simulation, the number of counterions to add must be chosen. The screening of charged side-chain interactions depends not only on the distribution of counterions, but also on the concentration of counterions used. Proteins and peptides have been simulated with explicit counterions at high (1 M) ionic strength (Marlow et al., 1993; Avbelj et al., 1990). For most simulations, however, counterions are only placed near the protein surface and the full ionic environment is not modeled. In the present paper we show that, for maintaining the stability of the WW domain of YAP (Yes kinase-associated protein), explicit modeling of full ionic strength (at 0.2 M) was necessary. The usual procedure of only placing charge-balancing counterions near charged groups on the protein surface did not result in a stable protein fold, nor did simulation without any counterions.

The WW domain is a small protein module of ~40 residues with two highly conserved tryptophans that binds proline-rich peptides (Chen and Sudol, 1995). The binding of ligands by the WW domain is thought to play a role in several human genetic disorders, such as Liddle's syndrome, muscular dystrophy, and Alzheimer's disease (Sudol, 1996), and in the budding process of HIV (human immunodeficiency virus) and RSV (Rous sarcoma virus) (Garnier et al., 1996). The first structure of a WW domain, that from YAP, was solved by NMR (Macias et al., 1996) and consists of a twisted three-stranded antiparallel -sheet. No secondary structure was found for the N- and C-termini (Fig. 1). A difficulty during the structure determination of this particular WW domain was the poor stability of the protein. Sufficient stability for structure determination could only be obtained in a construct of 50 amino acid residues with a disordered N-terminal region that the structure determination revealed formed a "buckle" covering a hydrophobic part of the surface of the -sheet. The small size, -sheet structure, and weak stability of the YAP WW domain make it very sensitive to conditions in molecular dynamics simulations, and thus a good system on which to test different protocols and models. It was necessary to investigate a number of different conditions to obtain a stable WW domain during the MD simulations. Under standard simulation conditions, the WW domain was unstable and its secondary structure was destroyed within 10 ps. We found that a more complete and realistic representation of salt conditions than is usually necessary was required to obtain a stable equilibrated system for this protein of marginal stability.

View larger version (22K): [in this window] [in a new window]   FIGURE 1   Secondary structure of the WW domain. This figure was prepared using the MOLSCRIPT program (Kraulis, 1991).

	MATERIALS AND METHODS

Top Abstract Introduction Materials & Methods Results & Discussion Conclusions References

The initial positional coordinates for MD simulation were taken from the NMR structure of the complex of the WW domain with a proline-rich peptide (Macias et al., 1996). As the experiment indicated very little conformational change of the WW domain upon binding the peptide, the peptide was removed. To model the charge distribution on the surface of the WW domain at pH 6.0, the pH at which the NMR structure was obtained, Arg and Lys residues were protonated and Asp and Glu residues were deprotonated. The single histidine present, His-32, was singly protonated at the N¹ position as this was the most favorable site according to the hydrogen placement algorithm (Hooft et al., 1996), based on optimization of the hydrogen bond network, in the WHATIF package (Vriend and Sander, 1993).

Both initial setup and dynamics runs were performed with the AMBER 4.1 program (Pearlman et al., 1995) using the all-atom force field (Cornell et al., 1995) and the PME method (Darden et al., 1993); van der Waals interactions were computed within a cutoff distance of 9 Å.

Two systems were generated for simulation. The first consisted of the WW domain with charge-balancing counterions at the protein surface surrounded by water molecules. The second consisted of the WW domain with charge-balancing counterions surrounded by ionic solution with water molecules and explicit Na⁺ and Cl ions; i.e., with full ionic strength modeled. For both systems, the 6 Na⁺ and 4 Cl charge-balancing counterions were added to neutralize charges on the protein surface using the standard AMBER 4.1 program CION. This explores the electrostatic energy surface near charged residues to find ion positions with favorable counterion-residue interactions. To set up the second system, so as to perform MD simulation of the WW domain under ionic strength conditions similar to those employed in the structure determination (0.1 M NaCl), the WW domain with charge-balancing counterions was immersed in the middle of a 56 × 56 × 56 Å³ box of equilibrated 0.23 M NaCl solution with water molecules removed (V. Lounnas, unpublished data). Explicit Na⁺ and Cl ions lying within 4 Å and 5 Å of protein atoms and charge-balancing counterions, respectively, were removed. After this, the system contained 54 counterions (28 Na⁺ and 26 Cl).

For both systems, the WW domain and counterions were surrounded by ~5000 water molecules represented by the TIP3P model (Jorgensen et al., 1983). For placement of water molecules, the minimum distance of a water oxygen atom from any of the solute atoms was 2.6 Å, and the minimum distance of a water hydrogen atom was 2.0 Å. The systems were both ~56 × 56 × 56 Å³ with a minimum distance from a protein atom to the box edge of 12 Å and contained ~16,000 atoms.

Each system was coupled to a thermal bath of temperature T₀ = 300 K and to a pressure bath of pressure P₀ = 1 atm by applying the algorithm of Berendsen (Berendsen et al., 1984) with temperature relaxation time t_T = 0.4 ps and pressure relaxation time t_P = 0.4 ps. The pressure scaling was anisotropic. A time step of 2 fs was used for all the simulations and covalent bonds were constrained using the SHAKE algorithm (Ryckaert et al., 1977).

All the MD runs were set up using the same protocol. First, all hydrogen atoms, ions, and water molecules were subjected to 200 steps of minimization by steepest descent to remove close van der Waals contacts and to allow hydrogen bonds between water molecules in the periodic box and the protein to form. Then hydrogen atoms, ions, and water molecules were heated up to 300 K gradually during 15 ps. Then, a second minimization of 200 steps was performed for these atoms only. Thereafter, all atoms in the system were allowed to move. They were first energy-minimized by steepest descent with 200 steps and then heated up to 300 K gradually during 15 ps. Then simulation at a constant temperature of 300 K under NPT conditions was started. A simulation of 260 ps was run for the system with only charge-balancing counterions, and a simulation of 1 ns was run for the system with a full salt model.

A modified version (W. Bitomsky, unpublished data) of the CARNAL program (B. Ross, San Francisco, CA) was used to calculate RMSD values and the radius of gyration of the WW domain. The QUANTA program (MSI, San Diego, CA) was used to analyze secondary structure changes during the MD simulations. The QUANTA program uses a modified DSSP algorithm (Kabsch and Sander, 1983) based on hydrogen bond formation to define secondary structure. The ANAL program (B. Ross, San Francisco, CA) was applied to recalculate and analyze energy values for separate parts of the systems simulated. Because the PME method is not implemented in the ANAL program, a nonbonded cutoff of 100 Å was applied to both van der Waals and electrostatic energy evaluations. This large cutoff value exceeded the longest diagonal of the periodic box, allowing all possible nonbonded interactions in the systems simulated to be taken into consideration.

	RESULTS AND DISCUSSION

Top Abstract Introduction Materials & Methods Results & Discussion Conclusions References

Protein structural features

Because the overall fold of the WW domain is represented by a three-stranded -sheet, the secondary structure was examined to characterize conformational changes of the WW domain during MD simulations. The RMSD of backbone atoms of the -sheet (15 residues) from their initial positions were evaluated during the MD runs (Fig. 2). The trajectory of the simulation of the WW domain with full ionic strength modeled appeared to be well equilibrated with an average RMSD value of 1.1 ± 0.09 Å over the last 800 ps of the trajectory. In contrast, when simulated with charge-balancing counterions only, the secondary structure of the WW domain experienced larger conformational changes. After rapid initial structural perturbations, a plateau in the RMSD of 1.7 Å was reached. However, after 180 ps of simulation time, the RMSD started to increase again, and continued to do so until the simulation was stopped at 260 ps.

View larger version (36K): [in this window] [in a new window]   FIGURE 2   RMSD values of the backbone atoms of the -sheet of the WW domain. (gray), With charge-balancing counterions only. (black), With full ionic strength.

The RMSD curves correspond well with visual observation of the maintenance of the secondary structure in the WW domain (Fig. 3). With full ionic strength modeled, the -sheet kept its structure throughout the simulation with the exception of a momentary disappearance of the third -strand at the beginning of the simulation. With only charge-balancing counterions modeled, the complete -sheet was observed only rarely during the beginning of the simulation. The first and third strands disappeared in the first 10-20 ps, coincident with the initial increase in RMSD values. The most stable part, the second -strand, disappeared after ~50 ps of simulation time.

View larger version (12K): [in this window] [in a new window]   FIGURE 3   Secondary structure of the WW domain throughout MD simulation. (A) With charge-balancing counterions only. (B) With full ionic strength.

Energetic features

The energies of nonbonded interactions within the WW domain became more negative and then leveled out during both simulations (Fig. 4). Energetically, the structure of the WW domain relaxed to a similar extent in both simulations. There was only one difference between these two MD runs, namely the number of salt ions present. To explain the possible role of the salt concentration in the periodic box on the stability of the WW domain, the energy of the interactions between the WW domain and the ions was evaluated. As Fig. 5 shows, the magnitude of the WW domain-ion interaction energy decreased steeply when simulations included only charge-balancing counterions. The reason for this is that the counterions diffused away from the protein surface during the simulation. Whereas the average distance of charge-balancing counterions from the closest nitrogen or oxygen atom of the side chains of the charged residues on the protein was 3.1 Å initially, after 260 ps of simulation time all distances from the charged side chains of the protein (nitrogen or oxygen atoms) to counterions are at least twice as large. The average distance from charged residues to the nearest counterion is 11.2 Å.

View larger version (37K): [in this window] [in a new window]   FIGURE 4   Energy of nonbonded interactions within the WW domain. (gray), With charge-balancing counterions only. (black), With full ionic strength.

View larger version (25K): [in this window] [in a new window]   FIGURE 5   Interaction energy between the WW domain and counterions. (gray), With charge-balancing counterions only. (black), With full ionic strength.

With full ionic strength modeled, the decrease of WW domain-ion interactions occurred over a longer time. After 260 ps of simulation time, most of the charge-balancing counterions still hold positions near the protein surface (the average distance between charged residues (nitrogen or oxygen atoms) and the nearest counterion is 5.4 Å). As charge-balancing counterions diffuse away from the protein surface, other counterions approach the charged residues. As the simulation continues, the magnitude of the WW domain-ion interactions decreases and then levels off. Thus, while the fully equilibrated system does not have as strong protein-ion interactions as initially when the charge-balancing counterions were placed close to charged surface residues, there is a continual significant protein-ion interaction. The net charge on the WW domain is relatively small (+2e) and it has a rather hydrophobic face where the proline-rich peptide binds. This protein does not have highly charged electrostatic features that would have strong interactions with counterions. Thus, it is appropriate that, when fully equilibrated, the protein-ion interactions are rather weak and nonspecific for this system. Nevertheless, the explicit ions present in the periodic box provided an essential electrostatic potential that prevents diffusion of counterions away from the protein at the beginning of the MD simulation. The maintenance of counterions near the WW domain surface may be thought of as navigating the protein on the potential energy surface along those pathways that keep its secondary structure stable.

The full ionic strength model is important despite the moderate ionic strength of the experimental system (0.1 M) and lack of highly polar or charged features in the protein. Indeed, it is the lack of strong electrostatic features on the protein that necessitates the full ionic strength model: the protein-ion interactions are insufficiently strong to keep the charge-balancing ions close to the protein surface unless there are other ions present in the surrounding solution.

A concern with the introduction of counterions in MD simulations is that they may lead to protein distortions because of the long times for equilibration of counterion positions. However, for YAP-WW, a simulation without counterions (but with peptide bound) (G. Ibragimova and R. Wade, unpublished data) actually results in loss of secondary structure. An additional factor is that recent evidence (Kollman et al., 1997) suggests that the stability of -structure may be underestimated in the force field used. This might increase the sensitivity of the protein to its salt environment.

	CONCLUSIONS

Top Abstract Introduction Materials & Methods Results & Discussion Conclusions References

Small proteins of marginal stability such as the YAP WW domain can explore significantly larger surfaces on the potential energy landscape during MD simulations than more stable proteins of the same size. Their structure is more sensitive to simulation conditions. Our simulations of the YAP WW domain demonstrate sensitivity to the ionic strength model used and show the importance of employing a complete model of ionic strength with explicit ions throughout the periodic system studied for maintenance of protein secondary structure. This is likely to be important for other proteins, particularly those that do not have strong electrostatic features.