INTRODUCTION
Cardiac Steroid (CS) is one of the promising anticancer.1-8 The previous study showed that CS exhibited the most potent antitumor activity among 9,000 screened compounds.2 In the last decade, there is an increasing number of studies that reported the anticancer activities of CS.1,3-6 CS inhibited the proliferation and migration of cells, and also sensitised the multidrug resistant (MDR) strain.7 Moreover, Wei and colleagues reported that CS showed significant inhibition tumor growth in vivo below their lethal dose.8 One of the CS derivatives is bufadienolide. It has an extended six-membered lactone ring at position 17 of the steroidal scaffold. Different with the structure of cholesterol which is all in trans-conformation, the ring A/B and C/D in bufadienolides are in cis-fused conformation Figure 1. The cis-conformation was suggested as the pharmacophore of bufadienolides.4 Previously, bufadienolide was used as a treatment for heart disorder disease. However, bufadienolide has also been known to control the growth of cancer cells.9,10 One of the proposed molecular mechanisms of its anticancer activity is through the inhibition of Na+/K+-ATPase (NaK-ATPase).9-14 Supratman et al. reported the sub-micromolar inhibitory activities of three isolated bufadienolides from Indonesian plant Kalanchoe pinnata, namely C1-C3 Figure 1, on the tumor-promoted Raji cell line.15 Computational methods, such as molecular docking, have been used to predict the molecular mechanism of bioactive compounds at the atomic-level.16-19 Therefore, the aims of this work were to study the structure-activity relationship of NaK-ATPase inhibition of three bufadienolides using molecular docking method.
METHODS
Preparation of Receptor and Ligands
The crystal structure of NAK-ATPase in complex with bufalin was retrieved from Protein Data Bank (rcsb.org/pdb/) with PDB ID 4RES.14 The structure of ligands (C1, C2, and C3) were obtained from PubChem NCBI.20 Furthermore, the NaKATPase in complex with ligand was minimised to remove the sterical clashes between ligands and receptor. Since NAK-ATPase is located in the membrane, then lipid bilayer system was assembled using CHARMM-GUI membrane generator21 with 0.07 M KCL in the water explicit solvent. Then, the complete structure was converted to amber PDB format using Charmmlipid2amber.py. Hydrogen atoms were added to the protein structure using Leap module of AMBER 14.22 The AM1-BCC partial atomic charges of ligands were calculated using antechamber module of AMBER 14. Minimisation was performed using Amber force field ff14SB,23 General Amber Force Field, and Amber Lipid force field 1424 using 250 steepest descents and 750 conjugate gradient algorithms. The minimised structure was converted to PDB format using ambpdb module of AMBER 14.
Protein-Ligand Docking and Psycho-chemical Property Calculations
Docking of all ligands was performed with AutoDock 4.2.25 The structures of protein and ligand were extracted from the membrane and solvent molecules. Furthermore, they were converted into the pdbqt format by adding the atom type and charges using AutoDock Tools 1.5.6.26 The K+ parameter was added to the AutoDock Parameter library, including its ionic charge. The size of the grid box of the receptor was 56 x 56 x 56 points with a grid spacing of 0.375 Å, centered in the binding site of bufalin. The Lamarckian Genetic algorithm was used with the following parameters: 100 number of runs, five million of energy evaluations, and 250 of populations. The energy breakdown analysis was performed using summary_result4.py program from AutoDockTools Utilities. The total binding energy was break down into hydrogen bond and electrostatic energies. In addition, the ligand efficiency (total binding energy divided by a number of atoms) of each ligand was calculated by this program. Psycho-chemical property such as PSA and A Log P have been computed by Biovia Draw 2016.27
RESULT
The binding mode of bufalin in the NaK-ATPase was investigated from its crystal structure (PDB ID 4RES). It is shown that the 3-OH and 14-OH groups of bufalin formed hydrogen bonds with Glu117 and Thr797, respectively. The steroidal core of bufalin formed CH-pi interaction with Phe783, while the rest of interactions between bufalin and NaK-ATPase was hydrophobic. Furthermore, all the bufadienolides (C1, C2, and C3) were superimposed to the crystal structure of bufalin inside the cavity of NaK-ATPase. It is shown that an oxygen atom of 1,3,5-orthoacetate moieties of C1 and C2 were in short distance with the carboxylate group of Glu117 (2.72 Å). Hence, a repulsive force was expected. Moreover, the 1-OH and 5-OH groups of C3 were within the hydrogen bond distance with Glu117 (3.03 Å and 2.48 Å, respectively), while the aldehyde group (10-CHO) of C1 and C3 were predicted to form a hydrogen bond with Gln111 and Asn122. It is noted that the 14-OH group of all bufadienolides potentially formed a hydrogen bond with Thr797 Figure 2.
Protein-Ligand Docking and Psycho-chemical Property Calculations
A repulsive force between the Glu117 and the 1,3,5-orthoacetate group of C1 and C2 was predicted due to the sterical hindrance from a loop which Glu117 is located. For this reason, this loop was objected to a short minimisation scheme to avoid the unfavorable interaction with C1 and C2. After minimisation, the interaction between 1,3,5-orthoacetate group and the loop was refined Figure 3. In general, the minimised structure of receptor has better interaction with all bufadienolides than the original crystal structure. Therefore, this structure was further used to dock all bufadienolides (C1, C2, and C3). As a result, there was a good correlation between the docking binding energy and and the experimental data15 (R2 = 0.99). Table 1 reveals that the C1 has the lowest binding energy, followed by C2 and C3. Regarding the polar surface area (PSA) and log P, all compounds were predicted to have good absorptivity properties Table 1. However, it is noted that the PSA and log P values of C1 were better than the C2 and C3. This ranking might also have explained the best activity of C1 as compared to the others, due to a fact that NaKATPase is located in the lipid bilayer membrane.
Further interaction analysis using BIOVIA Discovery Studio Visualizer showed that the C1 has four hydrogen bonds, four hydrophobic and seven CH-pi interactions, while the C2 has four hydrogen bonds, three hydrophobic and seven CH-pi interactions. Lastly, the C3 has the lowest number of interactions with only three hydrogen bonds, three hydrophobic and five CH-pi interactions with the receptor. The 10-CHO group of C1 and C3 formed hydrogen bonds with Asn122 and Gln111, respectively, while the 10-CH2OH group of C2 formed a hydrogen bond with Gln111. An oxygen atom at the position 1 of the 1,3,5-orthoacetate moiety in C1 and C2 and the 1-OH group of C3 formed hydrogen bonds with Gln111. The 11-OH group of C1 and C2 formed hydrogen bonds with Asn122 and the backbone of Ile315, respectively. The 14-OH group of all bufadienolides formed hydrogen bonds with Thr797.
Figure 2
(A) Superimposed ligands in the binding cavity of NaK. An unfavorable interaction between the orthoacetate moiety and GLU117 was observed. (B) 14-OH of all bufadienolides potentially form hydrogen bond with THR797.

Moreover, the energy breakdown analysis of docking result showed that C1 has the highest ligand efficiency with strongest hydrogen bond and electrostatic interaction energies, followed by C2 and C3. Table 2.
DISCUSSION
The previous study reported that a phorbol derivateive induced the growth of tumor cell by over-stimulating the NaK-ATPase.5 Since a phorbol derivative was also used in Supratman’s experiment,15 then it is suggested that the reduction of tumor cell growth by the bufadienolides was due to specific inhibition towards NaK-ATPase. The docking result showed good correlation with the experimental value. The 1,3,5-orthoacetate moiety, 10-CHO, 11-OH, and 14-OH were suggested as the important substituents for antitumor activity of bufadienolide.
Table 1
Binding Energy and Experimental IC50 of Bufadienolides from Kalanchoe Plants.
Compound | IC50 (µM) | Molecular Docking | Lipophilicity | ||
---|---|---|---|---|---|
RMSD (Å) | Binding Energy (kcal.mol-1) | PSA (Å2) | Log P | ||
C1 | 0.3 | 0.7 | -10.7 | 111.5 | 0.9 |
C2 | 1.6 | 0.9 | -9.5 | 114.7 | 0.7 |
C3 | 3.0 | 0.9 | -8.9 | 130.4 | 0.4 |
Table 2
Energy breakdown analysis and ligand efficiency of computed ligands.
Compound | Hydrogen bond energy (kcal/mol) | Electrostatic interaction energy (kcal/mol) | Ligand efficiency (Binding energy/Σatoms) |
---|---|---|---|
C1 | -1.9 | -0.3 | -0.32 |
C2 | -1.7 | -0.1 | -0.29 |
C3 | -1.1 | -0.1 | -0.26 |
Figure 3
Crystal structure (grey) and minimisation structure(cyan) binding cavity of NaK. The residues seem to move further outside the binding cavity after minimisation process.

Figure 4
Characteristic of hydrogen bond of CHO and CH2OH at the position 10, (A) 10-CHO of C1 form hydrogen bond with ASN122, (B) 10-CH2OH of C2 form hydrogen bond with GLN111, (C) 10-CHO of C3 form an intramolecular hydrogen bond with 1-OH and hydrogen bond with GLN111.

In this study, a crystal structure of NaK-ATPase which originated to genus Sus scrofa (PDB ID 4RES), with 98% identity with human NaK-ATPase, was used as a protein target. The molecular mechanics minimisation, besides removing the sterical hindrance from Glu117 to a 1,3,5-orthoacetate moiety, also relaxing the residues in the binding cavity and forming better interaction with bufadienolides.
The 1,3,5-orthoacetate moiety is a distinctive feature of C1 and C2 as compared to the structure of C3 and the other common bufadienolides. In C3, substituents at position 1 and 5 are the OH-groups, while at position 3 is an acetate group. Mijatovic et al. suggested that the substituents in positions 1, 3, and 5 were not very critical to the inhibition activity.1 However, when we introduced a 1,3,5-orthoacetate moiety to C3 and docked it to the NaK-ATPase, the binding affinity was improved.
Mijatovic1 also suggested that the polar contact of the 10-CHO group with the Gln111 and Asn122 was not essential to the optimum interaction with NaK-ATPase. However, Larsen28 showed that the 10-CHO group formed a hydrogen bond with Gln111 and Asn122. From our results, we agreed that the presence of 10-CHO group and the 1,3,5-orthoacetate moiety lowered the binding energy by forming a hydrogen bond with Asn122. It is noted that the 10-CHO group in C3 and 10-CH2OH group in C2 were not formed a hydrogen bond with Asn122, thus decreasing the inhibitory effect. The importance of asparagine at position 122 was tested by Canfield30 and Lingrel.31 It is shown that the mutations of asparagine to glutamic acid at the position 122 (N122E) in the NaK-ATPase decreased the inhibition activity of CS (ouabain). Therefore, this agrees with our proposal of the importance of hydrogen bond with Asn122, despite the simplified approach used in docking method.
The role of 11-OH group to the binding energy was observed in this study. We purposely removed the 11-OH from C1 compound Figure 5, which decreased the binding affinity by losing a hydrogen bond with Asn122. In contrast, Kamano et al. suggested that the 11-OH was not an important substituent for ligand binding.29 Interestingly, Supratman et al. found that the 11-OH group decreased the cytotoxicity of bufadienolides.15 Thus, it is indicated that the formation of hydrogen bonds with Ile315 (in C2) and Asn122 (in C1) would increase the ligand’s selectivity to NaK-ATPase. Therefore, a hydrogen bond with Asn122 might be one of the reasons for the better activity of C1 as compared to the C2. Furthermore, the energy breakdown analysis showed that C1 has the highest ligand efficiency among all, which means that the individual interaction in C1 was stronger than C2. Despite the similar number of hydrogen bond formed by C1 and C2, the 1,3,5-orthoacetate moiety and 10-CH2OH in C2 shared the hydrogen bonds with Gln111 Figure 4. This sharing interaction would result in a weaker hydrogen bond,32 and hence the total energy of hydrogen bond in C2 was less than that of C1.
Figure 5
The effect of 11-OH to the docking binding energy of C1 (red color) and C1mod (yellow color, the 11-OH is removed). Hydrogen bond is depicted in green-dashed line. C1mod loses a hydrogen bond with ASN122.

Regarding the physico-chemical properties of C1-C3, it is shown that their PSA and Log P values were correlated with their binding affinity, i.e., hydrophobicity increased the activity. As shown in Table 1, all of the Log P and PSA values were less than 5 and 140, respectively, indicating a good membrane permeability property.33,34 This result suggested that bufadienolides compounds, especially C1, would be able to enter the hydrophobic binding cavity of NaK-ATPase.
CONCLUSION
In summary, we have performed molecular docking to study the structure-activity relationship of bufadienolides derivatives from K. pinnata to the inhibition of NaK-ATPase. The crystal structure of NaK-ATPase from Sus scrofa bound to bufalin was used as protein target and as a control positive. The docking result showed a good agreement between the calculated binding energy and the experimental data. The presence of 1,3,5-orthoacetate moiety, 10-CHO, 11-OH, and 14-OH increased the activity of C1 than the other bufadienolides. Although the role of 11-OH group was still debatable, this substituent appeared to decrease the toxicity of the bufadienolide and hence better selectivity. The result is also supported by the psycho-chemical properties value (PSA and LogP). This study suggests that C1-C3 possess good selectivity to NaK-ATPase and low toxicity, thus interesting to be further studied for anticancer activity.