RESPEC Incorporates Residue Specificity and the Ligand Effect into

Subscriber access provided by UCL Library Services

Article

RESPEC Incorporates Residue Specificity and Ligand Effect into Elastic Network Model Burak T Kaynak, Doga Findik, and Pemra Doruker J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b10325 • Publication Date (Web): 22 Dec 2017 Downloaded from http://pubs.acs.org on December 26, 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.

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

The Journal of Physical Chemistry

1

RESPEC Incorporates Residue Specificity and Ligand Effect into Elastic Network Model Burak T. Kaynak∗†, Doga Findik‡, Pemra Doruker*‡

†Department of Physics, Bogazici University, 34342, Bebek, Istanbul, Turkey ‡Department of Chemical Engineering and Polymer Research Center, Bogazici University, 34342, Bebek, Istanbul, Turkey

*E-mail: [email protected]; [email protected]

ACS Paragon Plus Environment

The Journal of Physical Chemistry 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 32

2

Abstract RESPEC is a new framework that introduces residue specificity into elastic network modeling (ENM) to successfully render intact protein-ligand complexes as well as apo proteins. This framework establishes a broader application of coarse-graining idea via describing: (i) a coarsegrained residue/node through its heavy atoms as virtual nodes, (ii) an effective B-factor for such a node, directly obtained from the experimental data, and (iii) a node-node interaction by a cumulative distance-dependent force constant. RESPEC improves the level of correlations with B-factors after optimizing the parameters of the model. In the absence of ligands, the mean correlations excel 0.72, which is higher than the classical ENM results, based on a diverse set of proteins. Global modes satisfactorily describe the conformational transitions for apo structures. When the ligands are included at atomistic resolution in RESPEC calculations, mean correlation values exceed 0.9 over the same dataset.


Page 3 of 32 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 Ligand binding affects the conformations sampled by the protein and modifies its vibrational dynamics. The extent of conformational changes due to binding can be determined by comparing the apo and ligand-bound experimental structures. Assuming such structures represent local minima of the complex potential energy landscape of a protein, the vibrational dynamics of a specific conformational state (apo or ligand-bound) can be calculated by normal mode analysis (NMA). Classical NMA1–3 have been originally employed using empirical all-atom forcefields, which were followed by coarse-grained approaches named collectively as elastic network models (ENM). Such studies have shown that the lowest-frequency modes of the spectrum consist of global motions that dictate protein function and conformational changes4–7.

Originally introduced at the atomistic level8, ENM’s utility for large biological systems, such as supramolecular structures determined at low resolution, has become apparent, since the global motions can be captured at high levels of coarse-graining9. Gaussian10,11 and anisotropic12,13 network models (GNM and ANM) represent the simplest, one and three-dimensional versions of ENM, respectively, where the alpha-carbons are commonly chosen as the nodes of the network and close-neighboring node pairs are connected by springs with uniform force constants. Other versions of ENM have been developed with modified force constants, such as distancedependent functions14–17 to mimic short- and long-range interactions. Moreover, different types of nodes, either residue side chains18, or rigid blocks made of several residues19 can be used for representing nodes of the network. Common ENM versions have been recently evaluated based on collective dynamics from MD data20.



Page 4 of 32

4 Applications of ENM to protein-ligand complexes have generally considered only the protein, especially if the ligand is comparatively small. As such, the effect of conformational changes on vibrational dynamics is assessed based on apo and bound conformations of the protein part alone. However, the presence of the ligand may induce additional constraining or softening effect on protein dynamics. Recently, we have applied a mixed-resolution ENM approach21 to a protein-ligand system22, specifically for analysis of molecular dynamics trajectories of triosephosphate isomerase bound to an inhibitor. In such a mixed-resolution system, the protein is modeled by coarse-graining at the residue-level (alpha-carbon nodes), whereas the ligand is treated at atomistic detail excluding hydrogen atoms. Alternatively, the effect of certain entities in a complex or the environment on protein dynamics can be studied by dividing the overall system into subsystems and then integrating out the degrees of freedom corresponding to one of its subsystems, which can be performed in a controlled way despite being computationally more intensive23–26.

In order to perform vibrational analysis of proteins in apo or ligand-bound form in a systematic and efficient way, we introduce here a new model RESPEC, in which the harmonic interactions between node pairs are non-uniform across the network and distinguished depending on whether they belong to atoms of a ligand or a protein residue. This specialization is based on a broader application of the mixed-resolution approach regarding the types of the nodes, and is accomplished by the following steps: (i) local coarse-graining via appropriate mass-scaling of nodes’ coordinates and considering its heavy atoms as virtual nodes, (ii) determining the strength of pairwise interactions, either through a cumulative or a single power law function of atomatom distance, (iii) introducing a second cutoff that modifies the range of protein-ligand


Page 5 of 32 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 interaction, and (iv) coarse-graining the B-factors for protein residues, as it will be discussed in detail in Section titled RESPEC.

RESPEC is assessed using a large dataset of more than 300 proteins of various sizes, oligomerization states and topologies. Computational efficiency of RESPEC allows the study of a wide range of parameter values among which a relatively small set can be determined to assure high level of correlation with experimental data for both apo and complex structures.

Theoretical Methods In this section, we will study our model RESPEC, which is a novel extension of elastic network model (ENM). The main ideas of ENMs, related to our discussion on RESPEC, will be summarized before going into details of our model. As ENM is a broadly studied model, the pertinent literature is vast. Therefore, the readers are invited to consult the review27,28, and the references therein for extensive studies on this subject.

ENM Elastic network models briefly are all about studying a dynamical problem of molecular oscillations. Molecules are modeled as systems of sites interacting via simple harmonic potentials. In this picture, sites, which are widely treated as nodes of a network, are chosen to be -carbons of protein residues, and they oscillate about their equilibrium positions in a molecule.

The total potential energy of such a system consisting of nodes is defined by the sum of all

harmonic interactions between the nodes of unique , -pairs, where 1 ≤ , ≤ . Each



Page 6 of 32

6 potential is a quadratic function of the distance between the generalized coordinates of an , pair, and , respectively. Here each has a 3-dimensional spatial index as , , , albeit

being suppressed for the convenience. Let be the generalized coordinate of the node , and let

be the equilibrium position thereof, alternatively referred as a stable conformation in the context of protein structures. The potential energy of the , -pair is given by

=

− − − , 2

(1)

where and | ∙ | stand for the interaction strength and the Euclidean norm, respectively. The former is subjected to some conditions, depending on the flavor of ENM. Although the interaction strength is initially defined as a theoretical parameter, it turns out to be phenomenological as well. The total potential energy of the system in accord with the harmonic approximation near a stable equilibrium is, then, expressed as

1 = 2

! "

# # ,

(2)

where # = − is a small perturbation from the equilibrium position of the th node.

The term with the second order partial derivative in Eq. 2 is called the Hessian, $, which is a 3 × 3 matrix. As a standard practice, it is more convenient to work with a new set of

generalized coordinates, obtained by scaling 's in the following fashion: ' = ()/ , where (

is a diagonal mass matrix of dimension 3 × 3. These coordinates are called mass-weighted coordinates. Under this change of set of coordinates, the Hessian also becomes mass-weighted,

+ = ( ,)/ $(,)⁄ , as well. Now it contains all of the required mechanical parameters of the $ system from which the dynamics thereof can be inferred. Diagonalization of the Hessian matrix allows us not only to obtain the normal modes but also to have the normal coordinates of the


Page 7 of 32 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 system, in which the original system decouples into independent harmonic oscillators. This diagonalization is achieved by a similarity transformation, whose transformation matrix is a 3 × 3 special orthogonal matrix ., and each row thereof is an orthonormal eigenvector.

+ . = 0. Here 0 is an eigenvalue matrix, which is a diagonal The transformation is as follows: . / $

matrix of dimension 3 × 3, whose elements are 12 for 1 ≤ 3 ≤ 3. Since it is assumed that the system is not in an external field to interact with, not all of the eigenvalues are related to

the oscillations of the nodes. Specifically, six degrees of freedom out of 3 correspond to

translation and rotation of the molecule so that only 3 − 6 degrees of freedom are actually

nontrivial and determine the oscillatory motion of the system.

The set of nontrivial eigenvalues and corresponding eigenvectors are required for the calculation of the mean square fluctuation (MSF) of each node as well. They are of particular interest due to two reasons: firstly, MSFs allow one to assess how well the model renders the molecule in question by comparing those fluctuations with the experimental Debye-Waller factors (B-factors) coming from X-ray crystallography, and secondly, this comparison fixes/normalizes not only the free parameters of the model but also quantities with physical meaning, as it will be seen later in

this section. It is straightforward to calculate the MSF of the th node's coordinate, 〈# 〉 (not mass-weighted), through statistical mechanics, and it is given by

78 9 .2 〈# 〉 = , : ;2 2!)

(3)

where 78 is the Boltzmann constant, 9 is the absolute temperature, ;2 = ?12 is the frequency of the 3th normal mode, : is the corresponding mass of the th node, and .2 is its coordinate of



Page 8 of 32

8 the 3th normal mode in an orthonormal basis. B-factors indicate the relative vibrational motion of individual atoms around their stable conformations in molecules, thus linearly proportional to their mean square fluctuations. They take the following form for the anisotropic network model (ANM),

8B @ = C# D . 3

(4)

The factor 1/3 in the definition above takes place on the grounds that the experimental temperature factors given in Protein Data Bank (PDB)29 files are isotropic30.

Thus far, the interaction strength has not been defined due to being model specific. In ANM, it is

generally chosen to be uniform over the network, call it , and is only nonzero if the Euclidean distance between the nodes of a pair is less than a spatial cutoff. Since the characteristic equation of the Hessian matrix is not affected by an overall constant multiplier, can initially be set to 1.

Its value is later determined by matching the averaged values of both MSFs and B-factors over the whole network, hence being called experimental spring constant. That's also the reason why this parameter is also phenomenological by its very nature. Specifying in this manner leads to the normalization of the eigenvalues, 12 → 12 , as well.

The incorporation of coarse-graining procedure is another feature shared among all elastic network models. In ANM, G H atoms are generally taken as the nodes of the network, which will

be taken as reference for our model. Furthermore, the mass of each G H is also set to 1 for the sake of simplicity.


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

RESPEC Our model RESPEC is, in principle, based on the notion of the classical anisotropic network model. It is developed so as to represent protein-ligand complexes with high correlation between B-factors and MSFs. It turns out that RESPEC also conforms better to the experimental temperature factors, even if either apo or pseudo-apo states are of the main concern. Here “pseudo-apo” state is used for the protein part of a complex.

Definition of non-uniform coupling constant

RESPEC has a system wide spatial cutoff as ANM does, and we call it upper cutoff IJ . This cutoff is responsible for constructing the network. In other words, it determines unique

interacting pairs such that an , -pair is not formed unless − < IJ . On the other hand, RESPEC essentially differs from ANM in its two strands: firstly, the definition of its interaction strengths between nodes, and secondly, the way how the coarse-graining technique is applied to proteins and complexes. Coarse-graining technique can conceptually be viewed as an approximation to deal with systems with large degrees of freedom within the framework of statistical mechanics and renormalization group theory31–33. In this perspective, while the scale at which the system is studied evolves, the parameters of the system already defined at the previous scale should accordingly transform as well. In ANM this approach is uniformly applied. If the constituents of a system, however, differ in degrees of freedom as being in protein-ligand complexes, then uniform coarse-graining of such a system becomes ambiguous. Instead, we propose a new approach to this problem, in which interaction strengths between nodes are no longer uniform across the network. Henceforth, we prefer calling them coupling constants. Not



Page 10 of 32

10 only is each coupling constant exclusive to a unique pair, but also differs if the pair in question is a residue-residue or a residue-atom pair. Once the network is build up by means of IJ , we can

determine the value of each coupling constant of unique pairs. Since harmonic approximation

yields an effective model, we should assign a value to each in the spirit of coarse-graining method mentioned above with the knowledge of atomic interactions. The following function seems to be a natural choice regarding the discussion above,

=

O M∀R,S∈, M " M [\]^_`a

Q P T if , ∈ XY, YZ , QRS U

Q P T if , ∈ XY, bZ , QRS U

N M∀R,S∈, U M Q M P T if , ∈ Xb, bZ , Q L

(5)

where QRS is the Euclidean distance between the stable conformations of nodes and , c ∈ ℝe ,

and Ifg is another spatial cutoff, namely protein-ligand cutoff. Q is chosen to be 1 Å for the

future convenience. Here XY, bZ denotes the set of all pairs formed by a residue and a ligand

atom, and the other cases are self-explanatory. ∀7, h ∈ , in Eq. 5 also denotes that the sum

should be performed over all possible connections that can be built up between the nodes and .

Here what it is meant by “all possible connections” is that if at least one of the nodes in the , -

pair is a residue, then we take the effects of all individual heavy atoms in that residue into account to calculate the effective coupling constant of that pair as if those atoms were incorporated into the network. Therefore, one can consider these types of atoms in a residue as virtual nodes. If either an apo or a pseudo-apo structure is of interest, then only the first line in the definition above is relevant. However, if just one of the nodes in a pair is a residue, then the interaction between those nodes are assumed to represent an interaction between a protein


Page 11 of 32 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 residue and a ligand atom. And we constrain this type of interaction in such a way that an individual interaction between a node of ligand atoms and a virtual node of a residue can happen as long as the Euclidean distance between them is less than the protein-ligand cutoff Ifg . In this

way, one can control how a ligand can interact with surrounding protein residues. Furthermore, it allows us to represent protein-ligand complexes better since the model becomes more sensitive to the geometry of both parts, and thus making it more exclusive to the complex.

As we mention above, the role of coarse-graining technique can also be seen as a way to deal with the large degrees of freedom, and the sum in Eq. 5 specifically serves this purpose. In some sense, it integrates out the effect of fine details throughout the interactions occurring via the virtual nodes introduced above. We would like to remind that physical parameters in a model are also supposed to transform accordingly during a change of scale. It means that we have to reconsider not only masses but also B-factors of protein residues. In ANM, mass of a protein

residue is either set to 1 for simplicity or set to the mass of a G H . If the latter is chosen, then the mass-weighted Hessian is calculated as regards to it. However, none of these choices is consistent within the aforementioned framework discussed for RESPEC anymore. There should be a hierarchical difference between the scaling of the generalized coordinates of protein residues and the ones for ligand atoms. It is apparent that the best way to preserve this hierarchy is to scale each coordinate of nodes with its mass. Thus, one should scale the residue coordinates with their corresponding residue masses while scaling the coordinates of ligand atoms with their

atom masses. To achieve this, let us introduce the generalized mass matrix, (i . This time its

j m diagonal elements can be represented by the vector :)j , … , :

RESPEC Incorporates Residue Specificity and the Ligand Effect into

Recommend Documents