Exploring the effect and mechanism of Aloin A against cancer cachexia-induced muscle atrophy via network pharmacology, molecular docking, molecular dynamics and experimental validation

Systematic pharmacology and computational methods-based approach for the potential actions of Aloin A on CC-induced muscle atrophy

Construction of protein–protein interaction (PPI) network

PPI network was obtained from STRING database (https://string-db.org/, accessed on October 8, 2022). The species named “Homo sapiens”, and proteins corresponding to genes with high confidence were selected based on the minimum required interaction scores >0.7. The data were saved in TSV format and then imported into Cytoscape software (version 3.9.1) to visualize the PPI network.

GO and KEGG pathway enrichment analysis

To explore the biological function of potential targets of Aloin A on cancer cachexia, Gene Ontology (GO) analysis, including biological processes (BP), cellular components (CC) and molecular functions (MF), and the KEGG pathway were integrated Discovery (DAVID) (https://david.ncifcrf.gov/, accessed on October 10, 2022) database and conducted using bioinformatics.

Target identification based on the strategy of molecular docking

Three-dimensional (3D) structures of the key proteins and potential pocket were obtained from the RCSB PDB database (https://www.rcsb.org/, accessed on October 14, 2022) and DoGSiteScorer of Proteins Plus (https://proteins.plus/, accessed on October 14, 2022), respectively. Three-dimensional (3D) structure of Aloin A was downloaded from PubChem database (https://pubchem.ncbi.nlm.nih.gov/, accessed on October 14, 2022). Discovery Studio Visualizer 2019 has been used to hydrogenate proteins, remove water and ligand molecules. AutoDockTools ver. 1.5.7 was used to convert Aloin A and protein molecules into “pdbqt” format, and finally vina was used for molecular docking of HSP90AA1 (grid size: x = 0.847154, y = 29.73, z = 28.0159; box size: 58.392), AKT1 (grid size: x = 22.9599, y = −15.8276, z = −16.0423; box size: 72.98Å), CTNNB1 (grid size: x = 56.7304, y = 19.0907, z = 23.0016; box size: 93.012 Å), SRC (grid size: x = −10.4861, y = 45.1305, z = 23.0835 ; box size: 70 Å), AKT3 (grid size: x = 29.6112, y = 29.6112, z = −1.14428; box size: 80Å), TNF-α (grid size: x = 21.4268, y = 17.9164, z = 21.4864; box size: 76.992Å), RELA (grid size: x = 9.13265, y = 6.74423, z = 1.06598; box size: 124.31Å), BCL2L1 (grid size: x = 31.7336, y = 4.86445, z = 26.3051; box size: 110Å) with Aloin A. Ligand-macromolecule complexes with lower binding energy were regarded as more favorable potential target proteins and we selected binding energy less than −7.0 kcal/mol. The compound–target interactions and modes of binding were visualized using PyMoL and Discovery Studio 2019.

Target identification based on the strategy of molecular dynamic (MD) simulation

The related parameters of MD may better investigate the stability and conformational changes of a protein binding to the ligand and can be used to further explore its mechanism of action [20]. GROMACS 2021.5 software package was used to simulate the dynamics of the three systems. Proteins used OPLS-AA/L force field, and small molecular ligands were prepared with the LigParGen server (https://traken.chem.yale.edu/ligpargen). Each simulation system was placed in A SPCE water box model with a side length of 10 Å and randomly added sodium ions or chloride ions to maintain the system’s electrical neutrality. After the system was built, the energy optimization process was performed using the steepest descent method. Next, NVT and NPT balances of 1000 ps were performed, respectively. MD simulation production of each protein ligand complex runs up to 50 ns. Consequently, various geometric properties of the system, such as root mean square deviation (RMSD), root mean square fluctuation (RMSF), solvent accessible surface area (SASA), and hydrogen bond number, are calculated according to the trajectory for further analysis. The binding free energies (ΔGbind) including electrostatic interactions (ΔEelec), Vander Waals interactions (ΔEvdW), non-polar solvation energy (ΔGSASA) and polar solvation energy (ΔGpolar) were calculated using the molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) method implemented in GROMACS compatible tool “gmx_mmpbsa”, and the trajectory of the last 10 ns is selected for MM_PBSA calculation.

In vivo experiments to validate the anti-cachectic effect of Aloin A in mice bearing lung tumor

Tumor cell culture

Lung cancer cells (LLC) were provided by the Cell Resource Center, Institute of Basic Medicine, Chinese Academy of Medical Sciences (Beijing, China). The cells were cultured in Dulbecco’s Modified Eagle Medium (Biological Industries, Beit-Haemek, Israel) supplemented with 10% FBS (Biological Industries, Beit-Haemek, Israel) and 1% penicillin-streptomycin solution in digital incubator with 5% CO2 at 37°C. The cells were resuspended in phosphate-buffered saline (PBS) for establishment of cancer cachexia induced muscle atrophy animal model.

Model of cancer cachexia and drug treatment

For the establishment of cancer cachexia model, LLC cells (1 × 106 in 100 μL) or an equal volume of PBS was injected subcutaneously into the flanks of C57BL/6 mice as described previously [21]. One week after tumor inoculation, experimental mice were divided into three groups: control group, LLC model group, and LLC+Aloin A group. LLC+Aloin A group mice were orally administrated Aloin A (20 mg/kg/day) [22–24] every day for 16 days. The control and LLC group mice received the same volume of saline.

Rotarod test

Muscle function of tumor-bearing mice was evaluated by rotarod apparatus (ZB200-Taimeng, Chengdu, China) as previously reported [25]. On day 13, after acclimation for half an hour in the testing room, the mice (n = 8) were trained on a rotating rod (2–20 rpm) for 5 min. During habituation, the mice were immediately placed back on the rod if they fell off. After 0.5 h, muscle function was evaluated (3 trials, interval of 15 min) as the speed was accelerated from 5–30 rpm and latency to fall was recorded.

Inverted screen test

Muscle strength was evaluated using an inverted screen test [26]. Briefly, on day 14, the mice (n = 8) were placed on a wire mesh screen, and the screen was inverted. Muscle strength was calculated by measuring the hanging duration.

Wire grip test

The muscle strength of the mice was evaluated as previously described [27]. On day 15, all mice (n = 8) were allowed to grasp a metal wire (2 mm diameter), and the grasp of time was recorded until the mice fell off. The test score was calculated as the mean of at least three repetitions. Scores were assigned as follows: 0 points: fell off the wire within 10 seconds; 1 point: lifted one of the hind limbs; 2 points: tried to climb the wire; 3 points: grasped the wire with the front claw and at least one rear claw; 4 points: wrapped the legs and tail around the wire; and 5 points: tried to escape to the end of the wire.

Histological staining

Skeletal muscles (SM, n = 3) were minced, fixed in 10% formalin in 0.01 M PBS, dehydrated, and embedded in paraffin. 5 μm-thick slices of SM were cut. After deparaffinization, SM sections were rehydrated, stained with hematoxylin and eosin (H&E), and the slides were viewed, and photomicrographs were captured under a light microscope.

Real-time PCR

Real-time PCR (RT-PCR) was performed as described previously [28]. According to the manufacturer’s instructions, total RNA was extracted from muscle tissue (n = 6) with TRIzol (Sangon Biotech, Shanghai, China) and cDNA was prepared by using RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA). RT-qPCR was performed using StepOnePlus™ RT-PCR System (Qiagen, Hilden, Germany) with SYBR green (Qiagen, Hilden, Germany) and appropriate primers, and specific primer sequences are listed in Table 1.

Statistical analysis

The statistical data were analyzed with GraphPad Prism software (GraphPad Software, Inc., San Diego, CA, USA). One-way analysis of variance (ANOVA) was used to compare the differences between each group. P < 0.05 was considered statistically significant. All data are expressed as the mean ± SEM.

Data availability statement

All data generated or analyzed during this study are included in this article and are available from the corresponding author on reasonable request.

Network pharmacology analysis

PPI network analysis

To further clarify the correlations between Aloin A targets associated with cancer related muscle atrophy targets, a “target-target network” was built based on STRING results from 51 overlapping genes with high confidence (>0.7). Edges and nodes represent targets and interactions between targets in the PPI network, respectively. The color and size of the nodes demonstrate the magnitude of the degree. A larger sized node means a larger degree. A total of 52 nodes and 195 edges were involved in the PPI network, as shown in Figure 3A. The hub 10 targets included TP53, AKT1, SRC, CTNNB1, HSP90AA1, RELA, TNF, BCL2L1, PRKCZ and AKT3 (Figure 3B and Table 2).

Gene ontology enrichment analysis

For the 51 genes corresponding to potential therapeutic targets, we constructed a Gene Ontology (GO) enrichment analysis to identify relevant biological processes, cellular components, and molecular functions (Figure 4). We choose the top 10 terms for a brief demonstration. The results showed that a series of biological processes are involved in the treatment of cancer related muscle atrophy, including peptidyl-serine phosphorylation, negative regulation of apoptotic process, protein phosphorylation, inflammatory response, positive regulation of apoptotic process, positive regulation of smooth muscle cell proliferation, negative regulation of intrinsic apoptotic signaling pathway, and others.

KEGG pathway enrichment analysis and target-pathway network analysis

KEGG pathway enrichment analysis was performed, as shown in Figure 5A. The Y-axis indicates the KEGG pathways and the X-axis represents the p-value. The redder the color, and the smaller the value of p means more credibility and more importance. In this study, we selected the key 10 pathways from small to large according to the p-value for a brief demonstration. Among these pathways, the AGE-RAGE signaling pathway in diabetic complications, insulin resistance, chemokine signaling pathway, PI3K-Akt signaling pathway, MAPK signaling pathway, TNF signaling pathway, TLR signaling pathway, JAK-STAT signaling pathway and Wnt signaling pathway are considered the top priority. The top 10 signaling pathways with associated genes were mapped using Cytoscape3.9.0 software as shown in Figure 5B, indicating that Aloin A could alleviate cancer cachexia-induced muscle wasting by regulating multiple signaling pathways through multiple genes.

Molecular docking validation of Aloin A–key target interactions

Next, for validating the results of systematic pharmacology, we conducted molecular docking to evaluate the screened key targets. Based on the degrees of key targets in the PPI network and the KEGG results, Aloin A were analyzed by molecular docking with the hub genes HSP90AA1, AKT1, CTNNB1, SRC, AKT3, TNF-α, RELA, BCL2L1, respectively. The molecular docking results demonstrated that the binding energies (Vina scores) of Aloin A with the above key targets were in the range of −7.9 to −8.9 kcal/mol, (Table 3 and Figure 6A–6H), suggesting that Aloin A stably binds to those targets.

The bind energy of Aloin A with HSP90AA1 was −8.9 kcal/mol, and the binding affinity was contributed by the following: The conventional hydrogen bond formed with GLY137 and ASN51 residues; Pi-pi stacked bond formed with TRP162 residues; Carbon-hydrogen bond formed with GLY135 residues; The hydrophobic interaction (π-Alkyl) formed with LEU107 and VAL150 residues; Pi-pi T-shaped bond formed with PHE138 residues; Pi-sulfur bond formed with MET98 residues (Figure 6A).

The bind energy of Aloin A with AKT1 was −8.8 kcal/mol, and the binding affinity was contributed by the following: The conventional hydrogen bond formed with ALA226 and ALA295 residues; Carbon-hydrogen bond formed with ASN296 residues; Pi-pi stacked bond formed with TRP352 residues; The hydrophobic interaction (π-Alkyl) formed with ALA226 residues; Unfavorable Acceptor-Acceptor bond formed with HIS293 residues (Figure 6B).

The bind energy of Aloin A with CTNBB1 was −8.7 kcal/mol, and the binding affinity of was contributed by the following: The conventional hydrogen bond formed with ASP30 residues; Pi-pi T-shaped bond formed with PHE241 residues; Pi-pi stacked bond formed with TYR83 residues; The hydrophobic interaction (π-Alkyl) formed with PRO235 residues; Unfavorable Donor-Donor bond formed with THR233 residues (Figure 6C).

The bind energy of Aloin A with SRC was −8.5 kcal/mol, and the binding affinity was contributed by the following: The conventional hydrogen bond formed with ASP404, GLY406, GLY279 and GLU280 residues; Pi-sulfur bond formed with CYS277 residues; Pi-pi stacked bond formed with PHE278 residues; The hydrophobic interaction (π-Alkyl) formed with LYS295 residues (Figure 6D).

The bind energy of Aloin A with AKT3 was −8.3 kcal/mol, and the binding affinity was contributed by the following: The conventional hydrogen bond formed with GLU110, GLU94, GLN60, LYS58 and ASP90 residues; Pi-sigma bond formed with LYS58 residues; The hydrophobic interaction (π-Alkyl) formed with LYS58 residues (Figure 6E).

The bind energy of Aloin A with TNF-α was −8.3 kcal/mol, and the binding affinity was contributed by the following: The conventional hydrogen bond formed with SER154 residues; The hydrophobic interaction (π-Alkyl) formed with VAL38, LEU167, LYS53 and ALA51 residues (Figure 6F).

The bind energy of Aloin A with RELA was −8.3 kcal/mol, and the binding affinity was contributed by the following: The conventional hydrogen bond formed with GLU154 and GLY298 residues; Carbon Hydrogen Bond formed with GLY74 residues; The Pi-Sigma formed with ALA73 residues (Figure 6G).

The bind energy of Aloin A with BCL2L1 was −7.9 kcal/mol, and the binding affinity was contributed by the following: The conventional hydrogen bond formed with GLU129 and GLN125 residues; Pi-Alkyl bond formed with ALA119 residues; Unfavorable Donor-Donor bond formed with GLN111 residues; Pi-Donor Hydrogen bond formed with TRP169 and TYR120 residues;(Figure 6H). The above results implying that Aloin A may be exert anti-cachexia effect by interacting with these residues for stabling the key targets.

Validation of Aloin A-target interactions by molecular dynamics simulation

To better understand protein-ligand interactions, we selected representative molecular docking complexes for MD simulation and analysis. To investigate the stability of Aloin A binding to selected proteins, we calculated the RMSD of the skeleton atoms in the simulated trajectory. RMSD has been commonly used to estimate the degree of protein deviation from the original structure, which explains the stability of the complex [29]. As shown in Figure 7A, the plot of the HSP90AA1-Aloin A complex showed a stable equilibrium approximately after 15 ns, then a stable RMSD could be seen approximately at ~0.20 nm. The plot of AKT1-Aloin A complex showed a stable equilibrium approximately after 25 ns, then a stable RMSD could be seen approximately at ~0.3 nm. The RMSD value of the AKT3-Aloin A complex stabilized at about 0.2–0.28 nm between 18 and 50 ns. In general, the RMSD of the above three complexes remains stable after 15 ns. Among them, the stability of HSP90AA1-Aloin A complex is the best.

In addition, we analyzed the solvent accessible surface area (SASA) of the complexes by simulating the trajectory to evaluate the volume change of the complex. As shown in Figure 7B, the SASA of AKT1-Aloin A and HSP90AA1 complexes are similar and fluctuate less than the AKT3-Aloin A complex during entire simulation, indicating that although the SASA of three complexes is stable throughout the whole simulation process, the stability of AKT1-Aloin A and HSP90AA1 complexes is less affected by solvent and has better stability in aqueous solution.

The radius of gyration (Rg) reflects the binding tightness and degree of constraint of the system, and a higher Rg value is associated with a better chance of producing flexible ligands. Therefore, the higher the Rg value, the lower the stability [30]. As shown in Figure 7C, HSP90AA1-AloinA complex has the smallest Rg below 1.75 and is the most stable system compared to the Rg of AKT1-AloinA complex and AKT3-AloinA complex, their values fluctuate around 2.1 nm and 3.1 nm, respectively. In general, the Rg of the above three complexes is stable at 0–50 ns.

Hydrogen bond is one of the important parameters to reflect the stability of protein and ligand binding. Therefore, we further studied the variation of the number of hydrogen bonds in the complex during the MD simulation. As shown in Figure 7D, we found there was a maximum occupancy of 6 hydrogen bonds between HSP90AA1 and Aloin A, and the number of hydrogen bonds was stable between 1 and 4 during the MD simulation of 0–50 ns. The number of hydrogen bonds mostly concentrated between 1 and 3 for the complex of AKT1-Aloin A and the complex of AKT3-Aloin A.

RMSF mainly reflects the flexibility of residues of amino acids in the protein, which is important to deeply understanding the local conformational changes of protein side chains in MD simulated time [29]. The greater the fluctuation of RMSF of the system, the higher the flexibility of amino acids, resulting in unstable binding [31]. We further performed RMSF analysis to assess the positional fluctuation of each amino acid around its average position and stability of system. As shown in Figure 7E, there are large fluctuations in the protein terminal residues of the three complexes during the simulation process, it has little effect on the stability of the complexes because it is not the active site of the protein. Luckily, we found that there was little fluctuation in active sites of the protein in each complex. HSP90AA1-AloinA complexes fluctuated less at residues 101–107, 117–140, and AKT1-AloinA complexes fluctuated less at residues 122–128, 130–139, respectively, and AKT3-Aloin A complex fluctuated less at residues 77–78. In summary, the flexibility of these three compounds is low and has a good stability.

The analysis of Aloin A-target interactions by binding free energies

In order to better understand intermolecular interactions and the stability of complexes, we conduct the analysis of the binding free energy. The results showed that all binding free energy was less than zero, indicating that the reaction can proceed spontaneously (Table 4). The binding free energies of AKT3-Aloin A complex only has a van der Waals energy and electrostatic energy, but the other two complexes were decomposed into four energy components: Van der Waals energy, electrostatic energy, polar solvation energy and SASA non-polar solvation energy to get insights into their individual contributions. Among them, van der Waals force plays an important role in the binding process in this study. HSP90AA1-Aloin A complex displayed portrayed highest free energy (–79.450 kJ/mol), whereas AKT1-Aloin A complex possesses the least negative binding energy (−38.544 kJ/mol) (Table 4). To analyze the contribution of residues to protein ligand interaction, free energy decomposition per residue was employed (Figure 8). The results showed that MET-98, PHE-138 and LEU-107 make the greatest contribution to binding free energy for the HSP90AA1 system (Figure 8A). CYS-119, ILE-121 and TYR-122 make the greatest contribution to the binding free energy for the AKT1 system (Figure 8B), while GLU-115, ARG-116 and MET-117 contribute most to that of AKT3 system (Figure 8C).

Effect of Aloin A against cachexia-induced muscle atrophy in mice bearing lung tumor

Effect of Aloin A on body weight, muscle quality and adipose mass in mice bearing lung tumor

It is well known that cachexia is a complex metabolic disorder characterized by body weight loss, muscle wasting, and adipose tissue depletion. Therefore, we measured tumor free body weight, muscle mass and epididymis fat mass, and found that compared with the control group, tumor free body weight (Figure 10A), muscle masses, including gastrocnemius muscle (Figure 10B), quadriceps (Figure 10C), tibialis anterior muscle (Figure 10D) and soleus (Figure 10E), and epididymis fat mass (Figure 10F) were significantly reduced in LLC group. Luckily, we found that Aloin A (20 mg/kg) obviously recovered body weight and muscle mass reduction, but not epididymis fat mass in LLC tumor-bearing mice, suggesting that Aloin A could improve the cachexic condition of tumor-bearing mice.

Effect of Aloin A on muscle histology in mice bearing lung tumor

The results of HE staining revealed that compared with the control group, the cross-sectional area (CSA) of skeletal muscle in the LLC group was remarkably reduced (Figure 11A, 11B), and Aloin A (20 mg/kg) treatment obviously increased CAS in LLC tumor-bearing mice. Moreover, the muscle fibers in control group was more distributed between 1400-3400 μm2 while most of muscle fibers in LLC group are distributed between 200 and 1600 μm2. After treatment, we found that Aloin A increased the distribution and frequency of big fibers, most of them are between 1000 and 2600 μm2 in LLC tumor bearing mice (Figure 11C). These results further demonstrated that Aloin A can effectively rescue cachexia-induced muscle atrophy in tumor-bearing mice.

Effect of Aloin A on mRNA expression of key genes in skeletal muscle of mice

Based on the results of systems pharmacology, molecular docking, and MD, HSP90AA1, AKT1, and AKT3 may be key genes involved in the anti-cachectic effect of Aloin A against muscle atrophy in tumor-bearing mice, and we further validated these hub targets through RT-qPCR, and found that Aloin A (20 mg/kg) treatment significantly down-regulated the mRNA expression of HSP90AA1 (Figure 12A) while upregulated the mRNA expression of AKT1 and AKT3 (Figure 12B, 12C) in LLC tumor-bearing mice, suggesting that these targets are, indeed, contributing to muscle atrophy in tumor bearing mice, and Aloin A, at least in part, may alleviate the muscle atrophy via improving abnormal gene expression of HSP90AA1, AKT1 and AKT3.