Molecular docking and density functional theory studies on creatine,guanidinoacetic acid,and their phosphorylated analogues binding to muscle creatine kinase

Abstract

Comparative molecular docking studies on creatine and guanidinoacetic acid, as well as their phosphorylated analogues, creatine phosphate, and phosphorylated guanidinoacetic acid, are investigated. Docking and density functional theory studies are carried out for muscle creatine kinase. The changes in the geometries of the ligands before and after binding to the enzyme are investigated to explain the better binding of guanidinoacetic acid and phosphorylated guanidinoacetic acid compared to creatine and creatine phosphate.

Introduction

Creatine kinase (CK) is a crucial enzyme of cellular bioenergetics that regulates energy fluxes by the transfer of phosphates between different phosphorylated substrates and adenine nucleotides. ¹ Phosphocreatine (CrP) seems to be the primary CK substrate with the highest affinity for CK reported in both in vivo and in vitro studies.^2,3 However, many other guanidino compounds can act as substrate analogues for the enzyme, which can be particularly important in cases of CrP deficiency. Kan et al. ⁴ reported that phosphorylated guanidinoacetic acid (GAAP) can partly compensate for the lack of CrP, with GAAP used for high-energy phosphoryl transfer through reaction with CK in the skeletal muscle of mice. However, limited data are available concerning the comparison between GAAP and CrP in terms of CK kinetics. So far, there are no studies evaluating the docking environment of CK as well as differences between specific substrates (such as creatine and guanidinoacetic acid) and their phosphorylated analogues before and after docking. Molecular docking techniques are well-documented computational tools that suggest the most acceptable mechanism for intermolecular interactions of ligands with biomolecules. ⁵ Therefore, molecular docking studies on the above compounds are expected to provide valuable information about the mechanism of binding to CK and are highly justified.

Results and discussion

Density functional theory analysis

A molecule of phosphocreatine (CrP) forms a seven-membered ring (Figure 1(a)) through intramolecular hydrogen bonding between the H5-atom of the iminium group and the oxygen atom from the carboxylate group (length 1.70 Å). The angles between bonds from the guanidine groups C₂–N₃–C₄ (121.1°), N₃–C₄–N₆ (121.1°), and C₄–N₆–H₅ (113.8°) are close to the values obtained from X-ray analysis of guanidinium chloride ⁶ and guanidinoacetic acid ⁷ (121.5°, 120.8°, and 114°, respectively). Besides, the angle between the N₃–C₂–C₁ bonds (109.6°) is almost tetrahedral, which indicates that the formation of H-bonds does not disturb the geometry of the CrP molecule and contributes to its stability. However, steric hindrance between the methyl group on N₃ and the phosphate group restricts free rotation about the C₄–N₅ bond compared to the molecule of N-phosphoguanidinoacetic acid (GAAP). The values of the rotation barriers depend on the dihedral angles between the N₃–C₄–N₅–P bonds and are presented in Figure S1. As can be seen from Figure S1, the amount of energy needed for rotation around the C₄–N₅ bond in CrP is extremely high, even larger than the dissociation energy of a covalent bond. This result indicates that rotation in the range of 30°–70° will cause breaking of the C₄–N₅ bond, causing a restriction of free rotation. The existence of a methyl group on N₃ in CrP and its rotation in the stated range leads to overlapping of the Van der Waals volumes of the CH₃ group and the NH₂ group, as shown in Figure S2. The existence of the methyl group and steric hindrance along with overlapping of the Van der Waals volumes leads to a higher amount of energy being required for rotation in comparison with GAAP in which there is no methyl group in this position.

Figure 1.

Optimized structures of molecules of (a) CrP and (b) GAAP.

The molecule of GAAP forms six- and seven-membered rings through two intramolecular H-bonds (Figure 1(b)). The seven-membered ring is similar that of the CrP, but the H-bond in GAAP is shorter (1.66 Å) compared to that of CrP, and has large deviations of the bond angles from the preferred tetrahedral geometry. The second H-bond (length 1.58 Å), between H₃ and O₃ from the phosphate group, leads to the formation of a six-membered ring (–O₃–H₃–N₃–C₄–N₅–P–), which results in the almost planar geometry of this part of the molecule. From Table 1, it can be seen that the values of some of the dihedral angles are close to zero, indicating a high ring strain. Considering these results, the formation of this H-bond creates a rigid structure with restricted ability for free rotation and with relatively unfavorable geometry in the six-membered ring.

Table 1.

Dihedral angles before and after docking of CrP and GAAP.

Angle	Before	After
	CrP
C7-N3–C4-N6	147°	134.2°
C7-N3–C4-N5	−32.3°	−31.3°
N6-C4–N5-P	153.9°	148.4°
N3-C4–N5-P	−26.7°	−31.7°
	GAAP
H3-N4–C4-N5	0.7°	−66.3°
C4-N5–P-O3	−5.7°	72.0°
N3-C4–N5-P	4°	63.4°
C2-N3–C4-N6	14.9°	−5.5°
N6-C4–N5-P	−175.8°	−41.1°
C2-N3–C4-N5	−164.9°	−34.9°

Crp: creatine phosphate; GAAP: phosphorylated guanidinoacetic acid.

HOMO–LUMO orbitals play important roles in the stabilizing interactions between ligands and receptor proteins. Hence, the orbital energies of both the HOMO and the LUMOs were calculated to estimate the chemical reactivity of the selected compounds using density functional theory (DFT). HOMO–LUMO plots were generated to analyze the atomic contribution of these orbitals. The HOMO and LUMO plots show the positive electron density in red color and the negative electron density in blue (Figure 2). The HOMO and LUMO energy gap (ΔE = E_LUMO – E_HOMO) explains the eventual charge transfer interactions taking place within the molecule. As can be seen from Figure 2, a lower energy gap was exhibited by GAAP and CrP compared to non-phosphorylated analogues. The higher values of HOMO–LUMO gap indicate the higher kinetic stability of the molecules. On the phosphorylation of creatine and guanidinoacetic acid, the values of ΔE are lowered, indicating less kinetic stability of the phosphorylated analogues. Considering the chemical hardness, GAA and Cre molecules are harder than GAAP and CrP.

Figure 2.

Spatial distribution of the HOMO and LUMO molecular orbitals and the corresponding ΔE values of the investigated compounds. Ppositive electron density is shown in red, while negative electron density is shown in blue.

The molecular electrostatic potential (MEP) surfaces of the investigated ligands are shown in Figure 3. The color code is used to better visualize the density of the electrons, where blue indicates a minimal concentration of electrons and red indicates a high concentration of electrons. The regions with negative potential values were concentrated on oxygen and phosphorus, whereas regions with positive potentials were concentrated on the hydrogen atoms of NH₂, NH, CH₂, and CH₃ groups. A region with the highest positive potential is located over the hydrogens from the guanidino group, and this site of the molecule has an essential role in hydrogen-bonding interactions with the negative sites of CK. However, from Figure 3, it is evident that introduction of a phosphate group leads to a more pronounced negative potential, causing stronger H-bond donor ability in the case of GAAP and CrP compared to GAA and Cre.

Figure 3.

Molecular electrostatic potential (MEP) surfaces: (a) GAA, (b) Cr, (c) GAAP, and (d) CrP.

Molecular docking

Docking results for the enzyme CK with GAAP, CrP, GAA and Cre are shown in Figure 4, and the most significant interactions are tabulated in Table 2. As can be seen from Figure 4, the enzyme CK has strong interactions with the oxygen atoms of the carboxylate groups of GAAP and CrP, through hydrogen bonding with arginine side chains (Arg 96, Arg 320, Arg 341), glutamine (Gln 318), and serine (Ser 285), which are in a good agreement with the results of Bong et al., ⁸ Lim et al., ⁹ and Jourden et al. ¹⁰ These H-bond interactions decrease the electron density on the carboxylate anion, which result in breaking of intramolecular H-bonds between the oxygen atom of carboxylate group and a hydrogen atom from N₆ and allows for rotation around C₂–N₃ for 130° in the case of both molecules.

Figure 4.

Docking results between the enyzme creatine kinase and different ligands: (a) GAA, (b) Cre, (c) GAAP, and (d) CrP.

Table 2.

Most significant interactions between ligands and CK, along with lengths and type of interaction.

Interaction	Type	Bond length, r (Å)	Ebond (kJ∙mol−1)	Interaction	Type	Bond length, r (Å)	Ebond (kJ∙mol−1)
	CrP				GAAP
Ser285–O5	H-bond	1.728	25.00	Ser285–O4	H-bond	1.742	23.85
Gln318–O4	H-bond	1.729	25.00	Gln318-O5	H-bond	1.729	25.00
Glu232–H5	H-bond	1.704	23.45	Glu232–H6	electrostatic	4.683	–
Asn286–O2	H-bond	1.748	25.00	Glu232–H3	H-bond	1.939	24.43
Arg320–O3	electrostatic	4.425	–	Arg320–O3	Electrostatic	3.671	–
Arg320–O2	H-bond	1.748	25.00	Arg320–O2	Electrostatic	4.520	–
Arg320–O4	Electrostatic	4.875	–	Arg320–O5	H-bond	1.487	22.88
Arg96–O5	H-bond	1.63	25.00	Arg96–O4	H-bond	1.737	25.00
Arg132–O2	H-bond	1.889	25.00	Arg96–O5	Electrostatic	4.106	–
Arg132–O3	Electrostatic	3.631	–	Arg132–O2	Electrostatic	4.967	–
Arg341–O5	H-bond	1.752	11.33	Arg132–O5	Electrostatic	4.444	–
Arg341–O4	Electrostatic	4.064	–	Arg341–O5	H-bond	1.56	25.00
Arg341–O2	Electrostatic	4.425	–
	Cre				GAA
Glu232–N6	H-bond	1.803	24.43	Glu232–N6	H-bond	1.837	19.23
Glu232–N5	H-bond	1.470	25.00	Glu232–N5	H-bond	1.875	17.78
Arg320–O5	Electrostatic	2.992	–	Glu232–N3	Electrostatic	4.544	–
Arg96–O5	Electrostatic	3.420	–	Arg320–O5	Electrostatic	3.579	–
Arg96–O4	H-bond	1.562	25.00	Arg96–O5	H-bond	1.784	21.10
Arg341–O5	H-bond	1.541	25.00	Arg341–O5	H-bond	2.145	21.55
Gln318–O4	H-bond	1.716	25.00	Arg341–O4	H-bond	1.841	25.00
				Ser285–O4	H-bond	1.578	25.00
				Asn286–O4	H-bond	1.998	23.28

Crp: creatine phosphate; GAAP: phosphorylated guanidinoacetic acid.

However, the most important difference between CrP and GAAP from an aspect of geometry and availability to bind with CK is a methyl group on N₃ of the CrP molecule. As can be seen from Table 1, in the case of CrP, there are no significant changes in dihedral angles, except around the C₂–N₃ bond before and after binding with CK. However, the GAAP molecule has significant conformational changes after docking. Breaking the intramolecular H-bond between O₃ and H₃ in the GAAP molecule allows rotation around the N₃–C₄, C₄–N₅, and N₅–P bonds (changes in the dihedral angles are given in Table 1). As a consequence of these geometry changes, a new H-bond between the O₃ atom from phosphate group and H₅ from the iminium group in GAAP (length 1.74 Å) is formed. These conformational changes lead to the formation of a new six-membered ring between H₅–N₆–C₄–N₅–P–O₃, which is not planar and has a more stable conformation compared to GAAP before binding to CK. The formation of this new intramolecular H-bond in GAAP decreases the electron density on the oxygen atom from the phosphate group, resulting in a weaker interaction of this group with CK compared to a CrP molecule. Oxygen atoms from the $-{\mathrm{PO}}_3^{2-}$ group in CrP form three H-bonds with the $-{\mathrm{NH}}_2^{+}$ group of amino acid side chains (Arg 320, Arg 132, Asn 286), and strong electrostatic interactions with the ${\mathrm{NH}}_2^{+}$ group of amino acids (Arg 320, Arg 132, Arg 341). However, GAAP forms only electrostatic interactions between oxygen atoms (O₂ and O₃) and the $-{\mathrm{NH}}_2^{+}$ group of the arginine side chains (Arg 132, Arg 320).

It is a well-known fact that magnesium ions play a crucial role in the formation of the transition state of the phosphoryl transfer reaction. Thus, in our case, the electrostatic interactions between magnesium ions and phosphate groups from CrP and GAAP were confirmed. Furthermore, the side chain of glutamic acid (Glu 232) had a significant role in binding CrP and GAAP molecules. A CrP molecule forms a H-bond between the H₅ atom and the side-chain carboxyl group of Glu 232, while GAAP forms a H-bond through the H₃ atom and a electrostatic interaction through H₆ atom with the previously mentioned carboxyl group of Glu 232. These observations are in excellent agreement with the literature data, according to which the side-chain carboxyl group of Glu 232 plays a crucial role in orientating the substrate in the correct location for catalysis. ⁸ From these results, it can be concluded that the conformational geometry of docked GAAP changed more than the geometry of docked CrP, leading to relaxation in the molecular structure of GAAP. This relaxation is one of the main reasons for the more negative value of the binding energy of a GAAP molecule compared to CrP (Table 3).

Table 3.

Binding energies of selected ligands with the enzyme creatine kinase.

Ebind (kJ·mol−1)	Cre	GAA	CrP	GAAP
	−126	−244	−151	−279

Crp: creatine phosphate; GAAP: phosphorylated guanidinoacetic acid.

A molecule of ADP, in the presence of GAAP and CrP, binds almost identically in both cases for CK. The NH₂ group from adenine forms two H-bonds with an oxygen atom from the carbonyl group of serine (Ser 128) and glycine (Gly 294). Atom N₁ from the adenine ring also forms one H-bond with the OH group of the serine (Ser 128) side chain. Histidine (His 296) plays the most significant role in ADP binding to CK, forming a π–π interaction with the adenine ring and a hydrogen bond with the 2′-OH group of the ribose ring. These results agree with the results of fluorescence quenching measurements on ADP docking on CK performed by Borders et al., ¹¹ and with X-ray crystallography measurements of docked ADP on CK by Lahiri et al. ¹² and Bong et al. ⁸ Interactions between CK and oxygen atoms from phosphate groups occur via H-bonds and electrostatic interactions with positively charged side chains of arginine (Arg 130, Arg 132, Arg 292, Arg 320), as well as via H-bonds between the oxygen atom on the α-phosphate of ADP and the H-atom from the peptide bond between Val 325 and Gly 324. The crucial roles for optimal positioning of the molecule before the phosphoryl transfer reaction ¹⁰ have utilized ${\mathrm{Mg}}^{2+}$ ion and Arg 320, which simultaneously interact with the oxygen atoms from the ADP, GAAP, and CrP phosphate groups. After transfer of the phosphate group from CrP or GAAP to ADP, there are no significant conformational changes in formed Cre or GAA molecules. The key role for stabilization of the guanidinium group of these molecules involves the carboxyl group of the glutamic acid side chain (Glu 232). This carboxyl group forms two H-bonds with a molecule of creatine, while with the guanidinoacetic acid molecule, forms two H-bonds and one electrostatic interaction.

Stabilization of the carboxyl group of Cre, after transfer of the phosphate group, is less efficient compared to CrP. The Cre molecule forms five non-covalent interactions (three H-bonds and two electrostatic interactions), representing two less than CrP (five H-bonds and two electrostatic, Table 2). The difference in the efficiency of the stabilization of the carboxyl group is notably smaller in the case of GAA/GAAP. A molecule GAA has only one electrostatic interaction less with the side-chains of CK amino acids, compared to CrP. The total impact of all the interactions on the molecules of Cre and GAA and its phosphorylated forms is displayed via the electrostatic binding energies in Table 3, which show that the phosphorylated forms have stronger interactions with CK, and that the molecules of GAA and GAAP also are more strongly bond to the enzyme compared with Cre and CrP. Moreover, from the change of the atomic charges of the substrate before and after docking (Figure 5), it can be concluded that Cre and GAA have more significant changes in the electron structure after docking compared to CrP and GAAP. The existence of a phosphoryl group reduces the electron acceptor potential of the guanidino group, thereby decreasing the change in its electron structure.

Figure 5.

Atomic charges before (left) and after (right) docking for (a) GAA, (b) Cre, (c) GAAP, and (d) CrP.

Conclusion

In this work, molecular docking of GAA, Cre, GAAP, and CrP toward muscle CK was performed, supported by DFT calculations, to better understand the binding affinity of these molecules. From the results of molecular docking, it is noted that optimal positioning of the molecule before the phosphoryl transfer reaction involves ${\mathrm{Mg}}^{2+}$ ions and Arg 320, which simultaneously interact with oxygen atoms from the phosphate group of ADP and the phosphate groups of GAAP or CrP. However, the crucial role for stabilization of the guanidinium group of GAA and Cre involves a side-chain carboxyl group of Glu 232. Analyzing the conformational changes of the ligands, it was concluded that more significant changes in geometry before and after binding to CK occurred in the case of GAAP, leading to relaxation of the geometry. As a consequence, the GAAP molecule has lower binding energy compared to CrP. The obtained results from computational simulations are in good agreement with literature results on Cre and CrP binding to CK.

Computational methods

The initial three-dimensional structure of muscle CK was retrieved from the protein data bank (PDB), its identifier being 1I0E. The raw CK structure was prepared using the Protein Preparation Wizard and the Glide module in the Schrödinger suite 2015-2 package. Desolvation was carried out by deleting the crystallized free water molecules beyond 5 Å. The enzyme structure integrity was checked and adjusted and missing residues and loops were added using Prime. The protonation, orientation, and tautomeric states of Asp, Glu, Arg, Lys, and His were adjusted to match a pH of 7.4 using the PROPKA option. The network of hydrogen atoms in the CK model was optimized at pH 7.4 using the optimized potential for liquid simulations (OPLS) force field (2005). The active sites of the modeled CK were investigated using SiteMap program implemented in Schrödinger Suite Package. The active pocket in CK was fixed with the help of the receptor generation program from the Glide module and used as such in the docking studies. This is a grid-based method with energetics that gives scores based on favorable interactions between a ligand and a protein. The binding sites were computed with a grid resolution of 0.3 Å. A Van der Waals scaling factor of 0.8 and partial charge cutoff of 0.15 were used in the docking process for all ligands.

The prepared ligands (Cre, GAA, GAAP, CrP, Mg²⁺, ADP, and ATP) were subjected to docking into the grid enclosed by active site of the enzyme, with flexible ligand sampling and using the extra precision mode. The grid size was set to 60 Å along all axes with a grid spacing of 0.5 Å, so the grid box allowed sufficient room for both translational and rotational freedom of the ligands. The docking parameters were as follows: population size 100, maximum number of iterations 2000, and the maximum number of steps 750. Interactions between the enzyme and the ligands were calculated based on the quality of the geometric contacts and their energies. Each ligand structure was geometrically optimized employing the empirically dispersion-corrected B3LYP exchange-correlation functional (B3LYP-D3) with the 6-31+G(d,p) basis set using Macro Model/Conformational Search and prepared for further docking using LigPrep with the force field OPLS-2005. All models were subjected to post-docking minimization and the best-docked structures for each ligand were determined based on the model energy score, which combines the energy grid score, the binding affinity, the internal strain energy, and the Coulomb–Van der Waals interaction energy scores. The models with the lowest energy were visualized using Maestro Pose Viewer (Schrödinger suite). The lowest binding energy conformer was chosen among 10 different conformers for each docking simulation and the one with the best score was used for further analysis after the post-docking minimization process had been finished. To estimate the relative binding affinity of the ligands and for comparison with experimental binding affinities, the free energies were computed via the Prime/Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) approach. MM-GBSA considers contributions from OPLS molecular mechanics energies (E_MM), polar solvation through the surface-generalized Born solvation model (GSGB), and a nonpolar solvation term composed of the nonpolar solvent-accessible surface area and Van der Waals interactions. The methodology adopts the surface-generalized Born model using a Gaussian surface for better representation of the solvent-accessible surface area. Residues in binding pockets of the protein were treated as flexible and ligand partial charges were assigned from initial charges. For a more reliable representation of the solvation effect as a dielectric medium, the generalized Born (GB) model was used.

The electronic effects of the bioactive molecules play an important role in the pharmacological effects. Subsequently, surfaces (molecular orbital, density, potential) and atomic electrostatic potential charges were monitored to compute the HOMO and LUMO energies. The HOMO energy represents the region of the small molecules which can donate electrons during binding, while the LUMO energy signifies the capacity of the molecule to accept electrons from the receptor. The difference in the HOMO and LUMO energies, known as the HOMO–LUMO gap energy, indicates the electronic excitation energy that is necessary to compute the molecular reactivity and stability of the compounds.

Supplemental Material

sj-pdf-1-chl-10.1177_1747519820978583 – Supplemental material for Molecular docking and density functional theory studies on creatine, guanidinoacetic acid, and their phosphorylated analogues binding to muscle creatine kinase

Supplemental material, sj-pdf-1-chl-10.1177_1747519820978583 for Molecular docking and density functional theory studies on creatine, guanidinoacetic acid, and their phosphorylated analogues binding to muscle creatine kinase by M Vraneš, S Ostojić, Č Podlipnik and A Tot in Journal of Chemical Research