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

Protein-protein interactions are the key processes responsible for signaling and function in complex networks. Determining the correct binding partners and predicting the ligand binding sites in the absence of experimental data require predictive models. Hybrid models that combine quantitative atomistic calculations with statistical thermodynamics formulations are valuable tools for bioinformatics predictions. We present a hybrid prediction and analysis model for determining putative binding partners and interpreting the resulting correlations in the yet functionally uncharacterized interactions of the ryanodine RyR2 N-terminal domain. Using extensive docking calculations and libraries of hexameric peptides generated from regulator proteins of the RyR2 channel, we show that the residues 318-323 of protein kinase A, PKA, have a very high affinity for the N-terminal of RyR2. Using a coarse grained Elastic Net Model, we show that the binding site lies at the end of a pathway of evolutionarily conserved residues in RyR2. The two disease causing mutations are also on this path. The program for the prediction of the energetically responsive residues by the Elastic Net Model is freely available on request from the corresponding author.


Introduction
Ryanodine receptors are large protein complexes consisting of approximately 5000 residues that form calcium channels that mediate the release of calcium from the sarcoplasmic reticulum, SR, to the cytosol, which is essential for muscle and cardiac rhythm and contractility. There are three forms of ryanodine receptors, RyR1, RyR2 and RyR3. RyR1 is the channel in the skeletal muscle, RyR2 is the type expressed in the heart muscle, and RyR3 is found predominantly in the brain 1 . The present paper focuses on RyR2. Ca ++ release from the SR mediated by RyR2 is a fundamental event in cardiac muscle contraction. These receptors form a group of four homotetramers, with a large cytoplasmic assembly and a transmembrane domain called the pore region. The tridimensional structure of the full assembly is known from cryo-electron microscope studies 2 with limited precision. However, the crystal structures of the first 520 amino acids of the N-terminal domain of RyR1 and the first 217 amino acids of the N-terminal domain of the wild type RyR2 and its mutated form are determined with high precision by van Petegem and collaborators 3 . The main mass of the receptor with dimensions of ca. 280 × 280 × 120 Å is located in the cytoplasmic region, with a stalklike transmembrane region 2 . The full shape of the channel and the N-terminal are shown in Figure 1.
The cytoplasmic region consists of more than 10 sub-domains that are responsible for the functioning of the receptor through binding to several modulator proteins and ligands 4 . The modulators include cyclic AMP and protein kinase A (PKA) 4 , calmodulin 5 , FKBP12.6 (Calstabin2) 6 , phosphatases 1 and 2A (PP1 and PP2A) 7 , sorcin 8 , and triadin, junctin and calsequestrin 9 , and several others. Among these, cyclic AMP activates PKA, which in turn phosphorylates RyR2 at SER2809 and SER2815. Despite the important role of the channel, the binding sites of the modulators on the channel are known only approximately. Calmodulin binds to residues located between the positions 3611 and 3642, FKBP12.6 binds to residues around the positions 2361-2496, PP1 around 513 and 808, PP2A around 1451 and 1768, sorcin, triadin, junctin and calsequestrin bind to the vicinity of the transmembrane domain 7 .
FKBP12.6 binds to RyR2 with a stoichiometry of four FKBP12.6 molecules per single RyR2 channel complex. Binding of FKBP12.6 to RyR2 is required to keep the receptor closed during diastole. In addition to stabilizing individual RyR channels, FKBP12.6 is also required for coupled opening and closing between RyRs. Dissociation of FKBP12.6 from coupled RyR2 channels results in functional uncoupling of the channels leading to heart failure 4 . Overphosphorylation of RyR2 leads to dissociation of the regulatory protein FKBP12.6 from the channel, resulting in disease 7 exhibited as arrhythmias with abnormal diastolic SR Ca++ release. Uncontrolled Ca++ release during the diastole when cytosolic Ca++ concentrations are low can cause delayed after-depolarizations (DADs) which can then lead to fatal arrhythmias. These abnormalities are linked to mutations in the RyR2, located on chromosome 1q42.1-q43 10 , which lead to familial polymorphic ventricular tachycardia, CPVT, and arrhythmogenic right ventricular dysplasia type 2, ARVD/C. More than 300 point mutations have been identified in RyR2, some of which are associated with the disorders observed clinically 11 . In this respect, the N-terminal domain of RyR2, which is known to form an allosteric structure, contains several disease-causing mutations. However, there is yet no information on the mechanisms of the mutations that lead to disease and on the role of these mutations on modulator binding.
None of the modulators discussed above, except PKA, bind to the N-terminal domain. PKA phosphorylates Ser2809 and Ser2815, and it has to anchor to nearby regions of the two serines. PKAs are known to anchor to their hosts at points other than the catalytic domains 12 . In this study, we generated a hexameric peptide library from the PKA and docked these on several points on the surface of the RyR2 N-terminal. Calculations showed that the hexapeptide PHE LYS GLY PRO GLY ASP from the unstructured C-terminal region of PKA binds to RyR2 with very high affinity, with a dissociation constant of 5.5 nM. For brevity, we will refer to this hexapeptide as the 'ligand' and represent it in single letter convention as FKGPGD.
In the last part of the paper, using a coarse grained Elastic Network Model 13 , we show that the binding site of the ligand lies on a path of energy responsive residues. Energy responsiveness of a residue is defined in terms of correlated fluctuations of that residue with others in the protein. In allosteric proteins, a path of highly correlated residues exists and plays crucial role in energy and signal transfer 13a,14 . In RyR2 we identify such a path of highly correlated residues which contains most of the evolutionarily conserved residues. The path also contains the known two disease causing mutations, A77V and R176Q.

Docking predictions
We used the commercial software Gold 15 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 ΔR ̭ i and Δ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 = å j Γ ij ΔR j . Techniques of statistical mechanics allow us to derive several relationships between the fluctuations of residues 16 .
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 Bfactors, B i , through the equation The mean-square fluctuations 〈(ΔR ij ) 2 〉 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 17 . 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 18 .

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.

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  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.
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.
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 19 , 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 20 . 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 21 . 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.

Data availability
Data of predicted PKA binding sites on RyR2 23 .

Author contributions
Both authors contributed equally to the present study.

Competing interests
No competing interests were disclosed.

Grant information
The author(s) declared that no grants were involved in supporting this work.  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 I tend to agree with the other reviewer that such a large affinity between the PKA kinase and RyR2 is strange, since kinase interactions are expected to be transient. However, I think that it is unlikely that such a high affinity arises by chance. If we assume that it is correct, it may suggest a molecular mechanism for the disease, as follows. After phosphorylation, the kinase has to dissociate, which would be very difficult given the high affinity, unless a conformation change takes place upon phosphorylation and triggers the release of PKA. This mechanism may be consistent with the network of highly dynamic residues predicted by the ENM, which connects the binding site and the disease-associated site. If the mutations perturb the dynamics of RyR2, making the release of the kinase more difficult, they would result in overphosphorylation, which is the known signature of the disease. It would be interesting to investigate this hypothesis through ENM (and/or MD) computations that try to connect the phosphorylation site and the binding site of the kinase through a signaling pathway and test whether this pathway passes through the disease associated residues, as the computations presented here may suggest.
Concerning the presentation, the result that energy-responsive residues are strongly conserved is very interesting. The authors relate it with the influential proposal by Lockless and Raganathan (LR) that pairs of residues energetically coupled are evolutionarily coupled as well (actually this field of correlated mutations is quite older than the LR article and it recently experienced an explosion of interest), but the LR proposal is not exactly the same as what is found here, since in this paper conservation of individual residues rather than evolutionarily correlations are observed. Thus, it would be interesting to study whether this correlation between responsiveness and conservation is general, or it is a special property of RyR2 and maybe other strongly allosteric proteins. Of course I do not require that the authors do this study in this paper, but it would be useful that they share their thoughts on this point. Since the responsiveness of a residues must be correlated with its intrinsic fluctuations (predicted B factors), it would be interesting to known whether the former or the latter correlates more strongly with conservation.
Furthermore, the authors expect that responsive residues are hotspots of binding. This expectation does not seem so natural to me. Actually, I expected almost the opposite: two residues i and j that bind the same molecule must have small fluctuations of their interatomic distance <(DeltaR_ij)^2>, so that I expected that residues in a binding site must be characterized by <(DeltaR_ij)^2> smaller than generic residues. I tested this expectation some time ago, finding that it is true on the average but much weaker than I expected. Thus, the fact that responsive residues are close to binding sites can be interpreted either as a general property or as a special property of allosteric proteins in which binding and conformation changes are related. It would be interesting to discuss it.
Finally, although the paper is very well written, there are some small typos in the equations: 1) In Eq. (1) the factor gamma has to appear in the second line as well; 2) Temperature has disappeared from Eq. (5) ; 3) between Eq. (5) and (6), it is said that the five largest eigenvalues are enough for representing fluctuations, I think that the authors mean "smallest eigenvalues", and it would be useful to say that the number five holds for RyR2, other proteins may have a different number of relevant modes I have read this submission. I believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.
No competing interests were disclosed.