Mulberroside A could serve as a pan inhibitor for the tyrosine kinase domains of the HER family

1. Introduction

Human epidermal receptors (HERs) have been shown to play a crucial role in a variety of cancer forms, including non-small cell lung cancer (NSCL), which kills more people than prostate cancer, breast, colorectal combined and around 20% of breast cancers, which are characterized by overexpression of human HER2 protein and associated gene (Majeed et al., 2021).

There are four HER receptors: HER1, also known as EGFR, HER2, HER3, and HER4. The members of the HER family are transmembrane receptor tyrosine kinases, which are enzymes that catalyze tyrosine phosphorylation, specifically the transfer of adenosine triphosphate (ATP) to tyrosine residues on protein substrates and the subsequent recruitment and activation of downstream signaling proteins. This results in the activation of downstream signaling pathways that promote cell proliferation, differentiation, migration, survival, adhesion, and angiogenesis. Numerous malignancies have been associated with a mutation that enables some tyrosine kinases to be constitutively active and/or occasionally act improperly.

The extracellular domain (ECD), a transmembrane segment, and an intracellular portion all comprise the structure of HER receptors. The intracellular region consists of a juxtamembrane segment and a functional protein kinase domain (with the exception of HER3, which lacks tyrosine kinase activity and must be activated in conjunction with another family member). The catalytic domain or the tyrosine kinase domain (Figure 1) contains the conserved ATP binding groove, which is required for ATP binding (Jura et al., 2009). The tyrosine kinase domain (TKD) is made up of an activation loop, a C-lobe, and an N-lobe, which contains a C α-helix. The activation loop, which is a polypeptide region outside the binding site cleft that, upon phosphorylation, induces kinase activation and C α-helix must adopt a particular, “active” configuration in order for a TKD to function (Wintheiser & Silberstein, 2021; Adams, 2003).

Figure 1. Demonstration of the structural parts of the tyrosine kinase domain of HER family.

There are two lobes (C-lobe and N-lobe) that are spanned by a loop-like structure called the activation loop (purple in color). The N-lobe contains the C α-helix (purple in color) which is essential for the catalytic activity of the enzymes. The sequence alignment reveals the residues that are primarily responsible for forming the binding site (green-shaded) and that are predicted to interact with mulberroside A, along with their RMSD values, which show the deviation from each other.

The development of acquired resistance to chemotherapy is one of the most perplexing concerns confronting oncologists. Although cancer cells initially respond to chemotherapy, later as treatment progresses, mutations and epigenetic alterations may occur, resulting in a high degree of plasticity of these cells that leads to the emergence of resistance. Due to receptor mutations, anti-HER resistance has risen over the years, necessitating an ongoing search for new inhibitors.

Numerous studies on the antioxidative, anti-inflammatory, and anticarcinogenic properties of dietary polyphenolic substances have been carried out extensively. In this regard, dietary polyphenols have been shown in preclinical models to display chemo-preventive and therapeutic benefits against a variety of malignancies. However, available medications have a high level of unspecific cytotoxicity, which has a detrimental effect on normal cells as well. As a result, the search for anticancer agents derived from natural medicinal plants has progressively increased throughout the years (Mechchate et al., 2021). Phenolic compounds are one of the most abundant and diverse classes of plant metabolites, with over 8,000 compounds described. Stilbenes, benzoquinones, acetophenones, flavonoids, anthocyanins, phenolic acids, catechins/epicatechins, ellagitannins, and proanthocyanidins are the various polyphenolic classes. The most extensively studied of these are (-)-epigallocatechin-3-gallate, resveratrol (found in grape skin), and curcumin (found in turmeric) (Pandey et al., 2009). In this context, certain polyphenols have been shown to act as chemosensitizers, enhancing the efficacy of chemotherapeutic agents when used in combination. For instance, it has been demonstrated that the flavonoid silibinin considerably showed activity against the resistant lung cancer cells to the tyrosine kinase inhibitor (TKI) erlotinib which led to a suggestion of using silibinin in combination with TKIs as a treatment regimen (Verdura et al., 2021). Curcumin is another polyphenol which, when combined with 5-fluorouracil, demonstrated antimetastatic activity (colon carcinoma) (Yang et al., 2017). Additionally, it has been reported that resveratrol enhances the anticancer effects of gemcitabine by reducing the emergence of resistance in pancreatic cancer cells and sensitizing glioma stem cells to radiation (Zhou, et al., 2019). Stilbenes are a non-flavonoid polyphenolic phytochemical with a 1,2-diphenylethylene nucleus. They hold exceptional promise for the prevention and treatment of cancer, due to their antioxidant, cell death activation, and anti-inflammatory properties, all of which are associated with low toxicity in vivo (Sirerol et al., 2016). Numerous methoxy and hydroxy derivatives of trans-stilbene have been shown to inhibit the growth of varied human cancer cells by causing cell apoptosis (Parida et al., 2018). Mulberroside A is a stilbenoid found in Morus alba (the white mulberry) and a diglucoside of oxyresveratrol. Due to its successful docking results, mulberroside A will be the focus of this study, and its effect on the dynamics of the proteins will be investigated further.

2. Methodology

2.2 Molecular docking

Prior to docking, all water molecules, ions, and ligands were eliminated from the enzymes. These enzymes were then prepared for docking by adding hydrogen atoms and charges (Charges were modeled by AMBER ff14SB for the residues) using Dockprep's UCSF chimera-building feature (Pettersen et al., 2004). Afterward, the ligand was obtained from PubChem and prepared for docking by adding hydrogen atoms and charges. Charges were modeled by AM₁-BCC (Austin model with bond and charge correction) for the ligand (Maier et al., 2015; Jakalian et al., 2002). The ligand structure was minimized and optimized by MM2 force field that built in chemoffice software (Hocquet & Langgård, 1998).

The docking process and the estimated free binding energy were handled by the AutoDock Vina (Trott & Olson, 2010). Up to 10 poses were produced for the ligand within the binding site of the enzymes, and the selection was for the pose that was characterized by a high binding energy. The docking process was validated (Figure 2) by redocking the native crystallized ligands (HER1; TAK285, HER2; TAK285, HER3; bosutinib, HER4; lapatininb), and the webserver DockRMSD was used to calculate the RMSD between the redocked and crystallized ligands (Bell & Zhang, 2019). The box grid dimensions were determined based on the production of the best proximity to the conformation of the inhibitors from the original x-ray crystal structure. The dimensions were as follows: HER1 (Center |18.544|34.556|4.399, size|37.219|19.332|32.395), HER2 (Center |8.879|1.626|24.569, Size |39.994|37.149|22.890), HER3 (Center|-46.127| 19.229 |21.339, Size |26.413|19.229|21.339), HER4 (center|-39.845|49.482|21.561, Size|39.904|33.612|32.051).

Figure 2. Validation of the docking process by comparing the conformations of the docked ligands (black) and the original ligands (brown).

The RMSD values indicating the difference between crystallized and redocked ligands for HER1, HER2, HER3, and HER4 are 1.125, 2.036, 2.515, and 1.341, respectively. The redocked native ligands have docking scores of -44.35, -41, -35, and -43.5 kJ/mol for TAK285 (HER1 and HER2) bosutinib (HER3), and lapatininb (HER4), respectively.

3. Results and discussion

Concerning the modeling of the missing parts of the receptor kinases, the zDOPE scores for HER1, HER2, HER3, and HER4 are -1.84, -1.94, -1.87, and -1.58, respectively; negative values indicate superior models. Upon docking mulberroside A into these domains, all displayed a high binding affinity; however, HER1 and HER2 displayed the highest affinity, with scores of -40.54 and -40.12 kJ/mol, respectively (Table 1). Compared to the docking scores of the co-crystallized native ligands (mentioned in the caption of Figure 2), the binding energy of mulberroside is very close.

As shown in Figure 4, five H-bonds were formed between the residues Ser₂₀, Leu₇₇, Met₉₃, Asp₁₅₅, and Phe₁₅₆ of HER1 and the ligand. The bonds that were established with Met₉₃ and Phe₁₅₆ were the strongest of the five, and these bonds remained within the range of H-bond length (which is 2.5-3.5 Å) throughout the entirety of the simulation. HER2 formed four H-bonds by the residues Ser₇₄, Leu₇₆, Leu₈₇, and Asp₁₅₄. The H-bond formed by Leu₈₇ is the most stable of the four because it did not exceed the length threshold for a strong H-bond while the simulation was being run. The ligand was unsuccessful to build strong H-bonds with HER3; upon docking, the two H-bonds were perfectly with good length (2.57 Å with Leu₉₂ and 2.39 Å with Ser₉₆) but unfortunately, did not survive the MD simulation. The five H-bonds that were recognized within the binding site of HER4 upon the binding of the ligand were with Thr₈₉, Glu₉₉, Arg₁₄₀, and Asp₁₅₄ residues. The strongest were the bonds with Asp₁₅₄ and Glu₉₉, which remained stable during the simulation.

Molecular mechanics/Poisson–Boltzmann surface area (MM/PBSA) calculation was performed for the last 100 frames (the last 1 ns) and the results revealed some relevant outcomes to the docking score except for the interaction with HER2, which was scored -52.16 kJ/mol, the lowest due to low electrostatic interaction (energy) and a high wasting of binding energy in polar solvation energy (Table 2). The binding energy with HER1 is -119.1 kJ/mol, and it is based on strong electrostatic and VDW forces. Even though binding to HER4 required a high polar solvation energy, there was a noticeable offset from the energy gained from VDW, electrostatic, and SASA interactions.

Figure 3 depicts the binding poses of mulberroside A within the four kinases. In the case of HER1, mulberroside A was close enough to interact with both the activation loop (via VDW interaction with Thr₁₅₄ and two H-bonds with Asp₁₅₅ and Phe₁₅₆) and the C α-helix (via VDW interaction with Met₆₅). In HER2, the ligand has successfully interacted through VDW with Met₆₅ of the C α-helix while no direct interaction with the activation loop was observed. In HER3, the ligand was unable to successfully establish direct interaction with the C α-helix; however, Val₁₅₇ of the activation loop was close enough to the ligand to engage in direct VDW interaction. Concerning HER4, no interaction was found between the activation loop and the ligand. On the other hand, VDW interaction was found with the C α-helix, which was mediated by the residue Met₆₅.

Figure 3. Depiction of the binding site, ligand occupation, hydrogen bond acceptor and donor spots, and 2D map of ligand-protein interaction for the TKDs of a) HER1, b) HER2, c) HER3, and d) HER4.

There were roughly 10 H-bond donors and 11 H-bond acceptor spots in the binding pocket of HER1. 10 donors and 8 acceptors made up HER2. There were 10 donors and 12 acceptors for HER3. 8 donors and 10 acceptors comprise HER4. This led us to deduce that an increase in acceptors within the binding site would have an impact on the binding affinity of mulberroside A. By modifying the responsible moieties, the binding affinity of the ligand could be increased by avoiding clashes (acceptor-acceptor clashes) with Lys₄₄ and Leu₈₇ in the HER3 binding site and with Cys₉₆ and Met₉₂ in the HER4 binding site.

For MD simulation, the RMSD plot (Figure 5a) over the duration of the simulation showed that the conformation of the TKDs had undergone a very slight shift. When the ligand bound to the TKD of HER1, the protein became more rigid, resulting in less conformational variation than when the enzyme was in its apo state. This conformational variability is caused by one of the terminals in the apo form fluctuating more, as demonstrated by the region 1—10 in the RMSF plot (Figure 5b). Whereas in holo TKD of HER2, region 150—200 was the most variable segment, representing a portion of the binding site and extending into the activation loop. Both the apo and holo TKD of HER2 exhibited behavior that is, to some extent, comparable. Both of them fluctuated more at the beginning of the simulation, but by 90 ns, they had reached a point of convergence. The region 40—60, which included the binding site and a portion of the C -helix, shared the same overtones as the region 150—200. When compared to HER1 and HER2, the TKD of HER3 appeared to have a greater degree of instability in regard to regions 40—60 (C α-helix region) and 150—200 (activation loop region). The apo form and the holo form had the same tone for region 150—200, but the apo form has a different tone for region 40—60 due to the fact that the apo form was more flexible. The last few nanoseconds of simulation could reveal that holo-HER3 became more stable and less prone to conformational variations than apo-HER3. The high RMSD records were the direct result of the extreme fluctuation that occurs at one of the terminals in the TKD of HER4. In the apo form, the fluctuation began at 50 ns, whereas in the holo form, it began at 170 ns. In general, the apo form is very dynamically active in the region of the terminal and the region 150—200, whereas the holo form is characterized by a calm motion to some extent, in contrast to the apo form.

All of the TKD receptors showed more compactness in unoccupied form rather than occupied form as shown in the radius of gyration plot (Figure 5c), the most noticeable being HER4 which upon binding of the ligand led to a decrease in the compactness by 2 Å which was the largest. The solvent accessible surface area calculation (Figure 5d) showed an expansion (unfolding) in the structures of HER1 and HER4 in the presence of the ligand comparable of its absence, without any remarkable change in HER2 and HER3.

Figure 4. The H-bonds distance (length) between the ligand and the protein during the 200 ns MD simulation, with the Y-axis representing the distance in nanometers while the X-axis representing the simulation time; the upper left square represents the H-bonds with HER1, which are five bonds, two of which are strong (the green and the yellow); the upper right square represents the H-bonds with HER2, which are four bonds, one of which is strong (the purple); the lower left square represents the H-bonds with HER3, which have shown inconvenient H-bonds because they are unstable and exceeded the threshold length; the lower right square represents the H-bonds with HER4, where the green and the brown could be considered as good H-bonds.

Figure 5. Displaying the RMSD plot depicted by section "a" (Y-axis for the deviation in nanometers and X-axis for the time), the RMSF plot depicted by section "b" (Y-axis for the deviation in nanometers and X-axis for the residue's ID number), the Rg plot depicted by section "c" (Y-axis for the change in the compactness in nanometers and X-axis for the time), the SASA plot depicted by section "d" (Y-axis for the change in the area in square nanometer and X-axis for the time).

In all plots, the red and black lines represent holo-HER1 and apo-HER1, respectively. In all plots, the blue and green lines represent holo-HER2 and apo-HER2, respectively. In all plots, the brown and yellow lines represent holo-HER3 and apo-HER3, respectively. In all plots, the purple and grey lines represent holo-HER4 and apo-HER4, respectively.

PCA for the backbone of the receptors was performed and the first three eigenvectors were calculated and evaluated; the first two eigenvectors were taken for consideration as the most dynamical conformations present (Figure 6; Figure 8). Calculation and identification of the essential dynamic (PCA) is a method to monitor the motion of the protein to recognize its biological function. Also, PCA can help in making the trajectories shorter to spot the protein configurational spaces with irregular rhythm. To put it another way, PCA demonstrates, that a relatively small number of modes are sufficient to account for the majority of the fluctuations (Balsera et al., 1996). Fifty frames structures were isolated from each vector (one and two particularly) to facilitate spotting the change in the proteins structure (Figure 8). In terms of HER1, the PCA for the apo-protein showed that the region 164—180 contained the most observable irregular motion, and this region was duplicated in fluctuations caused by the binding of the ligand. Additionally, a portion of the 1—32 region was more pliable in the holo form (Figure 7). In contrast, the RMSF plot of vector one uncovered the fact that a portion of the region 199—232 was more flexible in its apo form than in its holo form. According to the results of the PCA plane analysis for the TKD of HER1, the apo form traveled through more configurational spaces than the holo form did. Whilst the PCA plane for the TKD of HER2 (Figure 6) postulated that both the holo and apo form visited different spaces at most, the significant difference in the dynamics was related to the region around residue 166, as revealed by the RMSF for the vector one, whereas the vector two had not shown any change in its behavior. Mulberroside A did not exert a significant influence on TKD of HER3 due to the fact that the principal plane demonstrated a stability and occupation almost for the same configurational spaces regarding holo and apo-form (Figure 6). As can be seen in the RMSFs of both vector one and vector two, the small change in the conformation that occurred as a result of the ligand's occupancy was primarily attributable to a reduction in the flexibility of the multiple parts that make up region 1—100. According to the principal plane of the apo form, which completely visited different spaces with a very wide transition in the conformation in comparison to the holo-form, TKD of HER4 dramatically changed its conformations upon interaction with the ligand. The change was primarily represented by the C-terminal and the portion between 160 and 180. When these receptors are activated, their C-terminal regions undergo autophosphorylation and act as a binding site for signaling proteins that transfer the signal downstream (Kovacs et al., 2015; Monsey et al., 2010).

Figure 6. Illustration of the PCA plane; the upper left square indicates occupied (red) and unoccupied (black) HER1; the upper right square indicates occupied (blue) and unoccupied (green) HER2; the lower left square indicates occupied (brown) and unoccupied (yellow) HER3; and the lower right square indicates occupied (purple) and unoccupied (grey) HER4.

Concerning HER4, the fifty frames of vector one for both the apo (blue) and holo (khaki) forms were matched to better visualize the major conformational difference between them. These changes appear to originate from the highly fluctuated region 160—180 and the C-terminus.

Figure 7. An illustration of the RMSF for vectors one, two, and three for all receptor kinases.

The holo proteins are shown in red, while the apo proteins are shown in black. The docked and undocked HER1 is shown in the upper left square, the docked and undocked HER2 is shown in the upper right square, the docked and undocked HER3 is shown in the lower right square, and the docked and undocked HER4 is shown in the lower left square.

Figure 8. Presentation of the PCA for the TKDs of the four proteins.

The C -helix and the activation loop are colored purple, as should be noted. As illustrated by the aligned frames, the degree and pattern of the activation loop's fluctuation differ slightly between the holo and apo forms of HER1, as does the orientation of the C α-helix. Upon the binding of the ligand to HER2, the decrease in HER2's flexibility is readily apparent. The flexibility of HER3's activation loop and neighboring regions appears to be greater than that of other proteins, both in the occupied and unoccupied state. Regarding HER4, the degree and pattern of fluctuation for the C -helix and the activation loop did not change significantly, with the exception of the regions shown in Figure 6.

4. Conclusion

The dearth of chemotherapeutic agents that are still effective against cancers through the tyrosine kinase receptors due to pervasive resistance is a significant issue for the scientific community. Mulberroside A a natural polyphenolic compound has been approved for its capability to bind tightly to these receptors leading to significant changes in the dynamic of the vital parts of their structures. RMSD, RMSF, Rg, and SASA calculations have revealed a change in the conformation, compactness, and the folding pattern of the TKDs upon the binding of mulberroside A. These findings were also supported by the PCA, which elucidated a variation in the dynamic of the backbone of the proteins. Due to the polar properties of the binding site of the human epidermal receptors, the polyphenolic compound mulberroside A perfectly interacted with them. The C α-helix changed in behavior upon the binding of the ligand to the TKD in HER1, HER2, and HER3, which may lead to a change in the catalytic activity of the enzyme. The activation loop is affected as well in HER1, HER2, and HER3 and could lead to a prediction of changing the activation level of these enzymes. Regarding HER4, the region that was affected by the ligand is different and may affect the autophosphorylation process. Broadly speaking, mulberroside A is a promising tyrosine kinase domain inhibitor, and the HERs overexpressed cell-based assay is recommended for further investigation of its in-vitro activity.

Data availability statement

Protein sequences and structures are taken from Protein Data Bank:

The following link contains all molecular dynamic parameter (mdp) files:

https://drive.google.com/drive/folders/1zdkRAnzupDKNjmojd7n-vR3EkdZiq5Rd?usp=sharing

Ethical approval

The scientific committee of Al-Rafidain University College has approved the publication of this study, which employs only computational aspects and does not involve any human or animal subjects.