Network pharmacology and in silico approaches to uncover multitargeted mechanism of action of Zingiber zerumbet rhizomes for the treatment of idiopathic pulmonary fibrosis

2.1 Identification of bioactive constituents of Zingiber zerumbet

The active components reported in the ZZR were identified from the Traditional Chinese Medicine Systems Pharmacology (TCMSP)³³ (https://tcmsp-e.com/tcmsp.php), database and by performing literature search.³⁴

2.2 Targets excavating

The 3D conformers of the bioactive phytomolecules were retrieved from PubChem in an SDF format. This was subsequently fed into “PharmMapper,”^35,36 (http://www.lilab-ecust.cn/pharmmapper/) which predicted human targets. The predicted targets were uploaded to the UniProt database to get standard gene names and official gene symbols.³⁷ The targets involved in IPF were found utilizing the “GeneCards”³⁸ (https://www.genecards.org/) and “DisGeNET”³⁹ (https://www.disgenet.org/home/) databases with efficient keywords “idiopathic pulmonary fibrosis” and “lung fibrosis” for finding the differentially expressed genes in the disease condition based on the disease score. The top targets were sorted to build a network. The overlapped gene symbols, referred to as targets of zingiber zerumbet in IPF were then obtained by uploading the component gene and the disease gene using the Venny 2.1 (https://bioinfogp.cnb.csic.es/tools/venny) duplication method mapping tool.⁴⁰

2.3 Protein-Protein interaction and identification of potential components and targets for idiopathic pulmonary fibrosis

Step 1: The target proteins of the phytomolecules in IPF were used to construct the Protein-Protein interaction networks using String app embedded in Cytoscape 3.9.1 software.⁴¹ Cytoscape software served as a platform with multiple inbuilt applications, that are needed for the interaction network and analysis through the application manager. Using the string protein query, 1.7.1^42,43 and the selection of the homo sapiens species, the targets of ZZR in IPF were imported to the Cytoscape, with the confidence score cut-off limited to a minimum of 0.40 and no further extra interactions. Each component was sorted based on the number of nodes and edges, and more than 300 edges and 40 nodes were considered for additional virtual screening.

Step 2: The hub genes were found using the Cytohubba 0.1 built-in programme.⁴⁴ The hub targets were built in a separate subnetwork, and the components with the highest degree, closeness centrality, and betweenness were chosen for further potential component screening.⁴⁵

Step 3: Screening was done with the new technique to bring the highly disease-specific targets. From the Cytoscape, by selecting a string disease query from the dropdown box and entering a specified confidence score, highly expressed disease-specific targets were imported into the Cytoscape. The top 30 hub genes based on the degree were identified by employing the analyse network tool.⁴⁶

Step 4: The core targets and potential components (phytomolecules) were identified by merging the top 30 disease query target and enriched targets of ZZR using a merge tool in Cytoscape. The potential phytomolecules which interact with all the identified core targets were selected for molecular docking.

2.4 Enrichment analysis

Enriched hub targets of 30 nodes obtained from the subnetwork and were uploaded to the database for Annotation, Visualization, and Integrated Discovery (DAVID) database for bioinformatic analysis of KEGG (Kyoto Encyclopedia of Genes and Genome) pathway enrichment and GO (Gene Ontology) functional annotation to elucidate biological functions (biological process, molecular function, and cellular components) of targets of ZZR components in treating IPF, with significant screening criteria of P-value and enrichment score.^47–49 The enrichment bubble chart and bar chart were obtained from the bioinformatics tool (https://www.bioinformatics.com.cn/).

2.5 Compound-Target-Pathway network

The Compound-Target-Pathway network was built by merging the Compound-Target and Pathway-Target networks in order to better understand the interactions among the phytomolecules, targets, and pathways.⁴⁰ Based on the pathway analysis, the top 20 key pathways were considered with the enriched hub targets to build and analyse the pathway-target network. The relationship between the primarily screened compounds with their top subnetwork targets was developed and analysed.

2.6 Molecular docking

In silico docking was carried out using the Glide module of Schrodinger software.⁵⁰ Alternatively, free software such as AutoDock and Gromacs could also be used. Protein information regarding the docking was considered from the UniProt with the best PDB (protein data bank) having high resolution and screening criteria.

2.6.1 Protein and Ligand preparation

The structural coordinates (Solved by X-Ray crystallography) for the core targets were identified and downloaded to Schrodinger Maestro from the RCSB Protein Data Bank. These proteins were prepared by a three-step process via the protein preparation wizard of Schrodinger software.^51–53 Alternatively, AutoDock 4 software can be used which is freely available. The imported structure was pre-processed, hetero atoms and crystallographic waters were removed from the structure, H-bond was optimized to obtain the low energy state of protein at pH 7.4 using OPLS3e force-field.^54,55 Using the LigPerp programme, the structures of phytomolecules retrieved previously for biological target prediction were optimised to yield 3D coordinates.⁵⁶ The ligands’ possible ionisation states were generated using the Epik tool at pH 7.4, generated tautomeric states for each ligand, and the output was saved in maestro format.⁵⁷

2.6.2 Extra precision (XP) glide dock and binding free energy

A site map tool was utilised to determine and evaluate the druggable pocket in proteins that did not have binding ligands.⁵⁸ Subsequently, the receptor grid generation tool was used to generate a receptor grid around the co-crystalized ligand by picking it from the workspace. For the proteins, without any co-crystalised ligand-receptor grid generation tool with an entry option was used to generate a grid around the site with the highest SiteScore and bearing necessary amino acid residue to modulate protein activity. To determine and rank the potential components, rigid with XP-docking was performed and presented in an XP-visualizer, where we interpret the pose viewer file in a table and 3D visualization of the pose to interpret changes in the complex internal geometry.⁵⁰ Pose with a better favourable interaction score when an H-bond forms, rewards of the electrostatic and hydrophobic enclosure, lipophilic chemscore, and penalties at various conditions were all taken into account while choosing an MD simulation. In rigid docking, the hydrogen bond interaction with the important amino acids was evaluated as an independent quality measure. The Prime MM-GBSA of all the docked protein-ligand complexes was performed to identify the binding-free energy.^59,60 Subsequently, the protein-ligand complex exhibited a good docking score, binding free energy, and XP-dock results with the interactions were considered for the molecular dynamics (MD) simulation.

2.6.3 Molecular dynamics (MD) simulation

MD simulations have become an effective method for identifying the crucial residues in biomolecular interaction and assessing the stability of protein-ligand complexes. In order to examine the conformational stability and steady state of the protein-ligand complexes, 100 ns MD simulations were run using the Desmond module.⁶¹ Alternatively, NAMD molecular simulation along with VMD software can be used. By using the SPC (simple point charge) solvent tool, the system builder tool was employed to develop an orthorhombic-shaped simulation box around the docked ligand complex, with a buffer distance of 10 Å from any side of the complex to keep water molecules around to prevent false interference from occurring on active sites.⁶² By adding the necessary quantity of Na⁺ or Cl^- salts, the water-solvated system was neutralised and iso-osmolarity was maintained by adding the 0.15M salt concentration. The system was relaxed and minimised to a 100 ps local energy minimum and achieved by a gradient threshold till it reaches 25 kcal/mol/Å using the steepest descent (SD) method and the limited memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) algorithm. The custom seed MD simulation was performed with the 300K temperature and 1.0315 bar pressure, using the Nose-Hoover chain thermostat and Martyna-Tobias-Klein Barostat method coupled to an isotropic style for 100 ns using NPT.^63,64 At the end of the MD simulation, a simulation interaction diagram (SID) was generated to analyse the MD simulation parameters like root mean square deviation (RMSD) of protein α-carbon atoms to measure the average change in the superimposition of the selected atoms at a particular frame on the first frame, and protein root mean square fluctuation (RMSF) of the ligand fragment fluctuating while binding with the protein, and the hydrogen bonding plays a crucial role in the biomolecular interactions of the protein-ligand complex, which are then followed by hydrophobic interactions and are necessary for the drug’s metabolization, specificity, and adsorption. The strong binding of molecules to proteins preserves the protein’s compactness, as assessed by the radius of gyration (rGyr), and Solvent Accessible Surface Area (SASA) was calculated to quantify the surface area of a molecule accessible by a water molecule.⁶⁵

2.6.4 Principal component analysis (PCA) and Domain cross-correlation map (DCCM)

The PCA was employed to evaluate the dynamic motion of the protein affected by the K3R and 4AFZ binding required for the biological function to calculate eigenvectors and eigenvalues. The Desmond simulation trajectory was converted to dnd format using the VMD v1.9.4 programme.⁶⁶ The PCA was carried out using the BioEdit R package^67,68 and the interaction pattern of the binding site of ligand was analysed by performing the DCCM during the overall equilibrium period.⁶⁹

3.1 Bioactive constituents of Zingiber zerumbet rhizomes

The TCMSP and literature search yielded a list of the 39 active components that were reported in ZZR (Table 1).

3.2 Disease targets and compound targets mining

Using the keywords “idiopathic pulmonary fibrosis” and “lung fibrosis” as input, the databases “GeneCards” and “DisGeNET” yielded a list of 6696 and 803 disease targets, respectively. Overall, 416 potential targets in IPF were identified by performing duplication analysis and combined based on the disease score from both databases with a score above 50 in GeneCards and a disease score above 400 in DisGeNET. Candidate targets were identified from Venny 2.1 by uploading disease and compound targets of PharmMapper. kaempferol derivatives like kaempferol-3-O-rhamnoside (K3R), kaempferol-3-O-methylether (K3ME), kaempferol-3,4,7-O-trimethylether (K374TME), kaempferol-3,4-O-dimethylether (K34DME), kaempferol-3-O-(2,4-di-O-acetyl-α-L-rhamnopyranoside) (K3(23DARP)), kaempferol-3-O-(2″-O-acetyl) rhamnoside (K32AR), 4-O-acetyl afzelin (4AFZ) and 3,4-O-diacetylafzelin (34DFZ) were shared more than 60 common targets as shown in Table S1 (Supplementary file 1, Extended data⁷⁰).

3.3 Protein-Protein interaction

Step 1: A network of 76 nodes and 1010 edges was obtained for a component K3R and represented in Figure 1; similarly, networks were built for all the components and screened primarily using the criteria to obtain 25 components (Table S1, Supplementary file 1, Extended data⁷⁰), based on the topological parameters of networks. Targets with more than a 50-degree score are represented in gradient green colour and gradient light blue between 20-50 degrees. Targets with thick borders (more than 15) represented good betweenness centrality. A circular shape represented good closeness centrality of more than 0.666, hexagon represents closeness centrality score between 0.595 to 0.666, and octagon represented lesser than 0.595 closeness centrality.

Figure 1. Protein-Protein interaction network and its legends of gradient colour and border width.

Step 2: Table 2 lists the hub targets for all components based on the degree, betweenness centrality, and closeness centrality using Cytohubba.

Step 3: To identify significant targets, the top 200 disease-specific targets for IPF (DOID: 0050156) were imported into the Cytoscape as a new screening method with a minimum confidence score of 0.90. After examining the network, the top 30 targets with the greatest degrees were chosen (Table 2).

Step 4: Overall, seven significant components were identified with a minimum of six nodes and six edges after merging the disease-specific network and enriched targets of ZZR for molecular docking, these components with common targets were considered. According to Table S2 (Supplementary file 1, Extended data⁷⁰) and Figure 2, compounds such as K3R, K3ME, K374TME, K34DME, K3(23DARP), K32AR, and 4AFZ have shared common targets of TGFB1, MMP9, MMP2, EGFR, AKT1, SRC, and TNF.

Figure 2. Final potential components (hexagon) and core targets network (V shape), constructed after the screening of disease-specific targets.

3.4 Pathway enrichment analysis and Gene ontology functional annotation

3.4.1 Pathway enrichment

DAVID database was the keen gene ontology (GO) and pathway enrichment analysis tool with a significant P-value <0.05. The KEGG pathway was enriched by performing the enrichment analysis. Targets are enriched in 147 pathways, and the top 20 entries were considered based on enrichment score (Table S3, Supplementary file 1, Extended data⁷⁰), as shown in Figure 3 and 4, which comprise the TGF-β signaling pathway, PI3K-Akt signaling pathway, MAPK signaling pathway, FoxO signaling pathway, and TNF signaling pathway in IPF (Figure S1-S5, Supplementary file 1, Extended data⁷⁰).

Figure 3. KEGG enrichment analysis of target genes.

Figure 4. Zingiber zerumbet inhibits FoxO signalling pathway which is affected by the TGF-β and other enriched signalling pathways.

3.4.2 Functional annotation

The enriched targets were uploaded to get specific functional annotations in every discipline. The biological process includes collagen degradation with the highest 51.8 enrichment scores, as represented in Figure 5 and Table S4 (Supplementary file 1, Extended data⁷⁰). Among the cellular components, the ECM was mostly investigated with a 6.9 enrichment score, and the molecular functions of the included entry targets consisted of the tyrosine-protein kinase with a 17.8 enrichment score.

Figure 5. Gene ontology-functional annotation (BP: Biological process, CC: Cellular components, and MF: Molecular function).

3.4.3 Compound-Target-Pathway network

Each gene interacts with several active components, while each component targets several genes. The Compound-Target-Pathway network in Figure 6 shows compounds such as K3R, K3ME, K374TME, K34DME, K3(23DARP), K32AR, and 4AFZ interact with multiple targets. These targets disrupt the pharmacological activities of many pathways involved in the progression of IPF. The Compound-Target-Pathway (C-T-P) network was built with multi-components, multi-targets, and multi-pathways interaction, which can be correlated with a clear understanding of the pharmacological activity of ZZR.⁷¹

Figure 6. Compound (light blue-hexagon)-Target (green-V shape)-Pathway (amber-triangle) network to understand the multi-compound and target effect.

3.5 Molecular docking

The identified components of ZZR were docked with 3KFD(TGFB1); 4XCT(MMP9); 3POZ(EGFR); 1QIB(MMP2); 2AZ5(TNF); 2UVM(AKT1) and 4HXJ(SRC) PDB IDs. Proteins were prepared by updating the missed side chains or loops to generate states and generated a receptor grid. The proteins without bound ligand SiteMap was used and generated a grid, with a site score of 0.976, a D-score of 1.024, and a volume of 302.18 ų, the binding site of the 3KFD was visualised and assessed as a druggable pocket at site 1. The amino acids implicated in latent TGF-β1 activation are found at this Site. The site score 0.995, D-score 0.987, and volume 395.48 ų were found in 1QIB site 1 with the essential amino acid residues to inhibit the protein (Table S5, Supplementary file 2, Extended data⁷⁰).

As illustrated in Table S6 (Supplementary file 2, Extended data⁷⁰) and Figure S6, K3R formed a hydrogen bond with GLY 46 and ALA75 residues of TGF-β1 dimer and residues play a major role in activating latent TGF-β1. It also includes hydrophobic interaction with the reported inhibitory amino acid residues.^72,73 K3R interacts with ASN842 and LYS745 by hydrogen bond, and MET793 interacted by hydrophobic with EGFR.⁷⁴ K3R interacts hydrophobically with MMP9 amino acids of LEU188 and ALA189 residues, creating hydrogen bonds with PRO246, GLU227, and GLY186 residues.⁷⁵ K3R interacts with TYR119 with pi-pi interaction, SER60 with hydrogen bond and other hydrophobic interactions with LEU57, and TYR59 residues of TNF-α dimer.⁷⁶ Hydrogen bond interaction was noted between the K3R and GLU202 and hydrophobic interaction with LEU164, ALA165, and HIS201 residue of MMP2.⁷⁷ K3R exhibited higher binding free energy (∆G); -34.22, -57.35, -66.37, -48.20, and -57.85 with MMP9, TGF-β1, EGFR, TNF-α, and MMP2 respectively (Table 3).

As demonstrated in Table 3 and Table S7 (Supplementary file 2, Extended data⁷⁰), 4AFZ provided a decent docking score and binding free energy, however specific amino acid residual interactions were only seen with the MMP9, MMP2, and TNF-α. Following glide docking, the reported 2D interactions of 4AFZ with MMP9, MMP2, and TNF-α were shown in Figure S7 (Supplementary file 2, Extended data⁷⁰). 4AFZ shared two hydrogen bonds with the residues GLU227 and ALA191, as well as other significant hydrophobic interactions with MMP9.^75,78 The MMP2 inhibitory amino acid residues of GLU202 and ALA165 have formed hydrogen bond interaction with the 4AFZ.⁷⁷ Ligand formed hydrogen bond interaction with SER60 and GLN61 as well as hydrophobic contacts with TYR119 residues in TNF.⁷⁹ 4AFZ displayed a decent binding free energy (∆G); -59.52, -45.62, and -49.94 with MMP9, TNF-α, and MMP2 respectively.

The structure of K3R and 4AFZ were illustrated in Figure S8 (Supplementary file 2, Extended data⁷⁰) and considered for dynamics simulation at 100 ns based on good docking score except for AKT1 and SRC, hydrogen bond interaction scores obtained from the XP-visualizer (Table S8, Supplementary file 2, Extended data⁷⁰), good binding free energy, and substantial interaction with the desired amino acid residues of all proteins to inhibit their activity.⁸⁰ The hydrogen bonds were formed within the geometric criteria of 2.5 Å in all the protein-ligand complexes as shown in Table S9 (Supplementary file 2, Extended data⁷⁰) except EGFR-K3R of GLY791 (2.68 Å) and MMP2-4AFZ of ALA165 (2.79 Å). The pi-pi interactions were observed within the geometric criteria of 4.5 Å between the TNF-K3R and TNF-4AFZ complexes.

3.6 Molecular dynamics simulation of kaempferol-3-O-rhamnoside and 4-O-acetyl afzelin

TGFB1: Moderate hydrogen bond interactions were formed from the K3R to GLY46 (chain A), ASN66 (chain A), and CYS44 (chain A and B) with occupancy of 44%, 46%, and 47%, respectively, as shown in Figure 7a.⁷² Figure 8a shows that ALA75 had hydrogen bond and hydrophobic interactions in the binding pocket. RMSD Figure 9a illustrates the stable interaction of K3R with the binding pocket of the protein within 3 Å up to 43 ns, thereafter, fluctuated for 2 ns due to conformational changes that occurred during the spatial ligand fitting and maintained continuous equilibrium till the end of the simulation period. Even during the 2 ns fluctuation, ALA75 (chain A) and CYS44 (chain B) maintained a consistent contact, but some interactions were shifted between the dimer and exposed to the solvent system. Protein RMSD is stably within 2 Å throughout the MD simulation. Inhibition of TGF-β1 suppresses mesenchymal transition, angiogenesis, ECM deposition, and myofibroblast differentiation.⁸¹

Figure 7. MD-2D interaction of K3R with (a) TGFB1, (b) EGFR, (c) MMP9, (d) TNF, (e) MMP2 and 4AFZ with (f) MMP9, (g) MMP2, and (h) TNF.

Figure 8. MD-contact histogram of K3R with (a) TGFB1, (b) EGFR, (c) MMP9, (d) TNF, (e) MMP2 and 4AFZ with (f) MMP9, (g) MMP2, and (h) TNF.

Figure 9. MD-RMSD of K3R with (a) TGFB1, (b) EGFR, (c) MMP9, (d) TNF, (e) MMP2 and 4AFZ with (f) MMP9, (g) MMP2, and (h) TNF.

EGFR: MD simulation retains the stability of protein and ligand interaction by succeeding to keep the molecular docking interactions within the RMSD limit. As represented in the Figures 7b & 8b, strong hydrogen bond interactions were retained in the MD simulation of K3R to residues such as LYS745, GLN791, ASP800, MET793, and ALA722 with occupancy of 90%, 99%, 87%, 64%, and 31% respectively and bridged hydrogen bonds were noted with ARG841, ASN842, and ASP855 (occupancy of 39%, 58%, and 61%).⁷⁴ Amino acid residues LEU844, CYS797, and ALA743 formed hydrophobic interaction, as illustrated in Figure 8b. The ligand fitting to the protein remained consistent throughout the simulation, as illustrated in Figure 9b, with 2Å RMSD being the same as protein RMSD. Inhibition of EGFR helps to decrease mucus secretion, collagen deposition, fibroblast proliferation, and disrupts lung morphogenesis.⁸²

MMP9: The protein and K3R interaction stabilised at the binding pocket. A strong and direct continuous hydrogen bond interaction was observed between K3R and GLU227 and ALA191 residues with 99% and 61% of occupancy (Figures 7c and 8c). The continuous Pi-Cation (Zn) interactions were observed with HIS226, HIS230, HIS236 and GLU227. While ALA189 and LEU188 interact via bridging interactions with water molecules to K3R with occupancy of more than 50%.⁷⁵ From Figure 9c, the RMSD of the ligand and protein resides within 2Å. Because MMP9 is a zinc-dependent, the compound interaction required a pi-cation binding to inhibit the protein. Figures 7f and 8f illustrated 4AFZ to PRO246 formed a continuous and strongest hydrogen bond interaction at an 83% occupancy, along with the 100% Pi-Cation interactions of HIS226, GLU227, HIS230, and HIS236 via a Zn coordination, which is required for a compound’s potency to inhibit MMP9.⁷⁵ Figure 9f illustrates the RMSD of the 4AFZ and MMP9 within the limit of 2.5 Å and maintained the phase of stabilization throughout the period of simulation. Inhibition of MMP9 decreases the activation of latent TGF-β1 and extracellular matrix deposition.⁸³

TNF: From the Figures 7d and 8d, K3R displayed continuous binding stability in the binding pocket of TNF due to strong hydrogen bond interactions with TYR151 (Occupancy of 53%), moderate H-bond with SER60, and GLN61 through a bridge contact (Occupancy of 54% and 34%) and a strong 93% Pi-Pi interaction with TYR119, as well as other important hydrophobic interactions involving LEU55 and TYR59.⁷⁶ RMSD of ligand-protein interaction resides within 2 Å, and interaction was stabilized during the simulation (Figure 9d).

Figures 7h and 8h show a moderate hydrogen bond formed from the 4AFZ to SER60 via a water bridge with a 36% occupancy and other key hydrophobic contacts to TYR59 and LEU57.⁷⁶ Figure 9h depicts the first 10 ns of fluctuation and maintains 1.5 Å RMSD for the remainder of the 90 ns simulation time to remain stable. Inhibition of TNF-α decreases inflammatory reactions and activation of macrophages.⁸⁴

MMP2: From Figure 7e, a strong hydrogen bond interaction was observed between ILE222 residues and K3R through the water bridge at an occupancy of 55%.⁷⁷ Hydrophobic interaction was observed in Figure 8e between the residue GLU202, GLY162, LEU164, TYR193, and TYR223 and the ligand. From Figure 9e, the RMSD was found to be within 2 Å for the Ligand and protein interaction.

Moderate water-bridged hydrogen bond interactions were formed from 4AFZ to ILE222 at 44%, ALA220 at 57%, HIS201 at 40%, and PRO221 at 35% occupancy as represented. The continuous Pi-cation interactions were seen to HIS201, GLU202, HIS205, and HIS211 through a Zn, which is required to inhibit MMP2 (Figure 7g and 8g). The ligand fitting to the MMP2 remains consistent throughout the simulation at 2.5 Å RMSD (Figure 9g) and maintains constant contact with the residues as depicted in the Figure 8g.⁷⁷ Inhibition of MMP2 suppresses β-catenin and lowers the TGF-β1 levels which favors the EMT. Inhibition promotes ECM degradation by affecting collagen degeneration.⁸⁵

With the exception of a few individual residues that fluctuated in the range of 40–50 and 150–160 of the EGFR–K3R complex, and 90-100 residues of the TNF-4AFZ complex, all other complexes exhibited good and lower RMSF values (Figure S13, Supplementary file 2, Extended data⁷⁰). The RMSF of K3R and 4AFZ ligands were stabilized at lower levels of less than 1-2 Å to their respective protein complexes (Figure S14, Supplementary file 2, Extended data⁷⁰), whereas the RMSF of K3R fluctuated considerably more with the TGF-β1 complex. In most of the complexes, an intramolecular hydrogen bond was observed. The rGyr of K3R and 4AFZ ranged between 4.00 ± 0.08 Å and 4.44 ± 0.10 Å respectively and maintained stability. The SASA of K3R was shown to be very less in TNF and MMP9, at 70 and 120 Ų respectively, while the other complexes with K3R and 4AFZ had an equilibrium of 240-260 Ų (Figure S15 & S16, Supplementary file 2, Extended data⁷⁰).

3.7 Principal component analysis (PCA) and Domain cross-correlation map (DCCM)

The three principal components are represented in the Figure 10 for all the complexes of K3R with TGFB1, EGFR, MMP9, TNF, and MMP2, and 4AFZ with MMP9, TNF, and MMP2. The complex conformational variations are shown in the 2D PCA scatter plot to indicate system fluctuations by eigenvectors of PC1, PC2, and PC3. The biological functions are governed by the motion mode of 20 PCA patterns and represented in three PCA modes. The density of the distribution of the dots demonstrated stable system conformations. According to the Figure 10, except for TGFB1-K3R, PCA of all complexes exerted denser conformation, but EGFR-K3R has a small variance in the eigenvalue rank. The first three eigenvalues (PC1, PC2, and PC3) of TGFB1-K3R contributed 62.15, 31.74, and 2.09% of variables, respectively, while the other systems generated lower eigenvalues of variance (Figure 10), indicating higher stability.

Figure 10. 2D-PCA scatterplot, K3R with (a) TGFB1, (b) EGFR, (c) MMP9, (d) TNF, (e) MMP2 and 4AFZ with (f) MMP9, (g) TNF, and (h) MMP2.

DCCM analysis was used to examine the movement correlation of the protein complex system. When there is a positive value, the same direction emerges (positive correlated movement), however, when there is a negative value, the opposite way appears (anti-correlated movement). Except for TGFB1-K3R, the protein-ligand complex systems were shifted in the same direction as they shifted towards the cetacean blue, and dark cornflower blue implies the positive values in the Figure 11.⁶⁹ TGFB1-K3R correlated negatively by moving in the opposite direction towards the negative value. The colour distribution of the TNF-K3R was demonstrated to deflect the residual strength.

Figure 11. Domain cross-correlation map (DCCM) analysis, K3R with (a) TGFB1, (b) EGFR, (c) MMP9, (d) TNF, (e) MMP2 and 4AFZ with (f) MMP9, (g) TNF, and (h) MMP2.