Article pubs.acs.org/JPCB
Over-Destabilization of Protein−Protein Interaction in Generalized Born Model and Utility of Energy Density Integration Cutoff Yukinobu Mizuhara,†,∥ Dan Parkin,‡,∥ Koji Umezawa,§,⊥ Jun Ohnuki,† and Mitsunori Takano*,† †
Department of Pure and Applied Physics, Waseda University, Okubo 3-4-1, Shinjuku-Ku, Tokyo 169-8555, Japan Department of Advanced Science and Engineering, Waseda University, Okubo 3-4-1, Shinjuku-Ku, Tokyo 169-8555, Japan § Department of Biomedical Engineering and ⊥Institute for Biomedical Sciences, Shinshu University, 8304 Minami-minowa, Kami-ina, Nagano 399-4598, Japan ‡
S Supporting Information *
ABSTRACT: The generalize Born (GB) model is frequently used in MD simulations of biomolecular systems in aqueous solution. The GB model is usually based on the so-called Coulomb field approximation (CFA) for the energy density integration. In this study, we report that the GB model with CFA overdestabilizes the long-range electrostatic attraction between oppositely charged molecules (ionic bond forming two-helix system and kinesin-tubulin system) when the energy density integration cutoff, rmax, which is used to calculate the Born energy, is set to a large value. We show that employing large rmax, which is usually expected to make simulation results more accurate, worsens the accuracy so that the attraction is changed into repulsion. It is demonstrated that the overdestabilization is caused by the overestimation of the desolvation penalty upon binding that originates from CFA. We point out that the overdestabilization can be corrected by employing a relatively small cutoff (rmax = 10−15 Å), affirming that the GB models, even with CFA, can be used as a powerful tool to theoretically study the protein−protein interaction, particularly on its dynamical aspect, such as binding and unbinding.
■
INTRODUCTION The association−dissociation dynamics of biomolecules is the physical basis of all of the activities of living organisms. In addition to the hydrophobic interaction, the electrostatic interaction has been recognized as an important factor in the assembly of proteins1 and known to play a vital role in fast association.2,3 Accurate modeling of the electrostatic interaction has been a challenge of computational studies because of the intrinsic difficulty in dealing with the heterogeneous dielectric environment in water at the intermediate length-scale between the microscopic and the macroscopic regimes.2−9 To calculate the electrostatic energies of such a heterogeneous dielectric system, the continuum model based on the Poisson equation (referred to as PE hereafter) has been widely used, where the solvent (water) is treated as a dielectric medium. While the techniques to numerically solve PE have been greatly developed and successfully applied to biomolecular systems in water,5,6,10,11 the generalized Born (GB) model12−14 has attracted much attention as an alternative of the PE approach because it can yield the approximate solution of PE with much less computational cost.14−16 Furthermore, due to the invention of the pairwise analytical and differentiable forms for the atomic self-energy (Born energy)15,17−19 and for the interatomic Coulomb energy,12 which can substitute for the numerical integration of the energy density of the electric field, the GB © XXXX American Chemical Society
model is now widely used in the molecular dynamics (MD) simulation of biomolecular systems. In both PE and GB models, a protein molecule in water is modeled by assigning a low dielectric constant to the interior of the protein (εin) and a high dielectric constant to the exterior of the protein (εex) that is filled with water, usually with εin = 1 and εex = εw where εw is the dielectric constant of the bulk water (εw ∼ 80). Once the boundary of the two dielectric constants and the spatial distribution of charges are set, the electrostatic potential is determined according to PE, and hence the electrostatic energy of the system is obtained by integrating the energy density. In the GB model, the electrostatic energy of the system is obtained by adding the solvation energy (ΔGsol) to the electrostatic energy in vacuum (Evac el ) (the state with εex = 1, i.e., in vacuum, is usually chosen as the reference). ΔGsol is composed of the self-energy term with the Born-like form and the interaction energy term with the Coulomb-like form, where the “effective Born radius”, Reff, plays the central role.12 Reff reflects the effective dielectric environment around the charged atom; Reff becomes larger as the effective dielectric constant around the atom is lowered due to the increase of surrounding Received: December 18, 2016 Revised: April 12, 2017 Published: April 20, 2017 A
DOI: 10.1021/acs.jpcb.7b01438 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B low dielectric (proteinous) region. The accuracy of the GB model (compared to PE) relies on the accurate calculation of Reff,20 which is improved by proper setting of the dielectric boundary and the consideration of the reaction field.15,20−30 For fast calculation, an approximate evaluation of the reaction field is frequently employed, which is known as the “Coulomb field approximation (CFA)”.12,15,17 Since Reff is obtained by the spatial integral of the energy density, it is practically necessary to truncate the integral (summation) at a certain distance from the atom for which Reff is calculated. While the lower limit of the integral, which is given by the intrinsic radius of the atom, has been shown to exert significant influence on the stability of hydrogen-bonds and salt-bridges,31−33 the effect of the upper limit, or the integration cutoff, on the evaluation of Reff has been considered to be little32,34 and therefore has been less studied. However, for systems exhibiting such large-scale molecular dynamics as protein folding or protein−protein binding, where the dielectric boundary is largely changed, the integration cutoff may have substantial influence on the evaluation of Reff and accordingly on the stability of the folded or the bound state. In this study, we show that the energy density integration cutoff does influence substantially the electrostatic interaction between proteins, in such a manner that the accuracy of the GB model is worsened and the attractive electrostatic interaction between oppositely charged molecules becomes unstable as the cutoff is increased, in spite of the expectation that more accurate results would be obtained with a larger cutoff. We show that the accuracy worsening with increasing cutoff originates from the overestimation of the unfavorable solvation energy (or desolvation penalty) upon binding, and that this overdestabilization is due to CFA employed in the GB models.18,19,27,28,30 We also point out that this overdestabilization is remedied by employing a relatively small cutoff (∼10 Å): note that this cutoff is usually set to a larger value in MD simulation using the GB models (for example, in AMBER,35 this cutoff, “RGBMAX”, is set to 25 Å by default). Therefore, from the practical standpoint, our study indicates that the GB models with CFA is still quite useful and powerful, because a small cutoff, which makes computation less demanding, can well reproduce the stability of the protein−protein interaction in water.
Gel(d) =
ΔGsol(d) =
+ ΔGsol
2ρ2
+
q1q2 d
+ ΔGsol
(2)
⎞ q2 ⎞ q2 1⎛ 1 1⎛ 1 − 1⎟ 1 + ⎜ − 1⎟ 2 ⎜ 2 ⎝ ε1(d) 2 ⎝ ε2(d) ⎠ ρ2 ⎠ ρ1 ⎛ 1 ⎞q q +⎜ − 1⎟ 1 2 ⎝ ε3(d) ⎠ d
(3)
The first and second terms of the right-hand side of eq 2 are the self-energies, the third term is the interaction energy, and ε’s in eq 3 represent effective dielectric constants in water: ε1, ε2, and ε3 take the dielectric constant of bulk water (εw) when d is large enough, whereas they are lowered as d becomes smaller because the regions with a low dielectric constant come closer (we assume εin = 1 inside the particles). At an infinite separation, eq 3 becomes, ΔGsol(∞) = ΔGsol,1(∞) + ΔGsol,2(∞)
(4)
where ⎞q2 1⎛ 1 − 1⎟ i (i = 1, 2) ⎜ 2 ⎝ εw ⎠ ρi
ΔGsol, i(∞) =
(5)
36
The Born energy appears in its original form in eq 5, which is derived by subtracting the spatial integral of the energy density of the electric field due to particle i in vacuum from that in water, i.e., 1 8π
ΔGsol, i(∞) =
−
=
1 8π
∫|r |≥ρ i
i
∫ Ei(ri)·Di(ri)dri 1 8π
∫ Eivac(ri)·Divac(ri)dri
⎛ q2 q2 ⎞ ⎜ i − i ⎟dri ⎜ ε |r | 4 |ri|4 ⎟⎠ ⎝ w i
(6)
(7)
⎞q2 1⎛ 1 = ⎜ − 1⎟ i 2 ⎝ εw ⎠ ρi
(8)
where Ei and Di are the electric field and the electric flux density due to the particle i, and ri = x−xi, with xi being the positional vector of particle i. In the GB model,12 an approximate functional form of eq 3 is introduced, i.e.,
THEORETICAL AND COMPUTATIONAL METHODS GB Model and Cutoff for Energy Density Integral. We first describe the basic framework of the GB model for general readers, so that those who are familiar with the GB model can skip to eq 11, after which the cutoff for the energy density integral is described. As mentioned in the introduction, the electrostatic energy of the system in water, denoted as Gel, can be obtained by adding the electrostatic component of the solvation energy ΔGsol to the electrostatic energy in vacuum Evac el , i.e., Gel =
2ρ1
+
q22
where
■
Eelvac
q12
⎛ ⎜ q2 ⎛ ⎞ q2 1 1 GB (d ) = ⎜ ΔGsol − 1⎟⎜ 1 + 2 2 ⎝ εw R2 ⎠⎜⎜ R1 ⎝
(1)
+
2q1q2 −d 2 4R1R 2
( )
d 2 + R1R 2exp
For simplicity, consider a system composed of two spherical particles solvated in water; the two particles possess charges q1 and q2 at their centers, the sphere radii are ρ1 and ρ2, respectively, and the distance between the centers of the particles is d. Then, the electrostatic energy of this system is given by,
⎞ ⎟ ⎟ ⎟ ⎟ ⎠
(9)
which takes into account the lowering of the effective dielectric constant with smaller d and approaches asymptotically to the solvation energy of the Onsager’s reaction field for a dipole (oppositely signed charge pair) at small d limit, and approaches B
DOI: 10.1021/acs.jpcb.7b01438 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B ⎛ ρj d − ρj ⎞ 1 ⎟ + f (d ; ρj ) = ⎜⎜ log 2 2 d + ρj ⎟⎠ 4d ⎝ 2(d − ρj )
asymptotically to eq 3 at large d limit. R1 and R2, called the effective Born radii, depend on d and act as the substitutes for the effective dielectric constants ε1 and ε2 in eq 3, as defined by, ⎞ qi2 1⎛ 1 − 1⎟ = ΔGsol, i(d)(i = 1, 2) ⎜ 2 ⎝ εw ⎠ Ri
(if d ≥ ρi + ρj )
The above framework for the two particle system can be extended to a system consisting of many particles (atoms) such as a macromolecule system. For a macromolecule system, the spatial integral of |r|−4 as in eq 14 runs over the interior of the solute macromolecule(s), which can be approximated by summation of the pairwise function in eq 16,17−19 i.e.,
(10)
ΔGsol,i(d), on the other hand, is obtained by calculating the spatial integral of the energy density as in eq 6, ΔGsol, i(d) =
1 8π
∫ Ei(ri)·Di(ri)dri
1 − 8π
∫
Eivac(ri) ·Divac (ri)dri
1 4π
(11)
=
1 4π
i
i
∫|r |≥ρ i
i
|ri − (xj − xi)|≤ ρj , j ≠ i
⎛ ⎞⎜ 4π qi2 ⎛ 1 = − ⎜ − 1⎟⎜ 8π ⎝ εw ⎠⎜⎜ ρi ⎝
∫|r |≥ρ i
i
|ri − (xj − xi)|≤ ρj , j ≠ i
⎞ ⎟ 1 ⎟ 4 dri ⎟ |ri| ⎟ ⎠
∫|r |≥ρ i
i
|ri − (xj − xi)|≤ ρj , j ≠ i
∑ f (dij ; sjρj ) (17)
i
i
1 1 dri ≃ 4π |ri|4
∑
∑ j(≠ i)
∫ρ ≤|r |≤r i
i
max
1 dri |ri|4
|ri − (xj − xi)|≤ sjρj
f (dij ; sjρj ) (18)
For example, in AMBER, the cutoff distance rmax, referred to as “RGBMAX”, is employed with the default value of 25 Å. MD Simulation and Intermolecular Energy Calculation. Model System To Investigate Protein−Protein Interaction. We first used a model system composed of two 15residue alanine-based α-helices that interact electrostatically with each other: one containing lysine at the center (Ala7-LysAla7) and the other containing glutamic acid at the center (Ala7-Glu-Ala7). Lysine was positively charged (+1e; e is the elementary charge) and glutamic acid was negatively charged (−1e). N and C termini were neutralized by acetylation and Nmethylamidation, respectively. The two helices (each of them was pre-equilibriated in explicit water) were aligned so that the helix axes, calculated by TWISTER,37 become parallel to each other, with Lys and Glu facing each other as shown in Figure 1a. The distance between Nζ of Lys and Cδ of Glu, denoted by d, was used for the reaction coordinate for the ionic bond formation between Lys and Glu. The potential of mean force (PMF) as a function of d was calculated using the umbrella sampling technique with the weighted histogram method.38 Four independent 20-ns MD runs at 300 K (using Langevin thermostat) were performed in 1 each window with the umbrella potential 2 k(d − d0)2 , where d0 was shifted from 3 to 40 Å at the intervals of 0.5 Å, resulting in 75 windows (6 μs-long sampling in total), and k was set to 2 kcal/mol/Å2. As the GB model, we used the one proposed by Onufriev et al.,27 which is one of the most popular GB models, with mBondi2 parameters for ρ and with εin = 1 and εex = 78.5. As for the energy density integration cutoff, four different cutoff
⎞ ⎟ 1 ⎟ 4 dri ⎟ |ri| ⎟ ⎠
(13)
1 dri |ri|4 (14)
and the integral in eq 14 can be further replaced by an analytical function,17 1 1 = − f (d ; ρj ) Ri ρi
i
35
By combining eq 10 with eq 13, the effective Born radius is given by, 1 1 1 = − ρi Ri 4π
i
1 dri |ri|4
|ri − (xj − xi)|≤ sjρj
{j | dij − sjρj ≤ rmax }
|ri − (xj − xi)|≥ ρj , j ≠ i
∫|r |≥ρ
j(≠ i)
∫|r |≥ρ
= (12)
∫
1 dri − |ri|4
∫|r |≥ρ
∑
interior
∫
⎛ ⎞⎜ qi2 ⎛ 1 = ⎜ − 1⎟⎜ 8π ⎝ εw ⎠⎜⎜ ⎝
1 4π
where sj (≤1) is a scaling factor to consider possible overlaps between atoms.18,19 In the GB models based on eqs 9, 14, and 17,18,19,27,28,30 a cutoff distance, rmax, is introduced here to truncate the integral (summation) and make the computation faster, i.e.,
1 1 Ei(ri)·Di (ri)dri − Eivac(ri)·Divac (ri)dri 8π 8π ⎞ q2 ⎛ 1 1 ≃ i ⎜ − 1⎟ dri |ri|≥ ρi 8π ⎝ εw |r |i4 ⎠
∫
i
1 dri |ri|4
j
With CFA, the integral over the region inside the particles is canceled in eq 11, i.e., ΔGsol, i(d) =
i
≃
qi
ri , Ei(ri) |ri|3 ⎧ qi r (for region inside particles) ⎪ 3 i ⎪ |ri| ≃⎨ ⎪ qi ⎪ ε |r |3 ri (for region outside particles) ⎩ w i
∫|r |≥ρ interior
but the spherical symmetry of Ei and Di employed in eq 7 is no longer valid as the two particles come closer, so that PE must be solved to obtain ΔGsol,i(d). Since solving PE is timeconsuming, the so-called Coulomb field approximation (CFA)12,15,17 is used in the GB models,18,19,27,28,30 where the reaction field due to the polarization charges that arise on the dielectric boundaries is approximately considered on the assumption of the spherical symmetry of Ei and Di, i.e, Di (ri) ≃
(16)
(15) C
DOI: 10.1021/acs.jpcb.7b01438 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
ionic bonds between side-chain analogues.43) During the MD simulations, angular restraints were applied to prevent the two helices from tilting and staggering relative to each other and to prevent each helix from rotating around its helix axis, and distance restraints were applied to the α helical hydrogenbonding O−H pairs to prevent helix melting. In addition to PMF, potential energies as a function of d was examined (by single point potential energy calculations at the intervals of 0.25 Å) using GB(OBC), GB(Neck2),30 GB(NSR6)29 models, and the Poisson equation (PE). We also conducted potential energy calculation using GB(OBC) model with the partial-charge parameter set employed in the FF99 family, which is identical with that in FF9444 but is different from that in FF03.41 Potential energy calculations for the GB models were carried out using AmberTools15,45 and PE calculation was carried out using PEP.46 Larger System To Investigate Protein−Protein Interaction. We then studied a larger system where the electrostatic interaction substantially contributes to the protein−protein interaction: a system composed of monomeric kinesin (KIF1A) and tubulin dimer (α and β tubulin heterodimer). KIF1A is a molecular motor that moves along a microtubule which is an assembly of tubulin dimers. Binding between KIF1A and tubulin is driven largely by electrostatic interaction,47 where positively charged residues in kinesin are involved in the electrostatic interaction with negatively charged residues in tubulin.48 We used the KIF1A-tubulin(dimer) complex structure49 (PDB ID 2HXH). KIF1A (with 5723 atoms and net charge of +3e including bound Mg-ADP) and tubulin (with 13843 atoms and net charge of −51e including bound Mg-GTP and GDP) were manually dissociated from each other by 20 Å in the direction perpendicular to the microtubule surface. We used this dissociated state as the initial structure of the MD simulation to observe binding process. Missing regions (e.g., tubulin-binding loops of KIF1A and KIF1A-binding C-terminal tails of tubulin) were complemented by using MODELLER50 and thermally equilibrated before binding MD simulation. As in the two helices system mentioned above, we used the GB(OBC) model27 and examined two different cutoff distances for the energy density integration, rmax = 10 and 25 Å, while the nonbonded interactions (Lennard-Jones and Coulomb interactions) were calculated without cutoff. The ionic strength was set to 50 mM considering that 50 mM K-acetate was used in the in vitro experiment47 and the temperature was set to 310 K using Lagevin dynamics with the frictional constant of 1 ps −1. To examine whether KIF1A binds to the tubulin, 100-ns MD runs (4 independent runs) were conducted with rmax = 25 Å, and 15-ns MD runs (4 independent runs) were conducted with rmax = 10 Å. During the binding MD simulation, the main-chain atoms of the tubulin dimer (except for the KIF1A-binding Cterminal tails) were positionally restrained by weak harmonic restraints with the force constant of 0.02 kcal/mol/Å2. To confirm the result obtained by the GB(OBC) model, we also conducted the same KIF1A-tubulin binding simulation using the GB(Neck2) model.30 For comparison, we also conducted the binding simulation in explicit water: KIF1A and the tubulin dimer were immersed in a truncated octahedron water box containing 83666 TIP3P water molecules and 75 K+ and 27 Cl− ([K+] ∼ 50 mM), and 100-ns MD runs (16 independent runs) were conducted at 310 K and 1 bar.
Figure 1. (a) Model system to study the effect of rmax on the intermolecular electrostatic interaction: 15-residue alanine-based α helices (A7KA7 and A7EA7) where Lys and Glu are facing each other. (b) The effect of rmax on PMF as a function of d in the GB(OBC) model.27 PMFs for rmax= 10 Å (magenta), 15 Å (green), 25 Å (red), and 999 Å (blue) are shown. For comparison, PMF obtained by using explicit water as the solvent is also shown (black line). PMFs are vertically aligned so as to conform to the conventional Coulomb potential, −332/78.5/d = −4.23/d (kcal/mol) (dotted line), at larger d region (≥30 Å). Error bars represent the 95% confidence interval. (c) The effect of rmax on the potential energy profile as a function of d in the GB(OBC) model. The potential energy includes all of the nonbonded interaction energies (i.e., Coulomb, GB solvation, and Lennard-Jones energies). d was changed by translating the whole molecule in the direction perpendicular to the helix axes (see (a)) at 0.25 Å intervals.
distances, rmax = 10 Å, 15 Å, 25 Å, and ∞ (999 Å), were examined, while the nonbonded interactions (Lennard-Jones and Coulomb interactions) were calculated without cutoff. For comparison, we also calculated PMF in explicit water: the two helices were immersed in a truncated octahedron water box containing 22501 TIP3P water molecules.39 The periodic boundary condition was applied with the particle mesh Ewald method40 for long-range Coulomb interactions. PMF was calculated in the same way as was done using the GB model except that the umbrella sampling was enhanced by conducting eight independent 20-ns MD runs for each window (total 12 μs). Both the GB and the explicit water simulations were carried out using AMBER1235 with FF03 force-field41 and time step of 2 fs using H-bond SHAKE.42 (Note that the partialcharge parameter set of FF03 presents a relatively good agreement with the experimental results for the strength of D
DOI: 10.1021/acs.jpcb.7b01438 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
■
RESULTS AND DISCUSSION To investigate the influence of rmax (the energy density integration cutoff introduced in eq 18) on the intermolecular electrostatic interaction between protein molecules, we first studied the interaction between oppositely charged alaninebased α-helices as shown in Figure 1a: one containing positively charged lysine at the center (Ala7-Lys-Ala7) and the other containing negatively charged glutamic acid at the center (Ala7Glu-Ala7). First of all, we can see in Figure 1b that PMF is significantly affected by rmax over a wide range of d, up to as far as d ∼ 20 Å. As rmax is increased, PMF becomes higher and hence the intermolecular interaction becomes unstable, which is most prominently seen in the energy barrier with a peak at d ≈ 5.5 Å and a long tail toward larger d. For comparison, we calculated PMF by conducting MD simulation of the same two helix system immersed in explicit water (black line in Figure 1b). PMF resulting from the explicit water simulation shows that the electrostatic interaction between the two helices (i.e, the ionic bond between Lys and Glu) is gradually stabilized as the two helices come closer, manifesting the existence of a longrange attractive force, with a major energy barrier at d ≈ 4.3 Å. This overall profile of PMF observed in the explicit water solvation is well captured by the GB(OBC) model when rmax is relatively small (rmax= 10 and 15 Å), even though the energy barrier in the GB(OBC) model appears in a coarse-grained manner (i.e., broad and unimodal) covering major and minor energy barriers observed in PMF of explicit water (minor energy barriers are recognized at d ≈ 6.8 and 9.3 Å). However, as rmax is increased up to 25 Å or greater, the GB(OBC) model obviously overdestabilizes the energy barrier, so that the longrange attractive interaction (d ≳ 10 Å) turns into repulsion (the slope of PMF is reversed). This overdestabilization of PMF is inherited from the potential energy as shown in Figure 1c, where the energy barrier becomes so high that the long-range attraction (d ≳ 10 Å) is disappeared as rmax is increased. Note that the same result as in Figure 1c was obtained when the partial-charge parameter set of the FF99 family44 was used (Figure S1). The same result was also obtained when direct numerical integration of 1/r4 in eq 14 was conducted without employing the approximation by the pairwise function as in eq 1718,19,27 (Figure S2). The destabilization of the intermolecular ionic bond with increasing rmax is also observed in a newer GB model,30 referred to as GB(Neck2), as shown in Figure 2. In the GB(Neck2) model, the Coulomb field approximation (CFA), as described
in eq 12, is also employed, but the dielectric boundary between solute and solvent is more accurately treated than in the GB(OBC) model and the key parameters used in the GB(OBC) model are reparametrized by reference to the numerical solution of PE. The destabilization with increasing rmax is observed in the GB(Neck2) model as well, but the extent of the destabilization is smaller than that observed in the GB(OBC) model. In addition, the GB(Neck2) model captures the overall potential energy profile of PE (Figure 3), preserving
Figure 3. Effect of the Coulomb field approximation (CFA) on the potential energy profile as a function of d (the system is identical with that in Figure 1). The potential energy resulting from the exact (numerical) solution of PE (black) and that calculated by using the GB(NSR6) model29(cyan) are shown. CFA is not employed in the numerical solution of PE nor in the GB(NSR6) model. Potential energy profiles of the GB(OBC) model already shown in Figure 1 are included for comparison [rmax= 10 Å (magenta), 25 Å (red), and 999 Å (blue)].
the long-range attraction (even with large rmax) and the location of the major energy barrier at (d ≈ 4.3 Å). Note that PE appears to overestimate the energy barrier and underestimate the stability of the ionic bond, indicating that there is room for PE to be improved by taking the variation of the intermolecular local dielectric environment into account (e.g., charges would remain partially solvated upon ionic bond formation). The reduction of the ionic bond stability is also seen in the GB(Neck2) model, compared to the result of the GB(OBC) model, which is anticipated because the mBondi3 parameter set was introduced in the GB(Neck2) model to reduce the ionic bond stability.30 The physical mechanism of how the intrinsic radius ρ (i.e., the lower limit of the energy density integration) affects the ionic bond stability will be examined in detail in our subsequent work (where an effective way to consider the variation of the local dielectric environment is offered to improve the ionic bond stability). Considering that the long-range attraction is seen in the numerical solution of PE where CFA is not employed, it is supposed that the overdestabilization that changes the longrange attraction into repulsion as found in the GB(OBC) model is due to the employment of CFA. It has been recognized that employing CFA leads to underestimation of the Born energy [or overestimation of the effective Born radius (see eq 10)],14,15,26 so there has been a lot of effort to correct this flaw of CFA.20,23−26,28−30,51 In particular, Onufriev and coworkers proposed a GB model, referred to as the GB(NSR6) model,29 that does not employ CFA but instead utilizes the Kirkwood-theory-based formulation by Grycuk.26 Here we examined the potential energy profile of the GB(NSR6) model. As shown in Figure 3, the GB(NSR6) model can well reproduce the potential energy profile of PE, particularly the
Figure 2. Effect of rmax on the potential energy profile as a function of d (the system is identical with that in Figure 1) in the GB(Neck2) model with the mBondi3 parameter set.30 E
DOI: 10.1021/acs.jpcb.7b01438 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
desolvation penalty exhibited by PE and the GB(NSR6) model in which the desolvation penalty almost disappears for d ≳ 10 Å. This disappearance of the desolvation penalty for d ≳ 10 Å is essential for the emergence of the long-range attraction in the GB(OBC) model (rmax = 10 Å) as well as in PE and the GB(NSR6) model (Figure 3). We note that the same physical mechanism as found for the ionic bond formation between oppositely charged residues underlies the hydrogen bond formation between polar residues where the electrostatic bond between oppositely signed polarized charges (partial atomic charges) is involved. Consider a general case of association of two solute molecules, A and B, which interact electrostatically with each other in water. Then the desolvation penalty of solute A upon association is determined by ∑i∈solute A[ΔGsol,i(d)−ΔGsol,i(∞)]. Eq 13 indicates that (0 ≥ ) ΔGsol,i(d) ≥ ΔGsol,i(∞) because of the occluded volume for the spatial integration of the energy density (the “solvation energy density”52) due to solute B (see the second line of eq 13). Thus, the overestimation of the desolvation penalty is expected to be reduced for the solutes with smaller volume. This expectation is proven to be the case as shown in Figure 5, where the overdestabilization (even with rmax = 999 Å) disappears for the smaller system composed of lysine and glutamic acid monomers.
long-range attraction. Thus, from the results of the GB(Neck2) and GB(NSR6) models as well as PE, we can confirm that the overdestabilization caused by setting rmax to large values (e.g., 25 Å) in the GB(OBC) model is due to the employment of CFA. Looking at this result from the opposite standpoint, we can point out that the flaw of CFA in the GB models that should become fatal in studying protein−protein interaction can be greatly remedied if we use relatively small rmax (e.g., 10 or 15 Å). Irrespective of the effect of rmax, the attractive electrostatic interaction becomes destabilized as the two molecules come closer to form ionic bond, because the electrostatic self-energy, which is always positive, is raised (or the Born energy, which is always negative, is reduced in magnitude) as the two molecules come closer and the effective dielectric constant around each charge becomes lower than εw. This destabilization is generally called the “desolvation penalty” that accompanies the ionic bond formation. Thus, the overdestabilization observed in the GB(OBC) model is supposed to be attributed to the overestimation of the desolvation penalty. So we calculated the desolvation penalty as a function of d, which is defined by ∑i[ΔGsol,i(d) − ΔGsol,i(∞)], where ΔGsol,i(d) is the Born energy for the ith atom in the molecules (see eq 10). As is clearly seen in Figure 4a, the desolvation penalty is actually
Figure 5. Effect of rmax on the potential energy profile as a function of d for the system composed of Lys and Glu (monomers) as shown in (a) (N and C termini of each monomer are neutralized by acetylation and N-methylamidation, respectively). (b) Potential energies are calculated using the GB(OBC) model27 with rmax= 10 Å (magenta), 15 Å (green), 25 Å (red), and 999 Å (blue).
Figure 4. Effect of rmax on the desolvation penalty (see the text) and the interatomic electrostatic interaction energy (i.e., the sum of interatomic Coulomb energy and interatomic GB solvation energy12), which were calculated by using the effective Born radii obtained by eq 10.
For the interaction of proteins with larger size, it is easily foreseeable that the GB models that employ CFA overdestabilize the protein−protein interaction if rmax is set to a relatively large value (e.g., 20 or 25 Å). Therefore, to examine how significantly rmax affects the protein−protein interaction, we conducted MD simulation using the GB(OBC) model for a larger system composed of molecular motor KIF1A (monomeric kinesin) and αβ tubulin heterodimer (building unit of microtubule). It has been experimentally demonstrated that KIF1A binds to the tubulin dimer at the junction region between the α and β tubulins49 and the binding is driven by the electrostatic interaction;47 the binding interface of KIF1A is
overestimated in the GB(OBC) model when rmax = 999 Å; furthermore, this overestimation remains substantial even at large intermolecular separation distances (up to d ∼ 30 Å). This overestimated desolvation penalty is largely compensated for by the interatomic electrostatic interaction energy (Figure 4b), but the compensation is not enough to recover the net long-range attraction. In contrast, the GB(OBC) model with small rmax (rmax = 10 Å) well captures the overall profile of the F
DOI: 10.1021/acs.jpcb.7b01438 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
ionic bonds to the stability of the folded state; in fact, the rootmean-square deviation (RMSD) from the initial structure for KIF1A with rmax = 10 Å remained the same (∼4 Å) as that with rmax = 25 Å. In addition, it is noteworthy to comment here that the binding dynamics can be made faster in the GB model than in the explicit water by using a low frictional constant. In fact, reflecting that we used a low frictional constant of 1 ps−1, which is about a magnitude lower than that of bulk water, the binding of KIF1A to the tubulin dimer observed in the MD simulation using the GB model (rmax = 10 Å) was about 10 times faster than that observed in the MD simulation using explicit water (Figure 6b). It should also be mentioned that we do not have to worry about the detrimental effect of the discontinuity arising from adopting a small cutoff (which is usually prominent in the Coulombic force calculation) because the energy density as a function of r rapidly deceases (in proportion to r−4 as in eq 18) and moreover the pairwise function to be used for the energy density integration is designed to avoid the discontinuity.17−19 Thus, our study suggests that once the overall thermodynamic aspect of the binding (e.g., the long-range electrostatic attraction) can be reproduced properly (at least semiquantitatively), the potential advantage of the GB model in studying the dynamical aspect of the protein−protein interaction (in addition to the obvious advantage in reducing the computational cost) would come to stand out.
positively charged and that of the tubulin dimer is negatively charged, as is seen in the characteristic electrostatic potentials on the surface (Figure 6a). Figure 6b clearly demonstrates that
■
CONCLUSIONS It is important to keep in mind that the cutoff for the energy density integral, rmax, exerts a significant effect on the result of MD simulation using a GB model where CFA is employed, in the way that larger rmax overdestabilizes the electrostatic interaction between protein molecules. We proposed an effective way to remedy this problem: setting rmax (e.g., the parameter RGBMAX in AMBER) to a relatively small value (10−15 Å). With this quite simple treatment, the GB model with CFA can be used as a powerful and useful means to study protein−protein interaction, particularly on the dynamical aspect, not only requiring less computational cost but also reproducing semiquantitatively the results of more accurate yet more computationally demanding MD simulation using explicit water.
Figure 6. Effect of rmax on the binding dynamics between KIF1A and the tubulin dimer. (a) The system is composed of KIF1A (right) and the tubulin dimer (left). Electrostatic potentials are mapped on the molecular surfaces (red: negative, blue: positive). (b) MD trajectories of the intermolecular distance between KIF1A and the tubulin dimer, d [the distance between the mass centers as shown in (a)], using the GB(OBC) model27 with rmax = 10 Å (magenta) and 25 Å (red). For comparison, MD trajectory using explicit water as the solvent is also shown (black). Each trajectory is the averaged trajectory of multiple independent MD runs. The error bars represent the 95% confidence interval.
rmax exerts a profound effect on the binding of KIF1A to the tubulin dimer; KIF1A does not bind to the tubulin dimer when rmax= 25 Å, while binding is observed when rmax= 10 Å. A similar result is obtained when using the GB(Neck2) model (see Figure S3), reinforcing the profound effect of rmax on the binding dynamics. Note that binding of KIF1A to the tubulin dimer is also observed in the MD simulation using the explicit water (black line in Figure 6b). (If we employed the GB(NSR6) model, though MD simulation using the GB(NSR6) model has not yet been available in AMBER, we would observe binding of KIF1A to the tubulin dimer, but the binding would be incomplete due to the overestimated energy barrier and the underestimated ionic bond strength as seen in Figure 3.) These results, taken together, have exposed the potential pitfall of the GB model with CFA, and at the same time, affirm that the GB models remain useful and powerful if the basic physics underlying the GB models is properly understood, which provides us with the way to avoid the pitfall. Here, the way is to set rmax to a relatively small value (especially for the GB(OBC) model). Note that reducing rmax also has an effect to increase the stability of intramolecular ionic bonds, but the effect would be limited for a relatively large protein, such as KIF1A, because of the limited contribution of intramolecular
■
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.7b01438. The effect of rmax on the potential energy profile using the GB(OBC) model with the partial-charge parameter set employed in the FF99 family; effect of rmax on the potential energy profile using the GB model without employing the approximation by the pairwise function; effect of rmax on the binding dynamics between KIF1A and the tubulin dimer using the GB(Neck2) model (PDF)
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. ORCID
Yukinobu Mizuhara: 0000-0001-9819-9617 Dan Parkin: 0000-0003-4897-0412 G
DOI: 10.1021/acs.jpcb.7b01438 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
(18) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. Pairwise Solute Descreening of Solute Charges from a Dielectric Medium. Chem. Phys. Lett. 1995, 246, 122−129. (19) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. Parametrized Models of Aqueous Free Energies of Solvation Based on Pairwise Descreening of Solute Atomic Charges from a Dielectric Medium. J. Phys. Chem. 1996, 100, 19824−19839. (20) Onufriev, A.; Case, D. A.; Bashford, D. Effective Born Radii in the Generalized Born Approximation: The Importance of Being Perfect. J. Comput. Chem. 2002, 23, 1297−1304. (21) Qiu, D.; Shenkin, P. S.; Hollinger, F. P.; Still, W. C. The GB/SA Continuum Model for Solvation. A Fast Analytical Method for the Calculation of Approximate Born Radii. J. Phys. Chem. A 1997, 101, 3005−3014. (22) Onufriev, A.; Bashford, D.; Case, D. A. Modification of the Generalized Born Model Suitable for Macromolecules. J. Phys. Chem. B 2000, 104, 3712−3720. (23) Lee, M. S.; Salsbury, F. R.; Brooks, C. L., III Novel Generalized Born Methods. J. Chem. Phys. 2002, 116, 10606−10614. (24) Lee, M. S.; Feig, M.; Salsbury, F. R.; Brooks, C. L., III New Analytic Approximation to the Standard Molecular Volume Definition and Its Application to Generalized Born Calculations. J. Comput. Chem. 2003, 24, 1348−1356. (25) Im, W.; Lee, M. S.; Brooks, C. L., III Generalized Born Model with a Simple Smoothing Function. J. Comput. Chem. 2003, 24, 1691− 1702. (26) Grycuk, T. Deficiency of the Coulomb-field Approximation in the Generalized Born Model: An Improved Formula for Born Radii Evaluation. J. Chem. Phys. 2003, 119, 4817−4826. (27) Onufriev, A.; Bashford, D.; Case, D. A. Exploring Protein Native States and Large-scale Conformational Changes with a Modified Generalized Born Model. Proteins: Struct., Funct., Bioinf. 2004, 55, 383−394. (28) Mongan, J.; Simmerling, C.; McCammon, J. A.; Case, D. A.; Onufriev, A. Generalized Born Model with a Simple, Robust Molecular Volume Correction. J. Chem. Theory Comput. 2007, 3, 156−169. (29) Aguilar, B.; Shadrach, R.; Onufriev, A. V. Reducing the Secondary Structure Bias in the Generalized Born Model via R6 Effective Radii. J. Chem. Theory Comput. 2010, 6, 3613−3630. (30) Nguyen, H.; Roe, D. R.; Simmerling, C. Improved Generalized Born Solvent Model Parameters for Protein Simulations. J. Chem. Theory Comput. 2013, 9, 2020−2034. (31) Tsui, V.; Case, D. A. Molecular Dynamics Simulations of Nucleic Acids with a Generalized Born Solvation Model. J. Am. Chem. Soc. 2000, 122, 2489−2498. (32) Tsui, V.; Case, D. A. Theory and Applications of The Generalized Born Solvation Model in Macromolecular Simulations. Biopolymers 2001, 56, 275−291. (33) Geney, R.; Layten, M.; Gomperts, R.; Hornak, V.; Simmerling, C. Investigation of Salt Bridge Stability in a Generalized Born Solvent Model. J. Chem. Theory Comput. 2006, 2, 115−127. (34) Dominy, B. N.; Brooks, C. L., III Development of a Generalized Born Model Parametrization for Proteins and Nucleic Acids. J. Phys. Chem. B 1999, 103, 3765−3773. (35) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M. et al. AMBER 12 Reference Manual; University of California: San Francisco, 2012. (36) Born, M. Volumen und hydratationswärme der Ionen. Z. Phys. 1920, 1, 45−48. (37) Strelkov, S. V.; Burkhard, P. Analysis of α-Helical Coiled Coils with the Program TWISTER Reveals a Structural Mechanism for Stutter Compensation. J. Struct. Biol. 2002, 137, 54−64. (38) Kumar, S.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A.; Rosenberg, J. M. The Weighted Histogram Analysis Method for Freeenergy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011−1021.
Jun Ohnuki: 0000-0002-7351-3427 Mitsunori Takano: 0000-0001-8855-008X Author Contributions ∥
Y.M. and D.P. contributed equally to this work.
Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS We thank Alexey Onufriev for valuable comments and discussion at Protein Electrostatics Berlin 2016. This work was supported by Grants-in-Aid for Scientific Research and Top Global University Project from MEXT (Japan). D.P. acknowledges Leading Graduate Program in Science and Engineering in Waseda University (“Energy Next” program) funded by MEXT (Japan).
■
REFERENCES
(1) Perutz, M. F. Electrostatic Effects in Proteins. Science 1978, 201, 1187−1191. (2) Sheinerman, F. B.; Norel, R.; Honig, B. Electrostatic Aspects of Protein-Protein Interactions. Curr. Opin. Struct. Biol. 2000, 10, 153− 159. (3) Schreiber, G.; Haran, G.; Zhou, H.-X. Fundamental Aspects of Protein-Protein Association Kinetics. Chem. Rev. 2009, 109, 839−860. (4) Warshel, A. Calculations of Enzymic Reactions: Calculations of pKa, Proton Transfer Reactions, and General Acid Catalysis Reactions in Enzymes. Biochemistry 1981, 20, 3167−3177. (5) Davis, M. E.; McCammon, J. A. Electrostatics in Biomolecular Structure and Dynamics. Chem. Rev. 1990, 90, 509−521. (6) Honig, B.; Nicholls, A. Classical Electrostatics in Biology and Chemistry. Science 1995, 268, 1144−1149. (7) Schaefer, M.; Bartels, C.; Karplus, M. Solution Conformations and Thermodynamics of Structured Peptides: Molecular Dynamics Simulation with an Implicit Solvation Model. J. Mol. Biol. 1998, 284, 835−848. (8) Simonson, T. Macromolecular Electrostatics: Continuum Models and Their Growing Pains. Curr. Opin. Struct. Biol. 2001, 11, 243−252. (9) Warshel, A.; Sharma, P. K.; Kato, M.; Parson, W. W. Modeling Electrostatic Effects in Proteins. Biochim. Biophys. Acta, Proteins Proteomics 2006, 1764, 1647−1676. (10) Nakamura, H. Roles of Electrostatic Interaction in Proteins. Q. Rev. Biophys. 1996, 29, 1−90. (11) Baker, N. A.; Sept, D.; Joseph, S.; Holst, M. J.; McCammon, J. A. Electrostatics of Nanosystems: Application to Microtubules and The Ribosome. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 10037−10041. (12) Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. Semianalytical Treatment of Solvation for Molecular Mechanics and Dynamics. J. Am. Chem. Soc. 1990, 112, 6127−6129. (13) Cramer, C. J.; Truhlar, D. G. Implicit Solvation Models: Equilibria, Structure, Spectra, and Dynamics. Chem. Rev. 1999, 99, 2161−2200. (14) Bashford, D.; Case, D. A. Generalized Born Models of Macromolecular Solvation Effects. Annu. Rev. Phys. Chem. 2000, 51, 129−152. (15) Schaefer, M.; Karplus, M. A Comprehensive Analytical Treatment of Continuum Electrostatics. J. Phys. Chem. 1996, 100, 1578−1599. (16) Feig, M.; Onufriev, A.; Lee, M. S.; Im, W.; Case, D. A.; Brooks, C. L., III Performance Comparison of Generalized Born and Poisson Methods in the Calculation of Electrostatic Solvation Energies for Protein Structures. J. Comput. Chem. 2004, 25, 265−284. (17) Schaefer, M.; Froemmel, C. A Precise Analytical Method for Calculating the Electrostatic Energy of Macromolecules in Aqueous Solution. J. Mol. Biol. 1990, 216, 1045−1066. H
DOI: 10.1021/acs.jpcb.7b01438 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B (39) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (40) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577−8593. (41) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; et al. 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. (42) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. 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. (43) Debiec, K. T.; Gronenborn, A. M.; Chong, L. T. Evaluating the Strength of Salt Bridges: A Comparison of Current Biomolecular Force Fields. J. Phys. Chem. B 2014, 118, 6561−6569. (44) 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. (45) Case, D. A.; Berryman, J. T.; Betz, R. M.; Cerutti, D. S.; Cheatham, T. E., III; Darden, T. A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W. et al. Amber 2015; University of California: San Francisco, 2015. (46) Beroza, P.; Fredkin, D. R. Calculation of Amino Acid Pkas in a Protein from a Continuum Electrostatic Model: Method and Sensitivity Analysis. J. Comput. Chem. 1996, 17, 1229−1244. (47) Okada, Y.; Hirokawa, N. A Processive Single-Headed Motor: Kinesin Superfamily Protein KIF1A. Science 1999, 283, 1152−1157. (48) Woehlke, G.; Ruby, A. K.; Hart, C. L.; Ly, B.; Hom-Booher, N.; Vale, R. D. Microtubule Interaction Site of the Kinesin Motor. Cell 1997, 90, 207−216. (49) Kikkawa, M.; Hirokawa, N. High-Resolution Cryo-EM Maps Show the Nucleotide Binding Pocket of KIF1A in Open and Closed Conformations. EMBO J. 2006, 25, 4187−4194. (50) Eswar, N.; Webb, B.; Marti-Renom, M. A.; Madhusudhan, M. S.; Eramian, D.; Shen, M.-Y.; Pieper, U.; Šali, A. Comparative Protein Structure Modeling Using MODELLER. Curr. Protoc. Protein Sci.; 2007; Chapter 2, Unit 2.9. (51) Mongan, J.; Svrcek-Seiler, W. A.; Onufriev, A. Analysis of Integral Expressions for Effective Born Radii. J. Chem. Phys. 2007, 127, 185101. (52) Arora, N.; Bashford, D. Solvation Energy Density Occlusion Approximation for Evaluation of Desolvation Penalties in Biomolecular Interactions. Proteins: Struct., Funct., Genet. 2001, 43, 12−27.
I
DOI: 10.1021/acs.jpcb.7b01438 J. Phys. Chem. B XXXX, XXX, XXX−XXX
