This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.104.058644v1 88/6/3829    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via Google Scholar Google Scholar Articles by Klähn, M. Articles by Gerwert, K. Search for Related Content PubMed PubMed Citation Articles by Klähn, M. Articles by Gerwert, K.

Theoretical IR Spectroscopy Based on QM/MM Calculations Provides Changes in Charge Distribution, Bond Lengths, and Bond Angles of the GTP Ligand Induced by the Ras-Protein

Ruhr-Universität Bochum, Lehrstuhl für Biophysik ND 04, 44780 Bochum, Germany

Correspondence: Address reprint requests to Klaus Gerwert or Jürgen Schlitter, E-mail: juergen.schlitter{at}rub.de.

	ABSTRACT

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
The GTPase Ras p21 is a crucial switch in cellular signal transduction. Fourier transform infrared (FTIR) spectra of the substrate guanosine triphosphate (GTP) show remarkable changes when it binds to the enzyme. The reduced band widths indicate that the flexible GTP molecule is guided by the protein into a preferred conformation. The delocalized phosphate vibrations of unbound GTP become localized. The frequency shifts show an electron movement toward ß-phosphate, which probably contributes to catalysis by reducing the free activation energy. To quantify these qualitative observations we performed QM/MM molecular dynamics simulations of Ras·GTP and GTP in water. The triphosphate part of GTP was treated quantum mechanically using density functional theory (DFT). Vibrational spectra were calculated in harmonic approximation with an average deviation of 3% from the experimental frequencies. This provides a high confidence in the computational results as vibrational spectra are highly sensitive to conformation and charge distribution. As compared to GTP in water, Ras-bound GTP shows a shift of negative charge of 0.2 e toward the ß-phosphate from -phosphate and from -phosphate due to the positive charge of the magnesium ion, to a lesser extent of Lys-16, and surprisingly without any effect of the P-loop backbone. Magnesium and Gly-13 twist and bend the -O-ß bonds such that the crucial bond is stretched before cleaving.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Proteins of the Ras superfamily play a central role in signal transduction including cell proliferation and differentiation. They switch between an ON state where guanosine triphosphate (GTP) is bound and an OFF state with bound guanosine diphosphate (GDP) (Barbacid, 1987). Ras itself is activated by a guanine exchange factor that catalyzes the exchange of GDP by GTP. Activated Ras interacts with effector molecules that transmit the corresponding signal to the cell nucleus (Wiesmuller and Wittinghofer, 1994). Ras is deactivated by GTP hydrolysis to GDP and released Pi. This reaction can be accelerated by five orders of magnitude by a GTPase-activating protein (GAP) and thereby regulated. The hydrolysis rate is decisive for the signal transduction: the longer Ras remains activated the longer the signal is transmitted and amplified. It has been shown that mutations at the positions 12, 13, or 61 in Ras are oncogenic because GTP can still be bound but GAP can no longer catalyze the hydrolysis (Wittinghofer et al., 1997). These malfunctioning Ras proteins are found in at least 30% of human tumors (Bos, 1989; Vetter and Wittinghofer, 2001) emphasizing its central role in signal transduction.

It has been shown by x-ray structure analysis that the structural basis for all members of the Ras family, the G domain, is highly conserved (Vetter and Wittinghofer, 2001). This implies a similar reaction mechanism for GTP hydrolysis of all members of the Ras family. On the other hand structural models provided by x-ray structural analysis of oncogenic Ras mutants have not resolved deviations in the structural models between oncogenic mutated and wild-type proteins that convincingly explain why the hydrolysis rate is dramatically slowed down by oncogenic mutations in Ras (Krengel et al., 1990). Therefore, only higher spatial resolution below 1 Å might reveal those details that govern the catalysis. Furthermore, factors that were not provided by the x-ray structure analysis like the charge distribution in GTP, seem most important. Based on time-resolved Fourier transform infrared (FTIR) experiments it was elucidated that a shift of electronic charge toward ß-phosphate due to Ras binding induces already in the educt GTP a more product-like charge distribution similar to GDP (Allin and Gerwert, 2001). This charge shift appears to be even enhanced by GAP binding (Allin et al., 2001), and was recently confirmed by another approach, isotopic exchange measurements (Du et al., 2004). Because this is the only effect of GAP binding on GTP, it is proposed that the charge shift is crucial for the catalysis (Allin et al., 2001). Recently, a charge shift toward ß-phosphate was also found for another Ras family member, the rap GAP, (Chakrabarti et al., 2004) as also revealed by time-resolved FTIR studies. This points to a general feature in the Ras superfamily.

However, the infrared (IR) data provide only qualitatively the charge shift. By combined quantum mechanical and molecular mechanical (QM/MM) calculations structural details are obtained with subangstrom resolution and charge distribution of GTP bound to Ras can be elucidated. IR frequencies, band widths, and intensities of the GTP vibrations provide the necessary link to the experimentally determined FTIR spectra, which is an extremely sensitive check. Here the first QM/MM calculated IR spectra of GTP bound to Ras are presented and compared with experiments. It will be shown that 10 cm^–1 frequency shifts represent bond length changes of 0.01 Å.

Hydrolysis takes place after nucleophilic attack of a nearby water molecule in principle via a dissociative or via an associative transition state (Fig. 1). Despite a lot of efforts to track the actual pathway, this issue still remains unresolved. In an associative transition state a penta-coordinated phosphorus atom would be formed, in a dissociative transition a trigonal metaphosphate.

View larger version (15K): [in this window] [in a new window]   FIGURE 1  Shown are the two possible transition states of the associative and the dissociative reaction pathway of the hydrolysis of GTP. The arrows are pointing in the direction of hydrolysis.

Theoretical approaches up to now have tried to determine the pathway of the hydrolysis reaction and used essentially only one parameter, the free activation energy, for comparison with the experiments. The measured rate constant for the intrinsic GTPase mechanism is 4.7 x 10^–4 s^–1, the corresponding activation free energy is G = 22.1 kcal/mol (Kötting and Gerwert, 2004; Temeles et al., 1985). GTP in solution is barely hydrolyzable. The rate constant is 3 x 10^–8 s^–1, with an activation free energy of G = 27.9 kcal/mol. Coordination of a magnesium ion to GTP in solution does not change this activation barrier (Kötting and Gerwert, 2004).

Futatsugi et al. (1999) set up a model system for the active site of Ras consisting of -phosphate and ß-phosphate, the magnesium ion, the amino group of Lys-16 and the attacking water molecule mimicking the most dominant interactions of GTP with its immediate environment. For this model the reaction path of the intrinsic hydrolysis was determined with a Hartree-Fock approach. Their analysis suggests that Lys-16 initiates the hydrolysis by mediating a proton transfer from -phosphate to ß-phosphate. The calculated free activation energy of 42 kcal/mol deviates considerably from the measured activation energy of 22 kcal/mol, suggesting that the model system should be extended to include larger parts of the GTP environment.

Cavalli and Carloni (2002) concentrated on the structural and qualitative aspects of the GAP-catalyzed reaction. Their model system (Cavalli and Carloni, 2002) of the binding pocket of Ras consists of the triphosphate, magnesium ion, surrounding water, and residues of Ras and GAP in the vicinity of the triphosphate cut out from the crystal structure of a transition state analog, where GTP was substituted with GDP and AlF₃. This model was described with a density functional theory approach (DFT) using the Car-Parrinello method for ab initio molecular dynamics (MD) simulations (Car and Parrinello, 1985) implemented in CPMD (Hutter et al., 1999). The hydrolyzing water molecule was successively moved toward the -phosphate by applying a constraint force. They observed that this water molecule established hydrogen bonds to Gln-61 and Thr-35 and that a proton transfer to Gln-61 lead to a stabilized OH^– ion. The hydroxyl ion eventually initiated the hydrolysis reaction in their simulation. This mechanism is in agreement with the proposal of Franken et al. (1993). However, substitution of Gln-61 by an artificial amino acid with a pK below zero that cannot be protonated, did not change the reaction rate at all, which makes the reaction pathway very unlikely (Chung et al., 1993; Schweins and Warshel, 1996).

Warshel and co-workers examined the reaction path of the hydrolysis of the Ras·GAP with a hybrid force field (Glennon et al., 2000). The hydrolyzing water molecule and -phosphate were treated with the empirical valence bond method (Aqvist and Warshel, 1993) and the remaining part with a classical force field. The free energy of the reaction paths was evaluated with umbrella sampling and the free energy perturbation technique (Warshel, 1997). The reaction barriers of the hydrolysis in Ras·GAP relative to Ras were compared with each other. The associative and the dissociative reaction paths were taken from the hydrolysis of phosphomonoester in solution treated with ab initio methods. As a result it was proposed that GAP stabilizes the transition state by direct interaction with GTP, especially via Arg-789 and by an allosteric effect on Ras caused by the binding of GAP. The calculated reduction of the free-energy activation barrier was G = 7 kcal/mol and comes close to the measured value of 5.8 kcal/mol. These findings applied to both possible reaction paths but the associative reaction path was suggested due to a little smaller energy barrier. Also the role of the oncogenic Gln-61 was studied and an influence of this residue on the hydrolysis via direct interaction was excluded. In a recent work, Shurki and Warshel showed the indirect influence of Gln-61 on the GAP-catalyzed hydrolysis (Shurki and Warshel, 2004) and proposed that oncogenic mutations of Gln-61 destabilize the transition state. Only the associative case was analyzed and free energy activation barriers for the hydrolysis of GTP in water, bound to Ras and bound to Ras·GAP including some mutations at position 61 were determined using the empirical valence bond method. The calculated energy barriers were in good agreement with measured ones. The calculations showed that a charge shift toward ß-phosphate is catalytic for both pathways.

In the very recent work of Topol et al. (2004) a small and a larger model system of the Ras-GAP were studied using a QM (DFT) and a QM/MM approach, respectively. Similar reaction pathways were found for both systems. The reaction starts with a single-step approach of water to the -phosphate leading to an intermediate with already stretched bond to the bridging oxygen, which is interpreted as a substrate-assisted dissociative mechanism. Accumulation of negative charge on the ß-phosphate is observed on the way to the transition state and in the intermediate, in nice agreement with the FTIR studies.

We conclude from the recent investigations that realistic consideration of the environment is indispensable, which requires application of a QM/MM method. Beside the ones quoted above other approaches successfully used for proteins were introduced by Tavan (Eichinger et al., 1999) and Rothlisberger (Laio et al., 2002). The ab initio description of GTP is demanding due to the high polarizability of GTP that leads to a high sensitiveness to its environment. Of equal importance will be the link to spectroscopic experiments. Recent QM/MM computations have shown that spectra are indeed accessible by these methods (Hayashi et al., 2001; Klähn et al., 2004; Rajamani and Gao, 2002; Rousseau et al., 2004).

In our approach presented here the complete Ras-GTP system within its natural water and ion environment is used. The triphosphate (TP) group of GTP (called TP in the following text), is treated with ab initio accuracy using DFT, whereas the remaining GTP part, Ras, and the surrounding water and ions are described with a classical molecular mechanical (MM) force field. In this first approach we calculate the GTP structure and charge distribution as bound to Ras and compare it with GTP solvated in water. The differences between the two systems, one with hydrolysis activity and the other without, should reveal which factors contribute to catalysis. The resulting structure of GTP is expected to be beyond the resolution limit of the x-ray structural model and will set the stage for the next step in which we calculate the hydrolysis reaction.

The largely fluctuating electrostatic field at the substrate caused by the flexible structure of the protein demands a statistical analysis. From the thermodynamically equilibrated simulated system determined with a long-time simulation with MM force field we perform consecutive short simulations with a DFT-MM force field starting from different structures of the long trajectory. The MM-trajectory samples the coarse structure of the conformational space of the system, whereas the DFT-MM trajectories sample the local fine structure. The final charge distribution is an average over the sample of starting coordinates and subsequent DFT-MM trajectories. Analogously the IR spectra of GTP bound to Ras are calculated and compared to FTIR and Raman spectra. They provide a very sensitive check of the theoretical approach as they depend on the second derivatives of the Born-Oppenheimer surface. With a deviation of only 3% from and measured IR spectra, the calculations reach the accuracy recently obtained for small solvated monophosphates (Klähn et al., 2004) and show that this simulation approach provides an adequate description of the Ras-bound GTP structure and charge distribution.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Applied simulation method and software
The dynamics of GTP bound to Ras and solvated in water-ion environment were simulated by using MD simulations. To sample the temporal coarse structure of the dynamics we used a pure MM force field to calculate trajectories of several hundred picoseconds length. To examine the fine structure of the occupied phase space of the system at points taken from the MM trajectory we used a hybrid force field where the electrons of a part of the system are treated with DFT and the other part with the MM force field as before. These shorter DFT-MM MD trajectories cover a simulation time of 125 fs each. The approach of MD simulation with DFT-MM force field taking into account the mutual interactions of the DFT and the MM part is implemented in the parallelized MD software EGO-MM II (Mathias et al., 2002). In EGO-MM II the long-range Coulomb interactions are treated by a multiple moving-boundary reaction field approach (Mathias et al., 2003). Within that boundary the electrostatics is treated explicitly using a structure-adapted multipole method (Niedermeier and Tavan, 1994, 1996). The electrostatic field, including the reaction field correction, generated by the environment of a respective DFT fragment is imported into the DFT Hamiltonian (Eichinger et al., 1999). After each integration step the coordinates of the DFT fragment are handed over to the parallelized ab initio program CPMD (Hutter et al., 1999). CPMD determines the charge distribution for each configuration of the fragment and the corresponding discrete partial charges by fitting the electrostatic surface potential (ESP) generated by the continuous charge distribution (ESP charges) (Singh and Kollman, 1984). In CPMD the Kohn-Sham wave functions and the charge distribution are expanded in plane waves making it necessary to use pseudopotentials that replace the explicit treatment of core electrons.

The suitability of our DFT-MM approach concerning the simulation of IR spectra was already shown for polar DFT fragments in similarly polar environments (Eichinger et al., 1999; Nonella et al., 2003). In our preceding work on the simulation of IR spectra of solvated monophosphates (Klähn et al., 2004) we studied the quality of the method on an ionic DFT fragment where the polarization of the environment is expected to be much stronger than in cases analyzed previously. We had suspected the missing description of polarization for the MM fragment and the lack of possible charge transfer from the ion to the environment that was not accounted for in the model. It turned out, however, that specific vibrations whose frequencies deviated strongly from the experiment always involved hydroxyl groups of the protonated monophosphates (for further discussions we refer to this work). GTP in Ras and aqueous solution, however, is known to be deprotonated and lacks any hydroxyl group at the phosphorus atoms. So we confidently applied the method to these systems keeping in mind that the vibrations of the pure PO groups that are chemically very similar to the GTP phosphates, had been described satisfactorily.

Setup of the Ras·GTP system
The initial positions of the atoms of the protein and its substrate were taken from the crystal structure of Ras·GTP at 100 K determined by Scheidig et al. (1999) (Protein Data Bank code, 1QRA). Due to the low temperature no GTP analog had to be used to crystallize the educt state of the hydrolysis, and the resolution is as high as 1.6 Å. The unresolved protons were added to the system. The protonation of ionizable amino acids was chosen according to their pK_a values in solution at neutral pH because they are almost all found on the protein surface with sufficiently large mutual distances. To determine their pK_a values elaborated calculations would be required (see, for example, Juffer, 1998; Spassov et al., 2001), which seems not be necessary here. The nearest carboxylate group of Asp-33 lies at >7 Å from the negatively charged GTP on the surface, points away from it, and is separated from it by two waters. Lys-16, which is directly coordinated to GTP is expected to be protonated, i.e., positively charged in agreement with Chung et al. (1993) who studied point mutations at that position. To simulate natural conditions found in cells and similarly in IR experiments, Ras·GTP was solvated by putting it into a filled water box representing elementary cell of the simulated system on which periodic boundary conditions were applied. The box was chosen to be a rhombic dodecahedron to minimize the number of water molecules. The diameter of the inscribing sphere of the rhombic dodecahedron was 8.34 nm. Interaction of the protein with its mirror images was avoided with a Coulomb interaction cutoff of half the diameter of the inscribing sphere of the simulation cell. The value of 4.17 nm is also the maximum distance of two atoms of the Ras protein.

To fill the simulation cell 12,730 water molecules were needed. As NaCl is ubiquitous in the natural cell and solutions used for IR experiments, 29 Na⁺ ions were added to simulate a salt concentration of 130 mmol/l. At –27 e from Asp/Glu, +19 e from Lys/Arg, and a GTP/Mg²⁺ with –2 e, the total charge of the complex amounts to –10 e. As this charge must be counterbalanced in the immediate vicinity, only 19 Cl^– ions were added. The total system contains 40,497 atoms.

DFT method
For the DFT Hamiltonian we used the functionals of Becke (1988) and Perdew (1986) (BP86) and the pseudopotential from Hartwigsen and co-workers (HGH) (Hartwigsen et al., 1998). The cutoff for the plane wave expansion was set to 80 Ry. In our previous work we have already demonstrated the quality of this HGH/BP Hamiltonian; i.e., it is appropriate to describe phosphate spectra and the corresponding charge distribution (Klähn et al., 2004). For DFT-MM calculations the triphosphate of the GTP substrate together with its adjacent CH₂ group were treated quantum mechanically whereas the remaining part of the system was described by an MM force field. Both fragments were connected via a link atom changing the CH₂ group to a methyl. The hydrogen link atom was dynamically positioned to mimic the actual bond of the CH₂ to the ribose ring. The partition of the GTP substrate into a DFT and an MM fragment is shown in Fig. 3. To allow charge transfer and a polarization of the environment we should have applied DFT to the immediate neighboring atoms of GTP as well. At least in the case of solvated GTP where no dominant interaction partners are present, the whole first hydration shell should have been included. The treatment of systems of this size would have been beyond our available computational resources. So we decided to describe all non-GTP atoms in both simulated systems with an MM force field to maximize the comparability of the partial charges derived from the charge distribution of GTP. The spatial size of the DFT box in which the charge distribution is calculated on a grid, was always dimensioned so that the minimum distance between an atom of the DFT fragment and the edge of the box was 2.7 Å like before (Klähn et al., 2004).

View larger version (9K): [in this window] [in a new window]   FIGURE 3  The triphosphate and the adjacent CH2 group, TP–CH2 (black) of GTP were treated quantum mechanically (DFT fragment) and the remaining part (gray) molecular mechanically (MM fragment).

View larger version (14K): [in this window] [in a new window]   FIGURE 2  (a) Closest internal contacts of GTP in the crystal structure 1QRA are to Mg2+ and H-bonds to residues and backbone amides of which 15–18 represent the P-loop. One oxygen of P is exposed to solution. (b) The snapshot taken after 350-ps production run shows conservation of the octahedral coordination of Mg2+ in the simulation system.

Despite the satisfactory description of vibrational frequencies, it seems possible that the magnesium ion not only shifts electron charge on the substrate, but also subtracts electron charge. This effect has been studied numerically (Stefanov and Tulub, 2002) on different configurations of the model system (Mg²⁺ [H₂O]_n Glu-DP-CH_3.) where the magnesium charge computed with the Lowdin approach turned out to vary only from 1.94 to 1.98 e. Mulliken charges for Mg²⁺ in the F₁-ATPase were found to be smaller and suggest that the ion may serve as an electron sink (Dittrich et al., 2004). In view of such differences the model could be improved by including a larger environment of GTP that is treated quantum mechanically to account more precisely for charge transfer and polarization of neighboring groups. Additional quantum mechanical treatment of the magnesium ion is expected to improve the model only if the magnesium's complete coordination sphere shown in Fig. 2 b is included as well (e.g., Hong et al., 2000).

Adapted force field
The MM force field we used for the MD simulations is based on the CHARMM22 force field (file "parallh22x.pro") (MacKerell et al., 1998). For the water we applied a variant of the TIP3P model (Jorgensen et al., 1983), the same force field we used in our previous work (Klähn et al., 2004). For the nonbonding forces of TP the van der Waals parameters were taken from the CHARMM22 force field (file "parallh22x.nuc") and the partial charges were initialized according to ESP charges of the optimized HGH/BP structure of TP–CH₂ in vacuum. These charges are shown in Fig. 4 in parentheses. Note, however, that these charges in vacuum are considerably higher than what we expect for TP in solution or in Ras, so the partial charges of TP had to be replaced later by ESP charges from (HGH/BP)-MM MD simulations, where the whole system was considered.

View larger version (35K): [in this window] [in a new window]   FIGURE 4  These are the ESP charges for TP–CH2 that were used for our MM force field. They were derived by a (HGH/BP)-MM MD simulation of Ras·GTP. The vacuum charges that were used for the force field as initializations are given in parentheses.

The beginning of the MD simulation, where our pure adapted MM force field was used, was divided into several successive simulation phases. In these phases the degrees of freedom of the molecular system were gradually increased to ensure a smooth equilibration and to avoid artificial deviations from the initial crystal structure of the protein: after construction of the model system, an MM-MD simulation of 200 ps simulation time was carried out with all atoms fixed except the noncrystal water molecules and a deactivated barostat to allow the solvent a first relaxation and dispersion around the protein. During the next 400 ps also the Na and Cl ions were allowed to move freely. During the next 100 ps harmonic restraints were applied to all C atoms of the protein backbone, the GTP, the Mg-ion, and the crystal water. Beyond this point neither constraints nor restraints were used anymore. In the next step the partial charges of TP for the MM force field had to be improved. For that, after further 50 ps, the ESP charges for TP–CH₂ were derived from a subsequent (HGH/BP)-MM trajectory (500 integration steps with t = 0.25 fs) as averages over simulation time and chemically equivalent atoms. The results are shown in Fig. 4.

With the adjusted MM force field we continued the MM-MD simulation for further 1000 ps. The first 500 ps of this phase contributed to the equilibration of the system and the last 500 ps were used for evaluation. The simulation preserved all hydrogen bonds in the binding pocket of Ras seen in the crystal structure except a cleaved hydrogen bond between Gly-60 and -phosphate. The conformation of GTP was left unchanged with the magnesium ion and Lys-16 coordinated tightly right between ß- and -phosphate like in the crystal structure.

Cation coordinated to -phosphate
During equilibration MD, a Na⁺ ion diffused quickly toward the active site where the negative -phosphate group is open for solution contacts, and remained fixed at this position for the rest of the simulation. Restarts with different starting coordinates confirm the result. Two adjacent water molecules stayed coordinated to the localized sodium ion during the whole simulation, whereas all other waters are very mobile and permanently exchange their positions. The stability of the configuration seems reasonable because the cation compensates the net charge of the active site (GTP, Mg²⁺, and Lys-16⁺), which is –1 e. More generally speaking, there is a high probability to find a cation here that may be a Na⁺ or hydronium H₃O⁺, which cannot be distinguished from H₂O by x-ray scattering. We also did a simulation without Na⁺ at the -phosphate for estimating the influence on the vibrational spectra.

Solvated GTP
It is well known that the hydrolysis of solvated GTP without catalysis is very unlikely, the half-life period of the reaction is 250 days (Klähn et al., 2004). To understand the hydrolysis of GTP catalyzed by Ras it is important to compare this simulated system with a simulation of solvated GTP. The analysis of the differences of both systems in structure and charge distribution of TP is crucial to understand the enzymatic capabilities of Ras. For the model system we assumed deprotonated TP according to measurements from Sigel et al. (2001). The inclusion of counterions was abandoned because the position of such an ion is not known and difficult to predict. An artificial coordination of the counterion to GTP would change the structure and charge distribution of TP substantially so that the comparability with TP in Ras would be lost. It is known that solvated GTP without coordinated counterion undergoes a very slow hydrolysis reaction and it was also measured that a counterion changes the structure of GTP but reduces the rate constant only slightly from 3.1 x 10^–8 s^–1 to 2.8 x 10^–8 s^–1 at room temperature (Kötting and Gerwert, 2004).

The GTP was simulated in a rhombic dodecahedron elementary cell with a diameter of the inscribing sphere of 5.2 nm. To solvate GTP in that cell 1375 water molecules were needed; no additional ions were added. The coordinates of the GTP atoms were initialized according to the geometry in the crystal structure of Scheidig et al. (1999). All simulation parameters were chosen to be identical with the parameters for the simulation of Ras to maximize the comparability of the two systems. Only the simulation time to equilibrate the system was reduced. During the first 300 ps the water molecules were allowed to equilibrate while GTP was constrained. Additional 300 ps were simulated to equilibrate the whole system and the following 500 ps of the trajectory were used for evaluation.

Generation and analysis of structure samples
To generate a representative and accurate structure sample for each simulated system we proceeded the same way as in our previous work. Starting from the equilibrated systems as described above an MD simulation with pure MM force field was performed to obtain six snapshots of the respective MM trajectory, with a simulation time interval of 100 ps between two snapshots that ensures that the snapshots are statistically independent, defining the sample. At each snapshot we then switched to the DFT-MM force field and computed short (HGH/BP)-MM trajectories of 500 integration steps with a step width of 0.25 fs subsequently. These trajectories were evaluated with respect to interesting observables, i.e., the charge distribution and the IR spectra of TP in the two different environments. Following our previous approach we obtained ESP charges for each TP–CH₂ atom as double averages. We averaged the ESP charges first over the last 400 integration steps of our single (HGH/BP)-MM trajectories. We neglected the first 100 steps during which the system underwent a considerable part of a structural relaxation after switching the force field. We then averaged the ESP charges a second time over the sample.

To determine the total IR spectra of GTP we applied the instantaneous normal mode analysis (INMA) using the protocol we have suggested in our previous work. The attractive alternative way is to compute spectra by Fourier transformation from sufficiently long trajectories, see, e.g., Mavri and Grdadolnik (2001a,b). It was rejected because it requires longer trajectories, i.e., more computing time. Starting from the (HGH/BP)-MM structures geometry optimizations of the DFT fragment were done until the maximum force on a single atom dropped below 8.4 x 10^–5 au, whereas the environment described by the MM force field was constrained to its 300-K structure. In the next step the Hessian matrix of the potential energy and the dipole gradient of TP–CH₂ were calculated using analytical first derivatives and numerical second derivatives. With these two values, normal mode analysis was performed resulting in an IR single spectrum for each structure of the sample. The absorption frequencies were left unscaled because it is known that the density functional BP86 is able to predict unscaled absorption frequencies in the harmonic approximation accurately (Neugebauer and Hess, 2003). The fluctuation of the absorption frequency of a single normal mode in the sample gives the corresponding band width and the average absorption frequency the band position. With these two measures we were able to build the total IR spectra by a superposition of Gaussian functions for each IR band, weighting each with the average calculated IR intensity of the corresponding normal mode.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS ACKNOWLEDGEMENTS REFERENCES

 
Simulated spectra of the GTP phosphate bands
The IR spectra of TP of the Ras substrate GTP obtained by our simulation approach offer the most sensitive probe for comparison with biophysical experiments. In particular FTIR spectroscopy has been applied successfully to the Ras protein for elucidating the intrinsic (Allin et al., 2001; Cepus et al., 1998) and the GAP-catalyzed hydrolysis mechanism of Ras (Allin and Gerwert, 2001). Substitution of GTP with caged GTP as a photolabile trigger compound, measuring time-dependent difference spectra of the hydrolysis and assigning bands by isotopic labeling of GTP are the most important tools for obtaining and interpreting the IR spectra (Kötting and Gerwert, 2005). With these techniques it was possible to identify and to assign vibrational modes to the TP bands of GTP bound to solvated Ras in the frequency region of 850–1300 cm^–1. Raman spectroscopy has also been successfully applied to Ras providing complementary band assignments for the TP modes (Wang et al., 1998). Combining the band assignments due to FTIR and Raman spectroscopy the characteristics of nine normal modes of TP can be compared with calculated TP modes.

For comparison with the experimental data we determined the total IR spectrum of GTP bound to Ras with our INMA protocol as described above, including averaging over the sample described in the Methods section. The assignment of the calculated normal modes of TP to measured bands was straightforward due to the fact that the vibrations of the single phosphates of TP are decoupled in Ras (Allin and Gerwert, 2001) as already mentioned above. With isotopic labeling it was possible to identify three -, two ß-, and two -phosphate bands. It can also be readily differentiated between asymmetric and symmetric vibrations due to their dominant absorption intensities in FTIR and Raman spectra, respectively. In our calculations the decoupling of the phosphate modes was also observed and allowed a relative simple evaluation of vibrations. So we assigned our calculated normal modes straightforwardly to the measured bands by visual inspection of the corresponding normal amplitudes, where we could easily identify the particular vibrating phosphate group and the phase of the vibration. Asymmetric and symmetric modes were denoted by _a/s(PO)_/ß/, respectively, where the last index designates the phosphate group on which the vibration is localized. All these modes consist basically of P-O stretch vibrations with nonbridging oxygens. We also found two normal modes of P-O stretch vibrations with the two bridging oxygens that were omitted in Fig. 5 and will be discussed elsewhere together with the experiments confirming the calculations. These and the other normal amplitudes of the modes for GTP bound to Ras can be viewed as animations on our web site www.bph.rub.de. Four modes are displayed in Fig. 5 as examples.

View larger version (16K): [in this window] [in a new window]   FIGURE 5  The normal amplitudes of the four a/s(PO2)/ß and three a/s(PO3) modes of one sample structure of GTP bound to Ras. The arrows are pointing away from the equilibrium positions of the nuclei. The display of marginal amplitudes of single atoms were suppressed to make the figure more concise.

View larger version (19K): [in this window] [in a new window]   FIGURE 6  IR absorption frequencies calculated with our INMA protocol (middle column) and measured with FTIR/Raman of GTP in Ras (right column) are compared in this figure. The measured frequencies (in cm–1) and the assignment of types of vibration (notation is explained in the text) are given on the right side. Additionally the left column shows calculated absorption frequencies of a single spectrum for the case where the sodium ion that coordinates to the -phosphate was replaced by water (explained in the text). For all calculated frequencies the deviation from measured ones are given.

To check the influence of the cation on the IR spectrum we also studied the change of the absorption frequencies of TP; we removed that ion in our simulation and subsequently equilibrated the system to a (metastable) state without ion at that site. Then we calculated one single spectrum with the INMA method as before and compared in Fig. 6 the result with the spectrum of the complete system. The spectrum of the system without sodium ion exhibits a lesser agreement with the experiment. Especially interesting are the two _a(PO) modes. These modes are degenerated in solution due to the symmetry of the -phosphate, which was confirmed by gas phase measurements and calculations. FTIR measurements show that in Ras the degeneracy is broken only slightly with an energy difference of 14 cm^–1 (Cepus et al., 1998). This implies that the -phosphate essentially retained a symmetric structure in Ras. According to Fig. 6, the almost-degeneracy is best reproduced in the model with cation at the -phosphate.

Given the asymmetric charge structure around the three -phosphate oxygens that participate in the two _a(PO) modes in the crystal, it is difficult to understand the observed degeneracy. One oxygen is coordinated to the Mg²⁺ ion, one to Lys-16 and one only to water. Apparently only an additional cation can counterbalance the electrostatics at the -phosphate, thus reducing the symmetry breaking and restoring the degeneracy of the two modes. We found in our calculations that the additional sodium ion reduces the energy difference from 125 to 69 cm^–1. Although still larger than the measured difference of 16 cm^–1, it provides an additional argument for its presence. Notwithstanding the most natural occurrence of the cation at the solvent-exposed end of GTP described in the Methods section, this finding makes clear that we had no reason to reject the presence of a sodium ion at -phosphate as a result of our simulations. The favorable steric arrangement of charges is best understood from Fig. 8.

View larger version (28K): [in this window] [in a new window]   FIGURE 8  Mean ESP partial charges on the TP–CH2 atoms for solvated Ras (top; selected snapshot structure) and Ras bound GTP (bottom). The values in parentheses (bottom) are single values for a configuration without cation at the -phosphate.

View larger version (23K): [in this window] [in a new window]   FIGURE 7  Photolysis (Ras·GTP–Ras·caged GTP) FTIR difference spectrum of GTP bound to Ras (solid line) and the corresponding calculated INMA spectrum (dashed line). For the INMA spectrum the band positions were shifted toward their measured positions to simplify the comparison of band widths and intensities. Given are also the assignments of modes to the bands.

To check the quality of the ESP charges we compare the charges on chemically equivalent atoms of solvated TP–CH₂. Due to the isotropy and homogeneity of the H₂O solvent the three -nonbridging oxygens and the two ß-nonbridging oxygens have naturally the same charge. Actually the two ß-nonbridging oxygens differ by only 0.02 e and the three -nonbridging oxygens do not differ at all. So all the mean calculated ESP charges of chemically equivalent atoms of the solvated GTP have equal charges, indeed, in consideration of the calculated error of 0.03 e, proving the reliability of these charges additionally. The charges of the two -nonbridging oxygens differ by 0.06 e due to the interaction with the adjacent guanosine that cancels the equivalence of the two atoms. This check is actually the only possible straight test for the quality of the ESP charges because partial charges are not observables; they are not uniquely defined. The charge distribution of - and ß-phosphate are similar due to their chemical equivalence and the same environment whereas negative charge is accumulated at the -phosphate as expected.

For GTP bound to Ras the environment exhibits neither isotropy nor homogeneity and the charge distribution is determined by dominant interactions with counterions and amino acids of Ras that constitute hydrogen bonds to TP. Note that our simulation system conserves the octahedral coordination of Mg²⁺ observed in the crystal structure, in particular the bidentate binding to - and ß-phosphate shown in Fig. 2. The calculated averaged ESP charges are given in Fig. 8 The most dominant interactions of TP with its environment are highlighted, i.e., the interaction with the Mg²⁺ counterion, with Lys-16 and the coordinated cation. As a consequence of the interactions the charge distribution of TP in Ras is less uniform than of solvated TP leading to the known decoupling of vibrations. The ß-phosphate is polarized by Ras, i.e., negative charge from the ß-phosphorus is moved toward the ß-nonbridging oxygens, whereas the -phosphate is depolarized by Ras. Consistently with these charge movements the average bond length of the ß-phosphorus to the ß-nonbridging oxygens is increased from 1.511 to 1.518 Å and vice versa the average bond length of the -phosphorus to the -nonbridging oxygens is decreased from 1.508 to 1.504 Å due to binding to Ras. These findings are supported by FTIR measurements that observed a red shift of 14 cm^–1 for the _a(PO)_ß mode of GTP due to binding to Ras corresponding to a decrease of bond strength of the participating bonds and a blue shift of 30 cm^–1 for the _a(PO) mode corresponding to an increase of bond strength (Cepus et al., 1998; Wang et al., 1998). The calculations show also that the bond lengths are quite insensitive to changes of the charge distribution and frequencies. A change of 10 cm^–1 in frequency involves a change of <0.01 Å in length. We like to note that ab initio calculations cannot determine absolute bond lengths at this accuracy. However, as the statistical error of the bond lengths is as small as 0.002 Å, it seems possible to derive changes of lengths at this scale by using our means.

The situation is more complicated in the case of -phosphate. Due to the binding to Ras the symmetry of the -phosphate is broken and the degeneracy of the two _a(PO) modes is canceled. This symmetry breaking is a result of different bond lengths and strengths of the -phosphorus to the three -nonbridging oxygens. In solution all three bond lengths are 1.527 Å but binding to Ras resulted in bond lengths of 1.520, 1.526, and 1.552 Å, respectively. The calculated _a(PO) mode with higher energy involves a vibration along the two short bonds, whereas the other _a(PO) mode consists mainly of a single stretch vibration along the remaining stretched bond. That means that the high-energy _a(PO) mode involves vibrations of strengthened bonds, which is consistent with a measured FTIR blue shift of 41 cm^–1 for this mode. The low-energy _a(PO) mode vibrates along a weakened bond, whereas for this mode a blue shift of 27 cm^–1 is observed. Considering that the calculated absorption frequency is erroneously 53 cm^–1 red-shifted compared to the measured frequency, we have to conclude that the stretching of the corresponding bond is an artifact of our calculations.

The binding to Ras also induced a charge movement from the - and ß-phosphorus atoms toward the -ß bridging oxygen. This important result means that the corresponding bonds are weakened including the bond that has to be cleaved during the hydrolysis reaction. The bond lengths from - and ß-phosphorus to the -ß bridging oxygen are increased from 1.704 to 1.740 Å and from 1.581 to 1.605 Å, respectively.

Combining all these results we find that, except for one -phosphate mode, all Ras-induced band shifts and charge shifts observed in our simulations are in good agreement with the band shifts seen in FTIR measurements and the charge shifts derived from them via changes in bond order (Allin and Gerwert, 2001). We also observe that the considerable charge shifts are accompanied by changes in the bond length that are too small to be resolved by x-ray spectroscopy, which is a well-known fact. FTIR difference spectroscopy and QM/MM calculations arrive at a consistent picture of these phenomena.

Fig. 9 summarizes the total charges of the three phosphates and the CH₂ group shown for solvated TP and Ras-bound TP and the shifts induced by Ras binding. The mean error of the average charge per group is 0.03 e for solvated GTP and 0.01 e for GTP bound to Ras.

View larger version (12K): [in this window] [in a new window]   FIGURE 9  The total charges of -, ß-, and -phosphate and the CH2 group are given in units of the elementary charge. The top scheme explains our split-up of TP–CH2 for that purpose. The charges for these groups are given in the rows below for solvated GTP and for GTP bound to Ras. The mean error of the average charge per group is ± SE 0.03 e for solvated GTP and ± SE 0.01 e for GTP bound to Ras. The arrows indicate the direction of Ras-induced shifts of negative charge between adjacent groups and the numbers in the bottom row specify the magnitude of these shifts.

Like the bond lengths, also the absolute values of the charges are certainly less reliable than the trend of charge shifts arising from the substrate-enzyme interaction. We believe that our treatment of neighboring ions (Mg²⁺, Na⁺) as smeared-out point charges is a good starting point. The next step would be a quantum mechanical treatment of the ions including their complete coordination sphere. Regarding only the bare ions as charge donors is expected to give an unrealistic picture of charge transfer (Hong et al., 2000).

Induced fit of the substrate GTP
The most palpable effect on GTP binding to Ras is the general reduction of flexibility and the preference of a particular configuration of the -phosphate relative to the ß-phosphate as shown in Fig. 10.

View larger version (27K): [in this window] [in a new window]   FIGURE 10  Comparison of the six (HGH/BP)-MM structures of solvated GTP (a) and of GTP (b) bound to Ras, respectively. Panels a and b obtained by fitting all atoms with equal weight demonstrate the enormous reduction of flexibility imposed to GTP when binding to Ras. The next figures were obtained by fitting only the ß-phosphorus and its nonbridging oxygens. As compared with the solution state (c) the bound state (d) shows a fixed bond angle between the phosphorus atoms and the favored eclipsed relative orientation of the phosphate groups.

Fig. 10 c shows that in solution already -phosphate and ß-phosphate take an approximately eclipsed relative conformation, the angle varying around 20°. The bound conformation Fig. 10 b, however, is restricted to a more eclipsed narrow distribution around 13°.

For the two nonbridging ß-oxygens atoms in the two samples we noticed a Ras-induced stretching of their bonds to phosphorus from 1.51 Å in solution to 1.52 Å (with a standard deviation in the sample of 0.005 Å). This barely significant change in bond length is accompanied in contrast by a considerable charge shift of 0.09 e as discussed before. Most important is the observed increase of the bond length of the -phosphorus to the -ß bridging oxygen from 1.70 Å in solution to 1.74 Å when bound to Ras. The correlation between the changes of bond length and charge distribution was already discussed above.

The increase of bond length is accompanied by a decrease of the bond angle at the -ß bridging oxygen from 145° in solution to 135°. Although the distance of the phosphorus atoms of - and ß-phosphate remains unchanged (3.1 Å), two nonbridging oxygens of these two phosphate groups, respectively, approach each other due to coordination to the magnesium ion. The distance between these two atoms is reduced from 3.7 Å in solution to 3.1 Å as shown in Fig. 11.

View larger version (31K): [in this window] [in a new window]   FIGURE 11  Some average geometric measures of TP bound to Ras, all distances in Å. The corresponding values for GTP in solution are given in parentheses. The arrow points out the bond that is cleaved during hydrolysis.

The magnesium-induced conformational change can be described as a mutual tilting of - and ß-phosphate toward each other with the two phosphorus atoms as centers of that movement. More precisely, the phosphate groups are twisted relative to the mean solution structure to take a slightly more eclipsed conformation and the P-O-P bond is bent. To retain the apparently stiff tetrahedron geometry of the single phosphates, the -ß bridging oxygen has to move away from the magnesium ion. This movement is supported and oriented by Gly-13 that establishes a backbone hydrogen bond to that oxygen. As a result the two bridging bonds become longer and weaker, negative charge being shifted from the phosphorus atoms to the bridging oxygen. This mechanism is depicted schematically in Fig. 12.

View larger version (10K): [in this window] [in a new window]   FIGURE 12  When bound to Ras the magnesium ion induces a tilting of - and ß-phosphate toward each other. The two phosphates retain their stiff tetrahedron structure, the bridging oxygen moving toward Gly-13.

Next we have to understand why binding to Ras induced a total negative charge shift toward ß-phosphate. For this shift the positive charges of Lys-16 and especially the magnesium ion are responsible. The backbone of the P-loop of Ras that establishes numerous hydrogen bonds to the ß-phosphate contributes only insignificantly to that total charge shift. This is the result of additional calculations that continued our (HGH/BP)-MM MD simulations, but with the partial charges of Mg²⁺, Lys-16, or the P-loop backbone switched off, respectively. A major contribution to the shifted electron charge comes from the -phosphate. This is simply due to the highly positive potential emerging from Mg²⁺ and Lys-16 that attracts electrons from the CH₂-group and the -phosphate into the direction of the ß-phosphate.

It might be surprising that Mg²⁺ and Lys-16 also shift negative charge from the - to the ß-phosphate although they are positioned right between them (see Fig. 11). Naturally in an isotropic and homogeneous environment there is more negative charge at the - than at the ß-phosphate (see Fig. 9) due to the asymmetric structure of TP. In the absence of internal forces, a counterion placed at equal distance from the two charge centers would, of course, induce a rearrangement toward a symmetrical distribution. Hence, the magnesium ion coordinated to the more negatively charged - and the less negatively charged ß-phosphate moves negative charge to the ß-phosphate, thus making the charges on both groups more equal. The situation in Ras is clearly different from that in aqueous solution where magnesium tends to bind to GTP in a tridentate manner (Wang et al., 1998). This type of interaction apparently has little or negligible catalytic effect (Kötting and Gerwert, 2004). Therefore it is not only the presence of Mg²⁺, but also its specific position between the ß- and -phosphate enforced by the enzyme that is a necessary condition for catalysis.

Accumulation of negative charge on the ß-phosphate means that the product of the hydrolysis reaction is approached, which requires a negative charge transfer from the - to the ß-phosphate. This would also be characteristic of a hypothetical product-like dissociative transition state. All findings suggest that the dissociative reaction pathway of the hydrolysis reaction is favored, whereas hydrolysis via the associative reaction pathway is inhibited due to the loss of negative charge on the -phosphate. The two corresponding transition states are shown in Fig. 1. We like to mention that a sharp distinction between the types of mechanism may be questionable and should be made with caution.

We come to the conclusion that there is a local structural excitation of the potential energy of GTP when it binds to Ras, which remains to be confirmed by an analysis of the further reaction pathway. Given the release of 15.3 kcal/mol substrate binding energy (Rensland et al., 1995), the required energy is available, probably resulting primarily from favorable electrostatic interactions. The increase of energy of GTP in the educt state should reduce the energy gap to the transition state, thereby making the hydrolysis more likely according to transition state theory if the transition state is not shifted upward to the same degree. The results of Warshel and co-workers (Glennon et al., 2000) seem indeed to prove that there is a further catalytic effect by stabilization of the transition state. It is also widely believed that hydrolysis in Ras is a substrate-assisted process where the -phosphate supports the acquisition of H₂O or HO^– in some way. The withdrawal of electron charge from the -phosphate being part of the observed charge shifts will certainly also facilitate the nucleophilic attack.

Ion coordination at the -phosphate
The predominant solution configuration at the -phosphate that is stable throughout the simulation is displayed in Fig. 13. The presence of the cation is due to the net charge –1 e of the active site (GTP^–4, Mg⁺², and Lys-16⁺). The only cations in our solution model are sodium ions, one of which diffuses quickly into the active site and keeps this position for the rest of the simulation period. Restarts from different starting positions of the ion all amounted to the same result. As a consequence of the fixation two adjacent water molecules stay coordinated to this site whereas all other waters are found to be mobile. Instead of Na⁺ also a protonated water complex like as observed in bacteriorhodopsin (Garczarek et al., 2005) seems a possible candidate for the cation that is coordinated to the -phosphate, and there is at present no experimental counterargument as x-ray crystallography cannot distinguish between any of these cations and a water molecule.

View larger version (31K): [in this window] [in a new window]   FIGURE 13  -Phosphate and ß-phosphate of GTP and their coordination to selected groups and ions as seen in a simulation snapshot. Displayed are also the two water molecules coordinated to the sodium ion. The positions a, b, and c correspond roughly to crystal waters 175, 189, and 190 of Scheidig et al. (1999) or the waters 1098, 1117, and 1107 of structure 1QRA, respectively. Average coordination distances are given in angstroms.

The presence of a cation ion in the binding pocket can also explain the fact that GTP bound to Ras is not protonated with a pH above 4.2 whereas solvated GTP is protonated above pH = 6.5 (Cheng et al., 2001; Sigel et al., 2001). The additional cation ion would repel protons. It also offers a possible explanation for another experimental finding, the acceleration of the intrinsic hydrolysis reaction by a factor of 20 (Frech et al., 1994) after the mutation Q61E. In our simulation we saw that Gln-61 was very mobile with an average distan