J. Phys. Chem. B 2007, 111, 13107-13115
13107
Brownian Dynamics Simulations of Binding mRNA Cap Analogues to eIF4E Protein Elz3 bieta Błachut-Okrasin´ ska and Jan M. Antosiewicz* Department of Biophysics, Warsaw UniVersity, Zwirki i Wigury 93 St., Warsaw 02-089, Poland ReceiVed: July 25, 2007; In Final Form: September 7, 2007
The binding of five analogues of the 5′-end mRNA cap, differing in their electrostatic and hydrodynamic properties, to the eukaryotic initiation factor eIF4E was simulated by means of Brownian dynamics methods. Electrostatic and hydrodynamic models of eIF4E protein and the ligands were prepared using established molecular electrostatics and hydrodynamics simulation methods for predicting ionization states of titratable groups, adequate for given experimental conditions, and for computing their translational and rotational diffusion tensors, respectively. The diffusional encounter rate constants obtained from simulations are compared with bimolecular association rate constants resulting from stopped-flow spectrofluorimeter measurements. A very good agreement between simulations and experiments was achieved, which indicates that the kinetics of binding 5′-mRNA caps can be satisfactory explained by referring to the Brownian motion of the particles with the electrostatic steering of the ligands toward the eIF4E binding site and electrostatic desolvation contributions upon complex formation.
1. Introduction The eukaryotic translation initiation factor eIF4E regulates gene expression.1 One of the essential steps in the translation initiation consists of binding of a particular structure at the 5′end of mRNA, called cap, by this protein.2 The cap has the general form m7GpppN1, where m7G is the 7-methylguanosine, and N can be any of the four nucleosides.3 Important insights into structural characteristics of cap-eIF4E complexes were obtained from three-dimensional structures by X-ray crystallography4-6 and multidimensional NMR study,7,8 as well as from spectroscopic investigations in vitro, both equilibrium5,9,10 and kinetic.11-16 The most crucial element of the recognition consists of sandwiching of the alkylated base between the side chains of two conserved tryptophans (Trp-56 and Trp-102 in the murine eIF4E). The interaction can be explained in terms of enhancement of π-π stacking enthalpy because of charge transfer between the electron-deficient 7-methylguanine (which carries a delocalized positive charge secondary to methylation) and the electron-rich indole groups. Moreover, there is a number of hydrogen bonds stabilizing the complex, N1 and N2 make hydrogen bonds with the carboxylate oxygen atoms of Glu103, and Arg157 and Lys162 make such contacts with oxygen atoms of the phosphate groups. Many of these hydrogen-bonding interactions include contributions due to direct electrostatic interactions because the interacting groups are ionized, and these interaction energies include contributions due to shifts in the appropriate pKa values upon binding.13 In our recent study, we measured kinetics of binding of five mRNA cap analogues to murine eIF4E protein using the stopped-flow technique.16 Kinetic transients obtained after mixing equal volumes of the ligand and protein solutions were analyzed using the DynaFit program developed by Kuzmic.17 With this program, one can fit a postulated theoretical model to a given experimental data, or one can simulate pseudo* Corresponding author. Fax: +48-22-55-40-771. E-mail: jantosi@ biogeo.uw.edu.pl.
experimental data by numerical integration methods. The program allows a discrimination analysis of various models that results in a statistically justified choice of the most appropriate reaction scheme for a given system. We analyzed our data in terms of two simple binding models. The first model was a onestep mechanism kR
P + L y\ zK k
(1)
D
And the second model was a two-step binding kinetics, described by k+1
k+2
-1
-2
z C y\ zK P + L y\ k k
(2)
In the two above equations, ki are the appropriate rate constants, C and K represent an intermediate and a final form, respectively, of the complex of protein P and ligand L. Application of the DynaFit program indicated that the binding of cap analogues by eIF4E protein is at least a two-step process; however, it appeared that reproducibility of the rate constants derived from fitting more than one-step model for independent experiments is not particularly good. On the other hand, neither the association rate constants nor the equilibrium binding constants derived from one-step and two-step binding are much different.16 The aim of the present study is to apply a Brownian dynamics (BD) simulation method to see if the diffusional encounter rate constants derived from these simulations correlate with kR or k+1 obtained from fitting stopped-flow data. For BD simulations, we use the SDA program developed by Gabdoulline and Wade.18,19 It appears that both sets of the bimolecular association rate constants correlate well with the diffusional encounter rate constants derived from the BD simulations. The slope of the regression line for the correlation plot having the experimental rate constants derived from the one-step binding model is 0.93 and that for the rates derived from the two-step model is 0.99, which might indicate some advantage of the two-step model.
10.1021/jp0758521 CCC: $37.00 © 2007 American Chemical Society Published on Web 10/20/2007
13108 J. Phys. Chem. B, Vol. 111, No. 45, 2007
Błachut-Okrasin´ska and Antosiewicz
TABLE 1: Bimolecular Association Rate Constants in Units of [µM-1s-1] Resulting from Stopped-Flow Measurements16a binding mechanism
a
cap analogue
one-step kR
two-step k+1
m7GMP m7GDP m7GTP m7GpppG m7Gpppm7G
70.1 ( 1.5 272 ( 23 694 ( 43 191 ( 3.0 139 ( 21
75.1 ( 4.8 235 ( 8.7 649 ( 122 197 ( 5.6 151 ( 13
The errors are standard deviations of the average value.
2. Materials and Methods 2.1. Experimentally Determined Association Rate Constants eIF4E-Cap. The experimental association rates for five cap analogues investigated in the present work are listed in Table 1. Details of these measurements are described elsewhere.16 It should be noted that the experimental rate constants are derived from measurements of the quenching of eIF4E fluorescence upon ligand binding and that it is not known at which proteincap distance this quenching becomes effective. Therefore, in the simulations described below several protein-ligand distances were considered as a “reaction distance” corresponding to the complex formation and were checked regarding consistency of the computed association rate constants with the experimental data. Because all the cap analogues contain the same domain responsible for the complex formation, it can be expected that one reaction distance for all analogues considered here should lead to agreement between computed and experimental association rate constants. 2.2. Structures Used in Calculations. Electrostatic and BD calculations were based on a crystal structure of truncated murine eIF4E (28-217) complexed with a dinucleotide cap analogue, m7GpppG, available from Protein Data Bank20 under accession code 1L8B.5 Monomer B of this structure was used. In our previous works12,13 another crystal structure of this protein was employed, that of eIF4E complexed with m7GDP, available under accession code 1EJ1.4 The coordinate file 1L8B does not contain coordinates for the second guanosine of the dinucleotide cap analogue, which is missing in the electron density because of high flexibility. Therefore, the coordinates present in the original PDB file represent just the complex of eIF4E with m7GTP. The lacking guanosine was model built-in using molecular dynamics program CHARMM.21 The resulting structure of the eIF4E protein complexed with the m7GpppG is presented in Figure 1. The protein is shown in a ribbon representation. Important proteins’ residues and the ligand are shown in the CPK representation. Both residues important for the ligand binding and directly involved in the fluorescence quenching, that is, Trp 56 and Trp 102, are shown in green, the Glu 103 is red, Arg 157 and Lys 159 are blue, and the ligand is colored according to the chemical elements. The models of the complexes with m7GMP and m7GDP were obtained simply by removing the appropriate number of phosphate groups from the original coordinate file. The model for the complex with dimethylated dinucleotide cap analogue, m7Gppp-m7G, was obtained by adding a methyl group at position N7 to the second guanosine in m7GpppG described above, using CHARMM. The above steps resulted in atomic structures for complexes of the eIF4E protein with five cap analogues investigated in the present work. The atomic models of the protein and the ligands were used to construct their electrostatic and hydrodynamic models employed in the Poisson-Boltzmann and BD simulations, as described below.
Figure 1. Complex of eIF4E protein with m7GpppG cap analogue. Drawing was done with the Chimera program.22 See text for further information.
2.3. Hydrodynamic Models of eIF4E and Ligands. The protein and ligand are modeled by sets of spherical beads and the appropriate diffusion tensors, which take into account hydrodynamic interactions23 between subunits of the given molecule, are calculated using methods described in details elsewhere.12,24-26 The average translational and rotational diffusion coefficients, to be used in the SDA simulations, were computed as one-third of the trace of the translational and rotational diffusion tensors, respectively. In the case of the eIF4E protein, each amino acid is substituted by one bead of radius 4.6 Å. Computed average translational diffusion coefficients of eIF4E is 9.38 × 10-7 cm2/s at 20 °C, which is equivalent to a sphere with the radius of 22.8 Å, which is close to the previous estimates, that is, 9.97 × 10-7 cm2/s at 20 °C and 21.5 Å,12 respectively, based on the structure available under accession code 1EJ1. In the case of cap analogues, hydrodynamic models were built by substituting each of the nonhydrogen atoms by a spherical bead. The common radius of the beads (3.1 Å) results from parametrization done on the basis of known values of translational diffusion constant of related molecules ATP and ADP,27,28 because appropriate data is lacking for the original ligands (see Table 2). The diffusion coefficients of the cap analogues resulting from these computations are shown in Table 3. Modeling of the hydrodynamic properties described above results in a reasonable estimation of the diffusion coefficients of isolated protein and ligand molecules. The relative translation diffusion coefficient required for simulation of molecular association is given by a sum of both diffusion coefficients, which means that hydrodynamic interactions (HI)29,30 between the protein and the ligand is neglected. It is difficult to estimate how large would be the decrease of the calculated rate constants for detailed hydrodynamic model of protein-ligand interactions because it depends on many details. Simulations31 for a dumbbell protein model and a dumbbell ligand model showed that computations neglecting the protein-ligand hydrodynamic
Binding Cap to eIF4E
J. Phys. Chem. B, Vol. 111, No. 45, 2007 13109
TABLE 2: Experimentally Determined Values of Translational Diffusion Constants, D, of ATP and ADP, and the Values Computed Assuming Hydrodynamic Radius of 3.1 Å for Each Heavy Atom of ATP or ADP Molecule molecule
temperature [K]
viscosity [cP]
Dexp [106 cm2s-1]
Dcalc [106 cm2s-1]
ATP ATP ATP ATP ATP ADP ADP
278 298 313 293 310 293 310
1.63 ( 0.01a 0.97 ( 0.01a 0.72 ( 0.01a 1.002b 0.6915b 1.002b 0.6915b
1.75 ( 0.09c 3.68 ( 0.14c 4.64 ( 0.13c 3.5 ( 0.5d 5.3 ( 0.6d 3.3 ( 0.4d 5.2 ( 0.5d
1.91 3.44 4.87 3.28 5.02 3.50 5.36
a Determined by authors of the original ref 27. b Values for water at appropriate temperature taken from Handbook of Chemistry and Physics, 60th edition; CRC PRESS, Inc.: Boca Raton, FL, 1979-1980; p F-51. c Taken from ref 27. d Taken from ref 28.
interactions for stick boundary conditions for the solvent velocity on the surfaces of moving particles lead to the rate constants, which are 27-32% too high in comparison to the simulations that include these interactions. In the present work, hydrodynamic interactions between the associating protein and its ligand are not taken into account. We are forced to do so because calculations with the hydrodynamic interactions between 190 subunits of the protein and more than 20 subunits required to model the ligand are too computationally demanding for the present day computers. One can obviously reduce the number of spherical elements describing associating molecules. However, reduction in the number of beads required to make inclusion of the HI between the associating molecules computationally feasible needs to be so drastic that the molecular details of the eIF4E protein would be completely lost. And as was shown previously, the complementarity of shape is also hydrodynamically important for the resulting association rate constants.31 2.4. Electrostatic Models of eIF4E and Ligands. Charge distributions in molecules can be described by the atomic partial charges. An additional very important factor is the existence of protonation equilibria of the titratable groups located on molecules that can result in the appearance of the (1e net charges at positions of location of these groups within the whole molecule. The net charges of the titratable groups depend on solution conditions, mainly the pH of the solution. The protonation probabilities of these groups, corresponding to experimental conditions of the stopped-flow measurements,16 were calculated as described in detail elsewhere34,35 using the University of Houston Brownian Dynamics program (UHBD).36,37 The titration curves were computed using the HYBRID program38 and a Monte Carlo (MC) approach described elsewhere.26 The MC program provides a list of a predefined number of the protonation states of the investigated molecules with the lowest free energies found during the MC search. The calculations were done for the apo-form of eIF4E protein, and for eIF4E complexes with the five cap analogues. In our previous study,13 probabilities of protonation in the case of the eIF4E-cap analogue complexes were estimated including two titratable groups on the ligands. The first group was the N1 position in the guanine ring, and the second position was the terminal phosphate group in the nucleotide analogue of the cap, that is, m7GMP, m7GDP, and m7GTP. From this study, it was concluded that in the complex the position N1 has a strong preference to be protonated at pH of 7.2 (i.e., the pH used in our experiments16), so the guanine derivative methylated at position N7 was assumed to have a fixed charge of +1e, whereas nonmethylated guanine moiety has a charge of 0. The most
probable protonation state of the terminal phosphate group of m7GMP and m7GDP at pH 7.2 is doubly ionized, whereas that for the m7GTP is singly ionized. Therefore in the present work, prediction of the protonation states of the eIF4E protein complexed with the cap analogues was done with the protonation states of the titratable groups of the cap analogues fixed at the most probable state predicted previously. Additionally for m7GTP, the doubly ionized state of its γ-phosphate group was also considered. It can be added that for isolated mononucleotide cap analogues in aqueous solutions, their terminal phosphate groups are doubly ionized at a pH of about 7. Therefore, only for m7GTP one can expect a change in the preferred protonation state of the terminal phosphate group upon formation of the complex. Details of the electrostatic properties of the ligands models used in this work are presented in Table 3. 2.5. Summary of Calculation of Association Rates by Brownian Dynamics Simulations. To calculate protein ligand association rates, one simulates many trajectories of the ligand molecule, initiated on the surface of a sphere of radius rb, the b-surface, around the center of coordinates of the protein. This sphere is made sufficiently large so that the electrostatic forces between the protein and the ligand are approximately centrosymmetric for r > rb. Each trajectory is continued until the ligand satisfies a predefined encounter criterion or reaches an outer spherical surface of radius rq, the quit-sphere. The bimolecular rate constant, kd, for diffusional encounter between a protein and its ligand can be calculated as36,37
kd ) kb
where, in general39,30
kb(q) ) 4π
βb
(3)
kb 1 - (1 - βb) kq
[
eU(r)/kT dr b(q) r2Dw
∫r∞
]
-1
(4)
βb is the fraction of trajectories that finish with encounters, U(r) is a centrosymmetric potential of interaction between molecules, and Dw is the relative translational diffusion coefficient independent of r, as we neglect HI between the associating molecules as discussed above. Simulations of the trajectories in the present work employed the BD method,40,41 implemented in the SDA program.18,19 This program was designed for simulation protein-protein association. It was adapted for simulation protein-nucleic acid subunits association simply by extending the list of residues and atoms recognized by the program. No other changes were required to use this program for eIF4E-cap systems. 2.6. Encounter Definition. There are three options to define the encounter criteria in the SDA program. The selection is done by setting a particular value of the icomm parameter in the SDA input file (1, 2, or 3). In the present work, intermediate option icomm ) 2 was used. For both options icomm ) 2 and icomm ) 3, a set of reactive contacts between two associating molecules is defined by tabulating all pairs of hydrogen bond donor on one encounter partner and hydrogen bond acceptor on the other, separated in the docked complexes and known, for example, from crystallographic studies, by no more than a predefined value (4.5 Å in the present work). During BD simulations, it is checked if the current interatomic distance for any of the predefined pairs of donor-acceptor is lower than the assumed reaction distance, and the number of pairs satisfying this requirement is counted. Usually the complex is considered
13110 J. Phys. Chem. B, Vol. 111, No. 45, 2007
Błachut-Okrasin´ska and Antosiewicz
TABLE 3: Hydrodynamic and Electrostatic Parameters of Models of Cap Analoguesa values and locations of the formal charges in SDA calc
avg diffusion coeff
cap analogue
C8
PA
PB
PC
C8
translational [106 cm2/s]
rotational [108/s]
m7GMP m7GDP m7GTP(a) m7GTP(b) m7GpppG m7Gpppm7G
+1 +1 +1 +1 +1 +1
-2 -1 -1 -1 -1 -1
-2 -1 -1 -1 -1
-1 -2 -1 -1
0 +1
3.72 3.41 3.29 3.29 2.77 2.77
4.38 3.62 3.39 3.39 1.84 1.83
a
m7GTP(a) and m7GTP(b) are models of m7GTP differing with respect to the assumed ionization state of its γ-phosphate group.
TABLE 4: Intrinsic and Apparent Ionization Constants and Charges of the Indicated Titratable Residues at Four Protonation Patterns of eIF4E Protein with the Lowest Free Energies Found during MC Simulations, Predicted at pH of 7.2a eiF4E-m7GDP
apo eIF4E
eiF4E-m7GpppG
eIF4E residue
pKi/pKa
1234
pKi/pKa
1234
pKi/pKa
1234
His 33 His 37 His 78 Glu 103 Glu 140 Arg 157 Lys 159 Lys 162 His 178 Lys 184 His 200 Lys 206
6.0/5.9 5.0/7.3 6.1/6.1 5.3/4.9 11/9.4 10/15 5.9/4.2 9.1/8.0 4.1/4.6 1.5/4.9 1.7/-3.1 9.0/8.7
0000 +000 0000 ---0000 ++++ 0000 +++0 0000 0000 0000 ++0+
6.0/5.9 5.0/7.3 6.1/6.2 -0.6/-1.8 11/9.5 17/26 6.9/4.8 14/14 4.2/4.6 1.7/5.2 1.2/-5.2 9.7/9.1
0000 +0+0 0000 ---0000 ++++ 0000 ++++ 0000 0000 0000 ++++
6.0/5.9 5.1/7.3 6.1/6.2 -0.7/-1.9 11/9.5 17/26 7.5/5.2 15/15 4.2/4.7 1.7/5.2 1.6/-4.9 8.1/7.5
0000 +0+0 0000 ---0000 ++++ 0000 ++++ 0000 0000 0000 ++00
a As examples, the data are shown for apo-eIF4E and its complexes with m7GDP and m7GpppG. The pK s and pK s are rounded to two significant i a digits.
to be formed if there are 2 or 3 such pairs found. The difference between icomm ) 2 and icomm ) 3 is that for the latter the contact pairs taken into account must be “independent”, which means separated by more than some predefined distance (usually 5 to 6 Å). Several values of the reaction distances were considered, ranging from 8 to 3.5 Å with step 0.25-0.10 Å, because of uncertainty in the eIF4E-cap distance required for the effective quenching of the protein’s tryptophan fluorescence used in experimental studies of the kinetics of the eIF4E-cap binding. For each of these reaction distances, satisfaction of the reaction criterion simultaneously by 2 donor-acceptor pairs was considered as the complex formation. As it was mentioned above, we seek a reaction distance for which theoretical association rate constants of all five cap analogues agree with their experimental counterparts. 2.7. Modeling of Forces. BD simulations of molecular associations usually assume that the particles move under influence of electrostatic intermolecular forces and the random forces mimicking influence of the bombardment by the solvent molecules. The latter forces are modeled by generating random numbers from Gaussian distribution and by using estimated diffusion coefficients of the molecules. Additionally, the SDA program includes electrostatic desolvation forces acting at small separations that reflect the energetically unfavorable exclusion of high-dielectric solvent that accompanies the approach of a low-dielectric protein. Finally, short-range repulsive forces are modeled by prohibiting overlap of the exclusion volumes of the associating molecules. The excluded volume of each molecule is determined by the exclusion grid with grid points assigned a value of 1 if they are within the probe accessible surface of the molecule and a value of 0 otherwise. Details of computations of the electrostatic, desolvation, and excluded volume effects are given in the original references for the SDA program.18,19,42 Figure 2 shows constant electrostatic potential surfaces of eIF4E protein (with the electric charge distribution correspond-
ing to the protonation state of the lowest free energy determined by the MC simulations) and of m7GTP, superimposed on each other to correspond to the mutual placement of these molecules in the complex, with red (protein) and light red (ligand) representing a potential of -0.5kT (0.29 kcal/mol(e)), and blue (protein) and light blue (ligand) representing potential of +0.5kT. It can be seen that the potential maps are quite complicated. In the complex, positive potential of the ligand faces both negative and positive potential regions around the protein, and the same is true for the negative potential of the ligand. 2.8. Technical Details of the Simulations. All simulations were performed at 293 K (RT ) 0.58 kcal/mol) at ionic strengths equivalent to 150 mM with a solvent dielectric constant of 78 and that for the protein of 4. The protein and cap analogues interiors were defined as the regions within the van der Waals surface. Dielectric boundary smoothing43 was used in all finite difference calculations. Electrostatic potentials around the protein and cap analogues were generated on a 1953 cubic grid with a 1.0 Å spacing. For each case, that is, apo-eIF4E and its complex with cap analogues, the potentials were generated for four protonation states with the lowest free energy as obtained from our MC simulations. All atomic partial charges (the partial charges for titratable groups, for their two protonation states, and the partial charges of nontitratable groups) and radii for the protein and cap analogues were taken from the CHARMM27 parameter set for the standard amino acids and nucleic acids.32,33 The effective charges were placed on the carboxylate oxygen atoms of Asp and Glu residues and the C terminus, the amine nitrogen atoms of Lys and Arg residues and the N terminus for proteins, and on atoms listed in Table 3 for the cap analogues. The grid used in determination of desolvation energies was 1503 with 1 Å spacing, and 1.5 Å was assumed as the mobile ion radius. Electrostatic desolvation terms were computed using
Binding Cap to eIF4E
Figure 2. Surfaces of constant electrostatic potential of eIF4E protein (blue +0.5kT; red -0.5kT) and 7mGTP (light blue +0.5kT; pink -0.5kT) molecules forming the complex for solution conditions pH 7.2, ionic strength of 150 mM, and temperature 293 K. Drawing was done with the Chimera program.22
a scaling factor R ) 1.62 as employed for simulations of protein-protein associations with van der Waals definition of the molecular interior in the example of barnase and barstar, provided with the SDA package. To model exclusion forces, 0.5 Å spacing grids were attached to eIF4E molecule, and the probe radius of 1.5 Å was used. BD simulations were done using the same grid as for the electrostatic potentials calculations. The b sphere radius was 60 Å, and the q radius was 300 Å. A variable time step was used. It was equal to 0.5 ps for distances for those where the molecules are close to overlap, which is up to 50 Å, Then, it was increased linearly where at the distance of 90 Å it was 20 ps, and then it decreased smoothly to 0 when the ligand molecule was crossing the q-surface. Diffusional encounter rate constants reported in the present work are based on 40 × 2000 Brownian trajectories for each pair of eIF4E protein-cap analogue. Each 2000 trajectories simulation resulted in an association rate constant. The final rate constant is an average of the 40 determinations, and the indicated errors are standard deviations of the averages. Because of the symmetry of m7Gpppm7G, the rate constants computed for this ligand were multiplied by a factor of 2 and the statistical error was that of the original calculations multiplied by a factor of x2. 3. Results and Discussion 3.1. Protonation States of eIF4E Protein. Table 4 shows three examples of protonation states of selected residues of eIF4E protein at 293 K, pH 7.2, and ionic strength of 150 mM obtained from our titration procedure for the four lowest freeenergy states. The examples obey the apo-eIF4E and its complexes with m7-GDP and m7GpppG. The selected residues include those important for the binding, as suggested by atomic structures of eIF4E-cap complexes (i.e., Glu 103, Arg 157, and Lys 162), all His residues, and all these residues for which protonation states are different than usually expected for pH of 7.2. All other eIF4E titratable residues have their most probable protonation states as expected, therefore it was not necessary to include them in Table 4.
J. Phys. Chem. B, Vol. 111, No. 45, 2007 13111 For calculations including a complexed cap analogue, the ligand was treated as part of the low dielectric region with charges fixed as established in our previous study.13 The N1 of the methylated guanine was considered permanently protonated; hence, methylated guanine was assumed to have net charge of +1e, and the phosphate groups had partial charge distributions corresponding to the net charges listed in Table 4. The pKi refers to the predicted pKa value when all other eIF4E residues are neutral; therefore, if this is the pKi for a eIF4E-cap complex it includes the influence of the fixed charge distribution of the ligand. The very high or very low (even negative) pKa values should be understood as reflecting correspondingly large freeenergy contribution through the standard thermodynamic equation K ) exp(-∆G/RT), rather than the values that can be determined experimentally, as they might be located outside of the pH range accessible for experiments. First consider the Glu 103 residue. Its pKi and pKa in the apo form differ not very much (0.4 pH unit) and the predicted protonation state for the all four lowest energy states is to be deprotonated. When a cap analogue is attached, both the pKis and pKas are shifted significantly down in comparison to the corresponding apo case. Also, the difference between the pKi and pKa is larger (-1.2 pH unit in comparison to -0.4 pH unit) indicating in the presence of the ligand, interactions between titratable residues of the eIF4E further stabilize the deprotonated state of Glu 103 in addition to stabilization due to the net positive charge of the methylated guanine. Regarding the two basic residues in close contact with phosphate groups of cap analogues (Arg 157 and Lys 162), we can see that presence of the ligand causes strong increase in the pKi values (7 pH units) and even stronger increase in the pKa values (11 pH units). Therefore, we can see here the direct effect of the negative charges of the phosphate groups of the ligands, and again we can see increased interactions between titratable residues of the eIF4E, stabilizing protonated states of Arg 157 and Lys 162. The neighboring Lys 159 residue has lower pKi and pKa than Arg 157 and Lys 162 and remains deprotonated in all 4 states. Although binding of a cap analogue increases both the pKi and pKa of Lys 159, apparently electrostatic interactions with charged Arg 157 and Lys 162 prevent any substantial increase in the pKa Lys 159 and it remains deprotonated at a pH of 7.2. His residues are included in Table 4 because their expected protonation states at pH of 7.2 cannot be so easily guessed as it is the case for Asp, Glu, and C-terminus on the one hand, and Arg, Lys, and N-terminus on the other. Looking further into Table 4 it is interesting to see that Glu 140 and Lys 184 residues, which in general can be expected to be ionized in the proteins at pH of about 7, are predicted to be neutral in the case of eIF4E protein. We can see that this is because of the pKi values, which are high in the case of Glu 140 and low in the case of Lys 184. Therefore, these are desolvation effects that make residues 140 and 184 neutral. The interactions with other titratable groups in the eIF4E increase the probability of ionization of these residues but not sufficiently to have them ionized in any of the four lowest energy states of the protein. Another interesting feature seen from Table 4 is that influence of cap analogues is to shift the pKas in such directions that the existing ionization states of apo eIF4E are in many cases more stable in the complexed form. For example, the pKa of Glu 103 is shifted to lower values and those of Arg 157 and Lys 162 are shifted to higher values upon a cap binding. Simultaneously, the ionization states, predicted to exist at pH of 7.2 for all four lowest energy states, are not changed. For His 37, Lys 162, and Lys 206 we see a different behavior: their predicted charge
13112 J. Phys. Chem. B, Vol. 111, No. 45, 2007
Błachut-Okrasin´ska and Antosiewicz
TABLE 5: Diffusional Encounter Rate Constants from BD Simulations Obtained Using icomm=2 Option and Assuming a Two-Contact Reaction Criterion, as Functions of the Reaction Distance (RD) Defined in Section 2.6a cap analogue RD [Å] m7GMP 4.1
4.2
4.3
4.4
4.5
48 ( 3 48 ( 4 43 ( 4 27 ( 3 59 ( 4 62 ( 4 58 ( 5 39 ( 3 77 ( 4 78 ( 5 73 ( 5 47 ( 4 95 ( 5 91 ( 6 93 ( 6 57 ( 5 119 ( 6 107 ( 6 113 ( 7 70 ( 5
m7GTP
m7GDP
q ) -2
m7GTP q ) -3
197 ( 9 196 ( 9 180 ( 8 193 ( 8 251 ( 10 238 ( 10 226 ( 11 235 ( 9 302 ( 11 282 ( 12 269 ( 13 281 ( 10 353 ( 12 334 ( 12 312 ( 13 327 ( 12 399 ( 13 379 ( 14 366 ( 14 380 ( 14
464 ( 15 442 ( 17 428 ( 13 447 ( 14 533 ( 16 516 ( 18 491 ( 15 517 ( 16 614 ( 18 587 ( 19 556 ( 17 588 ( 18 683 ( 19 656 ( 21 628 ( 19 669 ( 19 753 ( 21 734 ( 22 709 ( 20 751 ( 20
917 ( 25 902 ( 27 905 ( 28 870 ( 34 1035 ( 28 1014 ( 29 1023 ( 30 988 ( 36 1154 ( 29 1142 ( 31 1137 ( 33 1095 ( 38 1267 ( 30 1260 ( 35 1255 ( 34 1202 ( 41 1387 ( 33 1370 ( 37 1371 ( 35 1319 ( 44
m7GpppG m7Gpppm7G 151 ( 7 143 ( 7 60 ( 4 140 ( 6 186 ( 7 181 ( 9 82 ( 5 176 ( 7 228 ( 8 219 ( 9 102 ( 5 215 ( 8 270 ( 9 261 ( 10 120 ( 6 258 ( 9 305 ( 11 302 ( 10 143 ( 7 308 ( 9
64 ( 4 71 ( 4 52 ( 4 69 ( 5 84 ( 5 94 ( 6 72 ( 5 89 ( 5 112 ( 5 113 ( 7 89 ( 5 120 ( 7 132 ( 6 139 ( 8 112 ( 6 144 ( 6 158 ( 8 168 ( 8 134 ( 6 167 ( 7
The rate constants are given in units of µM-1s-1. The indicated errors correspond to standard deviations of averages of 40 results each coming from 2000 BD trajectories. a
at pH of 7.2 in the third and fourth protonation states of eIF4E is changed on going from the apo to the complexed form of the protein. The fact that the residues crucial for the ligand binding do not change the preferred protonation state upon complex formation might be of biological significance. As is known, changes in protonation patterns of active site residues in enzyme-substrate association might be necessary to facilitate the catalyzed reaction. The eIF4E is a nonenzymatic protein, and it is the stability of the electrostatic field near the binding site that may be useful. 3.2. Diffusional Encounter Rate Constants from Brownian Dynamics Simulations. The main purpose of the present study was to check if relatively the simple BD algorithm implemented in the SDA program is able to reproduce the experimentally determined bimolecular association rate constant for the five cap analogues described above. Previously, we have tried to do this job using the UHBD program with ligands modeled by a small number of larger beads (up to 9 bead models for the dinucleotide cap analogues with internal flexibility and HI between the beads), but we could not obtain such good correlation with the experimental results as we achieved using the SDA program. The experimental rate constants for both onestep and two-step binding mechanism are listed in Table 1. The agreement between these experimental rate constants and the results obtained from our simulations would mean that for one particular reaction distance (the same for all five cap analogues) we reproduce the experimental values. The expectation that this agreement should occur for the same reaction criterion is based on the fact that all the ligands have a common structural element, and that this element is crucial for the complex formation. The diffusional association rate constants obtained from our BD simulations are presented in Table 5 as functions of the assumed encounter protein-ligand distance. The presented values were obtained assuming a two-contact reaction criterion at an indicated distance from simulations including effective charges and electrostatic desolvation contributions. Because of space limits, only a narrow range of the assumed reaction
Figure 3. Comparison of the computed diffusional encounter rate constants as functions of the reaction distance with and without desolvation penalty. The rate constants are computed assuming a twocontact reaction criterion from simulations including effective charges.
distances is presented in Table 5. Actually, we monitored a bit wider range from 3.5 to 8.0 Å, but the range shown includes the value of 4.3 Å for which the best agreement between experimental and theoretical association rate constants was observed. It can be noticed that the rate constants obtained for 4.3 Å reaction distance constitute from 65 to 83% of the values obtained for 4.5 Å reaction distance. As it was written above, the HI between associating molecules lead to about 30 ( 10% decrease in the association rate; therefore, we could conclude that having the HI taken into account in our simulation, we would observe the best agreement between experimental and theoretical values of the rate constants for those computed for the reaction distance of 4.5 Å, which is very close to the value obtained when the HI are omitted. Before we discuss a comparison between experimental and theoretical rate constants, let us discuss some other issues. First, it can be noticed from Table 5 that the diffusional association rate constants for m7GTP with assumed double ionization of the γ-phosphate group are substantially larger than those obtained with assumed mono-ionized γ-phosphate group, and as it is discussed in the next section the value obtained for the mono-ionized γ-phosphate group of m7GTP corresponds to the experimental value. In all discussions below, we refer to the results for m7GTP mono-ionized at the γ-phosphate group until it is clearly stated that the second case is considered. Figure 3 presents a comparison of the association rate constants for the five cap analogues as functions of the assumed reaction distance obtained with and without desolvation penalty included in the simulations. Significant changes in the observed picture can be noticed. This indicates the importance of the desolvation penalty contribution for the predicted rate constants. Many sets of simulations without desolvation penalty contributions were done with a variation of the other simulation parameters. None of these resulted in the rate constants resembling experimental data. Similar observations regarding importance of the electrostatic desolvation contributions in protein-protein association were made previously by Elcock and co-workers.44 At pH of 7.2, the value at which our measurements were done, the predicted average charge of the protein is positive, ranging from 1.2e to 1.7e, for different cap analogues, whereas the net charges of the ligands vary from -1e for m7GMP and m7Gppp-m7G, -2e for m7GDP, and m7GpppG. The net charge of m7GTP, free in solution at pH 7.2, is close to -3e, but in
Binding Cap to eIF4E
J. Phys. Chem. B, Vol. 111, No. 45, 2007 13113
Figure 4. Comparison of the diffusional encounter rate constants as functions of the reaction distance with and without inclusion of the electrostatic interaction in the simulations. The rate constants are computed assuming a two-contact reaction criterion.
Figure 5. Computed diffusional association rate constants as functions of the reaction distance for the four lowest energy protonation patterns of eIF4E protein. Error bars are shown only for the lowest energy state. Error bars for the remaining states are similar, and they are now shown for clarity of the picture.
TABLE 6: Influence of the Electrostatic Steering on the Computed Rate Constants Expressed as the Ratio of the Association Rate Constants Computed with and without, Respectively, Electrostatic Interactions Includeda assumed reaction distance cap analogue m7
GMP m7GDP m7GTP(a) m7GTP(b) m7GpppG m7Gpppm7G
4Å
5Å
6Å
0.94 4.79 13.3 26.6 4.43 0.92
1.30 3.93 6.53 11.2 3.73 1.04
1.31 2.78 4.26 6.63 3.08 1.26
a m7GTP(a) and m7GTP(b) are models of m7GTP differing with respect to the assumed ionization state of its γ-phosphate group.
the complex with the eIF4E it is -2e, as discussed elsewhere.13 Therefore, for separated eIF4E and cap analogues there is some electrostatic attraction that is expected to steer the ligand toward the protein. It is thus interesting to evaluate the strength of this steering. Figure 4 compares results of simulations obtained for simulations with and without electrostatic forces included. It is not easy to make general conclusions regarding the significance of the electrostatic attraction between eIF4E and the ligands because for some of them the effects look different than for the others. To discuss this problem in more detail, we present the values of these ratios for selected assumed reaction distances: 4, 5, and 6 Å. They are presented in Table 6. These results are interesting indeed. We can see that the ratios strongly depend on the ligand. For m7GMP and m7Gppp-m7G, we observe that the electrostatic steering is similar for the separation of 6 Å, but it is relatively small at only about 1.3. For the distance of 5 Å, the ratios for these two ligands are different: for m7GMP it is again 1.3 but for m7Gppp-m7G it is equal 1. For a still shorter distance of 4 Å, both ratios are again close to each other but what is interesting is that they are below 1. These relations are qualitatively different for the remaining cap analogues, which can be considered as belonging to the second group. For these analogues we observe a gradual decrease in the ratio of the rate constants when the reaction distance is increased, and therefore the electrostatic steering is clearly visible and pretty strong. Moreover we can see that both factors, the net charge and its distribution, are important. Figure 5 shows the computed association rates for four lowest energy state protonation patterns of the eIF4E protein as functions of the reaction distance, for all five ligands. For the
Figure 6. Correlations between experimental and computed association rate constants for eIF4E protein and five analogues of cap. Experimental association rate constants are from one-step (A) and two-step (B) binding models. Straight lines represent results of fitting ksim ) (a)kexp, and r is the correlation coefficient. The point located clearly out of the correlation line in each part of the picture represents the rate constant of m7GTP computed with assumption that its γ-phosphate group is doubly ionized.
rates computed with the eIF4E protonation pattern corresponding to the lowest free energy statistical errors are included, for the remaining they are very similar, but they are omitted for clarity of picture. It can be seen that the rate constants for all four lowest energy protonation states of eIF4E are very similar. Only for m7GpppG, the binding rate constants for the third protonation state of the protein differ clearly from those for the three remaining states. 3.3. Comparison with Experimentally Determined eIF4ECap Bimolecular Association Rate Constants. Figure 6 presents a comparison of the theoretical association rate constants computed for the assumed reaction distance of 4.3 Å with both sets of the experimental association rate constants: one derived assuming a one-step binding and the other assuming a two-step binding. It can be seen that in both cases the correlation between the association rates is very good, close to a quantitative one. The correlation coefficients for both data sets are almost the same (0.99 for one-step model and 0.98 for two-step model), whereas the slope of the correlation line for two-step model is slightly closer to the value of 1 expected for ideal agreement than for the one-step binding model. The
13114 J. Phys. Chem. B, Vol. 111, No. 45, 2007 difference is not large: 0.99 versus 0.93 but nevertheless we can conclude that the two-step binding model is better. An additional interesting thing visible in Figure 6 is that for the result for m7GTP with its γ-phosphate group doubly ionized, the computed rate is clearly not in agreement with the experimental data. This supports the conclusions of our earlier study on pH dependence of association rate constants by stopped-flow technique13 that m7GTP has mono-ionized γ-phosphate group when bound to eIF4E protein. In the simulations, all potential hydrogen bond donors and acceptors on the cap models were included in the considered list of pair contacts with the protein polar atoms. As the cap analogues contain the same common part, that is, the one constituted by the m7GMP itself, one can ask if the potential contacts should be limited to those atoms on the ligands that are present in all analogues. Otherwise, there is some kind of entropic increase in the chances of making contacts when larger ligands have more atoms that can be included on the list. The effect of the number of possible contacts is reflected in the results of simulations. However it appeared that a good correlation with the experimental results was obtained by us when all polar atoms present in the cap analogues were included in the list of potential interatomic contacts between the eIF4E protein and its ligands. Obviously, it cannot be excluded that it is possible to assume such values of the parameters used in the SDA program that would lead to a good agreement with the experimental data for assumption of the possible protein-ligand contacts limited to those present in the m7GMP. However, because the number of possibilities is large one cannot predict how long it would take to find such set of the parameters. We have tried many of them but without success. 3.4. Concluding Remarks. We have shown in this work that it is possible to reproduce association rate constants of eIF4E cap system derived from experiments by BD simulations. The work presented here can be extended in many ways, but we have left them for some future studies as the presented material is already quite lengthy. The SDA program used in our simulations allows one to perform BD simulations and association rate determination in several other ways than that described in this work. Therefore, it might be useful to investigate influence of different modifications in the assumed SDA protocol and values of the parameters in its input file on the results. For example, one can change the dielectric boundary model from van der Waals to a Richards probe-accessible surface;45,46 however initial calculations with this option did not result in diffusional rate constants in agreement with the experimental results. The other possibilities are to use a different location of the effective charges on cap analogues, to take into account more protonation states of the protein, or to prepare electrostatic models of the protein based on the average protonations for given solution conditions, which would correspond to the situation that protonation equilibria are much faster in establishing than ligand diffusion in the vicinity of the binding site. Such approach was used by us in a previous BD study, however, describing only eIF4E-m7GpppG system and employing the UHBD program.12 It is also interesting to see how important for the predicted rate constants are the unusual protonation states of Glu 140 and Lys 184, that is, if assuming their ionized states in the simulations would significantly change these constants. One substantial extension would be the inclusion of protein and cap analogues flexibility but this requires reduction of structural elements just to make simulations feasible. Some insight into protein flexibility can be obtained by making
Błachut-Okrasin´ska and Antosiewicz simulations with true apo form of eIF4E, which was determined recently for the human eIF4E protein by NMR methods.8 Structural data obtained from X-ray crystallography for human eIF4E protein complexed with m7GpppA are also available,6 and therefore simulations with the complexed form of the same protein are possible. For completeness, one needs stopped-flow kinetic transient analysis by numerical methods like those implemented in the DynaFit program, which according to our best knowledge are not yet available. Nomenclature G ) guanosine. m7G ) 7-methyl-guanosine. ppp ) 5′-5′-tri-phospha-te bridge. GMP (GDP, GTP) ) guanosine 5′-mono-(di- tri-)-phosphate. ATP ) adenosine-tri-phosphate. eIF4E ) eukaryotic initiation factor 4E. Acknowledgment. We thank Dr. Razif R. Gabdoulline and Professor Rebecca C. Wade for their help with modifications of the SDA program and fruitful discussions. Partial support of the CoE BioExploratorium Project WKP_1/1.4.3/1/2004/44/44/ 115/2005 is acknowledged. All plots were done using the xmgrace program (Evgeny Stambulchik, Rehovot, Israel). References and Notes (1) Pestova, T. V.; Kolupaeva, V. G.; Lomakin, I. B.; Pilipenko, E. V.; Shatsky, I. N.; Agol, V. I.; Hellen, C. U. Proc. Natl. Acad. Sci. U.S.A. 2001, 98, 7029-7036. (2) Sonenberg, N. In Translation Control; Hershey, J., Mathews, M. B., Sonenberg, N., Eds.; Cold Spring Harbor Laboratory Press: New York, 1996; p 245-269. (3) Shatkin, A. J. Cell 1976, 9, 645-653. (4) Marcotrigiano, J.; Gingras, A.; Sonenberg, N.; Burley, S. K. Cell 1997, 89, 951-961. (5) Niedz´wiecka, A.; Marcotrigiano, J.; Ste¸pinski, J.; JankowskaAnyszka, M.; Wysłouch-Cieszyn´ska, A.; Dadlez, M.; Gingras, A.; Mak, P.; Darzynkiewicz, E.; Sonenberg, N.; Burley, S. K.; Stolarski, R. J. Mol. Biol. 2002, 319, 615-635. (6) Tomoo, K.; Shen, X.; Okabe, K.; Nozoe, Y.; Fukuhara, S.; Morino, S.; Sasahi, M.; Taniguchi, T.; Miygawa, H.; Kitamura, K.; Miura, K.; Ishida, T. J. Mol. Biol. 2003, 328, 365-383. (7) Matsuo, H.; Li, H.; McGuire, A. M.; Fletcher, C. M.; Gingras, A.; Sonenberg, N.; Wagner, G. Nature Struct. Biol. 1997, 4, 717-724. (8) Volpon, L.; Osborne, M. J.; Topisirovic, I.; Siddiqui, N.; Borden, K. L. B. EMBO J. 2006, 25, 5138-5149. (9) Carberry, S. E.; Rhoads, R. E.; Goss, D. J. Biochemistry 1989, 28, 8078-8082. (10) Cai, A.; Jankowska-Anyszka, M.; Centers, A.; Chlebicka, L.; Stepinski, J.; Stolarski, R.; Darzynkiewicz, E.; Rhoads, R. E. Biochemistry 1999, 38, 8538-8547. (11) Sha, M.; Wang, Y.; Xiang, T.; van Heerdens, A.; Browning, K. S.; Goss, D. J. J. Biol. Chem. 1995, 270, 29904-29909. (12) Błachut-Okrasin´ska, E.; Bojarska, E.; Niedz´wiecka, A.; Chlebicka, L.; Darzynkiewicz, E.; Stolarski, R.; Ste¸pinski, J.; Antosiewicz, J. M. Eur. Biophys. J. 2000, 29, 487-498. (13) Długosz, M.; Błachut-Okrasin´ska, E.; Bojarska, E.; Darzynkiewicz, E.; Antosiewicz, J. M. Eur. Biophys. J. 2003, 31, 608-616. (14) Khan, M. A.; Goss, D. J. Biochemistry 2005, 44, 4510-4516. (15) Slepenkov, S. V.; Darzynkiewicz, E.; Rhoads, R. E. J. Biol. Chem. 2006, 281, 14927-14938. (16) Błachut-Okrasin´ska E.; Bojarska E.; Ste¸pinski, J.; Antosiewicz J. M. Biophys. Chem. 2007, 129, 289-297. (17) Kuzmic, P. Anal. Biochem. 1996, 237, 260-273. (18) Gabdoulline R. R.; Wade, R. C. Biophys. J. 1997, 72, 1917-1929. (19) Gabdoulline R. R.; Wade, R. C. Biophys. J. 2001, 306, 11391155. (20) Bernstein, F. C.; Koettzle, T. F.; Williams, G. J. B.; Meyer, E. F.; Brice, M. D.; Rodgers, J. R.; Kennard, O.; Shimanouchi, T.; Tasumi, M. J. J. Mol. Biol. 1977, 123, 557-594. (21) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187-217. (22) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. J. Comput. Chem. 2004, 25, 1605-1612.
Binding Cap to eIF4E (23) Yamakawa, H. J. Chem. Phys. 1970, 53, 436-443. (24) Garcia de la Torre, J.; Bloomfield, V. A. Q. ReV. Biophys. 1981, 14, 81-139. (25) Antosiewicz, J.; Porschke, D. J. Phys. Chem. 1989 93, 5301-5305. (26) Antosiewicz, J. Biophys. J. 1995, 69, 1344-1354. (27) Hubley, M. J.; Locke, B. R.; Moerland, T. S. Biochim. Biophys. Acta 1996, 1291, 115-121. (28) de Graaf, R. A.; van Kranenburg, A.; Nicolay, K. Biophys. J. 2000, 78, 1657-1664. (29) Deutch, J. M.; Felderhof, B. U. J. Chem. Phys. 1973, 59, 16691671. (30) Friedman, H. L. J. Phys. Chem. 1966, 70, 3931-3933. (31) Antosiewicz, J.; McCammon, J. A. Biophys. J. 1995, 69, 57-65. (32) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E., III; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586-3616. (33) Foloppe, N.; MacKerell, A. D., Jr. J. Comput. Chem. 2000, 21, 86-104. (34) Antosiewicz, J.; Briggs, J. M.; Elcock, A. E.; Gilson, M. K.; McCammon, J. A. J. Comput. Chem. 1996, 17, 1633-1644.
J. Phys. Chem. B, Vol. 111, No. 45, 2007 13115 (35) Briggs, J. M.; Antosiewicz, J. ReV. Comput. Chem. 1999, 13, 249311. (36) Davis, M. E.; Madura, J. D.; Luty, B. A.; McCammon, J. A. Comput. Phys. Commun. 1991, 62, 187-197. (37) Madura, J. D.; Briggs, J. M.; Wade, R. C.; Davis, M. E.; Luty, B. A.; Ilin, A.; Antosiewicz, J.; Gilson, M. K.; Bagheri, B.; Scott, L. R.; McCammon, J. A. Comput. Phys. Commun. 1995, 91, 57-95. (38) Gilson, M. K. Proteins: Struct., Funct., Genet. 1993, 15, 266282. (39) Kramers, H. A. Physica 1940, 7, 284-304. (40) Ermak, D. L.; McCammon, J. A. J. Chem. Phys. 1978, 69, 13521360. (41) Northrup, S. H.; Allison, S. A.; McCammon, J. A. J. Chem. Phys. 1984, 80, 1517-1524. (42) Gabdoulline R. R.; Wade, R. C. J. Phys. Chem. 1996, 100, 38683878. (43) Davis, M. E.; McCammon, J. A. J. Comput. Chem. 1990, 11, 401409. (44) Elcock, A. H.; Gabdoulline, R. R.; Wade, R. C.; McCammon, J. A. J. Mol. Biol. 1999, 291, 149-162. (45) Richards, F. M. Ann. ReV. Biophys. Bioeng. 1977, 6, 151-176. (46) Gilson, M. K.; Sharp, K. A.; Honig, B. H. J. Comput. Chem. 1988, 9, 327-335.
