Study of the Affinity between the Protein Kinase PKA and Peptide Substrates Derived from Kemptide Using Molecular Dynamics Simulations and MM/GBSA

We have carried out a protocol in computational biochemistry including molecular dynamics (MD) simulations and MM/GBSA free energy calculations on the complex between the protein kinase A (PKA) and the specific peptide substrate Kemptide (LRRASLG). We made the same calculations on other PKA complexes that contain Kemptide derivatives (with mutations of the arginines, and with deletions of N and C-terminal amino acids). We predicted shifts in the free energy changes from the free PKA to PKA-substrate complex (ΔΔGE→ES) when Kemptide structure is modified (we consider that the calculated shifts correlate with the experimental shifts of the free energy changes from the free PKA to the transition states (ΔΔGE→TS) determined by the catalytic efficiency (kcat/KM) changes). Our results demonstrate that it is possible to predict the kinetic properties of protein kinases using simple computational biochemistry methods. As an additional benefit, these methods give detailed molecular information that permit the analysis of the atomic forces that contribute to the affinity between protein kinases and their substrates.


Introduction
The reversible phosphorylation of proteins, catalyzed by protein kinases (PKs), regulates the most important processes in the cell such as the gene transcription and metabolic pathways [1]. On the other hand, aberrant phosphorylation is associated with several disorders. For this reason, PKs have become important targets for rational drug design against cancer, Alzheimer, etc [2,3]. The phosphoryl-transfer mechanism catalyzed by PKs has been studied in several works using computational chemistry methods [4,5]; however, no work on computational chemistry studied PKs specificity and selectivity for substrates, which are very important topics since the differential recognition of substrates is essential for the adequate regulation of the cell. The specificity of PKs, referring to discrimination between substrates, arises from the three-dimensional (3D) structure of their active sites, that are sterically and electrostatically complementary to their substrates. Although the overall structures of PKs are very similar, they are also very selective, being capable of discriminating between closely related sequences that contain serine, threonine or tyrosine with different neighbor amino acids.
The catalytic subunit of protein kinase A (PKA) is the best characterized member of the large family of PKs, considering structural, biochemical, and physiological data. In this sense, PKA serves as a paradigm for the whole family [6,7]. Protein kinase A (PKA) or cAMP-dependent protein kinase (since its activity is dependent on cellular levels of cyclic AMP) is a serine/threonine protein kinase (PK) that regulates glycogen, sugar, and lipid metabolism. The selectivity of PKA is essential for cell integrity. It phosphorylates specific targets in the cell that contain a sequence pattern in the amino acid residues surrounding the phosphorylation site [8]. The optimal recognition motif for PKA is the sequence RRX(S/T)X, with a great importance of arginines at the 23 and 22 positions where position 0 is the primed phosphorylation site [9]. Kemp et al. designed the peptide Kemptide (LRRASLG) with kinetic constants comparable to native protein substrates [10]. They also found that substituting the arginine residues for other residues (alanine, lysine, or histidine), or reducing the chain length in Kemptide, negatively affect the kinetic constants [11]. In a general analysis, authors noted that these changes in Kemptide resulted in less affinities between PKA and substrates. The principal aim of this work is the use of simple computational biochemistry methods to reproduce the effects of changes in Kemptide on the kinetics of PKA-catalyzed phosphorylation. As a result we will create a predictive theoretical model and provide atomistic details about the causes of the differential affinities between PKs and their substrates.

Preparation of the Initial Structures
We used the initial structure of PKA in complex with a peptide inhibitor (PDB code 1ATP; this crystal also includes ATP and Mn 2+ atoms) [12]. This structure was modified in VMD software [13]. To construct our system, Mn 2+ atoms were changed by Mg 2+ . All missing PKA hydrogens were added, and protonation states were assigned for all ionizable residues to their default values at neutral pH. Starting from the peptide inhibitor we build the heptapeptide Kemptide (denoted as WT in this manuscript), which corresponds closely to the peptide sequence of the pig liver pyruvate kinase, a natural PKA substrate [11]. The WT sequence was obtained by mutation of the amino acids of the inhibitor peptide using Mutation plugin of VMD software package [13]. To prepare the mutants, we mutated arginines of the WT sequence inside the transformed PKA model and obtained the models that contain PKA and the peptides LARASLG (R18A), LRAASLG (R19A), LKRASLG (R18K), LRKASLG (R19K), LHRASLG (R18H), and LRHASLG (R19H). We also prepared the models including the shorter chain peptides RRASLG, RASLG, and LRRASL. For the mutants that contain histidine we prepared the models containing neutral histidine with the hydrogen atoms at position d or e (H d and H e ), and the models containing the protonated histidine H p (see the next section). All these systems (protein-peptide substrate complexes with ATP and Mg 2+ ions) were solvated in a octahedron periodic box of TIP3P [14] water molecules with a minimum solute-wall distance of 15 Å . Sodium and chloride counter ions (0.15 M NaCl) were placed on the grid to mimic physiological condition and to guarantee the electric neutrality. This procedure was made using VMD software package [13]. Thereafter, we will perform energy minimization on the models using conjugate gradient method (1000 steps) to reduce any close contacts resulting from the inclusion of new residues.
The experimental values of the variables related to the kinetics of phosphorylation for all the studied peptides are in the Table 1.

Special considerations for histidine mutants
The histidine contains the ionizable imidazole ring at the side chain. The histidine pKa value is approximately 7; therefore, both the acidic and basic forms are present at physiological pH. The acidic form (named H p in this manuscript) has the imidazolium ion, which is a resonance hybrid of two practically equivalent contributing forms ( Figure 1). The basic forms (named H d and H e in this manuscript) are produced when either of the two ring nitrogens release a proton. At pH values near 7, all three forms, are present in equilibrium, but differences in the molecular environment can shift this equilibrium.
In simulations based on molecular mechanics formalism it is not possible to consider the gain and loss of a proton. The most common practice is the assignment of a single protonation state and using it throughout the MD simulation [15]. Since our objective is the study of the effect of mutated amino acids in the affinity, we did a special treatment for histidines in the mutants R18H and R19H. We constructed three models for each mutant containing H d , H e , and H p and we accomplished MD simulations for each of them. After the analysis of our results, we suggest which are the most important protonation states for histidines at histidine mutants of Kemptide when this amino acid interacts with the residues at PKA.

Computational biochemistry protocol: Molecular Dynamic simulations and MM-GBSA calculations
MD simulations were carried out using the NAMD software package [16]. Topology and parameters of proteins included in CHARMM27 force field were used [17,18]. All the systems were subjected to 2.0 ns of equilibration, and then the production was accomplished during 10.0 ns. The equations of motion were integrated with a 1 fs time step. The van der Waals cutoff was set to 12 Å and the temperature was maintained at 298.15 K with the weak coupling algorithm [19]. Long range electrostatic forces were  taken into account by means of the Particle Mesh Ewald (PME) approach [20]. We used a computational protocol combining MD simulation and MM-GBSA (Molecular Mechanics-Generalized Born Surface Area) to study PKA-substrate systems. Figure 2 shows the free energy changes associated with the phosphorylation of Kemptide (WT) and mutants (Mt) catalyzed by PKA. As seen in Figure 2, the affinity of the PKA for the substrate is determined by the Gibbs free energy change DG ERES(i) of the reaction system from the reactants (E + S(i)) to the formation of the PKA-substrate complex ES(i), and the catalytic efficiency, that is, k cat (i)/K M (i), is determined by the Gibbs free energy change DG ERTS of the reaction system from the reactants (E + S(i)) to the corresponding rate-determining transition state TS(i).
We executed our study on the PKA-substrate complex, but not the transition state structure. This consideration makes the work much easier because standard MD approaches use classical 'molecular mechanics' force fields to simulate a stable molecular system corresponding to a local minimum on the potential energy surface, and the transition state is associated with a maximum potential energy on the potential energy surface. Our study is focused on the study of the effect on DG of the amino acids that are close to the serine that is phosphorylated (the free energy shift caused by a structural change in the sequence that contains this serine); therefore, we consider that there is no big change in the shift value when comparing formation of ES(i) with formation of TS(i). Equation 1 represents this consideration for a comparison between the WT and a mutant Mt: Considering the thermodynamic cycle in Figure 2 and the Equation (1), it is possible to obtain DDG ERTS using the MD simulation and a free energy calculation of the molecular systems that represent reactants and the PKA-substrate complex.
The DG ERES values were calculated using the MM-GBSA method [21][22][23] by the difference between the PKA, the substrate peptides, and their complexes. The molecular structures were taken from the MD simulations of the complexes. We extracted 900 snapshots from the 10 ns production MD trajectories. The explicit TIP3P water molecules and ions were removed. The free energy was calculated as follows: where DE MM(i) contains DE internal (bond, angle, and dihedral energies), DE elect (electrostatic energy), and DE VDW (van der Waals energy). DG solv(i) is the change of the solvation free energy upon binding, and includes the electrostatic solvation free energy (polar contribution calculated using generalized Born model [24], the dielectric constant of the solvent was set to 78.5), and the nonelectrostatic solvation component (nonpolar contribution estimated by solvent accessible surface area calculated with a probe radius of 1.4 Å ). MM-GBSA calculations were achieved in NAMD 2.8 [16]. To reduce computational time, the entropy term TDS (i) was not calculated since we consider that its contribution is similar for the compared molecular systems. In previous literature, many authors have been reported that the lack of the evaluation of the entropy is not critical for calculating the MM-PBSA or MM-GBSA relative free energies for similar systems [25][26][27][28][29].

Results and Discussion
Analysis of the Kemptide Arg18 and Arg19 mutants MD simulations give trajectories that contain structural data of PKA-substrate complexes. Initially, our endeavors were focused on the analysis of the effect of the mutations in the affinity. The dynamics of the interactions between the amino acids of the substrate peptides and residues at the PKA active site were studied. With this information, we give a deeper knowledge about the relevance of the molecular interactions in the complexes between PKA and its substrates.
Firstly, we analyzed the stability of the MD trajectories using the rmsd of the positions for all the backbone atoms as a function of simulation time. We observed that rmsd was almost constant after the equilibration process (2.0 ns) for all the systems (Figure 3). We identified the hydrogen bonds (HBs) formed between the substrates and PKA by measuring the donor-aceptor distances during MD simulations. We used a cutoff of 3.5 Å to consider that a donor-aceptor distance is equivalent to an HB. The numbers of HBs between the PKA and its substrates during MD simulation are shown in Figure 4. In general, the network of HB interactions is higher in the complex that contain the WT peptide, and there are less HBs when Arg18 or Arg19 are mutated. In addition, we found less intermolecular HBs when Arg19 is mutated. This result suggests that Arg19 has a higher contribution to the network of HB interactions when comparing with Arg18. Interestingly, the complex that contains the peptide R19A has one of the lower number of HBs. This could explain the lesser binding affinity of the peptide R19A for PKA, and the difference between its affinity and the affinity of the peptide R18A.
The results of the HBs occupancy during the MD simulations (Tables 2 and 3) allow a comparison of the HB interactions of the residues of Kemptide (WT) and mutants with the amino acids at the substrate PKA binding site. The analysis of these data shows that the side chain hydroxyl of the substrate Ser21 (the residue that is phosphorylated) forms stable HB interactions with the side chain carboxylate of the PKA residue Asp166, the side chain amino group of the PKA residue Lys168, and one of the oxygens of the c-PO 3 2 of the ATP. Meanwhile, the backbone NH group of the Ser21 forms an HB interaction with other of the oxygen atoms of the ATP c-PO 3 2 , and the backbone CO group forms an HB interaction with the side chain hydroxyl of the residue Ser53 in the glycine-rich loop of the PKA small lobe (Figure 5a). The backbone NH of the residue at position +1 from the phosphorylated serine (residue Leu22) has a stable HB interaction with the backbone CO of the PKA residue Gly200 in all the studied complexes.
The guanidine side chain group of the residue Arg18 of the peptide substrates (the residue at position -3 from the phosphorylated serine) has stable HB interactions with the side chain carboxylate group of the PKA residue Glu127, with the backbone CO of the residue Thr51 in the glycine-rich loop of the PKA small lobe, and with the hydroxyl at the position 3 of the D-ribofuranose ATP moiety ( Figure 5b). The analysis of the occupancies also shows that the side chain of Arg18 can form non-stable HB interactions with the side chain hydroxyl of the PKA residue Tyr330. Table 3 shows that the side chain groups of the mutated residues Lys and the protonated histidine at position 18 of the WT substrate (in mutants R18K and R18H p ) also form HB interactions with the side chain   Table 2. Hydrogen-Bond occupancies analysis from the results of MD simulations for peptide-PKA interactions for the common groups of the mutants derived from Kemptide (WT). carboxylate group of the PKA residue Glu127 (Figure 5c). These interactions are weaker for histidine e (mutant R18H e ), and even weaker for histidine d (mutant R18H d ). This interaction is not present in the mutant R18A. The guanidine side chain group of the residue Arg19 of the peptide substrates (the residue at position -2 from the phosphorylated serine) has stable HB interactions with the side chain carboxylate groups of the PKA residues Glu170, Glu203, and Glu230 (Table 2, Figure 5d). At the same time, the backbone NH of Arg19 forms an HB interaction with the side chain carboxylate of the PKA residue Glu170; this interaction is conserved when Arg19 is changed by other amino acids in all the studied mutants. Table 3 shows that the side chain groups of the mutated residues Lys and the protonated histidine at position 19 of the WT substrate (in mutants R19K and R19H p ) also form HB interactions with the side chain carboxylate group of the PKA residue Glu170, while the mutated residue histidine d at the same position 19 (in mutant R19H d ) also forms an HB interaction with the side chain carboxylate group of the PKA residue Glu203 (Figure 5e). The histidine e at this position has no stable HB interactions with PKA residues.
It is interesting to note, after the analysis of occupancies of the HB interactions, that the mutants that contain alanine at positions 18 and 19 of the substrate (R18A and R19A) do not contain groups able to form HBs with PKA; therefore, these mutants cannot form the interactions represented in the Figures 5b and d, respectively. It is also noteworthy the presence of a more stable and bigger HB network established by Arg19 with the residues at PKA, when comparing with the network established by Arg18. The substitution of the arginine at this position for lysine or histidine retains only a part of the original interactions; therefore, arginine looks optimal in the pocket that contains substrate residues at position 19, where three glumatates capture the guanidine group with a great molecular force forming a very stable framework.
The hydrophobic interactions were analyzed by accounting for stable C-C distances ,4.5 Å . We observed that the PKA residue Phe187 in the activation segment has a hydrophobic interaction with the substrate residue Ser21 and the PKA residue Pro202 has a hydrophobic interaction with the substrate residue Leu22 for all the complexes (Figure 5a). In previous results, Li et al. [30] demonstrated the positive cooperativity of ATP and the PKA residues Thr51, Glu170, and Phe187 using a series of MD simulations and free energy calculations. They also identified the interactions mentioned here for the complex between PKA and Kemptide.
The MM-GBSA free energy calculations were employed to get quantitative estimates for the binding free energies of the WT and the mutants inside the PKA active site. We calculated DDG ERES values and compared them with experimental DDG ERTS values. There is a high correlation between experimental DDG ERTS and DDG ERES values (R 2 = 0.988 including all the peptides in Table 1) that confirm that our data satisfy the Equation (1). Therefore, the calculation of DDG ERES instead of DDG ERTS is pertinent to get thermodynamic information of PKA-substrate systems. In other words, it is not necessary to describe the transition state; this is important because it is possible the use of a simple method such as MM-GBSA instead of other more complex method accounting for the quantum mechanical effects during transition state formation.
We have to consider that three classical histidines (H d , H e , and H p ) were used to model R18H and R19H mutants. To establish the correlation between experimental and theoretical DDG values we selected only one of them for each mutant; although, a thorough chemical analysis indicates that all of them should Table 3. Hydrogen-Bond occupancies analysis from the results of MD simulations for the peptide-PKA interactions between mutated amino acids. contribute to the real system. The better contribution to the correlation of one of the histidines should indicate which is the preponderant protonation state of histidine at positions 18 and 19 of the substrate when they interact with PKA active site. We evaluated the correlations for the mutants when different pairs of classical histidines are incorporated. It is noteworthy that there is a  linear correlation (R 2 .0.6 when different pairs of histidines are included, Figure 6) between the predicted DDG ERES and the experimental DDG ERTS . These correlations indicate that the parameters applied to calculate free energy changes are feasible and reliable. The best correlation with a coefficient R 2 = 0.710 was found when H e is at position 18 and H p is at position 19. The plot when mutants R18H e and R19H p are included is represented in Figure 6. Mutations try to mimic the original PKA-substrate interactions of arginines when Kemptide is present. Histidines in mutants R18H and R19H should be protonated to mimic the original arginines, but the histidine 18 does not need to be protonated according to our results. This result suggests that a protonated residue is more important at position 19 (position 22 from phosphorylated serine); therefore, the arginine at this position is more important for PKA-substrate affinity. The presence of a protonated residue at the position 18 (position 23 from phosphorylated serine) is less important for PKA-substrate affinity.
Since the affinity trend is successfully reproduced, we can analyze the free energy components such as van der Waals (VDW), electrostatic, and solvation contributions to give detailed molecular information about the studied systems. The results of the predicted DG ERES values and DE vdw , DE ele , and DG solv components for the complexes were summarized in Table 4. To get a better view on which energy terms have more impact on the change in PKA-substrate affinity, these individual energy components were carefully compared. From Table 4, it can be seen that DE vdw has the major favorable contribution to the total free energy, but there is no big difference among the DE vdw value for different PKA-substrate complexes. In this sense, the VDW energy term is not primarily responsible for differentiating the binding affinity of Kemptide and its mutants. The solvation term also has the same effect in all the PKA-substrate complexes with a favorable contribution to the global free energy. These terms neither influence the difference between the calculated DDG ERES values.
The most important term which dictates the difference in the binding affinity is DE ele . The better binding of Kemptide gains over 6.81 kcal/mol of DE ele value compared with the next favourable value for a mutant (R18H p ). In this sense, DE ele is the key factor for the more favorable free energy value for Kemptide. This is in agreement with the previous analysis of the importance of HB contributions. It is noteworthy that the DE ele value opposes affinity for the mutant models that contain the histidine residues H d , H e , and H p at position 19.

Analysis of the shorter chain peptides derived from Kemptide
MD simulations of the PKA-substrate systems including Kemptide-derived shorter chain peptides have stable rmsd values after the equilibration process (2.0 ns). After this time step, the effect of the removal of residues in PKA-substrate affinity was examined. The occupancies of the HBs formed between the substrates and PKA were identified by measuring the donoraceptor distances during MD simulations (cutoff of 3.5 Å was used to consider that a donor-aceptor distance is equivalent to an HB; Table 5). The analysis of Table 5 shows that the interactions described above between the substrate residues Arg18, Arg19, Ser21, and Leu22 and the PKA residues and ATP are conserved in shorter chain peptides when the complexes are formed.
In general, the analysis of the models allows explaining the experimental thermodynamic data, but we also accomplished MM-GBSA free energy calculations to get quantitative estimates for the binding free energies of the shorter chain peptides inside the PKA active site (Table 6). We observed that the lack of the substrate residues Leu17 (in the peptide RRASLG) and Gly23 (in the peptide LRRASL) does not lead to a lack of important PKAsubstrate interactions, which is reflected in small experimental and calculated DDG values (Tables 1 and 6). On the other hand, the lack of the substrate residue Arg18 (in the peptide RASLG) leads to the lack of the interactions represented in the Figure 5b, which has a great influence in experimental DDG ERES and DDG ERTS values ( Table 2). This effect was perceived in the theoretical calculation (calculated DDG ERES = 20.20 kcal/mol). We analyzed the VDW, electrostatic, and solvation MM-GBSA contributions for the shorter chain peptides. The results of the predicted DG ERES values and DE vdw , DE ele , and DG solv components for these complexes are summarized in Table 6. To get a better view on which energy terms have more impact on the change in PKA-substrate affinity, these individual energy components were carefully compared. As in the mutants, it can be seen that DE vdw has the major favorable contribution to the total free energy, but the contribution is lesser for the systems containing the shorter peptide RASLG. In this sense, the VDW energy term is responsible for differentiating the binding affinity of different chain length substrates. On the other hand, the solvation term has the same effect in all the PKA-substrate complexes with a favorable contribution to the global free energy.
The term DE ele has differences for the shorter chain peptides. The lack of the substrate residue Leu17 in the peptide RRASLG has no effect on DE ele , but the lack of Arg18 in the peptide RASLG harms the electrostatic interactions. Curiously, the electrostatic contribution for the complex that contains the peptide LRRASL is lower than expected (considering that this peptide contains both arginines at positions 18 and 19). We attribute this to the presence of the C-terminal carboxylate in the substrate residue Leu22; the pocket in PKA where Leu22 is located is very hydrophobic, and charged groups have non favorable interactions.
When the calculated DDG ERES values are inserted in the plot represented in the Figure 6 (open circles), the correlation coefficient R 2 = 0.714 (this is the correlation coefficient where all the PKA-substrate systems studied in this work are considered). In this sense, we can conclude that our calculations are able to reproduce the trend of the experimental DDG ERTS values, and we can propose the MM-GBSA model as a linear predictor.

Contribution of computational biochemistry to the study of PKs specificity and selectivity for substrates
Proteins play a main role on the cell functionality, and a very big part of the research today is concerned with the understanding of the interactions of such molecular components. The system biology and proteomics consider the study of the proteins and their functions, which is important in areas such as drug discovery, and understanding of the cell functionalities and the underlying causes behind diseases. Computational biochemistry methods can predict structural and functional details of molecular components. There are many reports involving PK-substrate systems, that can be studied using these methods.
High specificity and catalytic efficiency are the most important characteristics of PKs. Specificity arises from the 3D structure of the PK active site, which is able to include ATP and large substrates (other proteins), and is complementary to the transition state of the reaction. The substrate must fit perfectly within the large PK active site, which can establish many interactions with the phosphorylated residue and the nearest-neighbor residues in the sequence. Small changes in the surrounding amino acids in the sequence can dramatically affect the enzyme action. We demonstrated that our computational biochemistry protocol was able to capture the effect of the structural changes in the PK-substrate affinity for the case of the complexes between PKA and Kemptide derivatives. Despite Kemptide was reported over thirty years ago, the study of the interactions between PKs and their substrates is still a subject of study at present.
We would like to cite in this section other cases of PK-substrate systems, where small changes in the substrate sequence affect the PK catalytic action. We think that a computational protocol similar to our protocol presented here could contribute to explain the experimental findings exposed in the following.
Earlier studies show that several PKs can recognize similar phosphorylation motifs, but exhibit different steady-state kinetics for the same substrate. For instance, Thomas et al. [31] found that bovine lung cGMP-binding cGMP-specific phosphodiesterase (cG-BPDE) is a potent and relatively specific substrate for PKG (cGMP-dependent PK) as compared to PKA. Colbran et al. [32] found that the synthetic peptide RKISASEFDRPLR (BPDEtide), which has the sequence surrounding the phosphorylation site in cG-BPDE, retained the PKG/PKA specificity demonstrated by native cG-BPDE. They also studied peptide analogs of BPDEtide to determine the contribution of specific residues to PKG or PKA substrate specificity. They reported different substitutions that did not reduce PKG/PKA specificity. However, they found that a truncated BPDEtide (RKISASE) served equally well as substrate for both PKs, and a restoration of the phenylalanine, to yield RKISASEF, reproduced the original PKG/PKA specificity.
In other work, Glass and Krebs [33] evaluated the phosphorylation of histone H2B mediated by PKA and PKG. They used as substrates histone H2B, RKRS 32 RKE and RKES 36 YSV, corresponding to the amino acid sequences around Ser32 and Ser36 in histone H2B, and the undecapeptide RKRS 32 RKES 36 YSV containing both serines. They found that PKA preferentially phosphorylated the heptapeptide containing Ser36, rather than that corresponding to Ser32, mainly on the basis of a higher V max for the former substrate. Interestingly, they found that the undecapeptide containing both Ser32 and Ser36 was phosphorylated by PKA with kinetic constants that were identical with those reported for the PKA phosphorylation of Kemptide. On the other hand, based on a lower K M , the peptide containing Ser32 is a better phosphate acceptor for PKG than for PKA.
In general, the above mentioned examples show that sequence manipulation can lead to a change in specificity for different PKs. It is possible to create atomistic models of PKA and PKG forming complexes with the above mentioned peptide substrates and develop a protocol of MD and MM/GBSA to help in the elucidation of the interactions that are responsible of the results exposed. In general, the vast information about PK structures allows to create reliable models of the substrates and PKs that are not crystallized to proportionate novel ideas about the structure and function of PK-substrate complexes.
Other examples show the possible applicability of the protocol reported here in the prediction of phosphorylation sites. There are reports on the PK actions on protein substrates that cause some specific response associated with any physiological cellular change. For instance, a direct phosphorylation of Ca 2+ channel protein by PKA enhances the voltage-gated L-type calcium current by isoproterenol [34]. PKA also increases potassium current when it is stimulated in cardiac ventricular myocytes with a membranesoluble cAMP analog, indicating that PKs also can regulate K + channels, by phosphorylating some exposed amino acid residues of these proteins [35].
In a more recent report, Michard et al. [36] investigated the modulation of AKT2 (K + channel subunit encoded by the genome of the plant Arabidopsis thaliana) current by effectors of phosphatases and PKs. The authors demonstrated that transitions between AKT2 gating modes 1 and 2 involve phosphorylation/ dephosphorylation. They observed that dephosphorylation can result in AKT2 channel silencing; and phosphorylation by a PKAlike PK in the plant favors both recruitment of silenced AKT2 channels and transition from gating mode 1 to gating mode 2. They identified two PKA phosphorylation sites (Ser210 and Ser329) in the region of the pore inner mouth of AKT2. These sites are conserved in AKT2-related channels cloned from other plant species. The PK regulation of the dual gating mode of AKT2 has important consequences in plants since AKT2 is mainly expressed in the mesophyll and phloem tissues [37,38].
For the application of a computation biochemistry protocol in this example, it is necessary the construction of a 3D model of the PKA-like PK of Arabidopsis and the peptides which has the sequence surrounding the phosphorylation sites Ser210 and Ser329 in AKT2. With the protocol application, it could be possible to evaluate which AKT2 phosphorylation sites are more likely to be identified by the PKA-like PK of Arabidopsis. In addition, these models could help in the design of mutagenesis experiments to study the molecular aspects of AKT2 regulation by phosphorylation, and the possible physiological meaning of such regulation in the plant context.
Our current endeavors are oriented to the application of the computational biochemistry protocol reported here to study the exposed cases.

Conclusions
The understanding of the basis of PK-substrate specificity is essential for the identification of physiologically important substrates and also to manipulate the metabolic pathway with therapeutic purposes. In the present work, we report a computational biochemistry protocol that gives a theoretical model able to predict the differential affinities of PKA for the classic substrate Kemptide and some derivatives including mutants and short length derivatives.
In summary, the protocol consists of two steps: first, molecular structures representing the complexes in solvent media were subjected to MD simulations, and second, a binding affinity between PKA and substrate was predicted using MM-GBSA free energy calculations. A comparison between the calculated DDG values and experimental ones indicates that the model was able to describe differential affinities. As a result, the trend of the calculated affinities for the studied PKA peptide substrates was highly related to the trend of their experimental catalytic efficiencies (R 2 = 0.714). In this sense, our model demonstrates that it is possible to use theoretical atomistic models to propose specific changes in protein sequences to improve the interactions between PKs and their substrates. In addition, our theoretical model can be used to analyze the atomistic interactions that have positive or negative influence into the binding affinity of PKsubstrate complexes.
The success of this simple protocol also demonstrates that it is not necessary to investigate the activation free energies of the substrate molecules to get information about PK specificities. Those investigations require of sophisticated and complex calculations (such as QM/MM) using a higher level of theory. Instead, MM-GBSA is a popular method because of its high speed and low computational cost.