Ensemble- and Rigidity Theory-Based Perturbation Approach To

Subscriber access provided by READING UNIV

Article

Ensemble- and rigidity theory-based perturbation approach to analyze dynamic allostery Christopher Pfleger, Alexander Minges, Markus Boehm, Christopher L McClendon, Rubben Torella, and Holger Gohlke J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00529 • Publication Date (Web): 07 Nov 2017 Downloaded from http://pubs.acs.org on November 16, 2017

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 46

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

Ensemble- and rigidity theory-based perturbation approach to analyze dynamic allostery

Christopher Pfleger†, Alexander Minges†, Markus Boehm§, Christopher L. McClendon§, Rubben Torella§, Holger Gohlke†,* †

Mathematisch-Naturwissenschaftliche Fakultät, Institut für Pharmazeutische und

Medizinische Chemie, Heinrich-Heine-Universität Düsseldorf, 40225 Düsseldorf, Germany §

Medicinal Sciences, Pfizer Inc., 1 Portland Street, Cambridge, Massachusetts 02139, United States

Running title: Allosteric coupling deduced from network rigidity Keywords: allosteric signaling, allosteric cooperativity, rigidity analyses, network topology

*

Address: Universitätsstr. 1, 40225 Düsseldorf, Germany.

Phone: (+49) 211 81 13662; Fax: (+49) 211 81 13847 E-mail: [email protected].

ACS Paragon Plus Environment


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

Allosteric coupling deduced from network rigidity – Pfleger et al.

Page 2 of 46

2

Abstract Allostery describes the functional coupling between sites in biomolecules. Recently, the role of changes in protein dynamics for allosteric communication has been highlighted. A quantitative and predictive description of allostery is fundamental for understanding biological processes. Here, we integrate an ensemble-based perturbation approach with the analysis of biomolecular rigidity and flexibility to construct a model of dynamic allostery. Our model by definition excludes the possibility of conformational changes, evaluates static, not dynamic, properties of molecular systems, and describes allosteric effects due to ligand binding in terms of a novel free energy measure. We validated our model on three distinct biomolecular systems, eglin c, protein tyrosine phosphatase 1B, and the lymphocyte functionassociated antigen 1 domain. In all cases, it successfully identified key residues for signal transmission in very good agreement with experiment. It correctly and quantitatively discriminated between positively or negatively cooperative effects for one of the systems. Our model should be a promising tool for the rational discovery of novel allosteric drugs.


Page 3 of 46

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



3

Introduction Allostery is the process by which biomolecules transmit the effect of binding at one site to another, often distal, functional site.1 Allosteric regulation occurs in enzyme activation,2 metabolism regulation,3 transcription control,4 and signal transduction.5 The widespread exploitation of allostery by nature renders a quantitative and predictive description of allostery fundamental for understanding biological processes at the molecular and cellular level.6-7 Likewise, such a description is expected to significantly impact drug discovery.8 The first two models for allostery by Monod, Wyman, and Changeux (MWC)9 as well as Koshland, Nemethy, and Filmer (KNF)10 dominated the field for decades and centered on the importance of conformational change between two well-defined structural end states.11 However, both models are phenomenological.11-12 For gaining insights into how biomolecules enable allosteric communication, further model developments were required13-14 that led to the “structural view” of allosteric mechanisms.15-17 This view culminated in some approaches positing the existence of conserved pathways of allosteric signal transmission between the sites.12 Accordingly, a number of computational techniques, including mapping of residue networks by evolutionary, topology, and simulation analyses, have been developed for identifying such pathways at the atomistic level.18-19 Noteworthy, the MWC model already posited an equilibrium shift between the two structural end states upon ligand binding, and the KNF model involved an induced fit of a binding site20 enabled by the inherent flexibility of proteins.10 Thus, both models touched on the dynamic nature of biomolecules in the context of conformational changes. Later, the role of changes in protein dynamics for allosteric communication in the absence of conformational changes was proposed by Cooper and Dryden.21 Based on a statistical thermodynamics analysis of ligand binding, they showed that allosteric communication can arise from changes in frequencies and amplitudes of biomolecular thermal fluctuations, and that this “dynamic allostery” is primarily an entropic effect. This view of allostery has gained much attention,2225

fueled by the development of experimental techniques that allow observation of the role of

dynamics in transitions at the residue level25-28 and a theoretical underpinning of allostery in terms of the free energy landscape of a biomolecule.29-32 According to this framework, all possible conformations of a protein are sampled, and binding of an allosteric ligand shifts the population by stabilizing a certain conformation. Based on this, the ensemble allosteric model33 proposes that allosteric mechanisms may be more statistical, and less deterministic, than the classical models suggest,12 and that it is the relative stability of different



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


Page 4 of 46

4

biomolecular states that determines the type and magnitude of coupling between both sites.12, 33

Here, we integrate an ensemble-based perturbation approach with the analysis of biomolecular rigidity and flexibility to construct a model of dynamic allostery. We apply our model to two systems encompassing ligand-based K- (i.e., the allosteric effector alters the affinity of the substrate) and V- (i.e., the allosteric effector alters the enzyme activity but does not change the affinity of the substrate) type allostery, as well as to a system showing allosteric changes in dynamics from (protein surface) mutations that do not affect the overall structure. We introduce a free energy quantity based on the mechanical stability of the biomolecule and show that it predicts coupling in all three systems in a manner consistent with experiments. By definition, our model excludes any conformational changes of the biomolecule upon perturbation. Hence, the observed allosteric communication must arise from changes in the width of the distributions of biomolecular states only, and so is entropic in nature.


Page 5 of 46

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



5

Overall strategy and theory The aspect of dynamics as a carrier for allosteric signaling between two distant binding sites s1 and s2 was stressed by Cooper and Dryden.21 The cooperative free energy, ∆∆G (eq 1), defined as the difference between the free energy of binding to both binding sites, ∆Gs1/s2, and the free energies of binding to each binding site, ∆Gs1 and ∆Gs2, quantifies the non-independency between both binding processes.

∆∆G = ∆Gs1/s2 – (∆Gs1 + ∆Gs2)

(1)

By separating the free energies into their entropy and enthalpy components, Cooper and Dryden concluded that allosteric cooperativity effects could be primarily entropic, and, thus, happen in the absence of conformational changes.21 In entropy-driven allostery, it follows that, in positive cooperativity, binding of the first ligand results in a major loss of entropy, and thereby lowers the entropic costs for a subsequent binding process. In contrast, if binding of the first ligand does not quench all dynamics in the biomolecule, the major loss of entropy can occur during the second binding process, resulting in negative cooperativity.22 The binding of cAMP to the dimeric catabolite activator protein (CAP) is a prominent example for such a negative cooperativity: NMR experiments revealed that binding of the first cAMP increases the entropic penalty for subsequent cAMP binding without changing the overall structure of CAP.23 However, the pure entropic nature of this allosteric regulation has recently been questioned.34 By means of rigidity theory,35 biomolecules can be decomposed into flexible and rigid regions by determining the number and spatial distribution of the internal independent degrees of freedom (also termed floppy modes).36-37 Floppy modes are equivalent to zero-frequency vibrational modes of a system.36 Recently, we showed that changes of vibrational entropy upon binding to biomolecules can be estimated from changes in the variation of the number of floppy modes with respect to the number of non-covalent bonds in the molecular systems.38 Monitoring the change in the number of floppy modes upon ligand binding should thus provide a means to quantify dynamic allostery as introduced by Cooper and Dryden.21 This is illustrated in Figure 1A22: Rigidifying the network due to ligand binding causes a loss of floppy modes. At the rigidity percolation threshold, a phase transition occurs at which the network loses (almost) all of its floppy modes, and rigidity starts percolating through the network.36,

39

At the threshold, the network is most sensitive to changes, i.e., adding or

removing a few constraints, e.g., due to ligand binding, causes a rigid or flexible state. Away ACS Paragon Plus Environment


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

Page 6 of 46


6

from the threshold, additional or fewer constraints have no or only minor effects on the network stability.38 The binding of an allosteric activator rigidifies the network beyond the percolation threshold, losing most of its floppy modes, such that a second ligand binding process incurs a lower entropic cost than without the activator (Figure 1A). In the case of an allosteric inhibitor, the stabilization upon binding does not exceed the percolation threshold. However, the second ligand binding process crosses the threshold, which leads to a marked loss of floppy modes, incurring a higher entropic cost than without the inhibitor (Figure 1A). Therefore, by analyzing the change in biomolecular rigidity and flexibility40 upon ligand binding, we aim at detecting effects that are linked to entropy changes of the biomolecule and that may lead to dynamic allostery.41-42 For the rigidity analysis, a biomolecule is modeled as a constraint network. Atoms are represented as bodies, and covalent and non-covalent interactions (hydrogen bonds, salt bridges, and hydrophobic tethers) are defined as sets of bars (constraints).39,

41

A fast

combinatorial algorithm, the pebble game,43 then counts the bond rotational degrees of freedom and floppy modes in the constraint network. The theory underlying this approach is rigorous,44 and results from rigidity analyses have been successfully compared with those from experiments and other computational approaches.45-53 Ensemble-based perturbation approach. To detect changes in biomolecular rigidity and flexibility upon ligand binding, we analyze ensembles of the biomolecule’s bound and unbound states in terms of a perturbation approach (Figure 1B). First, an ensemble of network topologies of the biomolecule bound to the allosteric effector, and eventually the orthosteric ligand, is generated from a conformational ensemble obtained by molecular dynamics (MD) simulations,46, 51 constituting the ground state. Then, (a) perturbed state(s) is (are) obtained by removing the covalent and non-covalent constraints associated with the allosteric effector and/or the orthosteric ligand from each network topology of the ground state. Finally, results from rigidity analyses of the ground and perturbed states are post-processed to predict putative residues involved in allosteric signaling and/or the type and magnitude of the cooperativity in binding (see below). Note that, by definition, a change of the biomolecular conformation between the bound and unbound states is excluded in our perturbation approach. Therefore, any observed changes in the biomolecular rigidity and flexibility must arise solely from local changes in the network topology due to the uncoupling of ligand(s). Free energy quantity derived from mechanical stability. Biomolecules generally display a hierarchy of structural stability reflecting the modularity of their structure.49 In order to


Page 7 of 46

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



7

identify this hierarchy, a “constraint dilution trajectory” of network states {σ} is analyzed, obtained from an initial network topology by successively removing non-covalent constraints.39, 45, 49, 52, 54 Here, hydrogen bond constraints (including salt bridges) are removed in the order of increasing strength39, 48, 54 such that for network state σ only those hydrogen bonds are retained that have an energy EHB ≤ Ecut(σ). Figure 1C displays a succession of four network states of the lymphocyte function-associated antigen 1 (LFA-1), showing how the network switches from an overall rigid to a floppy state. The hydrogen bond energy EHB is determined from an empirical energy function55 used previously.39, 45, 47, 51-53, 56 In a next step, neighbor stability maps rcij, neighbor are derived that contain information accumulated over all network states {σ} along the trajectory as to which regions of the network are (locally) mechanically stable at a given state σ, and which are not. More specifically, rcij, neighbor is set to Ecut when the “rigid contact” between residue pairs R{i, j} is lost (Figure 1C); a rigid contact is present as long as the two residues belong to the same rigid cluster c of the set of rigid clusters C E .57 In addition, as performed previously,53 in order to cut

exclude pairs of structurally non-neighboring residues, we focus on short-range rigid contacts by demanding that at least one pair of heavy atoms of the residue pair R{i, j}, A{k ∈ i, l ∈ j}, is separated by d ≤ 4.5 Å58 (eq 2).

(

)

rcij ,neighbor = min Ecut ∃c ∈ C Ecut : (Ri ∧ R j ∈ c ) ∧ (d (A{∃k∈i ,∃l∈ j} ) ≤ 4.5 Å )

(2)

The stability of a contact rcij, neighbor thus relates to the local stability in the network. Taken together, the local stabilities of the residue–residue contacts display the hierarchy of stability in the biomolecule. The sum over all rigid contacts (eq 3)

E CNA =

n

n

i

j >i

∑ ∑ rc

(3)

ij , neighbor

represents the chemical potential energy due to non-covalent bonding, obtained from the coarse-grained, residue-wise network representation of the underlying protein structure.53 Recently, ECNA has been used as a proxy for the melting enthalpy of a protein and, as expected,59 correlates with the protein’s melting temperature.53 The difference in the chemical potential energy ∆ECNA = ECNA,perturbed - ECNA,ground reflects the change in biomolecular stability due to the removal of a ligand or a mutation. Note that contributions from removing non-



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


Page 8 of 46

8

covalent interactions between the ligand or mutated residue and the protein do not contribute to ∆ECNA. The perturbation of a ground state network by removing constraints between ligand and protein is localized and small. For the examples investigated in this study, less than 2.4% of the total number of constraints in a network topology are removed (Tables S1-3). Hence, the ensembles of network topologies of ground and perturbed states are highly similar. This warrants60 applying a one-step free energy perturbation approach61-62 as an approximation63-65 to compute the free energy ∆GCNA (ground → perturbed state) associated with the change in biomolecular stability due to the removal of a ligand or a mutated residue (eq 4)

 ∆E  ∆GCNA = −k BT ⋅ ln exp − CNA   k BT 

(4) ground

where kB is the Boltzmann constant, and the temperature T is set to 300 K. 〈…〉ground denotes averaging over the ground state ensemble (for a derivation of eq 4 see Supporting information). Identification of pathways for allosteric signaling. A per-residue decomposition of ∆GCNA allows identifying to what extent single residues contribute to the allosteric signaling. While the total free energy is a state function and as such is independent of the pathway used to calculate it, free energy components, in general, are not.66-68 The decomposition scheme proposed here, however, depends only on the endpoints (ground and perturbed states), i.e. individual components do not depend on the choice of a specific pathway.69 Furthermore, the residue-based description follows a “natural” choice and provides a decomposition, which helps understanding the influence of substructural elements on the allosteric signaling. We apply a linear response approximation70 to obtain the per-residue decomposition (eq 5)

(

∆Gi , CNA = α Eiperturbed − Eiground , CNA , CNA

)

(5)

where Ei, CNA is the chemical potential energy of residue i obtained by summation over all n short-range rigid contacts the residue is involved in (eq 6)

Ei , CNA =

1 n ∑ rcij , neighbor 2 j ≠i

(6),


Page 9 of 46

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



9

and 〈…〉 indicates averaging over the respective ensembles of network topologies. Approaches to compute free energies based on a linear response approximation have been successfully applied before in the context of binding of two molecules.70-77 α has been in general determined empirically there,72, 78-79 except for when purely electrostatic interactions were considered, for which a value of α = ½ is usually used80-81 according to Marcus’ electron transfer theory82.

Methods Constraint network analysis Rigidity analyses were performed using the CNA software package.83 CNA was applied on ensembles of network topologies generated from conformational ensembles provided as input (see Supporting information).46,

51

Average stability characteristics are calculated by

constraint counting on each topology in the ensemble.43 Each network of covalent and noncovalent (hydrogen bonds including salt bridges and hydrophobic tethers) constraints was constructed with the FIRST (Floppy Inclusions and Rigid Substructure Topography) software (version 6.2),41 to which CNA is a front and back end. The strength of hydrogen bonds (including salt bridges) were assigned by the energy EHB computed by FIRST.55 Hydrophobic interactions between carbon or sulfur atoms were taken into account if the distance between these atoms was less than the sum of their van der Waals radii (C: 1.7 Å, S: 1.8 Å) plus Dcut = 0.25 Å.39 Fluor atoms of CF3 groups in the allosteric inhibitor of PDB ID 3bqm84 were modeled as carbons to account for the increased bulky and hydrophobic character of CF3 groups.85

Allosteric signaling through constraint networks Stability maps57 as derived from rigidity analysis were used to generate an undirected graph G = (V, E). In such a graph, nodes V represent the residues, and edges E represent the pairwise perturbation free energy ∆Gij between residue i and j given by eq 7

(

∆Gij = β rcijperturbed − rcijground

)


(7),


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

Page 10 of 46


10

where 〈rcij 〉 is the average strength of a rigid contact (eq 2) between the two residues i and j, and β = ½ α, with α being the fitting parameter from eq 5.82 The graph embedding is optimized based on the pairwise Cα atom distances in the PDB structure used as input for the perturbation approach. For a graph representation, we calculated the current flow betweenness centrality CCB(e) of each edge e to quantify the relative importance of e with respect to the information flow through the graph, starting from a source s ∈ V to a sink t ∈ V. Compared to the widely used (shortest-path) betweenness centrality86-87, which assumes that information flows along shortest paths between s and t only and does not split, CCB(e) builds on the same intuition but lets information flow and split like current in an electrical network.88 Here, a supply bst of unit strength enters at s and leaves at t, leading to a st-current. Given bst and the conductance or strength of each edge, in our case represented by ∆Gij (eq. 7) with (i,j) ∈ E, a current flowing along the directed edge is computed via the Laplacian (also called Kirchhoff) matrix of G (note that the actual choice of orientation of is of no importance; see ref. 88 for further details on how to compute ). CCB(e) is then computed for all combinations of s and t from the throughputs =

| | due to a st-current (eq 8)88-89

CC B (e ) =

1 τ st (e) (n − 1)(n − 2 ) s∑ ≠ t∈V

(8);

= ∥ ∥ and the term in front of the summation sign is a normalizing constant with respect to the size of G. Edges are part of a pathway, if their CCB(e) values are significantly (pC

CB (e)

< 0.05)

different from zero. By doing so, more than one pathway can be detected. We define the pathway enclosing the maximal number of nodes as the most pronounced one, and discard all others. For the graph construction and calculation of the current flow betweenness centrality, the NetworkX extension for Python90 was used. The underlying algorithm to calculate the current flow betweenness centrality is described in refs.91-92. For visualization of the graph, we used the software Visone.93


Page 11 of 46

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



11

Results Allosteric effects of mutations in eglin c The small serine protease inhibitor eglin c is an ideal test case for validation of our model of dynamic allostery because experimental data is available from (I) chemical denaturation experiments providing free energies of unfolding for 13 single and double mutation variants and (II) NMR experiments for three single mutation variants revealing that these mutations do not influence the eglin c structure but rather change the dynamics of other residues.94-95 In all single and double mutation variants, the residues were mutated to alanine. Our perturbation approach was applied to three ensembles of 2,500 conformations each, extracted from independent MD trajectories starting from the wild type structure of eglin c (see Supporting information). In all cases, the eglin c structures are stable during the course of the MD simulations (Figure S1). Computed relative free energies associated with the change in biomolecular stability agree very well with experimental data. To validate our model, we calculated the perturbation free energies ∆GCNA (eq 4) of 13 single and double mutation variants. For generation of the perturbed state from the wild type conformations, side chain atoms of the respective residues beyond the Cβ atom were deleted, leaving the rest of the structure unchanged. As reference, we used differences of free energies of unfolding (∆∆G) between wild type and variants from chemical denaturation experiments94-95 ranging from 0.16 to 3.74 kcal mol-1, which relate to a destabilization of the variants compared to wildtype. The comparison of free energies computed for all pairs of the three independent ensembles shows perfect correlations (R2 ≥ 0.98, p ≤ 4.0·10-10), indicating the robustness of our approach. For each variant, the free energies were averaged over the three ensembles (Table S4), resulting in standard errors of the mean (SEM, n = 3, eq S11) of ≤ 0.02 kcal mol-1. The computed and experimental free energies yield a good and significant correlation (R2 = 0.80, 95% confidence interval: 0.67 < R2 < 0.97, p = 0.0001) (Fig 2a). The low slope (m = 0.12) of the correlation line may be due to eglin c being small; the system-size dependence of results from rigidity analysis has been noted before in the context of thermostability predictions for proteins.49, 51, 54, 56 Furthermore, the underprediction observed for ∆GCNA with respect to ∆∆G may result from our strict definition for introducing a perturbation, by which a change of the biomolecular conformation between the ground and perturbed states is excluded. That way, only local changes in the network topology due to the uncoupling of side chains contribute to ∆ECNA (eq 3), and thus ∆GCNA, but not (small) structural relaxations that may occur upon a mutation. The V63A ACS Paragon Plus Environment


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


Page 12 of 46

12

variant was not considered in the correlation analysis because NMR experiments demonstrated that this mutation results only in small changes in eglin c dynamics95 and that other factors must give rise to its destabilizing effect. Residues predicted to have a significantly altered structural stability accurately agree with dynamically coupled residues as identified by NMR. Based on the per-residue free energy (eq 5), it is possible to examine how much each residue’s structural stability is affected upon perturbation of the ground state, and how residues that are spatially separated influence each other. NMR studies on the V14A, V34A, and V54A variants showed contiguous pathways of dynamically coupled residues up to 15 Å away from the mutation site.94-95 None of the three mutations shows a significant effect on the tertiary structure of eglin c, hence the response appears to be purely dynamics-driven.95 We focused our attention on the V34A variant because it showed the most far-reaching dynamic response in the NMR experiment. Upon perturbation of wild type eglin c, the network topologies lose on average 1.5-1.7 (= 4.1-4.6% of all) non-covalent constraints (Table S1). The per-residue free energy ∆Gi,CNA (eq 5) was calculated from differences between neighbor stability maps of the ground (wild type) and perturbed (variant) states (Figure S2A). Across all three ensembles, the calculated ∆Gi,CNA are highly correlated (R2 ≥ 0.99, p ≤ 2.2·10-16), and, consistently, > 84% of the residues differ < 0.05 kcal mol-1 in ∆Gi,CNA. The average per-residue free energies over the three ensembles have a SEM < 0.07 kcal mol-1 (n = 3). The free energies associated with the change in biomolecular stability (eq 4) and the sum of the per-residue free energies (eq 5) yield a significant and very good linear correlation (R2 = 0.98, p < 0.0001, Figure S3) over all single and double mutation variants except V63A, indicating that the linear response approximation of eq 5 is applicable in this case. Figure S3 also reveals for the fitting parameter α in eq 5 a value of 0.02. Almost all residues close to the mutation site were affected by the perturbation but only a few residues showed larger ∆Gi,CNA (Figure 2B). The highest ∆Gi,CNA values were found for residues that form direct contacts with V34; these are residues Y35, F36, V52, R53, and V54 with ∆Gi,CNA ≥ 0.61 kcal mol-1. When considering only residues with ∆Gi,CNA > 0.2 kcal mol-1 (i.e., with a ∆Gi,CNA value of at least 24% of the maximal ∆Gi,CNA), a contiguous pathway of coupled residues was revealed (Figure 2C). Striking examples for the long-range character of altered rigidity upon V34A mutation are residues V62 and K16, which are 15.9 and 15.6 Å away from the perturbation site, respectively. Notably, both residues also show altered dynamics in NMR experiments (Figure 2D). On the opposite side of eglin c, the connection with the proteinase binding loop via residue V43 (13.4 Å away) is missed, and none of the


Page 13 of 46

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



13

loop residues show significant non-zero ∆Gi,CNA. However, a strong influence is observed for R53, which is in contact with V43, indicating that the mutation effect may be propagated towards the proteinase binding loop (Figure 2C). Overall, the identified pathway is in good agreement with results from NMR experiments (Figure 2C and D): When evaluated in a binary classification system (see Supporting information), our results show a discrimination between pathway and non-pathway residues with a sensitivity of 79 ± 7%, a specificity of 75 ± 2%, and an accuracy of 76 ± 1% (Figure 2E). Note that NMR data on eglin c dynamics is available only for certain residues.94 Thus, the pathway of dynamically coupled residues identified by NMR represents a subset of the true pathway. E.g., data on residue dynamics is not available for R53 and V63. However, both residues are considered to be on the pathway as they connect clusters of residues with dynamic response.95 Results based on altered structural stability do not suffer from this shortcoming as structural stability is assessed for each residue. Finally, information on sequence conservation or changes in the residue-wise solvent accessible surface area (SASA) for the V34A variant perform inferior in the pathway prediction compared to that based on ∆Gi,CNA (see Figure S4). Additionally, we analyzed the V54A and V14A variants. For V54A, we found a contiguous pathway similar to that of V34A (Figure S2B-E). While the sensitivity is lower (57 ± 3%), a specificity of 79 ± 1% and an accuracy of 76 ± 1% demonstrate that we can reliably identify residues whose dynamics are not influenced by the mutation. In the case of V14A, only small changes of structural stability are identified (Figure S2F-H); this variant has also been reported to have the lowest impact on dynamics in the NMR experiments.94 In summary, the results for eglin c demonstrate that computed relative free energies associated with the change in biomolecular stability agree very well with experimental data, and that residues predicted to have a significantly altered structural stability by our approach accurately agree with dynamically coupled residues as identified by NMR, even those ~15 Å away from the mutation site.

V-type allostery in PTP1B In a second case study, the perturbation approach was applied to protein tyrosine phosphatase 1B (PTP1B), which is controlled by a V-type allosteric mechanism, i.e. allosteric effectors antagonize the enzyme activity without changing the affinity for substrate binding.96 The orthosteric site, composed of an active site and a neighboring non-catalytic site,97 is ~20 Å away from the allosteric site.96 The protein shows only small structural changes between effector-bound96 and substrate-bound97 state (1.6 Å non-hydrogen atom rmsd overall). A



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


Page 14 of 46

14

comparison of crystal structures in the allosteric inhibitor-bound and substrate-bound states revealed that helix α7* next to the allosteric site and the catalytically important WPD loop (residues T177-P185) may be involved in allosteric signal transmission.97-99 However, exactly how these structural regions are coupled for transmission is not fully understood.96 Similar to eglin c, the ground state ensembles with 2,500 conformations each were extracted from three independent MD trajectories starting from the allosteric inhibitor-bound structure of PTP1B (see Supporting information). The PTP1B structure was generally stable along the trajectories of the MD simulations (Figure S5A). Network topologies of the perturbed state were generated by removing all constraints associated with the allosteric inhibitor from each network topology of the ground state. Perturbation of PTP1B results in altered structural stability up to 21 Å away from the allosteric site in functionally relevant regions. Upon perturbation, the network topologies lose on average 2.9-3.6 (= 0.9-1.1% of all) hydrogen bond constraints and 15.1-17.0 (= 8.9-10.4% of all) hydrophobic tether constraints (Table S2). The per-residue ∆Gi,CNA (eq 5) show a change of structural stability of almost all residues upon perturbation (Figure S6). Across all three ensembles, the calculated ∆Gi,CNA are highly correlated (R2 = 0.87-0.92, p < 0.001), and, consistently, 70% of the residues differ ≤ 0.05 kcal mol-1 in ∆Gi,CNA. Averaging the perresidue ∆Gi,CNA over the three ensembles results in a SEM ≤ 0.28 kcal mol-1 (n = 3; Figure 3A). Analogous to eglin c, we only considered residues with ∆Gi,CNA > 0.2 kcal mol-1 for further analysis. The largest per-residue ∆Gi,CNA upon removal of the allosteric effector were found for regions enclosing helices α3, α6, and α7 (Figure 3A), with maximal ∆Gi,CNA = 3.0 kcal mol-1 for residue S187, which connects residues of helix α3 and the WPD loop. By mapping residues with ∆Gi,CNA > 0.2 kcal mol-1 onto the structure of PTP1B (Figure 3B), we found that residues up to ~21 Å away from the allosteric site are influenced by removal of the allosteric inhibitor. This result clearly demonstrates the long-range character of rigidity percolation (Figure 3B).100 Remarkably, affected residues are not widely dispersed but rather show a distinct connection to active site residues (Y20, I215, G218-S222, R254, M258, G259, I261, Q262, Q266; marked in green in Figure 3B) and residues in the hinge region of the WPD loop (T177-P180, G183-P185; marked in orange in Figure 3B). The strongly altered structural stability in helices α3 and α6 upon removal of the effector indicate their function as connectors between the allosteric site involving helix α7 and the active site/WPD loop. Remarkably, while changes in the structural stability affect the orthosteric site, 57% of the PTB1B residues do not feel the perturbation at the allosteric site, indicating a highly specific effect. To ensure that the observation of altered structural stability of PTB1B residues is of


Page 15 of 46

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



15

significance, we also calculated the packing density for each residue over the ground state ensembles based on Voronoi volumes (see Supporting information). The latter procedure follows the assumption that the propagation of allosteric effects is more efficient in tightly packed atomic environments.101 No correlation (R2 ≤ 0.01) was observed between computed ∆Gi,CNA and per-residue packing densities (Figure S7A), demonstrating that changes in the structural stability upon perturbation of the allosteric site occur independently from densely packed regions in PTP1B. Identification of pathways of allosteric communication in PTP1B by graph analysis. Above analyses identify residues whose structural stability is affected by the removal of the allosteric inhibitor. However, possible pathways of the transmitted allosteric signal have not yet been revealed. To address this question, we first generated an undirected graph consisting of nodes representing PTP1B residues with ∆Gi,CNA > 0.2 kcal mol-1 and edges whose weights are defined by the pairwise free energy components ∆Gij (eq 7) (Figure 3C). We further stress-minimized the graph with respect to the distance between nodes so that nodes close in the PTB1B structure are also close in the graph. The generated graph shows that no residue of the allosteric site shares an edge with a residue of either the orthosteric site or the WPD loop: Each path from the allosteric site to the WPD loop and the orthosteric site has a length of at least two or three edges, respectively, indicating that the allosteric effects cannot arise from only local changes in the constraint network. In addition, residues in helix α7* show ∆Gi,CNA values of up to 1.68 kcal mol-1. While this finding might not appear surprising due to the small distance between helix α7* and the allosteric inhibitor, it stresses the effect of the inhibitor to prevent the allosteric coupling between helix α7* and the WPD loop by interrupting the interactions between helix α7* and the two helices associated with the closing of the loop, α3 and α6.96 Next, we computed the current flow betweenness centrality (refs.

91-92

, eq 8) to identify

edges through which most of the allosteric signal is transmitted. From this, we can identify pathways by which effects on structural stability are channeled through the protein that result from removing the allosteric inhibitor. Remarkably, only ≈ 5% of all edges have centrality values significantly (pCB(e) < 0.05) different from zero, indicating that only a small subset of edges is important for the allosteric information flow. Considering only those edges reveals a narrow pathway originating from the allosteric site (L192 and P188 in helix α3) and extending to the orthosteric site (R254, I261, and Q266) and the WPD loop (W179) (Figure 3C), with the signal transmission occurring via F269 (helix α6) and T224 (helix α4). Branching at F269, the edge to F191 shows a significant betweenness centrality value. In



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


Page 16 of 46

16

summary, our analysis reveals that removal of the allosteric inhibitor leads to a channeled structural destabilization of F191, which has been implicated in the closure of the WPD loop.96 The identified pathway of allosteric signaling in PTP1B forms a condensed cluster (“small narrow pathway”) and connects to residues in the orthosteric site and other functionally important residues.

K-type allostery in LFA-1 We applied the perturbation approach to the lymphocyte function-associated antigen 1 (LFA-1) domain, which is part of β2 integrin and binds to intercellular adhesion molecules (ICAM), including ICAM-1.102 LFA-1 is controlled by a K-type allosteric mechanism, i.e., allosteric effectors either inhibit103 or activate104 complex formation of LFA-1 and ICAM-1 by binding to the allosteric pocket. The binding site of ICAM-1 is referred to as orthosteric site. The crystal structures of both the inhibitor- and activator-bound state show minor differences in the overall structure and orthosteric site only (1.7 and 1.2 Å non-hydrogen atom rmsd, respectively). Accordingly, the exact details of how allosteric inhibitors and activators exert their opposite effects have remained elusive.103-104 As before, the ground state ensembles with 2,500 conformations each were extracted from three independent MD trajectories starting from the effector-bound structures of LFA-1 (see Supporting information). The LFA-1 structure was stable along the trajectories of the MD simulations (Figure S5B, C). Each network topology of the ground state ensemble was perturbed by removing the constraints associated with the allosteric effector. Upon perturbation of the inhibitor-bound state, the network topologies lose on average 1.8-2.2 (= 0.9-1.1% of all) hydrogen bond constraints and 11.1-13.0 (= 9.5-11.0% of all) hydrophobic tether constraints (Table S3). Upon perturbation of the activator-bound state, the network topologies lose on average 2.6-3.1 (1.3-1.5% of all) hydrogen bond constraints and 11.4-12.1 (9.7-10.4% of all) hydrophobic tether constraints (Table S3). Perturbation of LFA-1 results in long-range altered structural stability different in efficacy for the allosteric inhibitor- and activator-bound states. Across all three ensembles, the calculated ∆Gi,CNA are highly correlated for both the inhibitor (R2 = 0.88-0.97, p < 0.001) and the activator (R2 = 0.92, p < 0.001); consistently, 51% (inhibitor) and 53% (activator) of the residues differ ≤ 0.05 kcal mol-1 in ∆Gi,CNA (Figure S8A, C). Averaging ∆Gi,CNA over the three ensembles results in a SEM ≤ 0.32 kcal mol-1 (n = 3; Figure 4A and B). As before, we only considered residues with ∆Gi,CNA > 0.2 kcal mol-1 for further analysis. When mapped onto the LFA-1 structure, these residues are located in very similar regions for both allosteric inhibition and activation (Figure S8B, D). In particular, β-strands (except for β2’) in the core ACS Paragon Plus Environment

Page 17 of 46

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



17

region show the largest per-residue ∆Gi,CNA upon removal of the constraints associated with the allosteric effectors (Figure S8B, D). In contrast, the opposite region enclosing helices α1, α2, α3, α4, α5, and α6 is much less altered. The ∆Gi,CNA do not correlate with the packing density of residues, indicating that our observations are of significance (Figure S7B, C). Large ∆Gi,CNA have also been found for residues L289, D290, and L295, which are located close to the functionally important hinge region between sheet β5 and helix α7. L142 in the orthosteric site is among the most distant (~20 Å away from the allosteric site) residues with altered ∆Gi,CNA and is part of a network of residues important for protein-protein interactions with ICAM.105 Besides L142, residues D137, T238, and D239 in the orthosteric site are also influenced by binding of the allosteric modulator. These findings again demonstrate the longrange aspect of rigidity percolation. In particular, the identification of D137 and D239 in the orthosteric site is of relevance because both residues have been reported to be important in ICAM-1 binding.106 Although similar protein regions are affected by perturbing either the inhibitor- or activator-bound state of LFA-1, the magnitude of the ∆Gi,CNA generally differs between both cases, with larger ∆Gi,CNA observed for the allosteric activator than the inhibitor (Figure 4C). Mapping these differences onto the structure of LFA-1 shows that activator binding results in a more strongly altered structural stability particularly in helix α7, the hinge region next to strand β5, and the region of strands β1-β3 (Figure 4D). The larger change in structural stability can be explained by interactions of the charged aminoalkyl group in the activator with the protein, which is not present in the inhibitor.104 Graph analyses of LFA-1 reveal a uniform allosteric signaling across the structure. To identify by which routes the allosteric signal is transmitted through the protein, we computed the current flow betweenness centrality for each edge of the graphs, as described above for PTP1B (Figure S9). Of all edges, only 3.9% have centrality values significantly larger than zero. These edges primarily connect nodes in the vicinity of the allosteric site, which is in contrast to what we observed for PTP1B. Thus, compared to PTP1B, the allosteric signaling in LFA-1 is not as much channeled but is transmitted more uniformly across the structure. Computed cooperative free energies correctly characterize the type of allosteric modulation. The LFA-1 system provides an ideal test case for our approach whether it can correctly predict the type of allosteric modulation in terms of negative and positive cooperativity (eq 1, Figure 1A) between the binding of the allosteric effector and that of ICAM-1 to LFA-1. To do so, we performed MD simulations each starting from a modelled structure of LFA-1 in complex with ICAM-1 and either the allosteric inhibitor84 or the activator104 (see Supporting information). Both modeled complexes are stable over the course



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


Page 18 of 46

18

of the trajectories of the MD simulations (Figure S10). Structural ensembles of the respective LFA-1effector/ICAM-1 complexes were extracted from the trajectories, which were used to generate perturbed ensembles of apo LFA-1 (LFA-1apo), LFA-1 bound to the effector (LFA-1effector), and LFA-1 bound to ICAM-1 (LFA-1ICAM-1). The cooperative free energy (eq 1) between the respective two binding processes was computed from ∆GCNA values for LFA-1effector/ICAM-1, LFA-1effector, and LFA-1ICAM-1 versus LFA-1apo, respectively. For LFA-1inhibitor and LFA-1ICAM-1, the free energies ∆GCNA associated with a change in biomolecular stability due to the absence of a binding partner (eq 4, Table 1) reveal that inhibitor binding affects LFA-1 stability markedly less than ICAM-1 binding. Together with ∆GCNA associated with the binding of both inhibitor and ICAM-1, this results in a cooperative free energy ∆∆G (eq 1) of 3.5 kcal mol-1 (Table 1), indicating that the entropic burden of ICAM-1 binding to LFA-1 is higher if the inhibitor is already bound, hence resulting in a negative cooperativity in agreement with experiment.107-108 In turn, for LFA-1activator, LFA-1ICAM-1, and LFA-1activator/ICAM-1, the ∆GCNA values (Table 1) result in a cooperative free energy ∆∆G of -1.4 kcal mol-1, indicating that the entropic burden of ICAM-1 binding is lower if the activator is already bound, hence resulting in a positive cooperativity in agreement with experiment.104 As Kd values are available for the dissociation of ICAM-1 from the LFA-1/activator complex and the dissociation of the activator from LFA-1, we computed Gibbs free energies for the states ∆Gs1, ∆Gs2 and ∆Gs1/s2 (Table 1). According to eq 1, this yields a cooperative free energy of -5.9 kcal mol-1 (Table 1). The lower cooperative free energy from experiment might result from an additional enthalpic effect due to small conformational changes upon activation of LFA-1, which are neglected in our perturbation approach. In summary, the predicted free energies of cooperativity derived from the perturbation approach allow us to correctly distinguish between non-additive structural stabilizations of LFA-1 due to inhibitor/activator and ICAM-1 binding, in agreement with the underlying mechanisms of negative and positive cooperativity in LFA-1.

Discussion In this study, we integrated an ensemble-based perturbation approach with the analysis of biomolecular rigidity and flexibility to construct a model of dynamic allostery. The characteristics of our model are that I) it excludes any conformational changes of the biomolecule upon perturbation by construction, such that the observed allosteric communication is entropic in nature,

109

II) static rather than dynamic properties of the


Page 19 of 46

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



19

molecular systems are evaluated, circumventing the challenge to separate the signal from the noise 111

110

and the computational burden of intensive sampling of a biomolecule’s microstates

that occur when analyzing simultaneous dynamic processes at distant sites, III) the

influence on biomolecular stability due to removal of a ligand or introduction of a mutation, and how this effect percolates through the structure to influence a distant region, is evaluated in terms of the free energy measure ∆GCNA, introduced here for the first time, and IV) the model is able to distinguish between effects due to an allosteric activator or inhibitor from the computed cooperative free energy ∆∆GCNA. Our model was validated on three biomolecular systems with distinct properties. First, for eglin c, computed relative free energies associated with the change in biomolecular stability agreed very well with experimental data for 13 single- and double mutation variants. This result is remarkable considering that for all mutations either a valine or leucine residue is changed to alanine. Hence, our model is apparently sensitive enough to detect differences in the influence of closely related mutations on structural stability depending on the mutation’s location in the constraint network. Furthermore, residues predicted to have a significantly altered structural stability upon introduction of a mutation accurately agreed with dynamically coupled residues as identified by NMR,94-95 demonstrating the predictive power of the residue-wise ∆Gi,CNA based on a linear response approximation. Notably, this approach can identify residues involved in allosteric signal transmission all at once upon perturbation of the system, avoiding the computational burden to probe each residue separately.112 Second, for PTP1B with a V-type allosteric mechanism, our model identified altered structural stability upon removal of the allosteric inhibitor up to 21 Å away from the allosteric site in functionally relevant regions, highlighting the long-range character of rigidity percolation.100 In addition, we were able to detect a channeled structural destabilization of F191 upon inhibitor removal, which is remarkable because F191 in the inactive state is fixed at a site that W179 of the WPD loop must enter during the closure of this loop (Figure 3D).96 Furthermore, the identification of residues R254, I261, and Q266 in the orthosteric site upon inhibitor removal suggests an influence on a catalytically important network of residues, and, hence, an influence on the PTP1B-catalyzed reaction, even if the affinity for substrate binding is not affected by binding of the allosteric inhibitor.96 Support for this suggestion comes from the finding that Q266 needs to be involved in a water-mediated hydrogen bonding network as a prerequisite for formation of an intermediate state during catalysis (Figure 3E).113 Third, for LFA-1 with a K-type allosteric mechanism, our model revealed that structural stability is altered at long-range and that this alteration was similar with respect to affected



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


Page 20 of 46

20

protein regions, but showed a larger magnitude for the allosteric activator-bound state compared to the inhibitor-bound state (Figure 1A). Our findings extended those of previous studies based on the analysis of crystal structures that proposed a local effect by which complex formation of LFA-1 and ICAM molecules is inhibited108, 114: The studies suggested that binding of an allosteric inhibitor holds the structure of LFA-1 in the inactive state by freezing the C-terminal helix neighboring the allosteric site.115 Furthermore, computed cooperative free energies correctly characterized the type of allosteric modulation. Addressing the question whether statistical or deterministic effects dominate in allosteric mechanisms,12, 33, 116-117 we believe the detailed analysis of only three systems precludes us from making a general statement. In the case of PTP1B, our graph-based analysis of allosteric signal flow identified a rather narrow pathway, extending to a few residues likely involved in catalysis, thus supporting a “deterministic view”. In contrast, in the case of LFA-1, a rather broad influence in terms of multiple weakly coupled residues was found, supporting a “statistical view”. It will be interesting to investigate in future studies to what extent these findings are related to the occurrence of V- versus K-type allostery. In addition, such findings may be of relevance to the question of network stability and robustness of the allosteric coupling pathway in view of potential mutations,118 whereby less detrimental effects can be expected in the case of a broad pathway. From a structural point of view, while a number of studies concluded that the allostery they observe takes place through sheer dynamics based on Cooper and Dryden’s concept of allosteric communication in the absence of conformational changes,21 the lack of observing conformational changes does not necessarily mean that there is no conformational change.34 In this context, we note that for both the PTP1B and LFA-1 system investigated here, our analyses are clearly indicative for the influence of changes in the dynamics on allosteric signal transmission, as they exclude per se conformational changes, although for both systems also small but significant conformational changes have been detected experimentally between apo and effector-bound states. With regard to the importance of secondary structures for allosteric signal transmission, our analyses for LFA-1 revealed that β-strands in the core region generally showed the largest per-residue ∆Gi,CNA upon perturbing the system, in contrast to LFA-1’s helical region. This result is not only in good agreement with findings from NMR experiments103 but also highlights that the β-strand motif is apparently more sensitive to perturbation effects than the surrounding helices, which can be understood from analyzing the inherent structural stability of secondary structures by constraint counting.119


Page 21 of 46

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



21

From a methodological point of view, we applied ensemble-based rigidity analysis to overcome the challenge that rigidity analyses of biomolecules are in general sensitive to the structural information used as input.46, 51 A related ensemble-based mechanical perturbation method has been introduced in connection with the Distance Constraint Model (DCM), which simulates the presence of a ligand by adding extra constraints to a certain site in the constraint network.112 Although this approach reveals the mechanical coupling between distant sites, it does not quantify the cooperativity of the binding processes in terms of a cooperative free energy. In this work, we analyzed ensembles of constraint network topologies generated from structural ensembles obtained from MD simulations. In order to improve the computational efficiency of our approach, a method to generate ensembles of constraint network topologies from a single input structure can be used54, essentially modeling the fluctuations of noncovalent constraints rather than the motions of atoms.

Conclusion We presented a robust approach for quantifying the long-range effects of dynamic allostery in biomolecules based on rigidity theory. We correctly discriminated between positively or negatively cooperative effects due to allosteric effector binding, and identified key residues for signal transmission within the biomolecules. In a prospective manner, this approach could give valuable insights for tuning the allosteric effect of already existing effector molecules, exploiting the fact that small structural modifications can also be treated as perturbations. Furthermore, for biomolecules with yet unknown allosteric mechanisms, the response of potential allosteric sites can be probed by adding constraints to the biomolecule’s network representation, thus simulating binding of a hypothetical effector molecule. We believe that our approach should be a promising tool for the rational discovery of novel allosteric drugs.



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


Page 22 of 46

22

Figure captions Figure 1: Allosteric coupling deduced from flexibility and rigidity characteristics. (A) Illustration of the occurrence of positive and negative cooperativity of ligand (blue and black points) binding from entropy changes due to altered structural stability. The vertical dashed line indicates the rigidity percolation threshold where the biomolecule switches from a flexible to a rigid state. (B) Workflow of the perturbation approach. As input, structural ensembles extracted from MD trajectories of the ground state are used. The corresponding ensemble of the perturbed state is obtained by removing all constraints associated to the allosteric effector. Chemical potential energies53 are derived from the neighbor stability maps of both states. Free energy perturbation calculations are applied, and the free energies are related to (changes in) protein stability, allosteric cooperativity, and allosteric pathways. (C) Rigid cluster decomposition during the bond dilution process of LFA-1 (left). The four states show the decay of rigidity where at each state all hydrogen bonds with an energy EHB > Ecut are removed from the network. A rigid contact between two residues Ri and Rj ceases to exist when both residues start to belong to different rigid clusters along the bond dilution process (top right). The energy values Ecut at which the rigid contacts are lost yield the neighbor stability map (bottom right).

Figure 2: Stability changes in eglin c deduced from the perturbation approach. (A) Correlation between predicted ∆GCNA (eq 4) and ∆∆G from chemical denaturation experiments for 13 variants of eglin c.94-95 The SEM of ∆GCNA over three MD trajectories is ≤ 0.02 kcal mol-1. The SEM of experimentally determined free energies are exemplarily shown for V18A/V54A as horizontal error bars and is ≤ 0.13 kcal mol-1 for all ∆∆G. The V63A variant (red) is considered an outlier because it shows no altered structural stability upon in silico perturbation and is reported to show only small changes in the structural dynamics by NMR experiments.95 (B) The histogram shows the per-residue ∆Gi,CNA (eq 5) with the red bar highlighting the site of mutation. The dashed line at 0.2 kcal mol-1 indicates the threshold above which residues are considered perturbed. Residues with ∆Gi,CNA above the threshold are depicted as spheres on the structure of eglin c (C, PDB ID 1cse) and are in good agreement with a continuous dynamic network of residues derived from NMR experiments (D).95 The site of mutation (V34) is shown in red (C and D), blue colors reflect predicted ∆Gi,CNA values in (C) with darker colors indicating larger ∆Gi,CNA values, and white residues in (D) are those that lack NMR data on changes in the dynamics but are included in the network as suggested by the authors of ref.

95

. (E) Results from evaluating the predictions


Page 23 of 46

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



23

versus the experimental NMR data as derived from a binary classification problem: Average specificity (red), sensitivity (yellow), and accuracy (blue) of the predicted pathways from three independent ensembles. Error bars depict the SEM.

Figure 3: Probing the allosteric mechanism in PTP1B. (A) The histogram shows the perresidue ∆Gi,CNA (eq 5). Secondary structure as well as allosteric (red) and orthosteric (green) site information are depicted at the top. (B) Residues with ∆Gi,CNA ≥ 0.2 kcal mol-1 are depicted as spheres on the structure of PTP1B with bound allosteric inhibitor (yellow) (PDB ID 1t4j). The residue color indicates if a residue is part of the allosteric (red), orthosteric (green), WPD loop (orange), or any other site (blue). Darker colors indicate larger ∆Gi,CNA values. The red arrow points to the regions with largest ∆Gi,CNA values (enclosing helices α3, α6, α7, and residue S187). (C) The same information is presented as a stressminimized graph where nodes represent residues, with darker colors indicating a lager ∆Gi,CNA value, and adjacent edges the associated pairwise perturbation energies ∆Gij (eq 7). The graph embedding is optimized with respect to the pairwise Cα distances in PDB ID 1t4j such that the nodes in the 2D representation preserve their relative orientation in 3D. Highlighted residues show the dominant pathways for allosteric communication through the protein as derived from calculating the current flow betweenness centrality for each edge (eq 8). One pathway is mapped onto the PTP1B structure and branches at position F269 to connect the allosteric site (red) with the hinge region of the WPD loop (D, orange) and residues in the orthosteric site (E, green). Additionally, panel (D) shows the aligned closed structure of PTP1B (black, PDB ID 1pty) where the black arrow indicates the closure movement of the WPD loop.

Figure 4: Probing the allosteric mechanism in LFA. The histograms show the per-residue ∆Gi,CNA (eq 5) for the inhibitor- (A) and activator- (B) bound states. The secondary structure as well as allosteric (red) and orthosteric (green) site information are shown at the top, respectively. (C) Differences of per-residue ∆Gi,CNA values between inhibitor- and activatorbound states shown as histogram and (D) mapped onto the structure of LFA-1. Bluish colors indicate a stronger change in ∆Gi,CNA due to the activator, and reddish colors a stronger change in ∆Gi,CNA due to the inhibitor. For differences ≥ 0.2 kcal mol-1, the information is depicted as spheres for the Cα atoms of the LFA-1 structure.



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


Figures

Figure 1


Page 24 of 46

24

Page 25 of 46

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



Figure 2


25


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


Figure 3


Page 26 of 46

26

Page 27 of 46

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



Figure 4


27


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


Page 28 of 46

28

Tables Table 1: Cooperative free energies for characterizing the type of allosteric modulation in LFA-1. System

Inhibitor ∆GCNA1

Activator ∆GCNA1

Activator ∆G2

LFA-1effector/ICAM-1

9.8

4.3

-19.84

LFA-1effector

0.7

1.7

-10.0

LFA-1ICAM-1

5.6

4.0

-3.9

∆∆GCNA3

3.5

-1.4

-5.9

1

Perturbation free energies (eq 4) in kcal mol-1. Gibbs free energies calculated from Kd values reported in ref. 104, in kcal mol-1. 3 Cooperative free energies (eq 1) in kcal mol-1. 4 Sum of ∆G from Kd values for the dissociation of ICAM-1 from the LFA-1/activator complex and the dissociation of the activator from LFA-1. 2


Page 29 of 46

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



29

Associated content Supporting information Detailed information and protocols for preparing the protein structures and setting up the MD simulations, for calculating ∆GCNA, for calculating the sequence conservation and the Voronoi volumes of proteins, and for doing the statistical analysis and binary classification are provided. Additional figures are showing the conformational stability during the course of the MD simulations (S1, S5, and S10), changes upon perturbation for each individual ensemble (S2, S6, and S8), comparison of the pathway prediction in eglin c based on sequence conservation and SASA (S4), correlation between calculated Voronoi volumes and per-residue perturbation free energies (S7) as well as between the sum of per-residue ∆Gi,CNA and ∆GCNA for mutations in eglin c (S3), and the allosteric communication as derived from calculating the current flow betweenness centrality in LFA-1 (S9). Tables S1-S3 show the loss of constraints upon perturbing the systems and Table S4 the perturbation free energies for 13 single- and double mutant variants of eglin c.



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


Page 30 of 46

30

Author information Corresponding Author *E-mail: [email protected]. Phone: (+49) 211-81-13662. Fax: (+49) 211-8113847.

Author Contributions H.G. conceived and designed research; C.P. and A.M. performed research; C.P. and H.G. analyzed data; C.P. and H.G. wrote the paper; M.B., C.L.M., and R.T. contributed to the writing of the paper.

Funding Sources This study was supported by funds from Pfizer Inc and, in part, by the Deutsche Forschungsgemeinschaft (DFG) through the Research Unit FOR 2518. Financial support by DFG for funds (INST 208/704-1 FUGG) to purchase the hybrid computer cluster used in this study is gratefully acknowledged.

Notes The authors declare the following competing financial interest(s): M. Boehm, C. L. McClendon, and R. Torella are employees of Pfizer Inc..


Page 31 of 46

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



31

Acknowledgements We thank D. Mulnaes, D. Schmidt, and S.M.A. Hermans for fruitful discussions and critically reading the manuscript, D. Mulnaes for help with computing the sequence conservation, and B. Frieg for providing force field parameters for the allosteric effectors used in MD simulations. Computational support and infrastructure was provided by the “Centre for Information and Media Technology” (ZIM) at Heinrich Heine University Düsseldorf. We are grateful for computing time granted by the John von Neumann Institute for Computing (NIC) provided on the supercomputer JURECA at the Jülich Supercomputing Centre (JSC) (NIC project id HDD17).



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


Page 32 of 46

32

Abbreviations CNA, constraint network analysis; FIRST, floppy inclusions and rigid substructure topography; ICAM, intercellular adhesion molecules; LFA-1, lymphocyte function-associated antigen 1; MD, molecular dynamics; NMR, nuclear magnetic resonance; PDB, protein data bank; PTP1B, protein tyrosine phosphatase 1B; SASA, solvent accessible surface area; SEM, standard error of the mean.


Page 33 of 46

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



33

References (1) Changeux, J. P., 50 years of allosteric interactions: the twists and turns of the models. Nat. Rev. Mol. Cell Biol. 2013, 14, 819-829. (2) Jeffrey, P. D.; Ruso, A. A.; Polyak, K.; Gibbs, E.; Hurwitz, J.; Massague, J.; Pavletich, N. P., Mechanism of CDK activation revealed by the structure of a cyclina-CDK2 complex. Nature 1995, 376, 313-320. (3) Blackmore, N. J.; Reichau, S.; Jiao, W. T.; Hutton, R. D.; Baker, E. N.; Jameson, G. B.; Parker, E. J., Three sites and you are out: ternary synergistic allostery controls aromatic amino acid biosynthesis in Mycobacterium tuberculosis. J. Mol. Biol. 2013, 425, 1582-1592. (4) Srinivasan, S.; Nwachukwu, J. C.; Parent, A. A.; Cavett, V.; Nowak, J.; Hughes, T. S.; Kojetin, D. J.; Katzenellenbogen, J. A.; Nettles, K. W., Ligand-binding dynamics rewire cellular signaling via estrogen receptor-alpha. Nat. Chem. Biol. 2013, 9, 326-332. (5) Jura, N.; Zhang, X. W.; Endres, N. F.; Seeliger, M. A.; Schindler, T.; Kuriyan, J., Catalytic control in the EGF Receptor and Its connection to general kinase regulatory mechanisms. Mol. Cell 2011, 42, 9-22. (6) Freiburger, L. A.; Baettig, O. M.; Sprules, T.; Berghuis, A. M.; Auclair, K.; Mittermaier, A. K., Competing allosteric mechanisms modulate substrate binding in a dimeric enzyme. Nat. Struct. Mol. Biol. 2011, 18, 288-294. (7) Nussinov, R.; Tsai, C. J.; Ma, B., The underappreciated role of allostery in the cellular network. Annu. Rev. Biophys. 2013, 42, 169-189. (8) Wagner, J. R.; Lee, C. T.; Durrant, J. D.; Malmstrom, R. D.; Feher, V. A.; Amaro, R. E., Emerging computational methods for the rational discovery of allosteric drugs. Chem. Rev. 2016, 116, 6370-6390. (9) Monod, J.; Wyman, J.; Changeux, J. P., On the nature of allosteric transitions: a plausible model. J. Mol. Biol. 1965, 12, 88-118. (10) Koshland, D. E., Jr.; Nemethy, G.; Filmer, D., Comparison of experimental binding data and theoretical models in proteins containing subunits. Biochemistry 1966, 5, 365-385. (11) Cui, Q.; Karplus, M., Allostery and cooperativity revisited. Protein Sci. 2008, 17, 12951307. (12) Motlagh, H. N.; Wrabl, J. O.; Li, J.; Hilser, V. J., The ensemble nature of allostery. Nature 2014, 508, 331-339. (13) Perutz, M. F., Stereochemistry of cooperative effects in haemoglobin. Nature 1970, 228, 726-734. (14) Changeux, J. P.; Edelstein, S. J., Allosteric mechanisms of signal transduction. Science 2005, 308, 1424-1428. (15) Freire, E., The propagation of binding interactions to remote sites in proteins: analysis of the binding of the monoclonal antibody D1.3 to lysozyme. Proc. Natl. Acad. Sci. U. S. A. 1999, 96, 10118-10122. ACS Paragon Plus Environment


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


Page 34 of 46

34

(16) Lockless, S. W.; Ranganathan, R., Evolutionarily conserved pathways of energetic connectivity in protein families. Science 1999, 286, 295-299. (17) Tang, S.; Liao, J. C.; Dunn, A. R.; Altman, R. B.; Spudich, J. A.; Schmidt, J. P., Predicting allosteric communication in myosin via a pathway of conserved residues. J. Mol. Biol. 2007, 373, 1361-1373. (18) Schueler-Furman, O.; Wodak, S. J., Computational approaches to investigating allostery. Curr. Opin. Struct. Biol. 2016, 41, 159-171. (19) Dokholyan, N. V., Controlling allosteric networks in proteins. Chem. Rev. 2016, 116, 6463-6487. (20) Koshland, D. E., Enzyme flexibility and enzyme action. J. Cell. Comp. Physiol. 1959, 54, 245-258. (21) Cooper, A.; Dryden, D. T., Allostery without conformational change. A plausible model. Eur. Biophys. J. 1984, 11, 103-109. (22) Kern, D.; Zuiderweg, E. R., The role of dynamics in allosteric regulation. Curr. Opin. Struct. Biol. 2003, 13, 748-757. (23) Popovych, N.; Sun, S. J.; Ebright, R. H.; Kalodimos, C. G., Dynamically driven protein allostery. Nat. Struct. Mol. Biol. 2006, 13, 831-838. (24) Tzeng, S. R.; Kalodimos, C. G., Dynamic activation of an allosteric regulatory protein. Nature 2009, 462, 368-372. (25) Guo, J.; Zhou, H. X., Protein allostery and conformational dynamics. Chem. Rev. 2016, 116, 6503-6515. (26) Sekhar, A.; Kay, L. E., NMR paves the way for atomic level descriptions of sparsely populated, transiently formed biomolecular conformers. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 12867-12874. (27) Wand, A. J., The dark energy of proteins comes to light: conformational entropy and its role in protein function revealed by NMR relaxation. Curr. Opin. Struct. Biol. 2013, 23, 7581. (28) Manley, G.; Rivalta, I.; Loria, J. P., Solution NMR and computational methods for understanding protein allostery. J. Phys. Chem. B 2013, 117, 3063-3073. (29) Frauenfelder, H.; Sligar, S. G.; Wolynes, P. G., The energy landscapes and motions of proteins. Science 1991, 254, 1598-1603. (30) McMahon, B.; Frauenfelder, H.; Austin, R.; Chu, K.; Groves, J. T., The role of structure, energy landscape, dynamics, and allostery in the enzymatic function of myoglobin. Biophys. J. 2001, 80, 286a-286a. (31) Kar, G.; Keskin, O.; Gursoy, A.; Nussinov, R., Allostery and population shift in drug discovery. Curr. Opin. Pharmacol. 2010, 10, 715-722.


Page 35 of 46

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



35

(32) Cuendet, M. A.; Weinstein, H.; LeVine, M. V., The allostery landscape: Quantifying thermodynamic couplings in biomolecular systems. J. Chem. Theory Comput. 2016, 12, 57585767. (33) Hilser, V. J.; Thompson, E. B., Intrinsic disorder as a mechanism to optimize allosteric coupling in proteins. Proc. Natl. Acad. Sci. U. S. A. 2007, 104, 8311-8315. (34) Nussinov, R.; Tsai, C. J., Allostery without a conformational change? Revisiting the paradigm. Curr. Opin. Struct. Biol. 2015, 30, 17-24. (35) Hermans, S. M. A.; Pfleger, C.; Nutschel, C.; Hanke, C. A.; Gohlke, H., Rigidity theory for biomolecules: Concepts, software, and applications. WIREs Comput. Mol. Sci. 2017. (36) Thorpe, M. F., Continuous deformations in random networks. J. Non-Cryst. Solids 1983, 57, 355-370. (37) Duxbury, P. M.; Jacobs, D. J.; Thorpe, M. F.; Moukarzel, C., Floppy modes and the free energy: Rigidity and connectivity percolation on Bethe lattices. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 1999, 59, 2084-2092. (38) Gohlke, H.; Ben-Shalom, I. Y.; Kopitz, H.; Pfeiffer-Marek, S.; Baringhaus, K. H., Rigidity theory-based approximation of vibrational entropy changes upon binding to biomolecules. J. Chem. Theory Comput. 2017, 13, 1495-1502. (39) Rader, A. J.; Hespenheide, B. M.; Kuhn, L. A.; Thorpe, M. F., Protein unfolding: rigidity lost. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 3540-3545. (40) Thorpe, M. F.; Lei, M.; Rader, A. J.; Jacobs, D. J.; Kuhn, L. A., Protein flexibility and dynamics using constraint theory. J. Mol. Graph. Model. 2001, 19, 60-69. (41) Jacobs, D. J.; Rader, A. J.; Kuhn, L. A.; Thorpe, M. F., Protein flexibility predictions using graph theory. Proteins 2001, 44, 150-165. (42) Gohlke, H.; Kiel, C.; Case, D. A., Insights into protein-protein binding by binding free energy calculation and free energy decomposition using a generalized born model. J. Mol. Biol. 2003, 330, 891–913. (43) Jacobs, D. J.; Thorpe, M. F., Generic rigidity percolation: The pebble game. Phys. Rev. Lett. 1995, 75, 4051-4054. (44) Katoh, N.; Tanigawa, S., A proof of the molecular conjecture. Discrete Comput. Geom. 2011, 45, 647-700. (45) Hespenheide, B. M.; Rader, A. J.; Thorpe, M. F.; Kuhn, L. A., Identifying protein folding cores from the evolution of flexible regions during unfolding. J. Mol. Graph. Model. 2002, 21, 195-207. (46) Gohlke, H.; Kuhn, L. A.; Case, D. A., Change in protein flexibility upon complex formation: analysis of Ras-Raf using molecular dynamics and a molecular framework approach. Proteins 2004, 56, 322-337.



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


Page 36 of 46

36

(47) Rader, A. J.; Anderson, G.; Isin, B.; Khorana, H. G.; Bahar, I.; Klein-Seetharaman, J., Identification of core amino acids stabilizing rhodopsin. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 7246-7251. (48) Livesay, D. R.; Jacobs, D. J., Conserved quantitative stability/flexibility relationships (QSFR) in an orthologous RNase H pair. Proteins 2006, 62, 130-143. (49) Radestock, S.; Gohlke, H., Exploiting the link between protein rigidity and thermostability for data-driven protein engineering. Eng. Life Sci. 2008, 8, 507-522. (50) Fulle, S.; Gohlke, H., Statics of the ribosomal exit tunnel: implications for cotranslational peptide folding, elongation regulation, and antibiotics binding. J. Mol. Biol. 2009, 387, 50217. (51) Rathi, P. C.; Radestock, S.; Gohlke, H., Thermostabilizing mutations preferentially occur at structural weak spots with a high mutation ratio. J. Biotechnol. 2012, 159, 135-144. (52) Rader, A. J., Thermostability in rubredoxin and its relationship to mechanical rigidity. Phys. Biol. 2009, 7, 16002. (53) Rathi, P. C.; Jaeger, K. E.; Gohlke, H., Structural rigidity and protein thermostability in variants of lipase A from Bacillus subtilis. PLoS One 2015, 10, e0130289. (54) Radestock, S.; Gohlke, H., Protein rigidity and thermophilic adaptation. Proteins 2011, 79, 1089-1108. (55) Dahiyat, B. I.; Gordon, D. B.; Mayo, S. L., Automated design of the surface positions of protein helices. Protein Sci. 1997, 6, 1333-1337. (56) Pfleger, C.; Gohlke, H., Efficient and robust analysis of biomacromolecular flexibility using ensembles of network topologies based on fuzzy noncovalent constraints. Structure 2013, 21, 1725-1734. (57) Pfleger, C.; Radestock, S.; Schmidt, E.; Gohlke, H., Global and local indices for characterizing biomolecular flexibility and rigidity. J. Comput. Chem. 2013, 34, 220-233. (58) Skolnick, J.; Jaroszewski, L.; Kolinski, A.; Godzik, A., Derivation and testing of pair potentials for protein folding. When is the quasichemical approximation correct? Protein Sci. 1997, 6, 676-688. (59) Robertson, A. D.; Murphy, K. P., Protein structure and the energetics of protein stability. Chem. Rev. 1997, 97, 1251-1267. (60) Kollman, P., Free-energy calculations - Applications to chemical and biochemical phenomena. Chem. Rev. 1993, 93, 2395-2417. (61) Zwanzig, R. W., High‐temperature equation of state by a perturbation method. I. nonpolar gases. J. Chem. Phys. 1954, 22, 1420-1426. (62) Jorgensen, W. L.; Ravimohan, C., Monte-carlo simulation of differences in free-energies of hydration. J. Chem. Phys. 1985, 83, 3050-3054.


Page 37 of 46

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



37

(63) Rod, T. H.; Brooks, C. L., How dihydrofolate reductase facilitates protonation of dihydrofolate. J. Am. Chem. Soc. 2003, 125, 8718-8719. (64) Acevedo, O.; Ambrose, Z.; Flaherty, P. T.; Aamer, H.; Jain, P.; Sambasivarao, S. V., Identification of HIV inhibitors guided by free energy perturbation calculations. Curr. Pharm. Des. 2012, 18, 1199-1216. (65) Oostenbrink, C., Free energy calculations from one-step perturbations. Methods Mol. Biol. 2012, 819, 487-499. (66) Jorgensen, W. L., Free-energy calculations - A breakthrough for modeling organicchemistry in solution. Acc. Chem. Res. 1989, 22, 184-189. (67) Jorgensen, W. L.; Buckner, J. K.; Boudon, S.; Tiradorives, J., Efficient computation of absolute free-energies of binding by computer-simulations - Application to the methane dimer in water. J. Chem. Phys. 1988, 89, 3742-3746. (68) Hermans, J.; Wang, L., Inclusion of loss of translational and rotational freedom in theoretical estimates of free energies of binding. Application to a complex of benzene and mutant T4 lysozyme. J. Am. Chem. Soc. 1997, 119, 2707-2714. (69) Roux, B.; Nina, M.; Pomes, R.; Smith, J. C., Thermodynamic stability of water molecules in the bacteriorhodopsin proton channel: A molecular dynamics free energy perturbation study. Biophys. J. 1996, 71, 670-681. (70) Lee, F. S.; Chu, Z. T.; Bolger, M. B.; Warshel, A., Calculations of antibody antigen interactions - Microscopic and semimicroscopic evaluation of the free-energies of binding of phosphorylcholine analogs to Mcpc603. Protein Eng. 1992, 5, 215-228. (71) Aqvist, J.; Medina, C.; Samuelsson, J. E., A new method for predicting binding affinity in computer-aided drug design. Protein Eng. 1994, 7, 385-391. (72) Hansson, T.; Marelius, J.; Aqvist, J., Ligand binding affinity prediction by linear interaction energy methods. J. Comput.-Aided Mol. Des. 1998, 12, 27-35. (73) Lamb, M. L.; Tirado-Rives, J.; Jorgensen, W. L., Estimation of the binding affinities of FKBP12 inhibitors using a linear response method. Bioorg. Med. Chem. 1999, 7, 851-860. (74) Hulten, J.; Bonham, N. M.; Nillroth, U.; Hansson, T.; Zuccarello, G.; Bouzide, A.; Aqvist, J.; Classon, B.; Danielson, U. H.; Karlen, A.; Kvarnstrom, I.; Samuelsson, B.; Hallberg, A., Cyclic HIV-1 protease inhibitors derived from mannitol: Synthesis, inhibitory potencies, and computational predictions of binding affinities. J. Med. Chem. 1997, 40, 885897. (75) Aqvist, J., Calculation of absolute binding free energies for charged ligands and effects of long-range electrostatic interactions. J. Comput. Chem. 1996, 17, 1587-1597. (76) Aqvist, J.; Mowbray, S. L., Sugar recognition by a glucose/galactose receptor Evaluation of binding energetics from molecular-dynamics simulations. J. Biol. Chem. 1995, 270, 9978-9981.



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


Page 38 of 46

38

(77) Jones-Hertzog, D. K.; Jorgensen, W. L., Binding affinities for sulfonamide inhibitors with human thrombin using Monte Carlo simulations with a linear response method. J. Med. Chem. 1997, 40, 1539-1549. (78) Aqvist, J.; Marelius, J., The linear interaction energy method for predicting ligand binding free energies. Comb. Chem. High Throughput Screening 2001, 4, 613-626. (79) Lamb, M. L.; Jorgensen, W. L., Computational approaches to molecular recognition. Curr. Opin. Chem. Biol. 1997, 1, 449-457. (80) Warshel, A.; Russell, S. T., Calculations of Electrostatic Interactions in BiologicalSystems and in Solutions. Q. Rev. Biophys. 1984, 17, 283-422. (81) Roux, B.; Yu, H. A.; Karplus, M., Molecular-Basis for the Born Model of Ion Solvation. J. Phys. Chem. 1990, 94, 4683-4688. (82) Marcus, R. A., Chemical + Electrochemical Electron-Transfer Theory. Annu. Rev. Phys. Chem. 1964, 15, 155-&. (83) Pfleger, C.; Rathi, P. C.; Klein, D. L.; Radestock, S.; Gohlke, H., Constraint Network Analysis (CNA): a Python software package for efficiently linking biomacromolecular structure, flexibility, (thermo-)stability, and function. J. Chem. Inf. Model. 2013, 53, 1007-15. (84) Guckian, K. M.; Lin, E. Y. S.; Silvian, L.; Friedman, J. E.; Chin, D.; Scott, D. M., Design and synthesis of a series of meta aniline-based LFA-1 ICAM inhibitors. Bioorg. Med. Chem. Lett. 2008, 18, 5249-5251. (85) Bégué, J.-P.; Bonnet-Delpon, D., Bioorganic and medicinal chemistry of fluorine. John Wiley & Sons: Hoboken, N.J, 2008; p XVII, 365 S. (86) Anthonisse, J. M., The rush in a directed graph. Stichting Mathematisch Centrum. Mathematische Besliskunde: Amsterdam, 1971. (87) Freeman, L. C., A set of measures of centrality based on betweenness. Sociometry 1977, 35-41. (88) Brandes, U.; Fleischer, D., Centrality measures based on current flow. Stacs 2005, Proceedings 2005, 3404, 533-544. (89) Newman, M. E. J., A measure of betweenness centrality based on random walks. Soc. Networks 2005, 27, 39-54. (90) Hagberg, A. A.; Schult, D. A. S., P. J. In Exploring network structure, dynamics, and function using NetworkX, Proceedings of the 7th Python in Science conference, Pasadena, CA USA, Varoquaux, G.; Vaught, T.; Millman, J., Eds. Pasadena, CA USA, 2008; pp 11-15. (91) Brandes, U., On variants of shortest-path betweenness centrality and their generic computation. Soc. Networks 2008, 30, 136-145. (92) Brandes, U., A faster algorithm for betweenness centrality. J. Math. Sociol. 2001, 25, 163-177.


Page 39 of 46

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



39

(93) Brandes, U.; Wagner, D., Analysis and visualization of social networks. In Graph Drawing Software, Jünger, M.; Mutzel, P., Eds. Springer Berlin Heidelberg: Berlin ;Heidelberg, 2004; pp 321-340. (94) Clarkson, M. W.; Lee, A. L., Long-range dynamic effects of point mutations propagate through side chains in the serine protease inhibitor eglin c. Biochemistry 2004, 43, 1244812458. (95) Clarkson, M. W.; Gilmore, S. A.; Edgell, M. H.; Lee, A. L., Dynamic coupling and allosteric behavior in a nonallosteric protein. Biochemistry 2006, 45, 7693-7699. (96) Wiesmann, C.; Barr, K. J.; Kung, J.; Zhu, J.; Erlanson, D. A.; Shen, W.; Fahr, B. J.; Zhong, M.; Taylor, L.; Randal, M.; McDowell, R. S.; Hansen, S. K., Allosteric inhibition of protein tyrosine phosphatase 1B. Nat. Struct. Mol. Biol. 2004, 11, 730-737. (97) Puius, Y. A.; Zhao, Y.; Sullivan, M.; Lawrence, D. S.; Almo, S. C.; Zhang, Z. Y., Identification of a second aryl phosphate-binding site in protein-tyrosine phosphatase 1B: A paradigm for inhibitor design. Proc. Natl. Acad. Sci. U. S. A. 1997, 94, 13420-13425. (98) Jia, Z.; Barford, D.; Flint, A. J.; Tonks, N. K., Structural basis for phosphotyrosine peptide recognition by protein tyrosine phosphatase 1B. Science 1995, 268, 1754-1758. (99) Salmeen, A.; Andersen, J. N.; Myers, M. P.; Tonks, N. K.; Barford, D., Molecular basis for the dephosphorylation of the activation segment of the insulin receptor by protein tyrosine phosphatase 1B. Mol. Cell 2000, 6, 1401-1412. (100) Thorpe, M.; Jacobs, D.; Djordjevic, B., Generic rigidity percolation. 1996; Vol. 11, p 407-424. (101) Nussinov, R.; Tsai, C. J.; Csermely, P., Allo-network drugs: Harnessing allostery in cellular networks. Trends Pharmacol. Sci. 2011, 32, 686-693. (102) Qu, A.; Leahy, D. J., Crystal structure of the I-domain from the CD11a/CD18 (LFA-1, alpha L beta 2) integrin. Proc. Natl. Acad. Sci. U. S. A. 1995, 92, 10277-10281. (103) Huth, J. R.; Olejniczak, E. T.; Mendoza, R.; Liang, H.; Harris, E. A.; Lupher, M. L., Jr.; Wilson, A. E.; Fesik, S. W.; Staunton, D. E., NMR and mutagenesis evidence for an I domain allosteric site that regulates lymphocyte function-associated antigen 1 ligand binding. Proc. Natl. Acad. Sci. U. S. A. 2000, 97, 5231-5236. (104) Hintersteiner, M.; Kallen, J.; Schmied, M.; Graf, C.; Jung, T.; Mudd, G.; Shave, S.; Gstach, H.; Auer, M., Identification and X-ray co-crystal structure of a small-molecule activator of LFA-1-ICAM-1 binding. Angew. Chem., Int. Ed. 2014, 53, 4322-4326. (105) Shimaoka, M.; Xiao, T.; Liu, J. H.; Yang, Y. T.; Dong, Y. C.; Jun, C. D.; McCormack, A.; Zhang, R. G.; Joachimiak, A.; Takagi, J.; Wang, J. H.; Springer, T. A., Structures of the alpha L I domain and its complex with ICAM-1 reveal a shape-shifting pathway for integrin regulation. Cell 2003, 112, 99-111. (106) Edwards, C. P.; Fisher, K. L.; Presta, L. G.; Bodary, S. C., Mapping the intercellular adhesion molecule-1 and -2 binding site on the inserted domain of leukocyte functionassociated antigen-1. J. Biol. Chem. 1998, 273, 28937-28944.



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


Page 40 of 46

40

(107) Welzenbach, K.; Hommel, U.; Weitz-Schmidt, G., Small molecule inhibitors induce conformational changes in the I domain and the I-like domain of lymphocyte functionassociated antigen-1 - Molecular insights into integrin inhibition. J. Biol. Chem. 2002, 277, 10590-10598. (108) Weitz-Schmidt, G.; Welzenbach, K.; Dawson, J.; Kallen, J., Improved lymphocyte function-associated antigen-1 (LFA-1) inhibition by statin derivatives - Molecular basis determined by X-ray analysis and monitoring of LFA-1 conformational changes in vitro and ex vivo. J. Biol. Chem. 2004, 279, 46764-46771. (109) Tsai, C. J.; del Sol, A.; Nussinov, R., Allostery: absence of a change in shape does not imply that allostery is not at play. J. Mol. Biol. 2008, 378, 1-11. (110) McClendon, C. L.; Friedland, G.; Mobley, D. L.; Amirkhani, H.; Jacobson, M. P., Quantifying correlations between allosteric sites in thermodynamic ensembles. J. Chem. Theory Comput. 2009, 5, 2486-2502. (111) Wang, J. M.; Morin, P.; Wang, W.; Kollman, P. A., Use of MM-PBSA in reproducing the binding free energies to HIV-1 RT of TIBO derivatives and predicting the binding mode to HIV-1 RT of efavirenz by docking and MM-PBSA. J. Am. Chem. Soc. 2001, 123, 52215230. (112) Mottonen, J. M.; Jacobs, D. J.; Livesay, D. R., Allosteric response is both conserved and variable across three CheY orthologs. Biophys. J. 2010, 99, 2245-2254. (113) Pannifer, A. D.; Flint, A. J.; Tonks, N. K.; Barford, D., Visualization of the cysteinylphosphate intermediate of a protein-tyrosine phosphatase by x-ray crystallography. J. Biol. Chem. 1998, 273, 10454-10462. (114) Weitz-Schmidt, G.; Welzenbach, K.; Brinkmann, V.; Kamata, T.; Kallen, J.; Bruns, C.; Cottens, S.; Takada, Y.; Hommel, U., Statins selectively inhibit leukocyte function antigen-1 by binding to a novel regulatory integrin site. Nat Med 2001, 7, 687-692. (115) Kallen, J.; Welzenbach, K.; Ramage, P.; Geyl, D.; Kriwacki, R.; Legge, G.; Cottens, S.; Weitz-Schmidt, G.; Hommel, U., Structural basis for LFA-1 inhibition upon lovastatin binding to the CD11a I-domain. J. Mol. Biol. 1999, 292, 1-9. (116) Liu, J.; Nussinov, R., Allostery: An overview of its history, concepts, methods, and applications. PLoS Comput. Biol. 2016, 12, e1004966. (117) Tsai, C. J.; Nussinov, R., A unified view of "How allostery works". PLoS Comput. Biol. 2014, 10, e100339. (118) Stetz, G.; Verkhivker, G. M., Computational Analysis of Residue Interaction Networks and Coevolutionary Relationships in the Hsp70 Chaperones: A Community-Hopping Model of Allosteric Regulation and Communication. PLoS Comput Biol 2017, 13, e1005299. (119) Whiteley, W., Counting out to the flexibility of molecules. Phys. Biol. 2005, 2, S11626.


Page 41 of 46

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



TOC Figure


41


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

Figure 1 232x257mm (300 x 300 DPI)


Page 42 of 46

Page 43 of 46

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


Figure 2 171x140mm (300 x 300 DPI)



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

Figure 3 221x232mm (300 x 300 DPI)


Page 44 of 46

Page 45 of 46

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


Figure 4 167x133mm (300 x 300 DPI)



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

TOC graphic 88x35mm (300 x 300 DPI)


Page 46 of 46

Ensemble- and Rigidity Theory-Based Perturbation Approach To

Recommend Documents