Predicting Mutation-Induced Stark Shifts in the Active Site of a

Article pubs.acs.org/JPCA

Predicting Mutation-Induced Stark Shifts in the Active Site of a Protein with a Polarized Force Field Xianwei Wang,† Xiao He,*,† and John Z. H. Zhang*,†,‡ †

State Key Laboratory of Precision Spectroscopy and Department of Physics, Institute of Theoretical and Computational Science, East China Normal University, Shanghai 200062, China ‡ Department of Chemistry, New York University, New York, New York 10003, United States S Supporting Information *

ABSTRACT: The electric field inside a protein has a significant effect on the protein structure, function, and dynamics. Recent experimental developments have offered a direct approach to measure the electric field by utilizing a nitrile-containing inhibitor as a probe that can deliver a unique vibration to the specific site of interest in the protein. The observed frequency shift of the nitrile stretching vibration exhibits a linear dependence on the electric field at the nitrile site, thus providing a direct measurement of the relative electric field. In the present work, molecular dynamics simulations were carried out to compute the electric field shift in human aldose reductase (hALR2) using a polarized protein-specific charge (PPC) model derived from fragment-based quantum-chemistry calculations in implicit solvent. Calculated changes of electric field in the active site of hALR2 between the wild type and mutants were directly compared with measured vibrational frequency shifts (Stark shifts). Our study demonstrates that the Stark shifts calculated using the PPC model are in much better agreement with the experimental data than widely used nonpolarizable force fields, indicating that the electronic polarization effect is important for the accurate prediction of changes in the electric field inside proteins.

1. INTRODUCTION Electrostatic interactions play a significant role in protein structures and functions.1−5 Charged and polar groups in folded proteins produce large electric fields that drive processes such as protein−ligand binding,6 protein−protein interactions,7 electron transfer,8 proton binding and release,9 and enzyme catalysis.10 Recent studies on a number of enzymes have shown that electric fields in the active sites of proteins make a substantial contribution to the stabilization of transition states.10−12 The electric fields inside a protein differ between different states of the protein, for instance, two conformations, two protonation states, or two mutants.13 Different mutations in a protein can directly impact the rates of enzyme-catalyzed reactions and thus could provide an effective method for the design of more efficient enzymes. The accurate description of electric fields in proteins and their changes due to environmental alterations such as mutations can be used to help analyze structures and functions and make quantitative predictions of certain properties of proteins. © XXXX American Chemical Society

To provide a proper description of the electrostatic environment inside the protein, many experimental and theoretical works have been carried out over the past decade.13−18 Pioneering works by Boxer and co-workers utilized vibrational Stark effect (VSE) spectroscopy19,20 to measure the electric field shifts in the active site of the enzyme human aldose reductase (hALR2) using site-directed mutations.13,17,21 In such experiments, the nitrile group is often used as an ideal vibrational probe to monitor the change of electric fields.22−26 Many other experimental and theoretical works have been extended to quantify changes in the frequency or the electric field between different solvent environments for various nitrile probes.18,25,27−31 Special Issue: Prof. John C. Wright Festschrift Received: December 7, 2012 Revised: March 21, 2013

A

dx.doi.org/10.1021/jp312063h | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

hc = −0.69 cm−1/(MV/cm).13 Thus, the observed nitrile stretching frequency shift, ΔṽCN, between the mutated hALR2 and the wide-type (WT) protein (ṽmutant − ṽWT) and the change in the electrostatic field, ΔF⃗∥, along the nitrile axis between the mutant and wide type (F∥,mutant − F∥,WT) have the linear dependence13

The observed vibrational frequency shift of the nitrile group, ΔvC⃗ N, is proportional to the change in the electrostatic field, ΔF⃗protein, at the probe location between two different states of a protein13,17,21−23,26 ⃗ hc ΔvC̃ ≡ N = −Δμprobe ⃗ ·ΔFprotein

(1)

where Δμ⃗ probe is the Stark tuning rate, representing the change in dipole moment between the ground state and excited states of the vibrational transition for the probe; h is Planck’s constant; and c is the speed of light. As verified by previous works,32−34 the linear relationship in eq 1 is valid when the vibrational probe is not involved in hydrogen-bonding interactions. However, in case of direct hydrogen bonding, the response of the vibrational frequency shift of the probe to the change of the electric field may become nonlinear.35,36 In a previous work by Suydam et al.,13 the nitrile-containing inhibitor IDD743 (see Figure 1A) was utilized as a probe to

ΔvC̃ ≡ N = [0.69 cm−1/(MV/cm)]·ΔF

(2)

The measured vibrational Stark shift quantifies the magnitude of the electrostatic field change in the active site of hALR2, which can be directly compared with electrostatics calculations to validate the accuracy of different theoretical methodologies. Suydam et al. concluded that molecular dynamics (MD) simulations and ensemble averaging are needed to accurately predict the change in electric fields.13 Molecular mechanics (MM) force fields, such as CHARMM, Amber, and OPLS, can be utilized to calculate the electric fields at a given site in a protein by using partial atomic charges. Despite much success in the application of standard force fields, there is a fundamental limitation due to the lack of electronic polarization effects. The widely used charge models are mean-field-like and do not include polarization effect due to the specific protein environment,15,16 which often makes MM force fields incapable of reproducing accurate electrostatic potentials.37,38 Recently, the polarized protein-specific charge (PPC) model based on a fragmentation scheme39−41 for electronic structure calculations of biomolecules and the continuum dielectric model for the solvent was developed in a self-consistent manner.16 Because the PPC model correctly describes the polarized electrostatic state of a protein at a given structure, it is able to provide a more accurate description of electrostatic forces and electric fields in the vicinity of a given structure (native structure, in general). It has been demonstrated in a number of applications that the PPC model gives significantly better agreement with experimental data than standard nonpolarizable force fields in MD simulations.16,42,43 In this study, MD simulations employing different charge models (PPC and Amber ff99SB44) were carried out to calculate changes in the electrostatic field at the nitrile site between wild-type hALR2 and its eight mutants. The calculated infrared frequency shifts were compared to experimental values, and new physical insights into accurate descriptions of electric fields are discussed.

2. COMPUTATIONAL APPROACH 2.1. Wild-Type and Mutant Models of hALR2/NADP+/ IDD743 Complexes. Because the crystal structure of hALR2 bound with the inhibitor IDD743 is not currently available, we constructed a model for wild-type hALR2/NADP+/IDD743 based on the complex structure of hALR2 with a closely related inhibitor, IDD393 [Protein Data Bank (PDB) code 2PZN], as was done by Suydam et al.13 Panels A and B of Figure 1 show the molecular motifs of IDD393 and IDD743 and the location of IDD743 in the active site of the wide-type hALR2, respectively. Specifically, fluorine was substituted for chlorine in IDD743, and the bond length between fluorine and its connected carbon atom was adjusted to 1.32 Å. A nitrile group was added to the corresponding position of IDD393. The carboxylate of IDD743 was left deprotonated based on inhibition studies of the IDD594 crystal structure bound to hALR2.45 For direct comparison with the study by Suydam et al.13 using a nonpolarizable force field, we selected the same set

Figure 1. (A) Chemical structures of IDD393 and IDD743. (B) Structural model of the hALR2/NADP+/IDD743 (modified from IDD393) complex. IDD743 and hALR2 are shown in ball-and-stick and cartoon models, respectively. The side chains of eight mutated residues are shown as sticks. All the mutated residues are within 11 Å from the midpoint of the nitrile of IDD743, where the electric field was calculated.

deliver a unique vibration to the active site of hALR2. Point mutations of hALR2 that do not cause hydrogen-bonding interactions with the nitrile group were made. Vibrational Stark shifts of the nitrile stretch in IDD743 were measured between these mutations and wild-type hALR2. Δμ⃗ probe for the IDD743 nitrile stretch was fitted against the vibrational Stark spectrum and yielded a quantity of 0.041 D. By defining the nitrile axis pointing from carbon toward nitrogen, one can obtain Δμ⃗probe/ B



Article

Figure 2. Calculated electrostatic potentials for the wild-type hALR2/IDD743 complex using the PPC and Amber ff99SB models. Electrostatic potentials are plotted on the plane of the phenyl ring. (A) Calculated electrostatic potential based on the PPC model, with the interior dielectric constant set to 1. (B−D) Calculated electrostatic potentials based on the ff99SB charge model, with the interior dielectric constant set to (B) 1, (C) 2, and (D) 4. (1 kT/eÅ = 2.57 MV/cm at 298 K, where k is the Boltzmann constant, T is the temperature, and e is the charge of an electron.)

computational protocol: The Poisson−Boltzmann (PB) solver Delphi52 was first used to calculate the induced charges on the solute−solvent interface, which was defined using a probe radius of 1.4 Å. The protein was then divided into capped amino-acid fragments, and the rest of the protein fragments were employed as external point charges for quantum-chemical calculations using density functional theory at the B3LYP/631G* level. The electrostatic potentials were saved, and a standard two-stage RESP fitting procedure was used to fit effective point charges of the protein.53,54 The newly fitted atomic charges were passed to the next calculation of charge fitting. This process was iterated until the corrected reaction field energy calculated with Delphi converged and its variations were smaller than a certain criterion. Usually, the criterion was reached within five iterations. 2.3. MD Simulations and Calculation of Electric Fields. Two separate MD simulations were performed for the wild type and each of the eight mutants of hALR2 using the standard Amber ff99SB and PPC models. In the simulation using the PPC model, the atomic charges of the Amber ff99SB force field were simply replaced by PPCs, and the rest of the ff99SB parameters are retained. In the simulation, the protein was placed in a periodic rectangular box of TIP3P water molecules. The distance from the surface of the box to the closest atom of the solute was set to 12 Å. Counter ions were added to neutralize the system. Two minimization steps were employed to optimize the initial structure. In the first step, only the solvent molecules were relaxed, whereas the protein atoms were constrained to their initial structure. In the second step, the entire system was energy-minimized until convergence was reached, then the system was brought to room temperature

of eight mutated hALR2 structures (Y48F, H110A, W20Y, V47N, F121E, V47D, Q49R, and K77M) in this study. In each mutant, only one residue was mutated from the wild-type structure.13 All of the selected residues lie within 11 Å of the nitrile of IDD743 at the site where electric fields were calculated and positioned at the active-site cleft of the enzyme. The mutagenesis tool in the program PyMOL46 was used to create mutated hALR2 structures. In this program, the most frequently observed conformation reported in the Protein Data Bank was selected for the mutated residue, and it was ensured that this structure does not have any nonbonded atoms in close contact. 2.2. Force Field Parameters for NADP+ and IDD743 and Derivation of PPCs for Proteins. The initial structure of NADP+ was taken from the crystal structure of 2PZN. The geometries of NAPD+ and IDD743 were optimized at the HF/ 3-21g and HF/6-31G* levels, respectively. Subsequently, force field parameters for NADP+ and IDD743 were obtained using the ANTECHAMBER module47 based on the generalized Amber force field (GAFF)48 with the HF/6-31G* restrained electrostatic potential (RESP) charges. All ab initio calculations were carried out using the Gaussian 09 program.49 The polarized protein-specific charges (PPCs) were fitted to electrostatic potentials by fragment quantum-mechanical calculations using an iterative approach as described in ref 16. Specifically, for the nine protein structures (the wild-type and eight mutated hALR2 proteins), the missing hydrogen atoms were added using the LEaP module of the Amber program.50 A series of minimizations using the Amber ff99SB force field51 was carried out to remove close contacts. The minimized structure was then used to calculate PPCs by the MFCC-PB C



Article

Table 1. Calculated Electric Field Changes (ΔF∥, MV/cm) from Eight Mutations Relative to Wild-Type hALR2 Based on the PPC Model and Amber ff99SB Force Field, As Well As Values Calculated Using the OPLS-UA Force Field by Suydam et al.13 for Comparison OPLS-UA (ε = 2) ff99SB (ε = 1) ff99SB (ε = 2) ff99SB (ε = 4) PPC (ε = 1) expt a

13

W20Y

V47D

V47N

Y48F

Q49R

K77M

H110A

F121E

RMSDa

+0.12 −1.38 −0.95 −0.72 −0.53 −1.01

−4.47 +7.79 +3.30 +1.12 −3.01 −2.32

−3.27 −9.05 −4.88 −2.77 −4.73 −4.49

−0.38 +0.76 +0.26 +0.00 −2.58 −2.03

+0.50 −0.50 −0.38 −0.32 −1.82 −0.58

−1.71 −1.06 −0.80 −0.64 −3.64 −3.19

+0.47 −1.24 −0.74 −0.50 −1.17 +0.72

+0.52 −0.96 −0.42 −0.13 −0.31 −1.30

1.00 2.88 1.65 1.30 0.66

RMSDs are referred to the experimental values.

fluorine than that calculated using ff99SB. The electrostatic potential at the midpoint of nitrile calculated by the PPC model is also higher than that given by ff99SB. Furthermore, the positive electrostatic potential area around the coordinate of (−4, 0) based on PPCs is narrower than that based on ff99SB, but the positive potential area around the coordinate of (−4, −6) reverses the trend. The appearance of several positive electrostatic potential areas in the cross section suggests that positively charged residues might be nearby. By visually examining the X-ray structure, we found that LYS77 is located near the nitrile group, as indicated in Figure 1. To evaluate the contribution of LYS77 to the electrostatic potential in the cross section, we set the atomic charge of residue LYS77 to 0 and recalculated the electrostatic potential. In this case, most of the positive potential areas in the cross section disappeared except for the small area around the coordinate space (−4, −6) under both the PPC and ff99SB models. Through the same procedure, we found that the primary contribution to the negative electrostatic potential is given by TRP20, followed by TRY48. The range of deviation in the atomic charges of these three residues (LYS77, TRP20, TRY48) between the PPC and ff99SB models is from −0.129 to 0.126 e. As shown in Figure 2B−D, the absolute value of the electrostatic potential becomes smaller when the interior dielectric constant increases from 1 to 4, as expected. The positive potentials around both the fluorine and the nitrile decrease significantly as ε increases and tend to be neutral for ε = 4. Therefore, we conclude that the dielectric constant has a substantial effect on continuum electrostatic calculations. Duan et al. used ε = 4 to mimic the average protein environment when fitting atomic charges for each type of residue in the Amber ff03 force field.15 It is understood that the local electrostatic environment inside a protein is inhomogeneous, and thus, the atomic charges of each residue of a specific protein should be affected by other residues and solvents surrounding the protein. In the traditional nonpolarizable force field, a higher value for the interior dielectric constant is expected to remove local errors, which does not always work.58,59 The PPC model considers the electronic polarization effect explicitly and derives the charges from quantum calculations of the protein with the implicit solvent model, which is able to accurately reflect the polarization effect arising from the solvents and nearby residues. As a result, for the calculation of electric fields based on the PPC model, an interior dielectric constant of ε = 1 should be appropriate, which is no longer an undetermined parameter as opposed to calculations using traditional nonpolarizable force field. 3.2. Electric Field at the Midpoint of the Nitrile. Our previous studies42,43 demonstrated that the PPC model can maintain the ligand in a more stable location inside the protein

(300 K) in 100 ps, and another 1-ns simulation with periodic boundary conditions at 300 K and 1 atm was conducted. The integration time step was 2.0 fs, and the SHAKE algorithm55 was applied to maintain all hydrogen atoms in reasonable positions. The particle-mesh Ewald method56 was used to treat long-range electrostatic interactions, and a 10-Å cutoff for the van der Waals interactions was implemented. Langevin dynamics57 was applied to regulate the temperature with a collision frequency of 1.0 ps−1. Trajectories were generated with structures saved every 0.5 ps. All MD simulations were performed with the Amber 10 program.50 To obtain the time-averaged electric field, each saved snapshot along the 1-ns MD simulation (2000 configurations in total) was used to calculate the electric field with the Delphi program using an interior dielectric constant, ε, of 4, 2, or 1 under the ff99SB force field. Under the PPC model, the interior dielectric constant of 1 was used in the calculation of the electric field because the electronic polarization effect was explicitly included in the PPCs.58 The specific location where the electric field was calculated was chosen as the midpoint of the nitrile of IDD743, and the calculated electric field was projected onto the direction pointing from the carbon to the nitrogen along the nitrile bond. If the point charges were too close to the midpoint of the nitrile, large fluctuations in the calculated electric fields would be obtained. Thus, the charges of all atoms on the nitrile ring of IDD743 (including nitrile) were set to zero for Delphi calculations, whereas the charges for the remaining atoms of IDD743 were retained to account for the conformational flexibility of the inhibitor itself. This approach did not affect the calculated difference in electric field, because the main focus of this study was to predict the changes in the electric field brought about by mutations.

3. RESULTS AND DISCUSSION 3.1. Electrostatic Potentials Inside the Protein. To examine the difference between the PPC and ff99SB models, the cross sections of the electrostatic potential in the plane of phenyl ring of IDD743 for the energy-minimized wild-type hALR2/IDD743 complex are plotted in Figure 2. Panels A and B of Figure 2 show the calculated electrostatic potentials using the same interior dielectric constant of 1 (ε = 1) based on the PPC and ff99SB models, respectively. Owing to the mean-fieldlike feature of the ff99SB charge model, the dielectric constant for the protein is an uncertain parameter.58 We also plot the electrostatic potentials in the same plane using the ff99SB charge model with ε = 2 and 4, as shown in Figure 2C,D. Comparison of panels A and B of Figure 2 shows that the most significant difference in the electrostatic potentials appears around the fluorine. The calculated electrostatic potential based on PPCs gives a more positive potential area around the D



Article

Figure 3. Calculated electric fields as a function of MD simulation time for the wild type (black line) and two mutants (K77M and V47D, red line) of hALR2. The left and right columns represent the calculated results based on the Amber ff99SB and PPC models, respectively. The time-averaged electric fields for the wild type and mutants are shown in green and blue lines, respectively. In addition, a yellow line that shows the time-averaged electric field for the last 800 ps is plotted for the mutant V47D based on the Amber force field. The interior dielectric constant was set to 1 for the electric field calculations.

experimental value (−2.1 versus −1.6 cm−1 as shown in Table 1). To further demonstrate how the calculated electric fields are affected by the electronic polarization effect of protein, we plot the side-chain positions of mutated residue MET77 in K77M and ASP47 in V47D during the MD simulation (see Figure 4). The coordinates of side-chain atoms were projected onto the plane of the phenyl ring of IDD743. As can be seen from Figure 4, the populated position of residue MET77 during the PPC simulation is close to that based on the ff99SB simulation for K77M, and the relative position of MET77 with respect to the phenyl ring of IDD743 during the MD simulation is near the initial position. Therefore, the difference in calculated electric field shifts between the PPC and ff99SB models might not be caused mainly by the relative position distinction of MET77 during the MD simulation, but rather might originate from the difference between the PPC and ff99SB models, which directly affects the continuum electrostatics calculations. On the other hand, during the ff99SB simulation for V47D (see Figure 4), the side-chain atoms of ASP47 moved away from the initial position to a different geometric center relative to the phenyl ring of IDD743 in a few picoseconds, which led the calculated electric field to change dramatically from negative (less than −8 kT/eÅ) to positive as indicated in Figure 3. On the contrary, one can see from Figure 4 that the side chain of ASP47 remained close to the initial position during the PPC simulation. The stabilization of the relative position is due to the polarization effect embedded in the PPC model, which has been shown to be able to correctly describe the hydrogen-bond network in proteins as demonstrated in our recent studies.42,61 For V47D, the discrepancy between experiment and the electric field shift calculated from the ff99SB simulation might also arise from the incorrect sampling of the configurational space of the protein. Therefore, we conclude that the electronic polarization effect not only influences the direct calculation of the electric field for a given protein structure, but also offers a better

than traditional nonpolarizable force fields. In addition, the electronic polarization effect is important in stabilizing the hydrogen bonds, which is critical to preserving the native secondary structure of proteins,60,61 and thus the PPC model is able to provide better conformational sampling,16 which is always required for accurate calculations of the electric field.13 We plot the root-mean-square deviation (RMSD) of protein backbone atom positions as a function of MD simulation time for wild-type and eight mutated hALR2 structures using the ff99SB and PPC models, as shown in Figures S1 and S2, respectively, of the Supporting Information. The RMSDs are all smaller than 2 Å with respect to the initial hALR2/IDD743 complex structure, showing that the protein trajectory is close to the native structure during the 1-ns MD simulation. To investigate the electronic polarization effect, we compare the calculated electric field changes by both the ff99SB and PPC models on two mutants, namely, K77M and V47D, which caused large vibrational frequency shifts of −2.2 and −1.6 cm−1, respectively (see Table 1), relative to the wild type according to experimental observations. Figure 3 shows the calculated electric field as a function of MD simulation time for the wild type and mutants K77M and V47D. For direct comparison of the ff99SB and PPC models, we set the interior dielectric constant equal to 1 in both calculations. One can see from Figure 3 that the shift of electric field by K77M (relative to the wild type) based on the PPC model is larger than that calculated by ff99SB, so that the vibrational frequency shift calculated by the PPC model exhibits better agreement with experiment (−2.5 versus −2.2 cm−1, where the calculated value is the ensemble average over the final 1600 snapshots in the last 800 ps of the simulation) than ff99SB, which predicted a much smaller frequency shift of −0.7 cm−1. On the other hand, for mutant V47D, the electric field calculated using ff99SB has a positive frequency shift of 5.4 cm−1. In contrast, the PPC model predicted a negative frequency shift, which is in accord with the E



Article

Figure 4. Distribution of projected coordinates of selected atoms of MET77 (top) and ASP47 (bottom) onto the plane of the phenyl ring of IDD743 during MD simulations based on the Amber ff99SB (left column) and PPC (right column) models. The triangles in the black circles represent the initial projected positions of the corresponding atoms.

conformational sampling for proteins, which is an essential prerequisite for obtaining any ensemble-averaged property of proteins. 3.3. Calculation of the Vibrational Frequency Shifts. Figure 5 shows the electric field calculated at the midpoint of the nitrile of IDD743 based on the PPC model using ε = 1 as a function of MD simulation time. For comparison, the electric fields calculated from ff99SB with ε = 4, 2, and 1 are shown in Figure 6 and Figures S3 and S4 of the Supporting Information, respectively. As can be seen from these figures, the calculated electric field exhibits an approximately linear dependence on the interior dielectric constant. In addition, the calculated electric fields mostly become steady after 200 ps, assuring that 1 ns of MD simulation is sufficient for the time-averaged field calculation. By comparing Figure 5 with Figure S4 (Supporting Information), one can see that the fluctuations of the calculated electric fields for some mutants based on the Amber ff99SB force field are slightly larger than those calculated using the PPC model, indicating that the polarized force field sustains a more stable local chemical environment than the Amber ff99SB force field. The calculated electric fields of the wild-type and eight mutated hALR2 structures based on both the ff99SB and PPC models are able to converge, except that for V47N, in

Figure 5. Electric fields for the wild ype and eight mutants of hALR2 calculated using the PPC model as a function of MD simulation time. The red line represents the time-averaged electric field. The interior dielectric constant was set to 1.

which the polar side chain of residue ASN47 is very close to the location where the electric field was calculated. A subtle swing F



Article

Figure 6. Electric field of the wild type and eight mutants of hALR2 calculated using the Amber ff99SB force field as a function of MD simulation time. The red line represents the time-averaged electric field. The interior dielectric constant was set to 4.

Figure 7. Correlation between the observed changes in frequency and the calculated changes in electric field brought from eight mutations relative to wild-type hALR2 based on the Amber ff99SB force field. Triangles, solid circles, and open circles represent results calculated with interior dielectric constants of 1, 2, and 4, respectively. The symbols in order from top to bottom denote the mutants H110A, Q49R, W20Y, F121E, Y48F, V47D, K77M, and V47N. The black line represents the linear correlation between changes in frequency and changes in electric field observed from experiment (see eq 2).

of the side chain of ASN would lead to large electric field change. In this case, the change of the electric field brought from positional variation of the side chain would far outweigh the mutation. As suggested by Suydam et al.,13 eight different initial structures were first obtained from MD simulations, and then multiple MD trajectories starting from these eight structures were carried out for electric field calculations. The most stable area (the data between the two red vertical lines for V47N in Figure 5) was used to obtain the statistically meaningful value of the time-averaged electric field for V47N, whereas for the wild type and other seven mutants, 1600 snapshots in the final 800 ps of the simulation (one snapshot every 0.5 ps) were used to calculate the ensemble average of the electric fields. The calculated electric field changes between the wild type and mutants were compared with measured Stark shifts. To reduce the possible errors resulting from the uncertainty in the interior dielectric constant for the nonpolarizable force field, we calculated the electric field changes caused by the eight mutants based on the Amber ff99SB force field using different dielectric constants (1, 2, and 4), as shown in Figure 7. The black line is the experimental fitting curve. From Figure 7, one can see that the absolute changes in electric fields for six mutants (H110A, Q49R, W20Y, F121E, Y48F, and K77M) are all relatively small (less than 1.0 MV/cm) for both ε = 2 and 4. In contrast, mutants V47N and V47D have large electric field changes as calculated by ff99SB for ε = 2 and 4. V47D was predicted to have a positive electric field shift, which is contradictory to the experimental observation. The reason for this discrepancy was elucidated in the previous section. For ε = 1, all mutants give larger electric field changes than the wild type. However, the results for mutants V47D and V47N are farther from the experimental values. As a result, we conclude that the value of the interior dielectric constant has a significant impact on the calculated electric field, which increases approximately linearly with the decrease of the dielectric constant. When the dielectric constant equals 4, the calculated changes in electric fields have the smallest RMSD with respect to the experimental values, and the trend of the results agrees better with the experimental fitting curve as compared to those for ε = 2 and 1.

Figure 8 shows the calculated electric field changes based on the PPC model with ε = 1. For comparison, the results from

Figure 8. Similar to Figure 7, but the changes in electric field were calculated based on the PPC model with ε = 1 (red solid circles), the Amber ff99SB force field with ε = 4 (hollow circles), and the OPLSUA force field with ε = 4 (data from ref 13).

the ff99SB and OPLS united-atom (UA) force fields (computed by Suydam et al.13) with ε = 4 are also presented in Figure 8. Overall, the PPC model gives the smallest RMSD with respect to the experimental values. The calculated changes in electric fields based on PPCs have the best agreement with the observed Stark shifts, especially for the point mutations (Y48F, V47D, K77M, and V47N), which caused relatively larger electric field changes. For the results based on the Amber ff99SB force field with ε = 4, most of the calculated electric field changes have large deviations from experiment. The changes from mutants of F121E, Y48F, and K77M based on ff99SB are not as significant as the experimental values. The results for V47D and V47N also have considerable errors. On the other G



Article

hand, results calculated by Suydam et al. were13 using the OPLS-UA force field and the generalized Born implicit solvation model. For direct comparison with ff99SB, the electric fields calculated with ε = 4 were chosen for OPLS-UA. The discrepancy between the OPLS-UA-calculated results and our ff99SB calculations might arise from the different choice of nonpolarizable force fields and solvation models. Suydam et al.13 also showed that the electric fields calculated using the OPLS-UA force field exhibit better agreement with the experimental values for ε = 2. As shown in Table 1, the results based on the PPC model have the lowest RMSD of 0.66 MV/ cm with respect to the experimental values. The results calculated based on the OPLS-UA force field with ε = 2 have an RMSD of 1.00 MV/cm, which is smaller than the best result predicted by the Amber ff99SB force field (RMSD = 1.30 MV/ cm with ε = 4). Our study clearly emphasizes that the electronic polarization effect is indispensible for accurately reproducing the observed Stark shifts inside proteins.

(2) Honig, B.; Nicholls, A. Classical Electrostatics in Biology and Chemistry. Science 1995, 268, 1144−1149. (3) Matthew, J. B. Electrostatic Effects in Proteins. Annu. Rev. Biophys. Biophys. Chem. 1985, 14, 387−417. (4) Perutz, M. F. Electrostatic Effects in Proteins. Science 1978, 201, 1187−1191. (5) Warshel, A.; Russell, S. T. Calculations of Electrostatic Interactions in Biological Systems and in Solutions. Q. Rev. Biophys. 1984, 17, 283−422. (6) Simonson, T.; Archontis, G.; Karplus, M. Continuum Treatment of Long-Range Interactions in Free Energy Calculations. Application to Protein−Ligand Binding. J. Phys. Chem. B 1997, 101, 8349−8362. (7) Lee, L. P.; Tidor, B. Optimization of Binding Electrostatics: Charge Complementarity in the Barnase-Barstar Protein Complex. Protein Sci. 2001, 10, 362−377. (8) Steffen, M. A.; Lao, K.; Boxer, S. G. Dielectric Asymmetry in the Photosynthetic Reaction Center. Science 1994, 264, 810−816. (9) Okamura, M. Y.; Feher, G. Proton Transfer in Reaction Centers from Photosynthetic Bacteria. Annu. Rev. Biochem. 1992, 61, 861−896. (10) Hilvert, D. Critical Analysis of Antibody Catalysis. Annu. Rev. Biochem. 2000, 69, 751−793. (11) Warshel, A.; Sharma, P. K.; Kato, M.; Xiang, Y.; Liu, H.; Olsson, M. H. Electrostatic Basis for Enzyme Catalysis. Chem. Rev. 2006, 106, 3210−3235. (12) Xiang, Y.; Duan, L.; Zhang, J. Z. Protein’s Electronic Polarization Contributes Significantly to Its Catalytic Function. J. Chem. Phys. 2011, 134, 205101. (13) Suydam, I. T.; Snow, C. D.; Pande, V. S.; Boxer, S. G. Electric Fields at the Active Site of an Enzyme: Direct Comparison of Experiment with Theory. Science 2006, 313, 200−204. (14) Abbyad, P.; Shi, X.; Childs, W.; McAnaney, T. B.; Cohen, B. E.; Boxer, S. G. Measurement of Solvation Responses at Multiple Sites in a Globular Protein. J. Phys. Chem. B 2007, 111, 8269−8276. (15) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T. A Point-Charge Force Field for Molecular Mechanics Simulations of Proteins Based on Condensed-Phase Quantum Mechanical Calculations. J. Comput. Chem. 2003, 24, 1999−2012. (16) Ji, C. G.; Mei, Y.; Zhang, J. Z. H. Developing Polarized ProteinSpecific Charges for Protein Dynamics: MD Free Energy Calculation of pKa Shifts for Asp26/Asp20 in Thioredoxin. Biophys. J. 2008, 95, 1080−1088. (17) Xu, L.; Cohen, A. E.; Boxer, S. G. Electrostatic Fields near the Active Site of Human Aldose Reductase: 2. New Inhibitors and Complications Caused by Hydrogen Bonds. Biochemistry 2011, 50, 8311−8322. (18) Sandberg, D. J.; Rudnitskaya, A. N.; Gascón, J. A. QM/MM Prediction of the Stark Shift in the Active Site of a Protein. J. Chem. Theory Comput. 2012, 8, 2817−2823. (19) Bublitz, G. U.; Boxer, S. G. Stark Spectroscopy: Applications in Chemistry, Biology, And Materials Science. Annu. Rev. Phys. Chem. 1997, 48, 213−242. (20) Hush, N. S.; Reimers, J. R. Vibrational Stark Spectroscopy. 1. Basic Theory and Application to the CO Stretch. J. Phys. Chem. 1995, 99, 15798−15805. (21) Webb, L. J.; Boxer, S. G. Electrostatic Fields near the Active Site of Human Aldose Reductase: 1. New Inhibitors and Vibrational Stark Effect Measurements. Biochemistry 2008, 47, 1588−1598. (22) Andrews, S. S.; Boxer, S. G. Vibrational Stark Effects of Nitriles I. Methods and Experimental Results. J. Phys. Chem. A 2000, 104, 11853−11863. (23) Andrews, S. S.; Boxer, S. G. Vibrational Stark Effects of Nitriles II. Physical Origins of Stark Effects from Experiment and Perturbation Models. J. Phys. Chem. A 2002, 106, 469−477. (24) Park, E. S.; Andrews, S. S.; Hu, R. B.; Boxer, S. G. Vibrational Stark Spectroscopy in Proteins: A Probe and Calibration for Electrostatic Fields. J. Phys. Chem. B 1999, 103, 9813−9817.

4. CONCLUSIONS Molecular dynamics simulations and electrostatic field calculations using the PPC and Amber ff99SB models have been carried out for the wild type and eight mutants of hALR2. Comparison with experimentally observed Stark shifts shows that calculations based on the PPC model give significantly improved result (change of electric field due to mutation) over those obtained using nonpolarizable force fields, especially for mutations that lead to large electric field changes in hALR2. Because the PPC model explicitly incorporates the polarization effect of proteins and solvents, our study demonstrates that the electronic polarization effect is crucial to accurate predictions of the electric fields in proteins.

■

ASSOCIATED CONTENT

S Supporting Information *

RMSDs of the protein backbone atoms during MD simulations based on the Amber ff99SB force field and PPC model, respectively (Figures S1 and S2). Calculated electric fields of the wild-type and eight mutants of hALR2 based on the Amber ff99SB force field as a function of MD simulation time, where the interior dielectric constant was set to 2 (Figure S3) and 1 (Figure S4). This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected] (X.H.), [email protected] (J.Z.H.Z.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the National Natural Science Foundation of China (Grants 10974054 and 20933002) and Shanghai PuJiang program (09PJ1404000). We also thank the Computational Center of ECNU for providing computational time.

■

REFERENCES

(1) Davis, M. E.; McCammon, J. A. Electrostatics in Biomolecular Structure and Dynamics. Chem. Rev. 1990, 90, 509−521. H



Article

NADPH and Alrestatin-like Inhibitors. Biochemistry 1994, 33, 7157− 7165. (46) DeLano, W. L. The PyMOL Molecular Graphics System; Version 1.5.0.1, DeLano Scientific: San Carlos, CA; 2012. (47) Wang, J. M.; Wang, W.; Kollman, P. A.; Case, D. A. Automatic Atom Type and Bond Type Perception in Molecular Mechanical Calculations. J. Mol. Graphics Modell. 2006, 25, 247−260. (48) Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (49) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Keith, T.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, Ö .; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, revision B.01; Gaussian, Inc.: Wallingford, CT, 2010. (50) Case, D. A.; Cheatham, T. E.; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. The Amber Biomolecular Simulation Programs. J. Comput. Chem. 2005, 26, 1668−1688. (51) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (52) Rocchia, W.; Sridharan, S.; Nicholls, A.; Alexov, E.; Chiabrera, A.; Honig, B. Rapid Grid-Based Construction of the Molecular Surface and the Use of Induced Surface Charge to Calculate Reaction Field Energies: Applications to the Molecular Systems and Geometric Objects. J. Comput. Chem. 2002, 23, 128−137. (53) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. A WellBehaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The RESP Model. J. Phys. Chem. 1993, 97, 10269−10280. (54) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Kollman, P. A. Application of RESP Charges to Calculate Conformational Energies, Hydrogen Bond Energies, and Free Energies of Solvation. J. Am. Chem. Soc. 1993, 115, 9620−9631. (55) Jean-Paul, R.; Giovanni, C.; Herman, J. C. B. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of n-Alkanes. J. Comput. Phys. 1977, 23, 327−341. (56) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N·log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (57) Pastor, R. W.; Brooks, B. R.; Szabo, A. An Analysis of the Accuracy of Langevin and Molecular Dynamics Algorithms. Mol. Phys. 1988, 65, 1409−1419. (58) Schutz, C. N.; Warshel, A. What Are the Dielectric “Constants” of Proteins and How to Validate Electrostatic Models? Proteins: Struct., Funct., Genet. 2001, 44, 400−417. (59) Gilson, M. K.; Honig, B. H. The Dielectric Constant of a Folded Protein. Biopolymers 1986, 25, 2097−2119. (60) Ji, C. G.; Zhang, J. Z. H. Electronic Polarization Is Important in Stabilizing the Native Structures of Proteins. J. Phys. Chem. B 2009, 113, 16059−16064. (61) Tong, Y.; Ji, C. G.; Mei, Y.; Zhang, J. Z. H. Simulation of NMR Data Reveals That Proteins’ Local Structures Are Stabilized by Electronic Polarization. J. Am. Chem. Soc. 2009, 131, 8636−8641.

(25) Reimers, J. R.; Zeng, J.; Hush, N. S. Vibrational Stark Spectroscopy. 2. Application to the CN Stretch in HCN and Acetonitrile. J. Phys. Chem. 1996, 100, 1498−1504. (26) Suydam, I. T.; Boxer, S. G. Vibrational Stark Effects Calibrate the Sensitivity of Vibrational Probes for Electric Fields in Proteins. Biochemistry 2003, 42, 12050−12055. (27) Ringer, A. L.; MacKerell, A. Calculation of the Vibrational Stark Effect Using a First-Principles Quantum Mechanical/Molecular Mechanical Approach. J. Phys. Chem. Lett. 2011, 2, 553−556. (28) Lindquist, B. A.; Corcelli, S. A. Nitrile Groups As Vibrational Probes: Calculations of the CN Infrared Absorption Line Shape of Acetonitrile in Water and Tetrahydrofuran. J. Phys. Chem. B 2008, 112, 6301−6303. (29) Lindquist, B. A.; Furse, K. E.; Corcelli, S. A. Nitrile Groups As Vibrational Probes of Biomolecular Structure and Dynamics: An Overview. Phys. Chem. Chem. Phys. 2009, 11, 8119−8132. (30) Lindquist, B. A.; Haws, R. T.; Corcelli, S. A. Optimized Quantum Mechanics/Molecular Mechanics Strategies for Nitrile Vibrational Probes: Acetonitrile and para-Tolunitrile in Water and Tetrahydrofuran. J. Phys. Chem. B 2008, 112, 13991−4001. (31) Getahun, Z.; Huang, C. Y.; Wang, T.; De Leon, B.; DeGrado, W. F.; Gai, F. Using Nitrile-Derivatized Amino Acids As Infrared Probes of Local Environment. J. Am. Chem. Soc. 2003, 125, 405−411. (32) Franzen, S. An Electrostatic Model for the Frequency Shifts in the Carbonmonoxy Stretching Band of Myoglobin: Correlation of Hydrogen Bonding and the Stark Tuning Rate. J. Am. Chem. Soc. 2002, 124, 13271−13281. (33) Choi, J. H.; Cho, M. Vibrational Solvatochromism and Electrochromism of Infrared Probe Molecules Containing CO, CN, CO, or CF Vibrational Chromophore. J. Chem. Phys. 2011, 134, 154513. (34) Saggu, M.; Levinson, N. M.; Boxer, S. G. Direct Measurements of Electric Fields in Weak OH···π Hydrogen Bonds. J. Am. Chem. Soc. 2011, 133, 17414−17419. (35) Saggu, M.; Levinson, N. M.; Boxer, S. G. Experimental Quantification of Electrostatics in X−H···π Hydrogen Bonds. J. Am. Chem. Soc. 2012, 134, 18986−18997. (36) Layfield, J. P.; Hammes-Schiffer, S. Calculation of Vibrational Shifts of Nitrile Probes in the Active Site of Ketosteroid Isomerase upon Ligand Binding. J. Am. Chem. Soc. 2013, 135, 717−725. (37) Gascón, J. A.; Leung, S. S. F.; Batista, E. R.; Batista, V. S. A SelfConsistent Space-Domain Decomposition Method for QM/MM Computations of Protein Electrostatic Potentials. J. Chem. Theory Comput. 2006, 2, 175−186. (38) Menikarachchi, L. C.; Gascón, J. A. Optimization of Cutting Schemes for the Evaluation of Molecular Electrostatic Potentials in Proteins via Moving-Domain QM/MM. J. Mol. Model. 2008, 14, 479− 487. (39) Gao, A. M.; Zhang, D. W.; Zhang, J. Z. H.; Zhang, Y. K. An Efficient Linear Scaling Method for ab Initio Calculation of Electron Density of Proteins. Chem. Phys. Lett. 2004, 394, 293−297. (40) Mei, Y.; Ji, C.; Zhang, J. Z. A New Quantum Method for Electrostatic Solvation Energy of Protein. J. Chem. Phys. 2006, 125, 094906. (41) Zhang, D. W.; Zhang, J. Z. H. Molecular Fractionation with Conjugate Caps for Full Quantum Mechanical Calculation of Protein− Molecule Interaction Energy. J. Chem. Phys. 2003, 119, 3599−3605. (42) Ji, C. G.; Zhang, J. Z. H. Protein Polarization Is Critical to Stabilizing AF-2 and Helix-2′ Domains in Ligand Binding to PPAR-γ. J. Am. Chem. Soc. 2008, 130, 17129−17133. (43) Tong, Y.; Mei, Y.; Li, Y. L.; Ji, C. G.; Zhang, J. Z. H. Electrostatic Polarization Makes a Substantial Contribution to the Free Energy of Avidin−Biotin Binding. J. Am. Chem. Soc. 2010, 132, 5137−5142. (44) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins: Struct., Funct., Bioinf. 2006, 65, 712−725. (45) Ehrig, T.; Bohren, K. M.; Prendergast, F. G.; Gabbay, K. H. Mechanism of Aldose Reductase Inhibition: Binding of NADP+/ I


Predicting Mutation-Induced Stark Shifts in the Active Site of a

Recommend Documents