Long-Range Conformational Response of a PDZ Domain to Ligand

Dec 18, 2015 - Altis , A.; Nguyen , P. H.; Hegger , R.; Stock , G. Dihedral angle principal component analysis of molecular dynamics simulations J. Ch...
4 downloads 0 Views 3MB Size
Subscriber access provided by CMU Libraries - http://library.cmich.edu

Article

Long-range conformational response of a PDZ domain to ligand binding and release: a molecular dynamics study Cheng Lu, Volker Knecht, and Gerhard Stock J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.5b01009 • Publication Date (Web): 18 Dec 2015 Downloaded from http://pubs.acs.org on December 24, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Long-range conformational response of a PDZ domain to ligand binding and release: a molecular dynamics study Cheng Lu, Volker Knecht∗ and Gerhard Stock∗ Biomolecular Dynamics, Institute of Physics, Albert Ludwigs University, 79104 Freiburg, Germany E-mail: [email protected],[email protected]

Abstract The binding of a ligand to a protein may induce long-range structural or dynamical changes in the biomacromolecule even at sites physically well separated from the binding pocket. A system for which such behavior has been widely discussed is the PDZ2 domain of human tyrosine phosphatase 1E. Here we present results from equilibrium trajectories of the PDZ2 domain in the free and ligand-bound state, as well as nonequilibrium simulations of the relaxation of PDZ2 after removal of its peptide ligand. The study reveals changes in inter-residue contacts, backbone dihedral angles, and Cα positions upon ligand release. Our findings show a long-range conformational response of the PDZ2 domain to ligand release in the form of a collective shift of the secondary structure elements α2 , β2 , β3 , α1 -β4 , and the C terminal loop relative to the rest of the protein away from the N-terminus, and a shift of the loops β2 -β3 and β1 -β2 in opposite direction. The shifts lead to conformational changes in the backbone especially in the β2 -β3 loop but also in the β5 -α2 and the α2 -β6 loop, and are accompanied by changes of inter-residue contacts mainly within the β2 -β3 loop as well as between the α2 helix and other segments. The residues showing substantial changes of inter-residue contacts, backbone conformations, or Cα positions, are considered “key residues” for the long-range conformational response of PDZ2. By comparing these residues with ∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

various sets of residues highlighted by previous studies of PDZ2, we investigate the statistical correlation of the various approaches. Interestingly, we find a considerable correlation of our findings with several works considering structural changes, but no significant correlations with approaches considering energy flow or networks based on inter-residue energies.

Introduction The binding of a ligand to a protein may change the properties at remote regions of the molecule. This phenomenon called allostery is crucial to many biological processes. 1,2 Traditionally, allostery has been attributed to changes in structure 3 or shifts in the populations of alternative conformations. 4 More recently it has been proposed that allostery may be (i) of dynamic origin with only marginal alterations in average structure, 5–9 or (ii) based both on dynamic and structural changes at the same time. 10 Along with ongoing progress in experimental and theoretical approaches, tremendous effort has been devoted to understand allosteric mechanisms on a microscopic level. 10–45 Ligand-induced allostery is observed for a subset of a class of modular domains for proteinprotein interactions denoted as PSD95/Discs large/ZO-1 (PDZ) domains. 11–16,26,27,37,40 PDZ domains mediate the clustering of membrane ion channels by binding to their C-termini. 13–15 They exhibit roughly 90 amino acids which fold into globular shape with a βββαββαβ topology where the second β-strand and the second α-helix form the canonical binding groove, 11–13,16 see Figure 1 (a). The example shown in Figure 1 is the second PDZ domain (PDZ2) of human tyrosine phosphatase 1E (hPTP1E), a protein involved in the regulation of multiple receptor-coupled signal transduction pathways including the regulation of cell growth and apoptosis in breast cancer. 37,46–48 The PDZ2 domain of hPTP1E interacts with the human Fas/CD95 receptor, 46,49 the zyxin-related protein ZRP1, 46,50 the tumor suppressor adenomatous polyposis coli protein (APC), 46 and the guanine nucleotide exchange factor RA-GEF2. 46,51 A C-terminal peptide derived from RA-GEF2 bound to PDZ2/hPTP1E is shown in Figure 1a (yellow). Due to their abundance and small size, and owing to the influential study on evolutionarily conserved pathways of energetic connectivity by Lockless and Ranganathan, 35 PDZ domains have become a popular model system to study allosteric communication in a protein. To this end, a number of “communication networks” have been proposed, which aim to investigate which residue “talks” to which. 10,36–40,43–45 The construction of the network can be based on the correlation of structural or dynamical changes of various protein residues, 10,37–40 which may report on the conformational propagation of a structural perturbation in the system. Alternatively, the analysis has been based on the exchange of vibrational energy between

2

ACS Paragon Plus Environment

Page 2 of 33

Page 3 of 33

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

evolution of the structural response of PDZ2, we perform nonequilibrium MD simulations of the relaxation of PDZ2 after removal of the ligand, which show that the process occurs on a timescale of 100 ns. By considering the various sets of “key residues” identified by our and nine other studies of PDZ2, we investigate the statistical correlation of the various approaches.

Methods Simulation set-up Following previous work, 53 all MD simulations of PDZ2 were performed using GROMACS with a hybrid GPU-CPU acceleration scheme. 56 Here PDZ2 was simulated in aqueous NaCl solution. The protein was described using the all-atom Amber99sb∗ -ILDN 57–59 force field, the water by the TIP3P 60 model, and the ions with the model of Ref.. 61 The side chains of all four histidine residues (32, 53, 71, and 86) of PDZ2 were chosen electrostatically neutral (mono-protonated). The lengths of bonds involving hydrogens were constrained, 62 allowing for a 2 fs time step. Long-range electrostatic interactions were evaluated in reciprocal space using Particle-Mesh Ewald 63 with a maximum spacing for the FFT grid of 0.12 nm and an interpolation via a 6th order polynomial. The minimal cut-off distance for electrostatic and van der Waals interactions was set to 1.2 nm. The temperature of 300 K was maintained via the velocity rescaling algorithm 64 (0.1 ps relaxation time) and the pressure P = 1 bar was controlled using the weak coupling method of Berendsen. 65

Non-equilibrium simulations MD simulations of the PDZ2 domain in the bound and free state in explicit aqueous salt solution under periodic boundary conditions were conducted as shown in Figure 2. Six 1 µs simulations for PDZ2 in the bound (b,eq) and the free state at equilibrium (f,eq) were taken from our previous study. 53 From the ensemble b,eq, n = 20 snapshots were selected randomly, yielding the sub-ensemble b,sel. Hence, from each configuration, the ligand was removed, the resulting void filled with water molecules, and the system obtained was energy minimized using steepest descent, yielding the set of final configurations denoted as f∗∗ . From each of these configurations, a 10 ns N P T simulation with harmonic restraints on the protein atoms was conducted. The final configuration was taken with the positions and box dimensions scaled such that the volume was equal to the average volume of the second half of the N P T simulations, resulting in the set of configurations f∗ /0. From each of these

4

ACS Paragon Plus Environment

Page 4 of 33

Page 5 of 33

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the overall change in the (φ, ψ) angles for a given residue was quantified as ∆2 (φ, ψ) = ∆2φ + ∆2ψ . An alternative measure for the change in backbone conformation for a given residue was obtained via computing the root mean square deviation (rmsd) between the average conformation of the main chain (NH and CO group as well as Cα and C atoms) of this residue in the bound and the free state. Cα shifts Ligand-induced changes in the average positions of individual residues of PDZ2 were monitored by computing the Cα shifts ∆ri ≡ hri ib − hri if ,

i = 1, . . . , Nres .

Here, ri denotes the position of the Cα atom of residue i after fitting the corresponding configuration to the initial structure such as to minimize the rmsd of the Cα atoms of the protein from the reference structure. Furthermore, Nres = 94 is the number of residues, and h. . .ib or h. . .if denote averages over the bound or the free state, respectively. Standard errors for the average positions at equilibrium were obtained from the averages over individual trajectories. Hence, the statistical errors of the shifts were determined via error propagation. When the size of the shift in one component was below its statistical error, the shift was considered insignificant (standard error criterion) and set to zero. This approach yielded results similar to those from a more rigorous analysis using hypothesis tests as described in the SI (Table S4). Rotational fit As the rotation applied for the fitting to the reference structure depends on the configuration of the molecule via the moment of inertia, the separation of internal and rigid body motions is is only straightforward for relatively rigid systems. In the case of the PDZ2 domain, issues might thus arise due to flexibility of the loops. In order to address effects due to loop flexibility, we compared Cα shifts obtained employing two different sets of atoms for the fit, using either (i) all the Cα atoms or (ii) only those of two β-strands shown to be rigid (β1 and β6 ). 53 The Cα shifts obtained in conjunction with these two fitting schemes are shown in Figure 3 (a). The error bars of the two curves overlap with each other, indicating that

6

ACS Paragon Plus Environment

Page 6 of 33

Page 7 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

within the statistical accuracy the shifts obtained using the two fitting schemes agree with each other. A second check was performed by determining the ligand-induced change in the vectors V6,12 along the β1 and V90,84 along the β6 strand. Here the vector V6,12 was chosen to point from Ile-6 to Ala-12 and V90,84 from Glu-90 to Val-84. The angle between the corresponding average vector in the bound and in the free state obtained using the two fitting schemes was calculated. As evident from Table S1, each scheme yielded a rotation angle of zero. Also the shifts of the corresponding Cα atoms were small (0.004 nm to 0.011 nm). Hence, subsequent analyses of Cα shifts were performed by fitting on the β strands, β1 and β6 . Initial conditions of nonequilibrium simulations Figure 3 (b) shows the sizes of Cα shifts between the set of configurations selected from the bound state at equilibrium and the initial configurations of the nonequilibrium simulations after ligand removal, showing effects from the preparation of the nonequilibrium simulations (especially from the energy minimization after replacing the ligand by water). Large shifts are observed in the β2 -β3 loop and moderate shifts in the β5 -α2 loop as well as in the α2 helix. Figure S1 displays the corresponding individual components of the shifts. The large size of the shift in the β2 -β3 loop arises from a strong shift in y-direction (see coordinate system in Figure 1) which is opposite to that observed at equilibrium. Figure 3 (c) shows the sizes of Cα shifts between the set of configurations selected from the bound state at equilibrium and the final configurations of the nonequilibrium simulations, indicating that the nonequilibrium simulations ultimately using solely the experimental structure of the bound state largely reproduce the transition from PDZ2-bound to PDZ2free at equilibrium (observed from simulations using the experimental structures of both the bound and the free state).

Results In order to understand the long-range conformational response of the PDZ2 domain to ligand release or binding, MD simulations of the PDZ2 domain in the bound and free state in solution at equilibrium, as well as nonequilibrium simulations of PDZ2 in solution after ligand release, were conducted. In the following, first the equilibrium and then the nonequilibrium simulations are discussed.

7

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

Cα shift (nm)

0.4

all Cα β1,β6

(a)

0.2

0 β1 10

Cα shift (nm)

0.6

β2 20

β3 30

α1

β4 β5

40 50 60 residue index

α2 70

80

β6 90

∆rbf ∆rsel,f*0

(b)

0.4

0.2

0

β1 10

0.4

Cα shift (nm)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

β2 20

β3 30

α1

β4 β5

40 50 60 residue index

α2 70

80

β6 90

∆rbf,eq ∆rsel,f*100

(c)

0.2

0 β1 10

β2 20

β3 30

α1

β4 β5

40 50 60 residue index

α2 70

80

β6 90

Figure 3: (a) Rotational fit: Sizes of Cα shifts upon ligand release at equilibrium obtained from two different sets of Cα atoms for the least squares fit to remove overall translation and rotation (either all Cα atoms or only Cα atoms in the β1 and β6 strands). (b,c) Nonequilibrium simulation conditions: Sizes of Cα shifts between configurations selected from the bound state at equilibrium and (b) the initial or (c) final configurations of the nonequilibrium simulations after ligand removal. In (b,c), the sizes of Cα shifts upon ligand release at equilibrium are shown for comparison. In (a – c), the widths of the shaded regions indicate standard errors.

Comparison to crystal structure The average rmsd of the simulated Cα atom positions for PDZ2 in solution in our simulations from the corresponding crystal structure inferred from x-ray diffraction was 0.16 ± 0.01 nm 8

ACS Paragon Plus Environment

Page 8 of 33

Page 9 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

for PDZ2-f and 0.14 ± 0.01 nm for PDZ2-b (the standard errors were obtained from averages over individual trajectories). This structural difference is similar to that observed when comparing structural models for proteins in crystals derived from x-ray diffraction with corresponding structures of proteins in solution inferred from nuclear magnetic resonance (NMR) spectroscopy. Andrec et al. 67 compared 148 pairs of protein structure pairs (NMR vs crystals). For each pair, the backbone rmsd between the mean NMR structure and the crystal structure was determined. The rmsds obtained were averaged over all protein structure pairs, yielding an average rmsd of 0.14 nm. A detailed analysis suggested that some of the differences arise from attractive intermolecular interactions in the crystal.

Contacts The probabilities of contacts of individual residues of PDZ2 (residues 1 – 94) with the ligand (residues -5 – 0) when the ligand is bound to PDZ2 are given in Table S2 and indicated in Figure 1. It is found that all six residues of the peptide interact with the PDZ2 domain. The α2 helix and the β2 strand forming the so called binding groove (or binding pocket) are most frequently in contact with the ligand (red color). Whereas the β2 strand contacts all peptide residues, the α2 helix only contacts the C-terminal half of the peptide. Besides the canonical binding groove, also several loops contact the peptide. In particular, the β2 -β3 loop contacts the N-terminus (Glu(-5) and Gln(-4)) and the β1 -β2 loop the C-terminus of the peptide (Val(0)). Table 1 shows the inter-residue contacts whose probabilities change more than 0.2 upon ligand removal. Significant changes mainly involve residues or secondary structure segments contacting the ligand in the bound state. An exception is an increase in contact probability within the β4 strand.

Backbone conformational change The changes in backbone dihedral angles ∆(φ, ψ) along the amino acid sequence are shown in Figure 4 (a) and mapped to the three-dimensional structure in Figure 4 (b). The residues showing changes in backbone dihedral angles ∆(φ, ψ) ≥ 0.16 are listed in Table S3. Significant conformational changes occur in the loops β2 -β3 , β5 -α2 , and α2 -β6 . Especially in the former two loops which are in contact with residues -5 to -2 of the ligand as shown in Table 1, ligand release leads to a significant configurational rearrangement. Table S3 and Figure 4 and S2 show that ten out of the 11 residues of the β2 -β3 loop and two out of the seven residues of the β5 -α2 loop change their average backbone conformation by more than the above given threshold.

9

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 10 of 33

Page 11 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Table 1: Contacts between residues i in segments Si and residues j in segments Sj whose probabilities change significantly upon ligand removal at equilibrium (∆Peq > 0.2). The symbols ∆Peq and ∆Pneq denote the changes in contact probabilities at equilibrium and nonequilibrium, respectively. The latter are the changes in contact probabilities in the final 20 ns of the nonequilibrium simulations and the contact probabilities in the bound state at equilibrium. The superscripts “*” or “**” indicate that, in the bound state, the respective residue is sometimes or always respectively in contact with the ligand. Contact changes mainly occur within the β2 -β3 loop (residues 24 – 34) and between the α2 helix (residues 73 – 80) and other segments. i Ala-74* Leu-78** Leu-78** Leu-78** Arg-79** Leu-78** Ile-20** Gly-24* Gly-24* Gly-25 Val-26 Asn-27* Gly-25 Val-58

j Ala-69 Leu-66 Ile-20** Val-22** Lys-13* Leu-87 Val-22** Asn-27* Val-30* Ser-29* Ser-29* Ser-29* Thr-70 Val-61

∆Peq -0.25 -0.43 0.23 0.24 -0.39 -0.40 -0.26 0.24 -0.27 0.25 0.21 0.21 -0.24 0.29

∆Pneq -0.32 -0.45 0.34 0.26 -0.25 -0.29 -0.27 – -0.23 – – 0.35 – 0.30

Cα Shifts Changes in backbone conformations upon ligand release are expected to lead to shifts in residue positions. The latter were monitored in terms of shifts in the average positions of Cα atoms within an internal coordinate system of the protein defined in Figure 1 (a). The Cα shifts of the individual residues along the three directions are shown in Figure 5. The figure and Table 2 indicate that ligand release leads to a shift of the α2 helix, the β2 and β3 strands, as well as the α1 -β4 and the C-terminal loop in positive x-direction, and of residues 14 – 17 of the β1 -β2 loop as well as residues 25 – 28 of the β2 -β3 loop in negative x-direction. In particular, the α1 -β4 loop which acts as a steric bridge (see Table S6) transduces the shift of the β3 strand to the C-terminal loop. Furthermore, Figure 5 shows that the Cterminus of the α2 helix moves in negative y- and negative z-direction, and the N-terminus of the helix in positive y-direction. This suggests a rotation of the helix. This rotation is presumably related to the contact changes of the α2 helix with other secondary structure segments described above.

11

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 12 of 33

Page 13 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Table 2: Cα shifts in PDZ2 upon ligand release at equilibrium conditions. The shifts are given with respect to the internal coordinate system defined in Figure 1 (a). Only residues showing shifts with sizes above 0.025 nm are listed. The symbol “–” indicates shifts which are not significant according to Welch’s test or the standard error criterion. S

R

N β1 -β2 β1 -β2 β1 -β2 β1 -β2 β1 -β2 β2 β2 β2 β2 β2 -β3 β2 -β3 β2 -β3 β2 -β3 β2 -β3 β2 -β3 β2 -β3 β2 -β3 β3 β3 β3 β3 β3 β3 -α1 α1 -β4 β4 -β5 β5 -α2 β5 -α2 β5 -α2 β5 -α2 β5 -α2 β5 -α2 β5 -α2 α2 α2 α2 α2 α2 α2 C C

2 14 15 16 17 18 20 21 22 23 24 25 26 27 28 29 30 31 35 36 37 38 39 43 54 62 66 67 68 69 70 71 72 73 75 76 78 79 80 92 93

∆xeq (nm) -0.01(1) -0.04(1) -0.06(1) -0.04(1) -0.02(0) – 0.03(0) 0.06(1) 0.05(1) 0.04(1) 0.02(1) -0.07(5) -0.13(8) -0.19(11) -0.16(11) – – – 0.03(0) 0.03(1) 0.03(1) 0.06(1) 0.05(1) -0.02(0) 0.03(0) – -0.01(1) -0.05(2) – 0.03(2) 0.06(3) 0.06(2) 0.09(2) 0.04(2) 0.03(1) 0.04(1) – 0.06(0) 0.03(0) 0.02(0) 0.02(1)

∆yeq (nm) – 0.01(0) – -0.03(0) -0.01(0) – – – – – 0.01(1) 0.05(1) – – 0.05(4) 0.08(6) 0.12(7) 0.07(7) 0.01(1) – -0.01(0) -0.02(0) -0.02(0) 0.02(0) -0.02(0) -0.02(0) 0.03(1) 0.06(2) 0.09(3) 0.06(1) 0.06(1) 0.04(2) 0.06(1) 0.06(1) – 0.02(1) -0.05(1) -0.03(1) 0.01(0) – –

13

ACS Paragon Plus Environment

∆zeq (nm) -0.03(2) 0.03(1) 0.04(1) 0.04(0) 0.03(0) 0.03(0) 0.02(0) 0.02(0) 0.01(0) – – – -0.05(4) 0.03(2) – – – – – 0.01(0) – -0.01(0) -0.02(0) – -0.01(0) -0.01(0) – 0.04(1) 0.06(2) 0.06(2) 0.03(1) – – 0.03(1) -0.03(2) -0.02(1) -0.02(1) -0.04(1) – -0.02(0) -0.03(1)

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Non-equilibrium simulations The transient behavior during the relaxation of PDZ2 after ligand removal was investigated using nonequilibrium simulations. The simulation set-up is shown in Figure 2. Table 1 shows that the changes in inter-residue contacts in the nonequilibrium simulations are (in terms of the size and sign of the alterations in corresponding contact probabilities) very similar to the changes at equilibrium. Figure 5 indicates that also the sizes of the final shifts averaged over the final 20 ns of the 100 ns nonequilibrium simulations (|∆rsel,f ∗ 100 |) are similar to the shifts at equilibrium. Especially the shifts of the α2 helix, the β2 and β3 strands, as well as the α1 -β4 and the C-terminal loop in positive x-direction (see also Table S5), and of the β1 -β2 and β2 -β3 loop in negative x-direction, are largely reproduced. As noted in the Methods section and displayed in Figure 3 (b), the preparation of the initial configurations of the nonequilibrium simulations lead to large shifts, especially in the β2 -β3 loop and in the α2 helix. Hence, the initial configurations of the nonequilibrium simulations do not fully mimic the bound state. The β2 and β3 strands, as well as the α1 -β4 and the C-terminal loop, though, are not affected and their positions are similar to those in the bound state. Hence, the transient shifts in these regions from f∗ /0 to f∗ /100 might reflect features of the time-dependent relaxation behavior after ligand removal under more natural conditions. The time dependent Cα shifts in x-direction, relative to the set of structures selected from the bound state used to derive the initial configurations of the nonequilibrium simulations, ∆rsel,f ∗ /t , are shown in Figure 6. The shifts of the β2 and β3 strands, as well as the α1 -β4 and the C-terminal loop in positive x-direction are seen to built up gradually.

14

ACS Paragon Plus Environment

Page 14 of 33

Page 15 of 33

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Discussion and Conclusions We first want to mention several methodological issues raised by our study. As a general point, we note that distal conformational changes due to ligand binding or release are typically quite small (Table 2). An appropriate investigation therefore requires sufficient statistical sampling (i.e., several long trajectories) as well as a sound analysis of the statistical significance of small structural changes (see SI). Employing Cartesian coordinates (e.g., to calculate Cα shifts), we found it important to carefully check for possible uncertainties associated with the rotational fit of the trajectory (Figure 3a), which may easily obscure small structural changes. Concerning the nonequilibrium simulations of PDZ2, we note the ambiguous choice of initial conditions reflecting the ligand release (Figure 3c), which may affect the interpretation of the resulting structural dynamics. Dealing with the propagation of small conformational changes, it appears not surprising that relatively small modifications of the molecular system may result in quite different behavior. For example, a recent study monitoring cross correlations between inter-residue interaction energies and structural fluctuations 45 showed that ligand binding affects PDZ3PSD95 and PDZ2 in a different manner. Moreover, it is interesting to compare our results for wild-type PDZ2 to a previous computational study of a photoswitchable PDZ2 domain (PDZ2S), which mimics the free-bound transition by cis → trans photo-isomerization of an azobenzene photoswitch. 53 Generally speaking, the initial conformation change due to the photoswitch in PDZ2S is larger than the rearrangement in PDZ2 due to ligand binding, which affects larger conformational heterogeneity and higher residue fluctuations of PDZ2S. 53 As a consequence, the subsequent propagation of the structural perturbation of the two systems agrees only in the general trend but differs in atomistic detail. For example, large changes of the backbone conformation occur for 14 residues in PDZ2 and 24 residues in PDZ2S, where the two sets have nine residues in common (residues 25, 27 – 29, and 32 – 34 in the β2 -β3 loop as well as 67 – 68 in the β5 -α2 loop). However, we only find a single common residue pair (Ala-74, Ala-69), when we compare inter-residue contacts of PDZ2 and PDZ2S (Table 1 and Table 1 in Ref. 53 ). To summarize our findings of the MD simulation study of PDZ2, the above results reveal that the release of the ligand affects a complex and concerted response of the PDZ2 domain, which originates from direct interactions between ligand and protein. Apart from obvious contacts with the binding-pocket segments α2 and β2 , the ligand directly interacts via its N- and C-termini with the β1 -β2 and β2 -β3 loops of PDZ2, respectively (Figure 1). These interactions are important, because they affect a concerted conformational change of the binding-pocket with respect to these loops, i.e., α2 and β2 move in +x direction and β1 -β2 and β2 -β3 in −x direction in Figure 5(a). This conformational shift of α2 and β2 extends via several relatively rigidly connected secondary structure elements all the way 16

ACS Paragon Plus Environment

Page 16 of 33

Page 17 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

to the C-terminus of PDZ2. That is, β2 is hydrogen bonded to β3 , which through residue 54 of the α1 -β4 loop is linked to the C-terminus. As a consequence, we find an overall conformational transition of interconnected secondary-structure segments (α2 , β2 , β3 , α1 -β4 and the C-terminal loop) with respect to the β1 -β2 and β2 -β3 loops of PDZ2. Introducing considerable strain in the protein, this structural change subsequently affects a significant conformational rearrangement of the β2 -β3 loop (Figure 4). Due to the finite rigidity of the inter-residue connections, the sizes of the conformational shifts tend to decrease with distance from the binding pocket (Figure 5 and Table 2). Taken together, these findings suggest that the interconnected secondary-structure segments α2 , β2 , β3 , α1 -β4 and the Cterminal loop may provide a mechanism of long-range conformational propagation in PDZ2. Moreover, our study identifies a list of “key residues” (i.e., residues that are important for allosteric communication), which were shown to undergo significant changes of inter-residue contacts (Table 1), backbone dihedral angles (Table S3) or Cα positions (Table 2). The long-range conformational changes in the PDZ2 domain of hPTP1E induced by the C-terminal peptide ligand derived from the guanine nucleotide exchange factor RAGEF2 revealed here could change the affinity of PDZ2/hPTP1E to its other interaction partners, the human Fas/CD95 receptor, 46,49 the zyxin-related protein ZRP1, 46,50 and the tumor suppressor adenomatous polyposis coli protein (APC). 46 The modulation of the affinity of PDZ2/hPTP1E with these other interaction partners via binding of RA-GEF2 may be an important factor in the regulation of the receptor-coupled signal transduction pathways involving hPTP1E, such as the regulation of cell growth and apoptosis in breast cancer. 37,46–48 It is interesting to compare our results to related studies of the allosteric communication in the PDZ2 domain of hPTP1E. As a simple (yet not unambiguously defined) observable, we follow Ref. 39 and consider the key residues predicted by each method. (A full list of all residues and methods is found in Table S7). By using Fisher’s exact test (see SI), we calculate corresponding p-values to assess the statistical correlation of the various approaches with our study. An apparent agreement between both studies is unlikely to be by chance if p ≪ 1 but may be by chance of p ≤ 0.1. The p-values obtained are listed in Table 3. On the experimental side, Fuentes, Der, and Lee 37 used side-chain (2 H-methyl) and backbone (15 N) NMR spin relaxation to identify residues that show a significant change of pico- or nanosecond fluctuations upon ligand binding of PDZ2. Interestingly, the resulting p-value of 0.0008 (line “NMR” in Table 3) reveals a highly significant correlation of the experiment with our study. Fuentes, Gilmore, Mauldin, and Lee 27 proposed a set of key residues based on the effect of mutations on the binding affinity of the ligand as well as side chain 2 H-based dynamics probed by NMR spectroscopy. The resulting p-value of 0.0699 (line “Mut/NMR” in Table 3) indicates significant correlation of the experiment with our work. On the theoretical side, various computational studies have reported communication networks of the allosteric response in PDZ2, which are based either on the correlation of 17

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

structural changes of various protein residues 10,38–40 or on the analysis of energy-based correlations. 41–43,45 Table 3: Comparison of the present work to previous studies in terms of key residues identified (see text), as given from the p-values obtained from Fisher’s exact test and the sensitivities σ. Method NMR 37 PRS 38 RestrMD 10 PSN-ENM 39 Mut/NMR 27 MC 40 E-RMSF 45 IEcorr 43 PEN 44 ET 42

p 0.0008 0.0009 0.0025 0.0555 0.0699 0.2012 0.4806 0.8231 0.8938 0.9595

σ 0.4 0.6 0.5 0.5 0.5 0.4 0.2 0.2 0.1 0.0

The structure-based network approaches includes the study of Dhulesia et al., 10 who conducted MD simulations of PDZ2 in explicit solvent using ensemble-averaged NMR restraints (“RestrMD”), as well as the perturbation response scanning (“PRS”) approach by Gerek and Ozkan, 38 who probed the average response of the protein structure (Cα positions of an elastic network model) to random perturbations of the binding pocket. The obtained p-values of ≈ 10−3 (Table 3) report on significant correlation of these methods of our study with experiment. (We note in passing, though, that the good agreement with the PRS results indicated from the p-value refers to the size of Cα shifts, whereas a closer analysis of the PRS results revealed no correlations to our findings in terms of the relative orientations of the shift vectors, see the SI). Another method was proposed by Raimondi et al., 39 who combined 20 ns MD simulations with an elastic network model to construct a protein structure network (“PSN-ENM”) by computing average numbers of interatomic contacts between residues. Finally Cilia et al. 40 applied Monte Carlo sampling (“MC”) of the side chain conformational space to estimate a mutual information matrix, that reports on ligand-induced changes of the individual residue pairs. The predicted key residues of the latter two methods still correlate well with our results, however, with comparatively low significance (Table 3). Alternatively, communication networks of the allosteric response in PDZ2 have been constructed via the analysis of energy-based correlations. This includes nonequilibrium MD studies of the transport of vibrational energy between protein residues, 41,42 (“ET”) a MDbased network of the correlations of inter-residue interaction energies, 43 (“IEcorr”) a protein energy network (“PEN”) based inter-residue interaction energies, 44 as well as a recent study of the cross correlations between energetic and structural fluctuations. 45 Interestingly, none of these methods show significant correlation with our work. In addition, the sensitivities denoting the fraction of residues highlighted by our work also highlighted by these energy18

ACS Paragon Plus Environment

Page 18 of 33

Page 19 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

based approaches were lower (0 – 0.2, see Table 3) than those of the other methods (0.4 – 0.6, see Table 3). The absence of significant correlation with our work for the energy-based methods may be caused by the considerable timescale separation between local energy flow and long-range structural change. That is, vibrational energy in disordered condensed phase systems typically dissipates within picoseconds, while the propagation of structural change was found occur on a 100 ns timescale in our nonequilibrium simulations of PDZ2 (and will be certainly much longer for larger proteins). Nonetheless, it has been suggested that residues that provide a preferred pathway of vibrational energy at the same time may indicate possible allosteric connectivity in proteins. 32,41,42 The present findings for PDZ2 do not support this idea, though. To conclude, we have presented an allosteric mechanism for the response of PDZ2 to ligand binding based on collective shifts of secondary structure segments and changes in backbone conformations as well as inter-residue contacts. Our results show significant correlations with previous studies based on structural changes but not with energy-based analyses which can be explained from the timescales involved in these two classes of phenomena.

Acknowledgments The authors thank Sebastian Buchenberg and Markus Schumacher for inspiring and helpful discussions. This work has been supported by the Deutsche Forschungsgemeinschaft and the China Scholarship Council.

Supporting Information Directional changes of β1 and β6 strand as well as individual components of Cα shifts, both at equilibrium for two different schemes for the rotational fit, probabilities of inter-residue contacts between PDZ2 and the ligand, changes in backbone conformation of individual residues of PDZ2 in terms of shifts in average backbone dihedral angles upon ligand release at equilibrium, residues showing large changes in their main chain conformations upon ligand release at equilibrium and extent of changes, nomenclature by Hollingsworth and Karplus for regions in the Ramachandran plot, thorough analysis of statistical significance of Cα shifts upon ligand release at equilibrium, list of substantial Cα shifts, analysis of contacts revealing the residue 54 in α1 -β4 loop is a steric bridge between the β3 strand and the C-terminal loop, description of Fisher’s exact test, detailed description of previous studies of allostery in PDZ2, as well as key residues highlighted by previous studies and present work. This information is available free of charge via the Internet at http://pubs.acs.org.

19

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References (1) Nussinov, R.; Tsai, C.-J. Allostery in Disease and in Drug Discovery. Cell 2013, 153, 293 – 305. (2) Motlagh, H. N.; Wrabl, J. O.; Li, J.; Hilser, V. J. The ensemble nature of allostery. Nature 2014, 508, 331 – 339. (3) Monod, J.; Wyman, J.; Changeux, J.-P. On the Nature of Allosteric Transitions: A Plausible Model. J. Mol. Biol. 1965, 12, 88–118. (4) Cui, Q.; Karplus, M. Allostery and cooperativity revisited. Protein Sc. 2008, 17, 1295– 1307. (5) Cooper, A.; Dryden, D. T. F. Allostery without conformational change - a plausible model. Europ. Biophys. J. Biophys. Lett. 1984, 11, 103–109. (6) Kern, D.; Zuiderweg, E. R. P. The role of dynamics in allosteric regulation. Curr. Opin. Struct. Biol. 2003, 13, 748–757. (7) Popovych, N.; Sun, S.; Ebright, R.; Kalodimos, Dynamically driven protein allostery. Nat. Struct. Mol. Biol. 2006, 13, 831–838. (8) Frederick, K. K.; Marlow, M. S.; Valentine, K. G.; Wand, A. J. Conformational entropy in molecular recognition by proteins. Nature 2007, 448, 325 – 329. (9) Smock, R. G.; Gierasch, L. M. Sending Signal Dynamically. Science 2009, 324, 198– 203. (10) Dhulesia, A.; Gsponer, J.; Vendruscolo, M. Mapping of two networks of residues that exhibit structural and dynamical changes upon binding in a PDZ domain protein. J. Am. Chem. Soc. 2008, 130, 8931–8939. (11) Harris, B. Z.; Lim, W. A. Mechanism and role of PDZ domains in signaling complex assembly. J. Cell Sci. 2001, 114, 3219–3231. (12) van Ham, M.; Hendriks, W. PDZ domains – glue and guide. Mol. Biol. Rep. 2003, 30, 69–82. (13) Chi, C. N.; Bach, A.; Stromgaard, K.; Gianni, S.; Jemth, P. Ligand binding by PDZ domains. Biofactors 2012, 38, 338–348. (14) Doyle, D. A.; Lee, A.; Lewis, J.; Kim, E.; Sheng, M.; MacKinnon, R. Crystal structures of a complexed and peptide-free membrane protein-binding domain: Molecular basis of peptide recognition by PDZ. Cell 1996, 85, 1067–1076. 20

ACS Paragon Plus Environment

Page 20 of 33

Page 21 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(15) Zhang, J.; Sapienza, P. J.; Ke, H.; Chang, A.; Hengel, S. R.; Wang, H.; Jr., G. N. P.; Lee, A. L. Crystallographic and Nuclear Magnetic Resonance Evaluation of the Impact of Peptide Binding to the Second PDZ Domain of Protein Tyrosine Phosphatase 1E. Biochemistry 2010, 49, 9280–9291. (16) Lee, H.-J.; Zheng, J. J. PDZ domains and their binding partners: structure, specificity, and modification. Cell Commun. Signal. 2010, 8, 8. (17) Bowman, G. R.; Geissler, P. L. Equilibrium fluctuations of a single folded protein reveal a multitude of potential cryptic allosteric sites. Proc. Natl. Acad. Sci. U. S. A. 2012, 109, 11681–11686. (18) McLeish, T. C. B.; Rodgers, T. L.; Wilson, M. R. Allostery without conformational change: modelling protein dynamics at multiple scales. Phys. Biol. 2013, 10, 056004. (19) Rodgers, T. L.; Townsend, P. D.; Burnell, D.; Jones, M. L.; Richards, S. A.; McLeish, T. C. B.; Pohl, E.; Wilson, M. R.; Cann, M. J. Modulation of Global LowFrequency Motions Underlies Allosteric Regulation: Demonstration in CRP/FNR Family Transcription Factors. PLoS. Biol. 2013, 11, 1001651. (20) Hub, J. S.; Kubitzki, M.; de Groot, B. L. Spontaneous quaternary and tertiary TR transitions of human hemoglobin in molecular dynamics simulation. PLoS Comput. Biol. 2010, 6, e1000774. (21) Krukau, A.; Knecht, V.; Lipowsky, R. Allosteric control of kinesin’s motor domain by tubulin: a molecular dynamics study. Phys. Chem. Chem. Phys. 2014, 16, 6189–6198. (22) Vesper, M. D.; de Groot, B. L. Collective Dynamics Underlying Allosteric Transitions in Hemoglobin. PLoS Comp. Biol. 2013, 9, e1003232. (23) Pandini, A.; Fornili, A.; Fraternali, F.; Kleinjung, J. Detection of allosteric signal transmission by information-theoretic analysis of protein dynamics. FASEB J. 2012, 26, 868–881. (24) Elber, R. Simulations of allosteric transitions. Curr. Opin. Struct. Biol. 2011, 21, 167– 172. (25) LeVine, M. V.; Weinstein, H. NbIT - A New Information Theory-Based Analysis of Allosteric Mechanisms Reveals Residues that Underlie Function in the Leucine Transporter LeuT. PLoS Comput. Biol. 2014, 10, e1003603. (26) Petit, C. M.; Zhang, J.; Sapienza, P. J.; Fuentes, E. J.; Lee, A. L. Hidden dynamic allostery in a PDZ domain. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 18249–18254. 21

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(27) Fuentes, E. J.; Gilmore, S. A.; Mauldin, R. V.; Lee, A. L. Evaluation of energetic and dynamic coupling networks in a PDZ domain protein. J. Mol. Biol. 2006, 364, 337–351. (28) Chennubhotla, C.; Bahar, I. Signal propagation in proteins and relation to equilibrium fluctuations. Plos Comput. Biol. 2007, 9, 1716. (29) Stacklies, W.; Xia, F.; Gr¨ater, F. Dynamic allostery in the methionine repressor revealed by force distribution analysis. PLoS Comput. Biol. 2009, 5, e1000574, 11. (30) Sharp, K.; Skinner, J. J. Pump-probe molecular dynamics as a tool for studying protein motion and long range coupling. Proteins 2006, 65, 347–361. (31) Gerek, Z. N.; Keskin, O.; Ozkan, S. B. Identification of specificity and promiscuity of PDZ domain interactions through their dynamic behavior. Proteins 2009, 77, 796–811. (32) Nguyen, P. H.; Hamm, P.; Stock, G. In Proteins: Energy, Heat and Signal Flow ; Leitner, D., Straub, J., Eds.; Taylor and Francis/CRC Press: London, 2009; pp 151– 169. (33) Leitner, D.; Straub, J. Proteins: Energy, Heat and Signal Flow ; Taylor and Francis/CRC Press: London, 2009. (34) Gnanasekaran, R.; Agbo, J. K.; Leitner, D. M. Communication maps computed for homodimeric hemoglobin: computational study of water-mediated energy transport in proteins. J. Chem. Phys. 2011, 135, 065103. (35) Lockless, S. W.; Ranganathan, R. Evolutionarily Conserved pathways of energetic connectivity in protein families. Science 1999, 286, 295–299. (36) Chi, C. N.; Elfstrom, L.; Shi, Y.; Snall, T.; Engstrom, A.; Jemth, P. Reassessing a sparse energetic network within a single protein domain. Proc. Natl. Acad. Sci. U. S. A. 2008, 105, 4679–4684. (37) Fuentes, E. J.; Der, C. J.; Lee, A. L. Ligand-dependent dynamics and intramolecular signaling in a PDZ domain. J. Mol. Biol. 2004, 335, 1105–1115. (38) Gerek, Z. N.; Ozkan, S. B. Change in Allosteric Network Affects Binding Affinities of PDZ Domains: Analysis through Perturbation Response Scanning. PLoS Comput. Biol. 2011, 7, 1002154. (39) Raimondi, F.; Felline, A.; Seeber, M.; Mariani, S.; Fanelli, F. A Mixed Protein Structure Network and Elastic Network Model Approach to Predict the Structural Communication in Biomolecular Systems: The PDZ2 Domain from Tyrosine Phosphatase 1E As a Case Study. J. Chem. Theory Comput. 2013, 9, 2504–2518. 22

ACS Paragon Plus Environment

Page 22 of 33

Page 23 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(40) Cilia, E.; Vuister, G. W.; Lenaerts, T. Accurate Prediction of the Dynamical Changes within the Second PDZ Domain of PTP1e. PLoS Comput. Biol. 2012, 8, e1002794. (41) Ota, N.; Agard, D. A. Intramolecular signaling pathways revealed by molecular anisotropic thermal diffusion. J. Mol. Biol. 2005, 351, 345–354. (42) Ho, B. K.; Agard, D. A. Conserved tertiary couplings stabilize elements in the PDZ fold, leading to characteristic patterns of domain conformational flexibility. Protein Sci. 2010, 19, 398–411. (43) Kong, Y.; Karplus, M. Signaling pathways of PDZ2 domain: A molecular dynamics interaction correlation analysis. Proteins 2009, 74, 145–154. (44) Vijayabaskar, M. S.; Vishveshwara, S. Interaction Energy Based Protein Structure Networks. Biophys. J. 2010, 99, 3704 – 3715. (45) Morra, G.; Genoni, A.; Colombo, G. Mechanisms of Differential Allosteric Modulation in Homologous Proteins: Insights from the Analysis of Internal Dynamics and Energetics of PDZ Domains. J. Chem. Theory Comput. 2014, 10, 5677–5689. (46) Milev, S.; Bjelic, S.; Georgiev, O.; Jelesarov, I. Energetics of peptide recognition by the second PDZ domain of human protein tyrosine phosphatase 1E. Biochemistry 2007, 46, 1064–1078. (47) Freiss, G.; Vignon, F. Protein tyrosine phosphatases and breast cancer. Crit. Rev. Oncol./Hematol. 2004, 52, 9–17. (48) Palmer, A.; Zimmer, M.; Erdmann, K.; Eulenburg, V.; Porthin, A.; Heumann, R.; Deutsch, U.; Klein, R. EphrinB phosphorylation and reverse signaling: regulation by Src kinases and PTP-BL phosphatase. Mol. Cell 2002, 9, 725–737. (49) Sato, T.; Irie, S.; Kitada, S.; Reed, J. FAP-1 – a protein-tyrosine-phosphatase that associates with FAS. Science 1995, 268, 411–415. (50) Murthy, K.; Clark, K.; Fortin, Y.; Shen, S.; Banville, D. ZRP-1, a zyxin-related protein, interacts with the second PDZ domain of the cytosolic protein tyrosine phosphatase hPTP1E. J. Biol. Chem. 1999, 274, 20679–20687. (51) Gao, X.; Satoh, T.; Liao, Y.; Song, C.; Hu, C.; Kariya, K.; Kataoka, T. Identification and characterization of RA-GEF-2, a Rap guanine nucleotide exchange factor that serves as a downstream target of M-Ras. J. Biol. Chem. 2001, 276, 42219–42225.

23

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(52) Buchli, B.; Waldauer, S. A.; Walser, R.; Donten, M. L.; Pfister, R.; Bloechliger, N.; Steiner, S.; Caflisch, A.; Zerbe, O.; Hamm, P. Kinetic response of a photoperturbed allosteric protein. Proc. Nat. Ac. Sc. U.S.A. 2013, 110, 11725–11730. (53) Buchenberg, S.; Knecht, V.; Walser, R.; Hamm, P.; Stock, G. Long-range conformational transition of a photoswitchable allosteric protein: molecular dynamics simulation study. J. Phys. Chem. B 2014, 118, 13468–13476. (54) Ikeguchi, M.; Ueno, J.; Sato., M.; Kidera, A. Protein structural change upon ligand binding: linear response theory. Phys. Rev. Lett. 2005, 94, 078102. (55) Yang, L.-W.; Kitao, A.; Huang, B.-C.; G¯ o, N. Ligand-Induced Protein Responses and Mechanical Signal Propagation Described by Linear Response Theories. Biophys. J. 2014, 107, 1415 – 1425. (56) Pronk, S.; Pall, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 2013, 29, 845–854. (57) 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 2006, 65, 712–725. (58) Best, R. B.; Hummer, G. Optimized molecular dynamics force fields applied to the helix- coil transition of polypeptides. J. Phys. Chem. B 2009, 113, 9004–9015. (59) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved side-chain torsion potentials for the Amber ff99SB protein force field. Proteins 2010, 78, 1950–1958. (60) 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. (61) Joung, I. S.; III, T. E. C. Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations. J. Phys. Chem. B 2008, 112, 9020–9041. (62) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463–1472.

24

ACS Paragon Plus Environment

Page 24 of 33

Page 25 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(63) 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. (64) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 0141011–0141017. (65) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A. R. H. J.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684–3690. (66) Altis, A.; Nguyen, P. H.; Hegger, R.; Stock, G. Dihedral angle principal component analysis of molecular dynamics simulations. J. Chem. Phys. 2007, 126, 244111. (67) Andrec, M.; Snyder, D. A.; Zhou, Z.; Young, J.; Montellone, G. T.; Levy, R. M. A large data set comparison of protein structures determined by crystallography and NMR: Statistical test for structural differences and the effect of crystal packing. Proteins 2007, 69, 449–465.

25

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 26 of 33

Page 27 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

Journal of Chemical Theory and Computation

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47

b, eq

Δrsel Δrmin b,sel f**

Select n snapshots

Remove ligand, Add water

Δrbf,eq

Page 28 of 33

Δrpr MD: 10ns position restraints

Δrbf

f, eq

f*/0 MD: 100ns

f*/100ns ACS Paragon Plus Environment

Page 29 of 33

Cα shift (nm)

0.4

all Cα β1,β6

(a)

0.2

0 β1 10

β2 20

β3 30

40

α1 50

β4 β5 60

α2 70

80

β6 90

residue index

Cα shift (nm)

0.6

∆rbf ∆rsel,f*0

(b)

0.4

0.2

0

β1 10

β2 20

β3 30

40

α1 50

β4 β5 60

α2 70

80

β6 90

residue index 0.4

Cα shift (nm)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

∆rbf,eq ∆rsel,f*100

(c)

0.2

0 β1 10

β2 20

β3 30

40

α1 50

β4 β5 60

residue index ACS Paragon Plus Environment

α2 70

80

β6 90

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment

Page 30 of 33

Page 31 of 33

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

Journal of Chemical Theory and Computation

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48

ACS Paragon Plus Environment

Page 32 of 33

Page 33 of 33

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 ACS Paragon Plus Environment