9114
J. Phys. Chem. B 2007, 111, 9114-9124
In Silico Screening of Mutational Effects on Transmembrane Helix Dimerization: Insights from Rigid-Body Docking and Molecular Dynamics Simulations Daniele Dell’Orco,†,‡ Pier Giuseppe De Benedetti,† and Francesca Fanelli*,†,‡ Department of Chemistry and Dulbecco Telethon Institute, UniVersity of Modena and Reggio Emilia, Via Campi 183, 41100 Modena, Italy ReceiVed: February 18, 2007; In Final Form: April 29, 2007
In this study, a docking-based protocol has been probed for its ability to predict the effects of 32 single and double mutations on glycophorin A (GpA) homodimerization. Rigid-body docking simulations have been paralleled by molecular dynamics (MD) simulations in implicit membrane. The rigid-body docking-based approach proved effective in reconstituting the native architecture of the GpA dimer for the wild type and the wild-type-like mutants. The good correlative models between the intermolecular interaction descriptors derived both by rigid-body docking simulations and by MD simulations and experimental relative free energies support the assumption that the mutation-induced changes in the association free energy of GpA helices are essentially ascribed to differences in the packing interactions, whereas almost all the variations in the entropic contributions to the association free energy are constant and/or negligible. The MD-based models achieved provide insights into the structural determinants for disruptive and restoring mutational effects. The computational approaches presented in this study are fast and effective, and constitute simple and promising tools for in silico screening of mutational effects on the association properties of integral membrane proteins.
Introduction The association of helices within the apolar milieu of lipidic bilayers is an attractive issue for several reasons. One reason is that helix-helix interactions are considered the last crucial step of the membrane protein folding mechanism1-4 according to the widely accepted two-stage folding model.5,6 Another reason is that increasing evidence is emerging about the importance of membrane protein oligomerization in the regulation of a variety of biochemical processes. Homooligomeric transmembrane helices are particularly attractive systems for the study of membrane protein folding and stability. Indeed, the noncovalent assembly of membrane peptides is likely to mimic the folding of larger native proteins, as demonstrated by various fragment-complementation studies (see refs 1 and 4 and references therein). Among such systems, glycophorin A (GpA) homodimer from erythrocyte membrane represents one of the most studied cases to date. High-resolution structural data from solution NMR in dodecylphosphocholine (DPC) micelles7 and from solid-state NMR in dimyristoylphosphatidylcholine (DMPC) bilayers,8,9 together with the plenty of mutagenesis work aimed at characterizing the thermodynamics of association,10-16 make this system intriguing for computational studies of helix-helix associations. Recently, we set up a computational protocol based upon rigid-body docking for fast predictions of mutational effects on the thermodynamics of association of a number of proteinprotein complexes. So far, the method has been used on a total of 67 heterogeneous variants of water-soluble protein-protein complexes,17,18 and on 10 mutants of the MEF2A transcription factor in complex with DNA.19 The fundamental requirement of the approach is an accurate structural model of the complex
between the wild type forms of the interacting macromolecules. The variants (i.e., mutations or deletions) of either partner or both partners can be achieved by molecular modeling. The basic assumption is that the architecture of the mutated complexes is almost the same as that of the wild type and no major conformational changes occur upon binding. The good performance of the method is related to the scoring function of the protein-protein docking algorithm ZDOCK,20-22 which effectively accounts for electrostatic and shape complementarities as well as desolvation. We have also shown that neglecting the desolvation term in docking simulations leads to effective quaternary structure predictions of oligomeric membrane proteins, independent of the oligomeric order and of the extension of the intracellular and extracellular hydrophilic domains.23 In this study, we test the same docking-based protocol in its ability to predict the effects of 32 single and double mutations on GpA homodimerization, hence reproducing in silico a set of in vitro alanine-scanning mutagenesis experiments done in recent years.11-14 Rigid-body docking simulations have been paralleled by molecular dynamics (MD) simulations in implicit membrane.24 The input structures for such simulations were the wild type and 32 mutated forms of the solid-state NMR structure.8,9 The results of this study corroborate the effectiveness of our rigid-body docking-based protocol in fast estimation of mutational effects on the thermodynamics of protein-protein association, independent of the environment, i.e., water or membrane. Moreover, combining the results of rigid-body docking with MD simulations add interpretative power to the correlative models. Methods
* Corresponding author. Telephone: +39 059 2055114. Fax: +39 059 373543. E-mail:
[email protected]. † University of Modena and Reggio Emilia. ‡ Dulbecco Telethon Institute.
Rigid-Body Docking Simulations of GpA Dimerization. Rigid-body docking simulations of GpA homodimerization were carried out by means of the rigid-body docking algorithm
10.1021/jp071383r CCC: $37.00 © 2007 American Chemical Society Published on Web 06/28/2007
Mutational Effects on GpA Dimerization
J. Phys. Chem. B, Vol. 111, No. 30, 2007 9115
TABLE 1: GpA Dimerization: In Vitro Thermodynamic Data and Structural Parameters from Rigid-Body Docking and Molecular Dynamics Simulations variant
∆∆G°expa ∆G°coupb (kcal/mol) (kcal/mol)
wt 0.0 ( 0.1 L75A 1.3 ( 0.1 I76A 1.8 ( 0.2 I77A 0.0 ( 0.1 F78A 0.0 ( 0.1 G79A 1.7 ( 0.1 V80A 0.4 ( 0.1 M81A -0.1 ( 0.1 G83A 3.1 ( 0.1 V84A 1.0 ( 0.1 I85A -0.2 ( 0.1 G86A -0.4 ( 0.1 T87A 0.9 ( 0.1 L75A/I76A 2.4 ( 0.1 L75A/G79A 2.3 ( 0.1 L75A/V80A 1.3 ( 0.3 L75A/G83A 3.6 ( 0.2 L75A/V84A 0.1 ( 0.1 I76A/G79A 2.9 ( 0.2 I76A/V80A 2.5 ( 0.3 I76A/G83A 3.4 ( 0.2 I76A/V84A 1.2 ( 0.1 I76A/T87A 2.8 ( 0.2 G79A/V80A 2.4 ( 0.1 G79A/G83A 3.6 ( 0.3 G79A/V84A 3.0 ( 0.2 G79A/T87A 2.8 ( 0.2 V80A/G83A 3.5 ( 0.1 V80A/V84A 2.0 ( 0.2 V80A/T87A 3.7 ( 0.2 G83A/V84A 3.2 ( 0.1 G83A/T87A 3.3 ( 0.1 V84A/T87A 1.0 ( 0.1
0 0 0 0 0 0 0 0 0 0 0 0 0 -0.7 -0.7 -0.4 -0.8 -2.2 -0.6 0.3 -1.6 -1.6 0.1 0.3 -1.2 0.3 0.2 0 0.6 2.4 -0.9 -0.7 -0.9
∆ZDc
cluster populationd
best ranke
ZDbest f
0.00 ( 0.05 0.04 ( 0.05 -0.04 ( 0.05 -0.20 ( 0.05 -0.14 ( 0.05 0.35 ( 0.05 0.05 ( 0.04 -0.11 ( 0.05 0.60 ( 0.04 0.04 ( 0.04 -0.07 ( 0.05 -0.10 ( 0.04 0.21 ( 0.05 0.08 ( 0.04 0.36 ( 0.04 -0.01 ( 0.04 0.55 ( 0.05 0.12 ( 0.04 0.36 ( 0.04 0.04 ( 0.05 0.77 ( 0.04 0.07 ( 0.04 0.50 ( 0.04 0.56 ( 0.04 0.62 ( 0.04 0.48 ( 0.04 0.25 ( 0.06 0.73 ( 0.03 -0.06 ( 0.05 0.27 ( 0.04 0.55 ( 0.06 0.78 ( 0.07 0.25 ( 0.06
58, 61, 59 (179) 67, 70, 70 (207) 57, 52, 56 (165) 59, 57, 56 (172) 62, 59, 62 (183) 41, 43, 35 (119) 49, 47, 50 (146) 55, 66, 60 (181) 33, 30, 34 (97) 47, 46, 50 (143) 53, 59, 51 (163) 63, 66, 63 (192) 45, 36, 32 (113) 49, 49, 48 (146) 45, 45, 40 (130) 63, 52, 66 (181) 22, 21, 23 (66) 66, 63, 61 (190) 40, 42, 43 (125) 49, 45, 48 (142) 24, 21, 26 (71) 51, 47, 51 (149) 40, 36, 40 (116) 32, 34, 29 (95) 28, 25, 23 (76) 29, 34, 33 (96) 30, 27, 22 (79) 20, 18, 20 (58) 47, 45, 45 (137) 35, 34, 38 (107) 13, 16, 15 (44) 14, 9, 12 (35) 38, 34, 33 (105)
3, 21, 7 1, 1, 3 6, 4, 9 4, 8, 2 2, 3, 2 13, 41, 34 6, 17, 5 4, 4, 1 172, 51, 153 6, 2, 14 9, 24, 8 1, 2, 3 15, 17, 66 9, 15, 2 67, 49, 27 5, 5, 10 57, 57, 128 7, 11, 8 26, 19, 29 14, 24, 13 97, 132, 280 6, 11, 15 34, 75, 107 62, 108, 62 85, 101, 168 14, 29, 45 20, 15, 35 884, 413, 499 3, 27, 5 17, 61, 54 116, 115, 90 47, 240, 269 7, 52, 21
8.08, 7.52, 7.84 8.32, 8.46, 8.16 8.30, 8.04, 8.14 8.48, 7.96, 8.54 8.22, 8.02, 8.86 7.70, 7.14, 7.14 8.16, 7.74, 8.40 8.02, 8.14, 8.48 6.60, 7.08, 6.70 7.94, 8.24, 7.62 7.96, 7.68, 8.00 8.42, 8.00, 7.96 7.98, 7.64, 7.14 7.96, 7.70, 8.10 7.06, 7.26, 7.42 8.20, 8.28, 7.86 6.96, 6.92, 6.72 7.94, 7.86, 7.90 7.24, 7.30, 7.22 7.94, 7.80, 7.96 6.76, 6.58, 6.32 7.92, 7.60, 7.66 7.26, 6.90, 6.80 6.96, 6.70, 6.94 6.80, 6.66, 6.54 7.52, 7.18, 7.04 7.66, 7.64, 7.40 6.00, 6.30, 6.22 8.42, 7.62, 8.12 7.74, 7.24, 7.34 6.76, 6.78, 6.88 7.16, 6.46, 6.44 8.16, 7.40, 7.76
∆IE1g CR-rmsd1h ∆IE2i CR-rmsd2j (kcal/mol) (Å) (kcal/mol) (Å) 0.00 5.27 9.15 0.71 1.29 12.23 15.53 3.18 14.62 2.95 2.52 -0.71 0.00 14 18.15 14.78 21.69 13.05 17.2 15.16 18.9 13.59 13.34 17.19 14.37 17.7 16.13 14.96 11.45 9.64 14.52 15.69 6.06
0.82 1.51 1.5 1.66 0.71 1.09 1.08 1.48 2.22 1.57 1.5 1.65 1.26 1.6 1.58 1.07 1.75 1.21 1.53 1.33 2.5 1.5 1.29 1.47 2.02 1.39 1.55 1.58 1.57 1.27 1.61 3.16 1.38
0.00 3.20 4.44 -0.58 0.19 7.52 2.33 -1.47 15.42 -3.14 0.00 -1.58 0.80 5.71 9.42 6.88 22.49 1.08 8.63 7.89 19.70 4.01 14.14 9.07 15.17 5.85 16.93 15.76 0.11 10.44 13.41 16.49 6.86
0.66 0.66 0.97 1.03 1.39 1.14 1.62 1.48 1.20 1.33 1.05 1.06 1.53 1.46 1.46 1.22 1.49 1.36 0.93 1.36 1.57
a Experimental changes of binding free energy of GpA dimerization upon mutation, defined by the relationship ∆∆G° ) ∆G° 11-14 mut - ∆G°wt. Here and in the following definitions, “mut” refers to the mutated GpA form, while “wt” refers to the wild type dimer. b Free energy of coupling, as defined in eq 2 in the text. c Relative average docking scores for the nativelike complexes obtained from three independent docking runs, defined by the relationship ∆ZD ) ZD-smut - ZD-swt. The standard errors are reported. d Number of nativelike solutions obtained from each independent docking run, out of 4000 overall configurations saved for each run. Data referred to each run are separated by commas. In parentheses is reported the population of each nativelike ensemble used in correlation analyses, out of 12 000 overall retained configurations. e Rank order of the bestscored nativelike solution resulting from each run out of 4000 overall retained configurations. Commas separate data from each run. f ZDOCK score of the best-scored nativelike solution resulting from each run out of 4000 overall retained configurations. Commas separate data from each run. g Changes in the interaction energy with respect to the wild type defined by the relationship ∆IE1 ) IE1mut - IE1wt. The interaction energies have been computed on the average minimized structures of the wild type and mutated forms derived by MD simulations in the absence of interhelical distance restraints between the hydroxy oxygen atom of T87 from one helix and the backbone oxygen atom of V84 from the other helix. h RCarbon root-mean-square deviation between the average minimized structure obtained by 1 ns MD simulation without interhelical distance restraints and the experimental solid-state NMR structure.9 i Changes in the interaction energy with respect to the wild type defined by the relationship ∆IE2 ) IE2mut - IE2wt. For the wild type and all mutants except for the G83A and T87A series, the interaction energies have been computed on the average minimized structures derived by MD simulations in the presence of interhelical distance restraints between the hydroxy oxygen atom of T87 from one helix and the backbone oxygen atom of V84 from the other helix. j R-Carbon root-mean-square deviation between the average minimized structure obtained by 1 ns MD simulation with interhelical distance restraints and the experimental solid-state NMR structure.9
ZDOCK2.1.20,21 The grid-based pairwise scoring function22 implemented in this version of the software accounts for shape and electrostatic complementarities neglecting the desolvation term, which is suitable for aqueous environments. The algorithm performs an efficient fast Fourier transform sampling of the entire six degrees of freedom in the translational and rotational space, and finally provides a score for each docking solution, which is representative of the complementarity of the complex. Docking simulations were carried out on the wild type and 32 single and double mutant forms (Table 1). Truncated forms of segments A and B corresponding to the transmembrane portions (i.e., I76-I95 segments) were considered for docking simulations. Both the solid-state NMR structure in DMPC bilayers8,9 and the minimized structure averaged over the 20 structures resolved by solution NMR in DPC micelles7 were initially employed to extract the I76-I95 A and B segments. However, the solid-state NMR structure was finally selected
for the gross of rigid-body docking simulations and analysis as well as for MD simulations (see the following section). The rationale behind the choice of the solid-state NMR structure instead of the solution structure was the more extended interface in the former compared to the latter and the consequent best performance of the docking algorithm. A posteriori, this choice is supported by the results of MD simulations in DPC micelle and DMPC bilayer, which converged in a helix-helix packing architecture similar to that of the NMR structure determined in bilayer.25 Single and double Ala mutations (Table 1) in both monomers were performed by means of the Protein Design module within the QUANTA2005 package. Three independent sets of docking runs were performed for each complex, i.e., one starting from the experimental coordinates and the other two randomizing the initial positions of the probe, in a fashion previously described.18 A 128 × 128 × 128 point grid with a 1.2 Å spacing was used for digitalizing the
9116 J. Phys. Chem. B, Vol. 111, No. 30, 2007 interacting molecules. In each docking run, the GpA monomer corresponding to segment A was kept fixed (i.e., target), whereas the one corresponding to segment B was allowed to rotate and translate around the target (i.e., probe). A rotational sampling interval of 6° was employed, i.e., dense sampling, and the best 4000 solutions were retained for each run and ranked according to their score. As for docking analysis, according to previous work,17,18,26 we selected as nativelike structures all the docked complexes characterized by an R-carbon atom root-mean-square deviation (CR-rmsd) less than 1.0 Å from the native complex. For each molecular system, we considered the ZDOCK score averaged over the scores of all the nativelike complexes resulting from the three independent runs and constituting the nativelike ensemble. These average scores (ZD-s) were then employed in the correlation analysis with thermodynamic data. In detail, correlations were carried out between the relative ZD-s (∆ZD), i.e., the difference between the ZD-s of a given form and that of the wild type, and the relative binding affinity (∆∆G°), i.e., the difference between the ∆G° of a given form and that of the wild type (Table 1 and references therein). The wild type form has ∆ZD-s and ∆G° values equal to 0. Molecular Dynamics Simulations of GpA Dimers. MD simulations were carried out on the solid-state NMR structure in its wild type and 32 mutated forms. Alanine replacements were made on both chains A and B in the GpA dimer by means of the CHARMM program.27 Such input structures of the wild type and of the mutants that differ only for the mutated side chains were, thus, subjected to energy minimization and MD simulations, by using the GBSW implicit membrane/water model recently implemented in CHARMM (i.e., c32b1 version).28 The all-hydrogen parameter set PARAM22/CMAP of the CHARMM force field was used. No cutoff was used for the nonbonded energy evaluation. A smoothing length of 0.6 Å (w ) 0.3 Å) was employed in both PB and GB calculations, by setting a0 ) -0.180, a1 ) 1.817, and s ) 0.952. As for the physical parameters representing the membrane in the GB model, the surface tension coefficient (representing the nonpolar solvation energy) was set to 0.03 kcal/(mol‚Å2). Furthermore, the membrane thickness centered at Z ) 0 was set to 29.0 Å with a membrane smoothing length of 5.0 Å (wm ) 2.5 Å). Minimizations were carried out by using 500 steps of steepest descent followed by adopted basis Newton-Raphson (ABNR) minimization, until the root-mean-square gradient was less than 0.001 kcal/(mol‚Å). MD production phases consisted in 1 ns of Langevin dynamics at 300 K with a friction coefficient of 5.0 ps-1. The lengths of the bonds involving the hydrogen atoms were restrained by the SHAKE algorithm, allowing for an integration time step of 0.001 ps. The secondary structure of the two helices was preserved by assigning distance restraints (i.e., minimum and maximum allowed distances of 2.7 and 3.0 Å, respectively) between the backbone oxygen atom of residue i and the backbone nitrogen atom of residue i + 4. The scaling factor of such restraints was 10, and the force constant at 300 K was 10 kcal/(mol‚Å). An equivalent set of MD simulations was also carried out by imposing interhelical distance restraints between the hydroxy hydrogen atom of T87 (i.e., HG1) and the backbone oxygen atom of V84 (i.e., O). These additional simulations were carried out on the minimized coordinates of the wild type and of all the mutants holding the interhelical H-bonds involving T87. Such a set of mutants excludes all the alanine substitutions for T87 and the alanine substitutions for G83, taken singularly or paired with alanine replacements of L75, I76, G79, V80, and V84 (Table 1).
Dell’Orco et al.
Figure 1. Structural predictions of wild type GpA dimerization. (A) Cartoon representation of the best-scored docking solution, i.e., solution N.3 from the first docking run (magenta), superimposed on the native solid-state NMR structure (green). The CR-rmsd between the two structures is 0.16 Å. (B) Cartoon representation of the minimized structure averaged over the 2000 structures constituting the 1 ns MD trajectory (blue) superimposed on the solid-state NMR structure (green). The CR-rmsd between the two structures is 0.82 Å. (C) Same representation as in (B), in which the amino acid side chains are represented by sticks.
The MD setup reported above was selected from among a large number of trial runs, because it was the one that better maintained the structure of wild type GpA. Trial tests included also the probing of (a) different positions of the input structure relative to the membrane center, (b) different membrane thicknesses (i.e., 25, 30, and 32 Å), (c) different friction coefficients for the Langevin dynamics (i.e., 0.5 and 10 ps-1), (d) Nose-Hoover dynamics at 300 K with and without intrahelix distance restraints, (e) Langevin dynamics without intrahelix distance restraints, and (f) finite cutoffs for the nonbonded interactions. The structures of the wild type and of the mutants averaged over the 1 ns equilibrated trajectories were considered for the structural and correlative analyses. The interaction energies, IE (i.e., the total energy of the helix-helix complex minus the energy of each helix in the complex), were computed both on each frame and on the average minimized structures. The difference between the IE of a given GpA form and that of the wild type (i.e., relative IE or ∆IE) was considered for the correlative analyses. By definition, the ∆IE of the wild type is equal to 0. Results Rigid-Body Docking Simulations of GpA Dimerization: Output Analyses. Three sets of repeated and independent rigidbody docking simulations were run to reproduce GpA dimers in their native and 32 mutated forms. Data concerning both relative affinities as determined by analytical ultracentrifugation experiments11,13,14 and relative docking scores are reported in Table 1. The basic assumption of the approach is that mutated GpA dimers hold almost the same architecture as the wild type form. In this respect, for each variant, the docking algorithm was able to reproduce nativelike dimers with a ranking number strictly related to the binding affinity. In detail, for each independent run concerning the wild type form, a nativelike solution was found within the best-scored 21 solutions out of the 4000 saved configurations (Table 1). The best-scored solution, i.e., solution N.3 from the first run (Table 1), holds a CR-rmsd from the native solid-state NMR structure8,9 equal to 0.16 Å (Figure 1A). Similar results were achieved for all the wild-type-like mutants, as far as binding affinity is concerned (Table 1). Interestingly, the relative population of the nativelike ensemble, i.e., the ensemble of all the nativelike solutions from each independent docking run, correlates linearly both with the relative binding free energy, ∆∆G° (N ) 33, R ) 0.81, p
