ProteinâLigand Complexes as Constrained Dynamical Systems

Subscriber access provided by UNIV OF NEW ENGLAND ARMIDALE

Computational Biochemistry

Protein-ligand complexes as constrained dynamical systems Burak T. Kaynak, and Pemra Doruker J. Chem. Inf. Model., Just Accepted Manuscript • Publication Date (Web): 26 Mar 2019 Downloaded from http://pubs.acs.org on March 27, 2019

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 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 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.

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 23 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 Information and Modeling

1

Protein-ligand complexes as constrained dynamical systems Burak T. Kaynak1† and Pemra Doruker1,2‡ 1

Department of Computational and Systems Biology, School of Medicine, University of

Pittsburgh, Pittsburgh, PA 15261, USA 2 Department

of Chemical Engineering and Polymer Research Center, Bogazici University,

34342 Bebek, Istanbul, Turkey * [email protected]

Abstract This study focuses on how the low-frequency end of the vibrational spectrum related to the functional motions changes as a protein binds to a relatively small ligand(s). Our recently proposed residue specific elastic network model (RESPEC) provides a natural laboratory for this aim due to its systematic mixed coarse-graining approach and parametrization. Current analysis on a large data set of protein-ligand complexes reveals a universal curve enclosing the frequency distributions, which bears the features of previous computational and experimental studies. We mostly observe positive frequency shifts in the collective modes of the protein upon ligand binding. This observation conforming to the Rayleigh-Courant-Fisher theorem points to a constraining effect imposed by ligands on protein dynamics, which may be accompanied by a negative vibrational entropy difference. Positive frequency shifts in the global modes can thus be linked to the harmonic well getting steeper due to interactions with ligand(s).

ACS Paragon Plus Environment

Journal of Chemical Information and Modeling 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 2 of 23

2

Introduction Binding of ligands modifies the conformational dynamics and flexibility of proteins to various extents. The preexisting equilibrium conformations accessible to the apo protein may be modified1–3 during its complexation with other molecules, hence causing possible shifts between the conformational states on the dynamical energy landscape4. Complexation may also alter the protein dynamics without a conformational shift/transition by modifying the steepness of the minimum energy well. These effects are closely linked to allosteric drug design5,6 and enzyme function7–9. Universal features of the vibrational density of states have been established for proteins using computational means of normal mode analyses (NMA)10,11, which exhibit a major peak below 100 cm-1 (3 THz). The lowest end of the spectrum is populated by the collective (delocalized) modes, which are known to guide transitions between apo and bound states, especially for hinge-bending type proteins using elastic network models (ENM)12–14. The distribution related to more localized motions also bears common features/peaks corresponding to vibrations of secondary structural elements, side chains, atoms etc11. Experiments have discussed ligand binding effects on protein vibrations as either stiffening or softening. For example, inelastic neutron scattering (INS) has shown that the cancer drug methotrexate binding softens DHFR low-frequency modes15, which was later explained by atomistic NMA performed with molecular dynamics (MD) simulation conformers16. In contrast, femtosecond optical Kerr-effect spectroscopy (OKE) has identified a positive frequency shift in a terahertz underdamped delocalized vibrational mode upon lysozyme-inhibitor binding in solution17. NMR experiments on protein-ligand binding have been reviewed18 in terms of decreasing or increasing effects on backbone and side-chain flexibility that are determinants of the protein conformational entropy. Underdamped modes could be detected through the application of time series analysis on the essential modes obtained from MD simulations19, which rather display diffusive character. For two proteins, the ligand binding indicated a shift to higher frequencies of the detected underdamped


Page 3 of 23 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 modes. Recently, such underdamped modes have been shown to exist in the terahertz range by optical microscopy20 and their relevance for protein-ligand binding by OKE17. Thus, changes in these delocalized, collective motions, which are known to promote allosteric control of protein function, are of interest especially regarding computational current drug design. In the present work, we will analyze how the vibrational spectrum of a protein changes during formation of a protein-ligand complex, based on a diverse data set of protein complexes. Recently, we have introduced a new framework RESPEC21 to perform vibrational analysis of proteins both in apo and in ligand-bound forms. This model is a novel extension of the classical ENM13,22, and is primarily based on its anisotropic version (ANM)23. RESPEC is developed so as to render protein-ligand complexes by introducing residue specificity into the classical ENM. Unlike the classical ENM, the interaction strengths between the nodes of unique pairs in RESPEC are nonuniform across the network, and their values are determined on the basis of whether they belong to atoms of a residue or a ligand. RESPEC allows for considerable improvement in the level of correlations between the experimental Debye-Waller factors (B-factors) and the theoretical mean square fluctuations (MSFs) for proteins and even higher for protein-ligand complexes.

Materials and Method RESPEC and its parameters The implementation of RESPEC is accomplished in the following stages once the network is constructed by an upper cutoff 𝛿𝑈: (i) local coarse-graining of protein residues by treating them as packs of their heavy atoms as virtual nodes localized around their 𝐶𝛼-atoms, (ii) scaling the nodes' coordinates either by residue masses or ligand atoms' masses, (iii) determining the coupling constants of unique pairs, either via a single or through a sum of a power law function of atomatom distance, (iv) introducing a secondary cutoff modifying the range of the interaction between a protein residue and a ligand atom, and (v) coarse-graining of B-factors so as to assign an effective B-factor for each protein residue by using its heavy atoms' B-factors. The second and the fifth steps also allow us not only to preserve the dynamical hierarchy between protein residues and



Page 4 of 23

4 ligand atoms consistently, but also to apply coarse-graining technique unambiguously within a mixed-resolution scheme24. As in the classical ENM, RESPEC models a protein molecule, which is a highly complicated many-body system, as a system of nodes of an elastic network, where unique pairs of nodes are interacting via simple harmonic potentials. Let 𝑁𝑅, 𝑁𝐴 and 𝑁 be the number of protein residues, the number of ligand atoms, and the total number of nodes, respectively. Then the potential energy of (𝑖,𝑗)-pair for a protein-ligand complex, where 1 ≤ 𝑖,𝑗 ≤ 𝑁, is given by

𝑉𝑖𝑗 =

𝛾𝑖𝑗

(|𝑞𝑖 ― 𝑞𝑗| ― |𝑞0𝑖 ― 𝑞0𝑗 |)2 ,

2

(1)

where 𝑞𝑖 and 𝑞0𝑖 stand for the generalized coordinates of the 𝑖𝑡ℎ node and the equilibrium position thereof, respectively. Here each 𝑞𝑖 has spatial indices (x, y, z), albeit being suppressed, so | ⋅ | stands for the Euclidean norm. The total potential energy is given by a summation over all unique pairs' potential energy terms. In Eq. (1) 𝛾𝑖𝑗 is the coupling constant between the nodes of the (𝑖,𝑗)-pair. Once the unique pairs are formed under the condition|𝑞0𝑖 ― 𝑞0𝑗 | < 𝛿𝑈, the coupling constant for each (𝑖,𝑗)-pair is defined by

𝛾𝑖𝑗 =

{

𝑝

( ) ∑ ( ) ( ) ∑

∀(𝑘,𝑙) ∈ (𝑖,𝑗) 𝑑0𝑘𝑙 < 𝛿𝑃𝐿

𝑑0

𝑑0𝑘𝑙 𝑑0

𝑑0𝑘𝑙

if (𝑖,𝑗) ∈ {(𝑅,𝑅)} ,

𝑝

if (𝑖,𝑗) ∈ {(𝑅,𝐴)} ,

∀(𝑘,𝑙) ∈ (𝑖,𝑗) 0 𝑝

𝑑

𝑑0𝑘𝑙

(2)

if (𝑖,𝑗) ∈ {(𝐴,𝐴)} ,

where 𝑑0𝑖𝑗 is the Euclidean distance within the (𝑖,𝑗)-pair, 𝑝 ∈ ℝ ≥ 0, and 𝛿𝑃𝐿 is the protein-ligand cutoff. 𝑑0 is set to be 1 Å for future convenience. {(𝑅,𝐴)} in Eq. (2) denotes the set of all pairs formed by a residue and a ligand atom, and the other cases are self-explanatory. ∀(𝑘,𝑙) ∈ (𝑖,𝑗) also indicates that the sums should be performed over all possible connections that can be built up between the nodes 𝑖 and 𝑗. By “all possible pairs” we mean that if at least one of the nodes in the


Page 5 of 23 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 (𝑖,𝑗)-pair is a residue, then the contribution of all individual heavy atoms in that residue should be taken into account so as to determine the effective coupling constant of that pair as if those atoms were incorporated into the network. Therefore, we call these atoms “virtual nodes”. The proteinligand cutoff by itself allows us to control how a ligand atom interacts with surrounding protein residues via restricting the extent of the interaction. In essence, the summation over all possible interactions occurring through the virtual nodes mimics integrating out the effect of fine details (within a cutoff distance in the case of a residue-atom pair) in the sense of renormalization group theory25,26. As explained in Ref.21, once we calculate the MSFs via the eigenvalues and eigenvectors of the Hessian of the system, we assess how well the model renders the protein-ligand complex via comparing them with the B-factors. This comparison yields a system wide uniform constant, which is determined in terms of kB T. This constant is called experimental spring constant in the classical ENMs, whereas we call it the overall normalization of coupling constants in RESPEC and denote it by γ0. This constant leads to the normalization of the eigenvalues λν as γ0λν. Such a normalization is important since it fixes/normalizes not only the free parameters of the model but also the quantities with physical meaning as regards to the scale at which the X-ray crystallography experiment is conducted. We give a rough estimate of coarse-grained B-factors21 by the following equation 𝐵𝑅𝑖 ≈

1

𝑁𝑅𝑖

∑𝐵 ,

(𝑁𝑅𝑖)2𝑗 = 1

𝑗

(3)

where the sum is taken over the individual B-factors, Bj, of heavy atoms in the ith residue. 𝑁𝑅𝑖 stands for the number of heavy atoms in the ith residue. In the case of ligand binding, a generalized Bfactor vector is needed, and it is defined as 𝐵𝐺 = (𝐵𝑅1 ,…,𝐵𝑅𝑁𝑅,𝐵𝐴1 ,…,𝐵𝐴𝑁𝐴). Comparison with B-factors indicates highest correlation for the parameter set of (𝛿𝑈,𝛿𝑃𝐿,𝑝) = (22 Å,12 Å,1), which will be used in this study. It is always possible to determine the best parameter set for a specific protein.



Page 6 of 23

6

Dataset We use the same dataset as in our previous work21, which consists of 315 protein-ligand complexes including 119 monomers, 152 dimers, and 44 tetramers. Each protein chain has at most 30% identity with others in the set and contains a maximum of 5% missing residues. Protein complexes are bound to one or more ligands, which have at least 10 heavy atoms and exclude RNA or DNA.

Results and Discussion The prediction of the functional motions of proteins via vibrational spectrum analysis with a low computational cost is widely acclaimed as the most powerful feature of coarse-grained ENMs. Therefore, we would like to understand what type of changes may occur in the vibrational spectrum when proteins and ligands form complexes. Henceforth, we follow the conventions to obtain frequencies in cm ―1 as explained in Ref.27.

Universal frequency spectrum of protein-ligand complexes The normalized frequency distribution for all complexes in our data set is plotted in Figure 1 based on the parameter set (𝛿𝑈,𝛿𝑃𝐿,𝑝) = (22 Å,12 Å,1). This parameter set is shown to have the highest mean correlation coefficient (0.92) with very small standard error of mean over the whole data set21. We would like to emphasize that the characteristic frequencies are arranged in increasing order. The full frequency distribution has its mean at 72.6 cm-1. A kernel-density estimate using the Gaussian kernel allows us to obtain a universal curve enclosing this distribution with its peak at 64 cm-1. The universal feature of the enclosing curve has been discussed including atomistic details with NMA11. Interestingly, our coarse-grained approach produces a quite consistent picture with the atomistic distribution. We should note that the extra nodes located at ligand atoms in RESPEC lead to additional high frequency modes that are present in Figure 1. Another issue of interest is how the vibrational spectrum belonging to pseudo-apo states are affected due to ligand binding. Here, the “pseudo-apo” state is used for a protein entity extracted


Page 7 of 23 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 from a complex, which disregards any conformational change that the protein would undergo during complexation. If we plot a frequency distribution for pseudo-apo states as in Figure 1, then the mean of that distribution is found to be at 65.8 cm-1 without a change in the form of the enclosing curve with a peak located at 60 cm-1. Comparison of the two distributions indicates that a significant shift occurs in the spectrum due to binding.

Figure 1. Normalized histogram of frequencies including all proteins in the data set for the parameter set (𝛿𝑈,𝛿𝑃𝐿,𝑝) = (22 Å,12 Å,1). A universal curve emerges for proteins of diverse folds, sizes and oligomerization states. The binding of ligands modifies the universal profile by shifting it toward higher frequencies, as indicated by the mean of the distribution in parentheses.

In comparison with the frequency distributions obtained from atomistic-detail NMA10,11, the left tail of the distribution in Figure 1 has fewer modes. Such a feature is also found in the frequency distributions from inelastic neutron scattering experiments on dihydrofolate reductase (DHFR)15. As such, RESPEC allows us to determine the collective motions of a protein/protein-ligand complex with fewer data since the functional, collective modes are found in this region. This is highly beneficial for conformational sampling when quite large globular proteins are of interest.



Page 8 of 23

8

Frequency shifts in collective modes due to ligand binding The extent of frequency shift observed in the presence of the ligand(s) highlights the dynamics of the binding process with regard to the collective modes. We should point out at this stage that we are not considering the experimental apo structures, which may reflect conformational changes with respect to the complex. Here, we will observe how much the steepness of the harmonic well around the complex structure is affected due to interactions with the ligand. In other words, we are interested in the changes occurring in the spectrum during the coupling-decoupling phase of a protein and its ligand(s). In order to measure such a change qualitatively, we can define the frequency shift for the 𝑖𝑡ℎ mode on the basis of the frequencies both of the protein complex 𝜔𝑐𝑖 and of its pseudo-apo 𝜔𝑎𝑖, which is as follows %𝑠ℎ𝑖𝑓𝑡(𝑖) =

(𝜔𝑐𝑖 ― 𝜔𝑎𝑖) 𝜔𝑎𝑖

× 100 .

(4)

Here, we would like to determine how the collective modes are affected by binding. Applying Eq. (4) to the slowest 100 modes allows us to plot their shifts as in Figure 2. In this work, mode index zero corresponds to the first mode with nonzero frequency. We provide both mean and median of the percent shifts over the whole data set to show the possible skew-symmetry in the data due to varying sizes and oligomerizations of the complexes. It is apparent that the positive shifts are highest for the most collective modes, having an average value > 8% for the first few modes. The shifts considerably level off after the first 20 modes, which effectually span the functional motions of proteins. Therefore, we conclude that ligand binding affects the very nature of the collective modes.


Page 9 of 23 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

Figure 2. Mean and median of the percentage shift in frequencies over the whole data set for the slowest nonzero 100 modes at (𝛿𝑈,𝛿𝑃𝐿,𝑝) = (22 Å,12 Å,1).

We would like to give a couple of specific examples to demonstrate specific trends before elaborating the shifts in the collective modes observed over the whole data set. In fact, 87% of the proteins in our data set exhibit positive shifts, i.e. increase in low frequencies, for the specific parameter set of (𝛿𝑈,𝛿𝑃𝐿,𝑝) = (22 Å,12 Å,1). In Figure 3, we display two examples, namely adenylate kinase (AK) in complex with the inhibitor AP5 (1ake28) and tyrosine kinase c-Src in complex with an allosteric inhibitor RL45 (3f3v29). While the frequency curves of the pseudo-apo structures are plotted in red, their counterparts for complexes at different 𝛿𝑃𝐿 values are plotted in black (4 Å) and blue (12 Å) for AK, and in cyan (10 Å) and blue (12 Å) for c-Src. For both proteins, the frequencies of the protein-ligand complex (blue or cyan curves) are higher than those of the corresponding pseudo-apo form, which is obtained by removing the ligand. In the next section we will further interpret these positive shifts as a constraint effect imposed by the ligands. For c-Src, a lower 𝛿𝑃𝐿 = 10 Å that corresponds to more localized protein-ligands interactions is shown to cause higher frequency shifts compared to 𝛿𝑃𝐿 = 12 Å. Decreasing the protein-ligand cutoff further to 4 Å indicates an unusual inflation in the frequencies for the AK complex. Such extremely localized interactions do not seem suitable for our model. Moreover, we emphasize that



Page 10 of 23

10 highest correlations with experiments have been observed for the zone with 𝑝 ≤ 1, and 10 Å ≤ 𝛿𝑃𝐿 ≤ 16 Å in our previous study on RESPEC21. Thus, the dynamics of protein complexes in the crystalline state is represented more realistically within this parameter zone.

Figure 3. Positive frequency shifts observed for AK and c-Src. (a) Frequencies of AK in its pseudoapo form (red) and in its complex with AP5, modeled with 𝛿𝑃𝐿 = 4, 12 Å (dashed-black, blue), (b) the corresponding frequency shifts for AK, (c) frequencies of c-Src in its pseudo-apo form (red) and in its complex with RL45, modeled with 𝛿𝑃𝐿 = 10, 12 Å (cyan, blue), and (d) the corresponding frequency shifts for c-Src. All curves are at p=1.


Page 11 of 23 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

On the other hand, 13% of the complexes in our data set display negative shifts in the collective modes, using the same parameter set. To introduce such a case, we apply our method to the tertiary complex of dihydrofolate reductase (DHFR) with NADPH and MTX, which is a folate antagonist used for cancer treatment. Neutron scattering experiments15 have shown that binding of MTX to the holo form (DHFR+NAPH) causes a shift towards lower frequencies, which has been reproduced by normal mode analysis at atomistic detail16. To study this so-called softening of the low-frequency vibrations by our mixed-resolution RESPEC analysis, we used the tertiary complex (1rx330) and the holo form (DHFR + NADPH) based on the same crystal structure. Figure 4(a) exhibits that the frequencies of the tertiary complex are lower than those of the holo form, indicating that the negative shifts (panel b) and the softening effect are in line with experiments. In contrast, panels (c) and (d) demonstrate the effect of sole NADPH binding to the pseudo-apo DHFR, which results in positive shifts in low-frequency modes.



Page 12 of 23

12

Figure 4. Vibrational softening of DHFR upon MTX and NADPH binding. (a) Frequencies of the first 100 non-zero modes for DHFR in its holo form (DHFR + NADPH) and for its complex (DHFR + NADPH + MTX), (b) the negative frequency shifts between the complex and holo forms of DHFR, implying a softening effect, (c) frequencies of the holo form in comparison with the pseudo-apo DHFR, (d) the positive shifts due to sole binding of NADPH.


Page 13 of 23 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

Positive frequency shifts as a constraint effect The important question here is how to interpret the effect of ligand binding on the vibrational spectrum. The case of positive shifts that are predominant in our data set is actually in accord with the behavior of characteristic frequencies of a dynamical system upon which a set of constraints is imposed. When we think of our system as a constrained dynamical system, the addition of ligand nodes into the network of the pseudo-apo structure can be viewed as a genuine effect from a set of constraints. In this picture, we can treat the original spectrum, i.e. the one before the inclusion of the ligand nodes, as the spectrum of the unconstrained system, whereas the spectrum of the proteinligand complex can be treated as that of the constrained system. This approach allows us to partially use the Rayleigh-Courant-Fisher Theorem31, one of the prominent theorems in classical mechanics, to compare the two spectrums and explain the difference between them. This theorem states that the characteristic frequencies of a system subjected to a set of constraints separate the frequencies of the original system so that the frequencies of the constrained system lie between the frequencies of the unconstrained system. In other words, the new and old characteristic frequencies interlace from the point of view of variational characterization of eigenvalues of Hermitian matrices32. The reason why these theorems can partially be applied in the case of protein-ligand complexes is that the dimensions of the vector spaces for the unconstrained and constrained systems do not match. However, interpreting these effects from this perspective still provides a dynamical insight to the collective motions. In this respect, we will analyze these shifts occurring in the spectrum over a set of (𝛿𝑃𝐿,𝑝)-pairs as well. Figure 5 shows the resolution of such shifts as a function of 𝛿𝑃𝐿 and 𝑝 for the first and fifth modes, where negative shifts are masked. The positive shifts become more significant for both lower 𝛿𝑃𝐿 and lower 𝑝 values, effectively forming a triangular region. In fact, this region has been shown to have both high correlation and low standard error of mean in our previous work21. Furthermore, the positive shifts in higher modes (not shown here) also populate the same region and fade off, consistent with the trend in Figures 2 and 3.



Page 14 of 23

14

Figure 5. Heat maps for different (𝛿𝑃𝐿,𝑝)-pairs showing the median of the percentage shift over the whole data set in (a) the first and (b) fifth modes.

Figure 5 provides crucial information about when a system behaves like a constrained dynamical system. Remember that the protein-ligand cutoff 𝛿𝑃𝐿 is responsible for the extent of the interaction between the ligand atoms and the protein residues. Positive shifts can be observed when 𝛿𝑃𝐿 is statistically within the range of [8 Å,16 Å]. As 𝛿𝑃𝐿 gets closer to the upper cutoff 𝛿𝑈 forming the protein network, the constraint effect and the accompanying positive shifts die out. In this regime, the ligands do not behave as distinct entities and coalesce with the protein molecule into a bigger


Page 15 of 23 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 network, thus only increasing the degrees of freedom of the total system, which corresponds to softening. As evident from Figure 5, the power is another determining parameter for the transition between constraining and softening regimes. For a fixed 𝛿𝑃𝐿, as the power increases, the interaction strength decreases so that the ligand atoms move more independently relative to protein residues, hence they do not impose any constraints on the global dynamics. We can elaborate more on the constraining effect in view of specific complexes by revisiting Figure 3. Both the frequencies of the complex and the pseudo-apo state for AK and c-Src monotonically progress at a similar rate when the protein-ligand cutoff is set to 10 or 12 Å. This profile has to be the case if an interlacing between the frequencies of those states is to be expected. Moreover, this trend is observed for other parameter pairs in the triangular region as well. On the other hand, extremely localized interactions at 𝛿𝑃𝐿 = 4 Å (Figure 3(a)) violates the conditions of the interlacing theorem: the most collective modes shift to the left whereas the rest shifts in the opposite direction.

Vibrational Entropy Change In accord with the discussion above, the loss in the degrees of freedom regarding the dynamics of the total system upon ligand binding may be expected to manifest itself as a negative vibrational entropy change. The vibrational entropy of a system over a set of (3𝑁 ― 6) nonzero frequencies of a specific conformation in the harmonic approximation is given in Ref.27. We will study such a change to calculate Δ𝑆𝑣𝑖𝑏 according to the following equation: Δ𝑆𝑣𝑖𝑏 = 𝑆𝐶 ― 𝑆𝑃 ― 𝑆𝐿 ,

(5)

where 𝑆𝐶 stands for the vibrational entropy of the complex. The ligand counterpart, 𝑆𝐿, is calculated by using the same parameter set, namely the upper cutoff and power, and its conformation is extracted from the corresponding complex structure. The entropy of the protein part may correspond to that of either the pseudo-apo form or the real apo crystal structure.



Page 16 of 23

16

Here, we will concentrate on two cases, namely AK and DHFR. AK undergoes a large conformational change of 7.2 Å upon binding to the ligand AP5 (complex: 1ake28, apo: 4ake33, 214 residues). The entropy changes of AK complex with respect to its pseudo-apo form and the apo crystal structure are displayed in panels (a) and (b) of Figure 6, respectively. For DHFR, the difference between its holo (DHFR+NADPH) and pseudo-apo forms is displayed in panel (c). Only the negative entropy changes are shown by masking the positive entities in panels (a) to (c). The complementary triangular regions in Figures 5 and 6 suggest that the positive shifts in the collective modes may be associated with a negative entropy change in line with a constraint effect on the dynamics of the system.


Page 17 of 23 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

Figure 6. Heat maps for different (𝛿𝑃𝐿,𝑝)-pairs showing entropy difference (kcal ⋅ mol ―1 ⋅ K ―1) for (a) AK between its complex with AP5 and its pseudo-apo structure, (b) AK between its complex and apo crystal structures, (c) DHFR between its holo form (DHFR + NADPH) and its pseudo-apo structure, (d) DHFR between its tertiary complex (DHFR + NADPH + MTX) and its holo form (DHFR + NADPH).

In Figure 6(d), we demonstrate the entropy change of the tertiary complex of DHFR with respect to its holo form. The net entropy change in this case is positive throughout all parameter pairs, which implies a softening of the holo DHFR upon MTX binding. If we go back to its frequency profile in Figure 4(a), we observe that this entropy difference is concordant with the decrease in its low frequencies. Finally, we should stress that constraining and softening effects cannot be distinguished based on sole entropy differences as the vibrational entropy is a summation over all frequencies. The dynamics picture based on frequency shifts and collective modes that is presented here is more central in terms of differentiating such effects and providing insight in terms of function.

Conclusions We presented a systematic analysis of the frequency shifts observed in the spectrum of proteins upon complexation with small ligands using a diverse dataset. We suggest that such a process can be interpreted within the perspective of constrained dynamical systems if positive shifts in the collective modes exist. As such, it is reasonable to infer that the positive frequency shifts in the collective modes could be related to the narrowing of the harmonic well associated with the bound structure, i.e. a steeper minimum. This picture conforms to both the dynamical energy landscape4 and population shift scenarios1–3. The tunable parameters of RESPEC model have been adjusted to mimic the crystal environment based on the highest correlation values with B-factors21. More specifically, the protein-ligand



Page 18 of 23

18 cutoff is lower than the upper cutoff that forms the protein network, and the power values fall in the range describing relatively strong atom-atom interactions, such as those of electrostatic ones. In this study, we show that this subset of parameters leads to significant positive shifts in the frequency spectrum of the collective modes. We interpret this as a constraining effect of the ligand binding on the protein dynamics. Another observation is the interlacing of the low-frequencies of complex and pseudo-apo structures, which can be elegantly explained by the Rayleigh-CourantFisher theorem31,32. Furthermore, we show that a negative vibrational entropy change could stem from the loss in the degrees of freedom of the complex due to constraints applied by the ligand. On the other hand, there are complexes in our dataset that present the reverse situation of softening (13% of the cases for the specific parameters chosen). Our model could also simulate an experimentally-verified softening case of holo DHFR upon MTX binding15. Another interesting finding is that the frequency distributions obtained by our coarse-grained approach exhibits a universal enclosing curve in line with the features of previous computational and experimental studies10,11,15. On the other hand, the left tail of the distribution from RESPEC contains less number of modes, making it easier to detect functional motions especially for large systems. Our scheme could be used to compare and distinguish the effects of allosteric and orthosteric ligands on the functional dynamics of enzymes. Allosteric ligands provide alternative and possibly more selective ways to regulate protein activity34, by modifying conformations and/or dynamics of proteins. In our recent computational study on triosphosphate isomerase35, binding of an inhibitor at the dimeric enzyme’s interface was shown to pose an allosteric effect by modifying the collective modes that are coupled to catalytic site dynamics. More importantly, a simple ENM-based scanning technique detected pronounced positive frequency shifts upon perturbation of residues located at this allosteric site35. As we have demonstrated in the current study, RESPEC presents a sophisticated approach to render proteinligand complexes and their frequency shifts in the collective modes upon complexation. Thus, the incorporation of a residue-based scanning into RESPEC may reveal important sites and/or ligands that perturb global modes to a large extent, thereby leading to allosteric effects that are crucial for drug design.


Page 19 of 23 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

The authors declare no competing financial interest.

Acknowledgements We would like to thank O. T. Turgut for valuable discussion at the initial stages of this work.

References (1)

Monod, J.; Wyman, J.; Changeux, J.-P. On the Nature of Allosteric Transitions: A Plausible Model. J. Mol. Biol. 1965, 12 (1), 88–118. https://doi.org/10.1016/S00222836(65)80285-6.

(2)

Tsai, C.-J.; Ma, B.; Nussinov, R. Folding and Binding Cascades: Shifts in Energy Landscapes. Proc. Natl. Acad. Sci. 1999, 96 (18), 9970–9972. https://doi.org/10.1073/pnas.96.18.9970.

(3)

Okazaki, K. -i.; Takada, S. Dynamic Energy Landscape View of Coupled Binding and Protein Conformational Change: Induced-Fit versus Population-Shift Mechanisms. Proc. Natl. Acad. Sci. 2008, 105 (32), 11182–11187. https://doi.org/10.1073/pnas.0802524105.

(4)

Frauenfelder, H.; Sligar, S.; Wolynes, P. The Energy Landscapes and Motions of Proteins. Science (80-. ). 1991, 254 (5038), 1598–1603. https://doi.org/10.1126/science.1749933.

(5)

Kar, G.; Keskin, O.; Gursoy, A.; Nussinov, R. Allostery and Population Shift in Drug Discovery. Curr. Opin. Pharmacol. 2010, 10 (6), 715–722. https://doi.org/10.1016/j.coph.2010.09.002.

(6)

Nussinov, R.; Tsai, C.-J. Allostery in Disease and in Drug Discovery. Cell 2013, 153 (2), 293–305. https://doi.org/10.1016/j.cell.2013.03.034.

(7)

Benkovic, S. J. A Perspective on Enzyme Catalysis. Science (80-. ). 2003, 301 (5637), 1196–1202. https://doi.org/10.1126/science.1085515.

(8)

Agarwal, P. K. Role of Protein Dynamics in Reaction Rate Enhancement by Enzymes. J. Am. Chem. Soc. 2005, 127 (43), 15248–15256. https://doi.org/10.1021/ja055251s.

(9)

Hay, S.; Scrutton, N. S. Good Vibrations in Enzyme-Catalysed Reactions. Nat. Chem.



Page 20 of 23

20 2012, 4 (3), 161–168. https://doi.org/10.1038/nchem.1223. (10)

ben-Avraham, D. Vibrational Normal-Mode Spectrum of Globular Proteins. Phys. Rev. B 1993, 47 (21), 14559–14560. https://doi.org/10.1103/PhysRevB.47.14559.

(11)

Na, H.; Song, G.; Ben-Avraham, D. Universality of Vibrational Spectra of Globular Proteins. Phys. Biol. 2016, 13 (1), 016008. https://doi.org/10.1088/14783975/13/1/016008.

(12)

Tama, F.; Sanejouand, Y.-H. Conformational Change of Proteins Arising from Normal Mode Calculations. Protein Eng. Des. Sel. 2001, 14 (1), 1–6. https://doi.org/10.1093/protein/14.1.1.

(13)

Bahar, I.; Lezon, T. R.; Yang, L.-W.; Eyal, E. Global Dynamics of Proteins: Bridging Between Structure and Function. Annu. Rev. Biophys. 2010, 39 (1), 23–42. https://doi.org/10.1146/annurev.biophys.093008.131258.

(14)

Uyar, A.; Kantarci-Carsibasi, N.; Haliloglu, T.; Doruker, P. Features of Large HingeBending Conformational Transitions. Prediction of Closed Structure from Open State. Biophys. J. 2014, 106 (12), 2656–2666. https://doi.org/10.1016/j.bpj.2014.05.017.

(15)

Balog, E.; Becker, T.; Oettl, M.; Lechner, R.; Daniel, R.; Finney, J.; Smith, J. C. Direct Determination of Vibrational Density of States Change on Ligand Binding to a Protein. Phys. Rev. Lett. 2004, 93 (2), 028103. https://doi.org/10.1103/PhysRevLett.93.028103.

(16)

Balog, E.; Perahia, D.; Smith, J. C.; Merzel, F. Vibrational Softening of a Protein on Ligand Binding. J. Phys. Chem. B 2011, 115 (21), 6811–6817. https://doi.org/10.1021/jp108493g.

(17)

Turton, D. A.; Senn, H. M.; Harwood, T.; Lapthorn, A. J.; Ellis, E. M.; Wynne, K. Terahertz Underdamped Vibrational Motion Governs Protein-Ligand Binding in Solution. Nat. Commun. 2014, 5 (1), 3999. https://doi.org/10.1038/ncomms4999.

(18)

Stone, M. J. NMR Relaxation Studies of the Role of Conformational Entropy in Protein Stability and Ligand Binding. Acc. Chem. Res. 2001, 34 (5), 379–388. https://doi.org/10.1021/ar000079c.

(19)

Alakent, B.; Baskan, S.; Doruker, P. Effect of Ligand Binding on the Intraminimum Dynamics of Proteins. J. Comput. Chem. 2011, 32 (3), 483–496. https://doi.org/10.1002/jcc.21636.

(20)

Acbas, G.; Niessen, K. A.; Snell, E. H.; Markelz, A. G. Optical Measurements of Long-


Page 21 of 23 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 Range Protein Vibrations. Nat. Commun. 2014, 5 (1), 3076. https://doi.org/10.1038/ncomms4076. (21)

Kaynak, B. T.; Findik, D.; Doruker, P. RESPEC Incorporates Residue Specificity and the Ligand Effect into the Elastic Network Model. J. Phys. Chem. B 2018, 122 (21), 5347– 5355. https://doi.org/10.1021/acs.jpcb.7b10325.

(22)

Dykeman, E. C.; Sankey, O. F. Normal Mode Analysis and Applications in Biological Physics. J. Phys. Condens. Matter 2010, 22 (42), 423202. https://doi.org/10.1088/09538984/22/42/423202.

(23)

Atilgan, A. R.; Durell, S. R.; Jernigan, R. L.; Demirel, M. C.; Keskin, O.; Bahar, I. Anisotropy of Fluctuation Dynamics of Proteins with an Elastic Network Model. Biophys. J. 2001, 80 (1), 505–515. https://doi.org/10.1016/S0006-3495(01)76033-X.

(24)

Kurkcuoglu, O.; Turgut, O. T.; Cansu, S.; Jernigan, R. L.; Doruker, P. Focused Functional Dynamics of Supramolecules by Use of a Mixed-Resolution Elastic Network Model. Biophys. J. 2009, 97 (4), 1178–1187. https://doi.org/10.1016/j.bpj.2009.06.009.

(25)

Wilson, K. The Renormalization Group and the ε Expansion. Phys. Rep. 1974, 12 (2), 75– 199. https://doi.org/10.1016/0370-1573(74)90023-4.

(26)

Fisher, M. E. Renormalization Group Theory: Its Basis and Formulation in Statistical Physics. Rev. Mod. Phys. 1998, 70 (2), 653–681. https://doi.org/10.1103/RevModPhys.70.653.

(27)

Levitt, M.; Sander, C.; Stern, P. S. Protein Normal-Mode Dynamics: Trypsin Inhibitor, Crambin, Ribonuclease and Lysozyme. J. Mol. Biol. 1985, 181 (3), 423–447. https://doi.org/10.1016/0022-2836(85)90230-X.

(28)

Müller, C. W.; Schulz, G. E. Structure of the Complex between Adenylate Kinase from Escherichia Coli and the Inhibitor Ap5A Refined at 1.9 Å Resolution. J. Mol. Biol. 1992, 224 (1), 159–177. https://doi.org/10.1016/0022-2836(92)90582-5.

(29)

Getlik, M.; Grütter, C.; Simard, J. R.; Klüter, S.; Rabiller, M.; Rode, H. B.; Robubi, A.; Rauh, D. Hybrid Compound Design To Overcome the Gatekeeper T338M Mutation in CSrc #. J. Med. Chem. 2009, 52 (13), 3915–3926. https://doi.org/10.1021/jm9002928.

(30)

Sawaya, M. R.; Kraut, J. Loop and Subdomain Movements in the Mechanism of Escherichia Coli Dihydrofolate Reductase: Crystallographic Evidence † , ‡. Biochemistry 1997, 36 (3), 586–603. https://doi.org/10.1021/bi962337c.



Page 22 of 23

22 (31)

Arnold, V. I. Mathematical Methods of Classical Mechanics; Graduate Texts in Mathematics; Springer New York: New York, NY, 1989; Vol. 60. https://doi.org/10.1007/978-1-4757-2063-1.

(32)

Horn, R. A.; Johnson, C. R. Matrix Analysis; Cambridge University Press: Cambridge, 2012. https://doi.org/10.1017/CBO9781139020411.

(33)

Müller, C.; Schlauderer, G.; Reinstein, J.; Schulz, G. Adenylate Kinase Motions during Catalysis: An Energetic Counterweight Balancing Substrate Binding. Structure 1996, 4 (2), 147–156. https://doi.org/10.1016/S0969-2126(96)00018-4.

(34)

Lee, G. M.; Craik, C. S. Trapping Moving Targets with Small Molecules. Science (80-. ). 2009, 324 (5924), 213–215. https://doi.org/10.1126/science.1169378.

(35)

Kurkcuoglu, Z.; Findik, D.; Akten, E. D.; Doruker, P. How an Inhibitor Bound to Subunit Interface Alters Triosephosphate Isomerase Dynamics. Biophys. J. 2015, 109 (6), 1169– 1178. https://doi.org/10.1016/j.bpj.2015.06.031.


Page 23 of 23 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

TOC Graphic


ProteinâLigand Complexes as Constrained Dynamical Systems

Recommend Documents

ProteinâLigand Complexes as Constrained Dynamical Systems