Regulation of ryanodine receptor RyR2 by protein-protein interactions: prediction of a PKA binding site on the N-terminal domain of RyR2 and its relation to disease causing mutations

Docking predictions

We used the commercial software Gold¹⁵ for docking the peptides to the surface of RyR2. The PKA chain (PDB code 2JDV) of 336 amino acids is partitioned into a library of 331 overlapping hexapeptides, such that the first peptide consists of the first six residues 1–6, the second of 2–7, and so on. Several binding sites are chosen on the surface of RyR2 as discussed below. A radius of 20 Å is used for docking. The GoldScore force field is used with rescoring on ChemScore. Flexible docking is used in the first round of calculations. Peptides with reasonable docking energies are chosen after the first run, and a more thorough and extensive docking is performed over this smaller subset. Additional calculations are made with hexapeptide libraries obtained from modulators of RyR2 that are known not to bind at the N-terminal domain as a partial check of the reliability of the method. Optimum binding is obtained for the hexapeptide FKGPGD from the residues 318–323 of 2JDV. The binding energy is obtained as -49 kJ/mol, which is significantly stronger than those of all other investigated hexapeptides. This binding energy corresponds to a dissociation constant k_D = exp(ΔA/kT) of 5.5 nM.

Our algorithm for the Elastic Net Model uses C^α based coarse graining which evaluates correlations between thermal fluctuations Δ$\stackrel{̭}{R}$_i and Δ$\stackrel{̭}{R}$_j in the position of residues i and j. On average, a residue has about eight to 12 neighboring residues to directly interact with. These fluctuation-based interactions are assumed harmonic as if the residues are connected by linear springs. Fluctuations in the distance between two neighboring residues induce changes in their interaction energy. Two residues are assumed neighbors in space if they are closer to each other than a given cutoff distance. This distance corresponds to the radius of the first coordination shell around a given residue, and is usually thought to be between 6.5–7.0 Å. Every pair of residues closer to each other than the cutoff distance is assumed to be connected by a linear spring. The knowledge of the tridimensional structure of the protein that has n residues allows us to write a connectivity matrix, C, where the rows and the columns identify the residue indices, from 1 to n, where the amino-end is the starting and the carboxyl-end is the terminating-end of the protein. If two residues i and j are within the cutoff distance, then C_ij = 1, otherwise it is zero. Another matrix, Γ_ij, is obtained from the connectivity matrix as

Where γ is the spring constant of the harmonic interactions. The relationship of the forces to the displacements is given by the equation ΔF_i = $\underset{j}{\Sigma }$ Γ_ijΔR_j. Techniques of statistical mechanics allow us to derive several relationships between the fluctuations of residues¹⁶. The correlation between the fluctuations of residues i and j is related, for example, to the inverse of the matrix Γ_ij as

Here, the angular brackets denote the time average of the product of fluctuations of residues i and j, k_B is the Boltzmann constant, T is the physiological temperature expressed in Kelvin scale, Γ^-1 is the inverse matrix Γ, and its subscripts i and j acknowledge the residue indices of interest. If i = j, then Equation 2 becomes

The left hand side of Equation 3 is the mean-square fluctuations of the i’th residue which is related to experimentally available B-factors, B_i, through the equation

The mean-square fluctuations 〈(ΔR_ij)²〉 of the distance ΔR_ij between residues i and j are obtained from Equation 3 as

The derivation of Equation 5 is given in the reference¹⁷. Calculations presented therein showed that the largest few eigenvalues of Γ give information at the residue level. Smaller eigenvalues are associated with global motions. Our calculations showed that the largest five eigenvalues and the corresponding eigenvectors are satisfactory for representing fluctuations at the residue level.

Fluctuations of the harmonic energy between two residues are proportional to the mean square fluctuations of the distance between the two. Thus, Equation 5 is representative of energy fluctuations, and summing over all the neighbors of the residue i shows the energy response ΔU_i of residue i with its surroundings:

This is a thermodynamically meaningful quantity showing the mean energy response of residue i to all fluctuations of its surroundings. These correlations extend throughout the protein, leading to specific paths along which the fluctuations propagate. Recent work shows that these paths are evolutionarily conserved^14a.

The N-terminal domain of RyR2 is a signal protein of 217 amino acids. The crystal structure of the N-terminal domain of physiological RyR2 (PDB code 3IM5) and the A77V mutated crystal structure (PDB code 3IM7) have been determined by x-ray with resolutions of 2.5 and 2.2 Å, respectively, by Van Petegem and Lobo^3a. The protein consists of a β-trefoil of 12 β strands held together by hydrophobic forces. A 10-residue α helix is packed against strands β4 and β5. A 3 residue 3–10 helix is present in the loop containing β3 and β4. The N-terminal contains two MIR domains, similar to the inositol 1,4,5-triphosphate receptor (IP3R), for which ligand-induced conformational changes have been studied more extensively¹⁸.

Docking results

The binding free energy of FKGPGD to the surface shown in Figure 2 is obtained as -49 kJ/mol by the ChemScore potential, which corresponds to a dissociation constant of 5.5 nM. The 42% of the binding energy comes from hydrogen bonds and 39% from lipophilic interactions. The dissociation constant of 5.5 nM is at least two orders of magnitude better than the values obtained for the other hexapeptides of the library. It is therefore highly likely that PKA anchors itself on RyR2 at the position shown.

Figure 2. The bound conformation of FKGPGD, shown in yellow ball and stick.

Residues with which it forms hydrogen bonds are shown in yellow wire, and labeled. The two disease causing mutation residues, ALA77 and ARG176 are shown in yellow CPK.

The energy conduction path of RyR2

In order to interpret the binding of the PKA on RyR2, we performed elastic net analysis of energetically responsive residues of RyR2. The residues that yield high values of the energy response defined by Equation 6 are calculated according to the scheme outlined in the Methods section. In Figure 3, the mean energy response ΔU_i of residue i is presented along the ordinate as a function of residue index. The circles indicate the highest conserved residues of 3IM5, obtained from the work of Goldenberg et al. (See also the PDBSum web site²²).

Figure 3. Energetically responsive residues (solid line) obtained with the Elastic Net Model, and the conserved residues (circles) obtained from Reference 22.

In Reference 20, conservation levels are ordered from 1 to 8, the latter being the highest degree of conservation. The filled circles correspond to residues with level 8. The ordinate values are in arbitrary un-normalized units.

Comparison of the solid curve peaks and the circles shows that there is a strong correlation between the energy responsive and conserved residues, in agreement with the recent suggestion of Lockless and Ranganathan^14a. The set of conserved residues, with the highest level of conservation according to Reference 20 of the protein, all lie within the set of energetically responsive residues and are located along or in the neighborhood of the path obtained from the energetically responsive residues. On the three-dimensional structure of the protein, the peaks shown in Figure 3 constitute a path of residues that are spatial neighbors.

A residue or set of residues at the surface of the protein which are energy responsive are expected to be the hotspots for binding, because these residues can exchange energy with the surroundings, and distribute the energy taken from the surroundings to the other residues of the protein. According to this conjecture, one needs to dock ligands only to the hotspots identified with the peaks in Figure 3. In our calculations, we adopted five such hotspot regions for docking. These hotspot regions are centered at: (1) VAL21, (2) VAL68, (3) ARG122, (4) SER185, and (5) ALA205. In the complex structure of the channel, some of these five surface regions may not be exposed to ligands but may be facing the other domains of the channel. However, a residue that neighbors another domain may become exposed to a ligand upon opening of the channel. We carried out the calculations for the five regions stated above, irrespective of their neighborhood.

In Figure 4, we show, in stick form, the evolutionarily highly conserved residues that lie along a path between ALA77, ARG176 and the ligand FKGPGD of PKA. The peak residues clearly form a path between the ligand and the mutation residues. The path shown in the figure contains the energetically responsive residues predicted by the GNM as may be seen from Figure 3. Using extensive docking calculations and libraries of residues obtained from regulator proteins of the RyR2 channel, we showed that residues 318–323 of PKA have a very high affinity for the N-terminal of RyR2. The location of binding is a pocket bordered by GLU171 and GLU189. GLU171 is a conserved residue and participates in calcium binding in inositol 3 receptors, IP3R. However, a ligand for RyR2 at GLU171 is not yet known. We also showed that the disease causing mutations ALA77VAL and ARG176GLN are joined by an energy interaction pathway to the ligand binding surface. Although these two mutations are responsible for arrhythmias, their exact mechanism is not known. The present model directs attention to the relationship between the residues at the binding site, the predicted path of energy responsive residues and the two disease causing mutation sites. Since binding of PKA to RyR2 results in phosphorylation of the latter, and since hyperphosphorylation leads to disease, one may indirectly conjecture that mutations in the two residues modify the binding characteristics of PKA.

Figure 4. Energy interaction paths from ALA77 and ARG176 to the ligand.

The residues shown in stick form are conserved residues which are also identified by the peaks in Figure 3. The hexamer ligand is shown in ball and stick form.

Relative orientations of RyR2 and PKA in bound form

Superposition of the three dimensional PDB structures of PKA and RyR2 in such a way that the residues FKGPGD of PKA are kept in the bound state gives the relative orientations of the two proteins. This is shown in Figure 5.

Figure 5. Relative orientations of RyR2, shown in surface, and PKA, shown in solid ribbon.

The sequence FKGPGD of PKA is shown in ball and stick form.

Using the Elastic Net Model, we identified the energy conduction pathway for the wild type RyR2. This path whose residues are shown in Figure 3 has several features of interest. Firstly, it contains most of the evolutionarily conserved residues. The remaining conserved residues are in the close neighborhood of the path, all making hydrogen bonds with the residues of the path. This important feature has recently been shown by Lockless and Ranganathan^14a implying that evolutionary conservation is driven by energy conduction in proteins. Although no ligands for the RyR2 N-terminal have been observed until now¹⁹, the three glutamic acids, GLU171, GLU173 and GLU189 at the pocket may potentially form a binding site. This suggestion is also based on the observation that in IP3R a potential calcium binding site is formed by GLU283 and GLU285 whose location on the path coincides exactly with that of RyR2. The residue GLU173 is exposed to water and is a candidate for possible binding.

The underlying determinant of the allosteric pathway, defined as the path of energy responsive residues in the present paper, is the graph structure of the protein²⁰. The view that proteins relay signals by energy fluctuations in response to inputs, have been recently discussed in an elegant paper by Smock and Gierasch^14b. In the present study, we showed that evolutionarily conserved residues lie on the pathway of energy responsive residues. RyR2 contains two interspersed MIR domains, residues 124–178 and 164–217²¹. Elastic net calculations show that the residues that lie on the energy conduction pathway are closely associated with these MIR domains: the energy responsive residues either lie on the MIR domains, or they are hydrogen bonded to the residues of these domains. There is no energetically responsive residue that is not closely associated with the MIR domain. We therefore conclude that the MIR domains of RyR2 play an active role in energy transport through the protein.