Palbociclib (PD 0332991) Interaction with Kinases. Theoretical and Comparative Molecular Docking Study

The optimized geometry of palbociclib, (PD 0332991) (8‐cyclopentyl‐6‐ethanoyl‐5‐methyl‐2‐(5‐(piperazin‐1‐yl)pyridin‐2‐ylamino)pyrido[2,3‐d]pyrimidin‐7(8H)‐one), electrostatic potential map, molecular orbitals were calculated using the density functional theory. The geometry was used in a molecular docking study of palbociclib‐kinase complexes, results could be explained by the charge of the nitrogen and oxygen atoms within the palbociclib. Energy gap of HOMO‐LUMO surfaces, could help to explain the reactivity of the ligand and the hydrogen bonding with three different kinases, two of CDK6 and one of CDK4 type. Docking results are similar and complementary with literature reports using molecular dynamics, were hydrogen bonding was obtained and analyzed. The promiscuity of three kinases with palbociclib was detected by the docking results, thus, palbociclib could be used in other types of cancer besides myeloid leukemia. Some similarities are found with CDK4/CDK6 kinases which allow us to determine that palbociclib could be used to control other resistant inhibitor types of cancer.


Introduction
Cancer is today one of the leading causes of death worldwide, with around 10.2 million cancer-related deaths in 2017 and new cases expected to rise up to 25 million in the following 2 decades. [1] In spite of the important progress in cancer therapy, controlling and effectively fighting this disease still remains a major challenge. Indeed, more effective therapeutic drugs are required to increase the number of drugs at disposal worldwide. Cohen (2002) has defined protein kinases as the 'targets of the twenty-first century' due to their role in most aspects of cell function and aberrant growth. [2] In fact, Kinases represent one of the largest protein families in the entire human genome, [3] for they are involved in intracellular signaling processes, catalyzing the transfer of the [7 -9] Nevertheless, despite the huge investments and the large arsenal of drugs, fighting cancer with kinase inhibitors is still severely hampered since most tumors can escape inhibition by a single kinase. [10] Protein kinases are enzymes that phosphorylate other proteins (typical) and other times different molecules (atypical) and thus play a vital role in signal transduction and cell differentiation. [11,12] Abundant clinical evidence from anticancer studies suggested that many cancers involve mutations harbored in the kinase, thus leading to signal disruption. [13,14] Similar is true for a variety of different disease states, and therefore, the design of safe and efficient kinase inhibitors is a promising area of research in the pharmaceutical industry, ever since the first kinasetargeted drug, imatinib, was designed, tested and approved in 2001. [15] As for today, 30 kinase inhibitors have been approved by the U.S. Food and Drug Administration (Nov 2015,FDA). [16 -18] Notwithstanding this success, developing a kinase-targeted inhibitor with a desired selectivity profile across the human kinome remains a challenge, for all kinases share a common catalytic domain and folding scaffold that binds ATP. [19] Thus, it is relevant to clarify the basic principles governing the kinase-ligand binding specificity across human kinome, if rational design of safe and efficient kinase-targeted (poly)-pharmacological therapeutics. [20] is to be achieved. With the availability of an increasing number of diverse kinase-ligand complexes structures released by the Protein Data Bank (PDB), [21] more efforts have been devoted to establishing subtle binding modes by utilizing solved kinase-ligand complex structures. Kinase-ligand interaction patterns have been analyzed by directly comparing the binding similarity and the structural-activity relationship of kinase-ligands. [22][23][24][25][26][27] These studies have allowed rational drug design and have clarified the inhibition mechanisms for certain drug molecules. [28,29] However, most studies are limited to the report of the binding mode of one kinase or a family of closely related kinases when a specific kinase is the focused target for drug design and discovery. Accordingly, the initial purpose of the present report is to characterize the ligand palbociclib by DFT, to obtain the optimized molecular geometry. Other aim is to elucidate the interaction of palbociclib optimized structure with CDK receptors and thus understand the inhibition mechanism of the ligand. Recently, molecular docking studies have been extensively used to understand these mechanisms and other techniques, such as molecular dynamics, are used as complements. Other aim is to compare the molecular docking results of the ligand with three different kinases, one of which was reported in the literature in a molecular dynamics study of palbociclib, and to explain the inhibitory action of the ligand.
Anticancer drug discovery has been strongly focused on the development of drugs focused to act against a specific target with high potency and selectivity. Clinical experience, including the discoveries of drug resistance in cancer chemotherapy, has revealed that single targeting might not always produce the desired biological effect, even if the target is inactivated or inhibited. [30][31][32] The reason of this is related to the development of resistance either by self-modification of the target through mutation or by the adoption of new pathways by a cancer cell, for the growth and multiplication. However, the approach of identifying and targeting a single oncoprotein has not produced a successful treatment and may not be sufficient as to achieve durable remission in patients. [33] Therefore, modulation of the biological network is known to be beneficial. Nowadays, there are two contrasting strategies to design the multi-targeting therapeutics. Combination drug therapy is the first strategy by creating an additive or synergistic effect of multiple drugs acting on separate targets. There are many successful treatments with the combination therapies, for example, preclinical evidence of increased apoptosis and delayed resistance to serine/ threonine-protein kinase B-Raf [34,35] has convinced the FDA to approve the combination of dabrafenib (BRAF inhibitor) plus trametinib (MEK inhibitor) for the treatment of metastatic melanoma with BRAF mutations. [36,37] Combination therapy with both RAF inhibitor (vemurafenib) and MEK inhibitor (cobimetinib) has been found promising in phase III clinical trials against BRAF mutated melanoma. [38] Other example of successful combined therapy is the use of palbociclib and letrozole in the treatment of advanced breast cancer. [39] Pfizer has recently introduced palbociclib (PD 0332991), a small-molecule CDK4/6 inhibitor, to inhibit cancer growth, specifically breast cancer. [40] Palbociclib preferentially inhibits proliferation of luminal estrogen receptor-positive human breast cancer cells in vitro. [41] According to Reuters, palbociclib could generate annual sales in excess of $5 billion in the next 10 years. [42] Thus, the ultimate goal of this present article is to study optimized geometry, molecular orbitals and the map of electrostatic potential by DFT studies, as well as molecular docking to demonstrate palbociclib's potential as a polytargeted inhibitor, using the values of the neutral molecule of palbociclib in the binding site of the kinases CDK4 and CDK6. We also analyze the hydrogen bonding, which is important evidence of the ligand-kinase interaction.

Results and Discussion
Quantum Chemistry

Structural Parameters
Gaussian, is an electronic structure program based the quantum mechanics. Density functional theory (DFT) was used to prepare the optimized palbociclib geometry. DFT predicts the potential energies, molecular structures; geometry optimization, and reactivity through the gap between LUMO (Lowest Unoccupied Molecular Orbital) and HOMO (Highest Occupied Molecular Orbital). The optimized geometry of palbociclib was selected as the lower energy conformer out of 43 conformers ( Figure 1). Distances, angles and dihedrals of the optimized palbociclib structure are summarized in Table 1.

Molecular Orbitals
The frontier molecular orbital energies HOMO and LUMO, are important for qualitatively predicting a chemical species' electronic properties and reactivity. The B3LYP/6-311 + + G(d,p) method was used to calculate this theoretically. HOMO is a parameter that is correlated with the electron donating ability of the molecule. High values of HOMO indicate the tendency of the molecule to donate electrons to the acceptor molecule with an empty molecular orbital and low energy, LUMO ( Figure 2). [44] The energy gap (ΔE = LUMO-HOMO) is a function of the reactivity of palbociclib towards its inhibitory cytotoxic activity. As ΔE decreases, the reactivity of the molecule increases, thus leading to an increase in the electron donating efficiency of the molecule. [45] Lower values of the energy difference will render good inhibitory efficiency because the energy to remove an electron from the last occupied orbital will be lower.
Using the values from LUMO and HOMO, the energy gap of palbociclib is ΔE = (À 2.28 eV)-(À 5.57 eV) = 2.77 eV = 116.797 kcal/mol. HOMO and LUMO analysis were used to determine charge transfer within the molecule Palbociclib, which has a low energy gap of 2.77 eV, which is associated with the ligand's high chemical reactivity and low kinetic. Experimental validation of our predictions is certainly encouraged, given the particular significance of kinase inhibitor targets and their inhibitory properties. The value of the energy gap can be used to obtain information by comparing it with that of similar compounds when the theoretical calculations were carried out at the same level of theory. The lower values of the energy difference will show more reactivity because the energy required to remove an electron from the last occupied orbital will be low. [46] Figure 1. Optimized structure and the numbered atoms of palbociclib. Oxygen in red, nitrogen in blue and carbon in grey. Other authors [43] reported a conformer with a higher value of energy (less stable), to a similar level of theory.

Electrostatic Potential
The electrostatic potential is a representation of an electric charge distribution. [47] A positive electric potential means that a positive charge will be repelled in that section of the space. A negative electric potential means that a positive charge will be attracted. The map of electrostatic potential (MEP) of the ground state is shown in Figure 3. The colors are representative of the values of the MEP energy at the electron density surface. Red color indicates a higher electron density around the oxygen of one keto group representing the most negative regions of the electrostatic potential for a positive test charge.
The present work reveals that the best conformation energy value of palbociclib is at À 9.2 × 10 5 kcal/ mol, which is the minimum energy value obtained by using Gaussian 09 software. Furthermore, the charge  Chem. Biodiversity 2023, 20, e202200554 values of O2, N2 and N4 atoms are À 0.563, À 0.606 and À 0.566, respectively. They are responsible for the various interactions between the ligand and the proteins, particularly the observed hydrogen bonding between the palbociclib and the kinases. Quantum chemical parameters such as HOMO, LUMO, energy gap (ΔE) and maps of the electrostatic potential of palbociclib have been obtained to gain a deeper insight into their molecular properties.

Docking
The docking analysis of optimized palbociclib and human cyclin dependent kinase receptor was carried by Molegro docking software. Molecular docking study allows virtually screening a database compounds and predicting the strongest binders based on their scoring functions. It explores ways in which molecules, such as drugs and a kinase receptor, fit together and dock well with each other. The molecules binding to a receptor, inhibit its function, and thus act as a drug. [48] Palbociclib and its receptor were identified via docking, their relative stabilities were evaluated and their binding affinities estimated by using molecular docking simulations. Palbociclib and the protein were docked using default parameters in the Molegro software, after the water molecules, the cofactors and the ligand used for crystallization were eliminated of the original structure. A 15 Å radius was selected as search area, then five potential docking surfaces of the kinase were located and finally the docking wizard of the Molegro software made 1500 interactions before the lowest RMSD pose value was obtained and selected if the RMSD was lowest than 2 Å, from the ten values, only the lower pose was compared with the original ligand position. If the RMSD value was below 2 Å and the MolScore has the lowest value, we can conclude that the ligand and the receptor were properly docked. It is worth mentioning that one of the important findings is the presence of hydrogen bonding in the docking results. It is very important to identify protein residues interacting with the ligand and forming hydrogen bonding. Then, the distance and the name of the residues could be obtained. Those data are of great importance for determining if a ligand could be a potential inhibitor of different kinases. Furthermore, these results could be useful to explain the action mechanism of the inhibitor. CDK protein receptor was docked against the ligand obtained structure by using the Molegro software package to find the hydrogen bindings and explore the protein ligand interactions. Molegro use a docking algorithm called MolDock which is based on a new hybrid search algorithm, called guided differential evolution, which combines the differential evolution optimization technique with a cavity prediction algorithm. Differential Evolution (DE) was introduced by Storn and Price [49] in 1995 and has previously been successfully applied to molecular docking. [50,51] The use of predicted cavities during the search process, allows a fast and accurate identification of potential binding modes (poses). Docking of the protein ligand complex was mainly targeted only on to the predicted active sites. Docking simulations were performed by selecting 'Docking wizard' as the docking engine. The selected residues of the receptor at least 3 Å away of the ligand were defined to be a part of the binding site.
It is worth mentioning that the only way to have a validated docking study is to have RMSD value less than 2 Å. RMSD values bigger than 2 Å, are indication that not docking was obtained between ligand and receptor. In the case of the complex palbociclib-5L2I, a value of RMSD = 0.3909 Å was obtained, for palbociclib-2EUF complex was 0.0734 Å, while for the complex Palbociclib-1BI8 RMSD = 0.1036, was obtained, thus indicating that the docking process was successfully carried out.
Comparison of the molecular docking study with residues interacting, distances and the number of hydrogen bonding are presented in Table 2. Also, a molecular dynamics report in the literature, of palbociclib and kinase PDB id: 2EUF is presented for comparison purposes.

Interaction Studies
The ultimate goal of ligand-protein docking is to predict the predominant binding model(s) of a ligand with a protein of known three dimensions structure. To study the binding modes of bioactive compounds in the binding site of human protein, intermolecular flexible docking simulations were performed and energy values were calculated from the docked conformations of the inhibitor complexes. Docking studies yielded crucial information concerning the orientation of the inhibitors in the binding pocket of the target protein. Docking results of the ligand and its protein complex via docking software reveals that the relative energy values are À 147.12 a.u. À 125.78 and À 115.25 a.u. for 5L2I, 2EUF and 1BI8, respectively. Amino acid residues of the protein 5 L2I interact with the ligand, as shown in Figure 4. Docking of palbociclib with this protein yield 3 hydrogen bonding. Palboci-clib-2EUF complex present also 3 hydrogen bonding, the same than that reported in the literature with molecular dynamics. Besides, in the literature, hydrogen bonding interactions are with different residues, however, distances were not reported. Differences in the results of both methods, could be explained by the fact that in our report, molecular docking of the complexes were carried out using no restrictions in ligand flexibility or in the receptor, after curating some properties of the proteins during the preparation prior to the docking process.
It is important to mention that the docking process with each protein-ligand complex is different, as expected and seen in Figure 4, for the complex palbociclib-5L2I, where the hydrogen bonding and the distance are shown as blue lines.
Hydrogen bonding in the palbociclib-2EUF complex, Figure 5, could be explained by the oxygen2 charge (with two hydrogen bonding), of (À 0.563), and nitrogen7 charge of (À 0.480). No other hydrogen bonding was detected, meanwhile for the other two complexes three hydrogen bonding interactions were present. The complex palbociclib-2EUF, with 3 H bonding and the palbociclib optimized structure, are also shown in Figure 5, with the residue names within a distance of 3 Å.
The docking position of the ligand in the 2EUFkinase complex, is presented in Figure 6.

Binding Site of the Protein
The active site of 1BI8 comprises of amino acid residues forming hydrogen bonding (4) such as Leu79 at (3.53 Å and 3.33 Å), Arg78 at (3.41 Å), Val77 at (2.71 Å), however there are residues within a distance of less than 3 Å with no H bonding. These residues are part of the interaction pocket, such as Val76, Val78, Glu92, Phe135 and His73. While the active site of the protein 5 L2I presents hydrogen bonding (3) with Val101 (2.72 Å), Val101(2.90 Å) and Asp163 (2.73 Å). As   Chem. Biodiversity 2023, 20, e202200554 most of the amino acid residues in the active site are hydrophobic so they are the main contributors to the receptor ligand interaction. Amino acid residues of the receptors 5L2I, 2EUF and 1BI8, undergo interactions with the ligand, thus modifying the respective electrostatic surface as shown in Figure 7.

Conclusions
Palbociclib optimized structure was obtained by using DFT, and allowed to characterize and obtain the orbital surface of LUMO and HOMO for this drug. In the present work we have chosen two human cyclin dependent kinases 4/6 and other type of kinase, to study their interactions with palbocilcib. From the docking results, palbociclib could be identified as potential kinase inhibitor. When the receptor 5 L2I was docked with palbociclib the energy value obtained was (À 147.12 a.u.), for receptor 2EUF was (À 125.78 a.u.) and for the receptor 1BI8 the energy value was (À 115.25 a.u.) using Molegro docking software. These results indicate the relative stability of the complexes; thus, it is possible to conclude that the complex with more stability is palbociclib-5L2I, and that with less stability is palbociclib-1BI8 complex. The last kinase was chosen since it was crystallized with a ligand different to palbociclib. The oxygen O2 with charge of (À 0.563), and nitrogen N7 charge of (À 0.480) could explain the hydrogen bonding (3) in the palbociclib-2EUF complex. No other hydrogen bonding was detected. However, on the other two complexes four hydrogen bonding interactions were present. From the low RMSD values of the comparative docking study we conclude that when interacting with the receptors, palbociclib could be used as a poly pharmacological agent to inhibit the biological action of cyclin dependent or other types of kinases. Furthermore, the energy gap from the LUMO and HOMO, is of a convenient and appropriate value to allow palbociclib to react and get the chemical stability of the complex ligand kinase.

Quantum Chemistry
To characterize the ligand palbociclib, a quantum chemistry study was carried out. The geometry of the molecules has been fully optimized ay Density Functional Theory (DFT) [53][54][55] particularly at the B3LYP/6-311 + + G(d,p) level, [56,57] after the conformer with the energy minimum was obtained ( Figure 1). Frequency calculations have been carried out at the same computational level to verify that the structures obtained correspond to energetic minima (0 imaginary frequencies). For MEP calculations, molecular orbitals surfaces and energy gap of the orbitals used to obtain the ligand reactivity. [58] All the calculations have been carried out with the Gaussian-09 package, using the supercomputer Miztli. [59] Docking Bioinformatics is seen as an emerging field with the potential to significantly improve how drugs are found and brought to the clinical trials. Protein target were downloaded from database Protein Data Bank, PDB: http://www.rcsb.org/pdb/home/home.do).
In the present study biological database PDB (Protein Data  Bank) and software like Molegro, AutoDock Vina, and ChemDraw were used. PDB contains structural information of the macromolecules determined by X-ray crystallographic, NMR methods, etc. and the protein structures were obtained from PDB database. For the molecular docking study Molegro and AutoDock Vina were employed. Three proteins were selected for the molecular docking study, PDB id: 5L2I, 2EUF and 1BI8 and used as the target proteins. The two first proteins were crystallized using palbociclib as ligand, the last protein was selected as blank for comparison, since palbociclib not was employed for its crystallization All water molecules were removed and on final stage, hydrogen atoms were added if needed to receptor molecules. Palbociclib was selected as ligand considering its probed biological activity as commercial kinase inhibitor.