Abstract
Keywords
Introduction
COVID-19 is an infectious disease caused by the SARS-CoV-2 virus. Globally, as of 10 February 2022, there have been 402,044,502 confirmed cases of COVID-19, including 5,770,023 deaths, reported to the World Health Organization (WHO). 1 COVID-19 first emerged in Wuhan, China, in December 2019; the most common symptoms are fever, fatigue, chest pain, and some other organ dysfunction.2,3 The SARS-CoV-2 has been identified as a new type of coronavirus (CoV), which is an enveloped RNA virus with a diameter of approximately 80–160 nm. 4 The spike (S), envelope (E) and membrane (M) proteins which are some main ingredients of the SARS-CoV-2 virus, have been reported to form the envelope of the CoV, and the nucleocapsid (N) proteins form the capsid to pack the genomic RNA. The spike protein can bind to angiotensin-converting enzyme 2 (ACE2) on the cell membrane to allow the virus to enter the cell. 5 According to some studies, the transmembrane spike (S) glycoprotein consists of two functional subunits, and the combination of the S1 subunit with ACE2 proven to be key in coronavirus entry into susceptible cells. 6 In addition, SARS-CoV-2-3-chymotrypsin-like protease (3CL.pro) was found to be mainly responsible for viral replication, which may also be a therapeutic target. 7 COVID-19 cases have been reported to affect all seven continents, although a total of 10,095,615,243 vaccines of COVID-19 doses have been administered globally, the high transmissibility and the constant mutation of the virus still pose a large burden and challenge for prevention and treatment in the future.1,8
B.1.617.2 (Delta) variant was first identified in India in late 2020, 9 and the Delta variant had a higher speed of transmission than the previous. 10 On 14 July 2021, the WHO declared that the third pandemic of COVID-19 had begun and that the Delta variant played the dominant role in it. 11 However, the Omicron variant (B.1.1.529), a novel variant named on 26 November 2021, has become the most prominent variant in the world now, which was declared from the B.1.617.2 (Delta) variant. 12 The Omicron variant possessed 59 mutations in its genome within 36 mutations in the spike protein. 13 A high-throughput pseudovirus neutralization assay proved that Omicron pseudovirus could still target the ACE2 receptor for entry but was 2-fold more efficient than the Delta pseudovirus which hinted at its stronger infecting efficiency and the ability to escape neutralizing immunity induced by vaccines. 14 Previous studies pointed the expression difference of human transmembrane serine protease 2 (TMPRSS2) between the Omicron variant and Delta variant, which may cause the difference in virus replication and explain why the clinical symptoms caused by the Omicron variant seemed milder than the clinical symptoms caused by the Delta. 15 However, Omicron is still infecting large numbers of people in the world and has caused more deaths than the Delta variant in the United States due to the rapid speed of transmission. 16 In the UK, the situation was similar. According to the news, the UK's worst daily death peaked at 1359 on 19 January 2021. 17 On 12 January 2022, the total COVID-19 cases reached 3,60,70,510 in India and the daily positivity rate of COVID-19 was 11.05%. 18
Some traditional medicines such as Ayurveda and traditional Chinese medicine (TCM), have been widely used for treating infectious diseases and for any other medical use. A large number of active compounds derived from medical plants hold the potential to target virus replication, promote cytotoxic activity against SARS-CoV and influenza virus, inhibit virus cell attachment and anchor other key targets for antiviruses.
19
Bharat Krushna Khuntia
Chuankezhi injection (CKZI) is an effective TCM injection that is popularly used in asthma. CKZI consists of two main Chinese herbs, Yingyanghuo (YYH, dried leaves of
Network pharmacology is widely used to study the "complex protein/gene-disease" path, which can describe the complexity between biological systems, drugs, and diseases from the perspective of networks and explore the study of "herb-disease" related mechanisms. 27 In this paper, the mechanisms between CKZ and COVID-19 were explored by the network pharmacology method of traditional Chinese medicine, and the related targets were analysed to provide ideas for further clinical application. Molecular docking and dynamics simulation studies have both been widely used in drug discovery, which could help predict ligand-target interactions at the molecular level.28,29 We used ACE2, spike protein S1, and SARS-CoV-2-3CL pro as receptors and to identify effective compounds to treat COVID-19.
Methods
Obtain the Candidate Compounds and Related Targets
The Pharmacology Database and Analysis Platform of Traditional Chinese Medicine System (TCMSP, http://lsp.nwu.edu.cn/tcmsp.php) is a systematic Chinese herbal medicine systems pharmacology platform that provides abundant information about drugs, targets, diseases, and their interactions. 30 It is common to make drug analyses and discoveries based on ADME (absorption, distribution, metabolism, and excretion). We searched the herbs “Yinyanghuo” and “BaJiTian” as keywords in the TCMSP database, setting bioavailability (OB) ≥30% and drug-like (DL) ≥0.18 as screening conditions to select candidate compounds. 31 Then we screened the related potential therapeutic target proteins of candidate compounds via the TCMSP database. All the target names were paired to standard gene symbols via UniProt (http://www.UniProt.org/). The species was selected as “Homo sapiens (Human)”. Eliminating repeated drug targets, we obtained the candidate treatment target set of CKZ.
Obtaining the Candidate Targets of COVID-19
We used two databases to search for COVID-19-related targets, the GeneCards database (https://www.genecards.org/) and the NCBI (https://www.ncbi.nlm.nih.gov). Finally, taking a union of the results, we obtained a COVID-19–associated gene set. To explore deeply, the obtained CKZ target set and COVID-19–associated gene set were imported into a Venn Diagram (https://bioinformatics.psb.ugent.be/webtools/Venn/), and the overlapping gene targets were identified as the main potential therapeutic targets.
Construction of the Protein-Protein Interaction (PPI) Network
The above key targets were used to construct a PPI network on the STRING 11.0 platform (https://string-db.org). As shown in the provided image, the colour and size of the nodes represent the degree value, with larger and darker nodes often meaning a greater degree value. 32 In further exploration, we used Cytoscape 3.8.0 software to analyse the PPT network and to obtain the key treatment targets.
GO and KEGG Pathway Enrichment Analysis
To further reveal the gene interactions and related signalling pathways, we input the main therapeutic targets into the Metascape database (www.metascape.org/), and we set a
Component-Target Molecular Docking
Preparation of Compound Structures
The constructs of docking compounds MOL004386, MOL000098, MOL000422, MOL004373, MOL000358, MOL004391, MOL003542, and MOL003044 were all from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/). After they were imported into Chem3D software (v 14.0) for optimization as well as energy minimization using the MM2 module, we saved them as SDF format files as ligand molecules for molecular docking.
Preparation of Target Protein Structures
ACE2 (PDB ID: 1R42), Spike1 (PDB ID: 3BGF) and 3CLh (PDB ID: 6LU7) protein structures were obtained from the Research Collaboratory for Structural Bioinformation (RCSB) database (https://www.rcsb.org/).34,35 The protein structure was processed on the Maestro 11.9 platform. The protein was processed with Schrodinger's Protein Preparation Wizard, crystal water was removed, missing hydrogen atoms were added, and missing bond information was repaired, missing peptides were repaired, and finally the protein was subjected to energy minimization, and optimization of geometry.36,37
Processing and Optimization of Molecular Docking
Molecular docking was performed using the Glide module in Schrödinger Maestro software. The Protein Preparation Wizard module was used for protein processing. The receptors were preprocessed, optimized, and minimized (restrained minimization using the OPLS3e force field) through it. Compound structures were prepared according to the default settings of the LigPrep module. For screening in the Glide module, the prepared receptor is introduced, and the natural ligand of the protein or the predicted pocket according to the structural characteristics of the protein is selected as the molecular docking site. The site box size was set to 20 Å × 20 Å × 20 Å, the lattice parameter was set to 0.375 Å, and the docking produced pose was set to 50. Finally, molecular docking and screening were performed by the standard docking method (Standard Precision, SP).
Screening and Analysis of Docking Results
By analysing the action mode of compounds and target proteins, we obtained the interaction between compounds and protein residues, such as hydrogen bonding, π-π interaction, hydrophobic interaction, etc, and then referred to the docking scoring of compounds to speculate whether the compounds to be screened have a certain active effect.
Molecular Dynamics Simulations
MD simulations of protein complexes with compounds were performed using Desmond version 2020.Here, the molecular force field for MD simulations was chosen as OPLS3e and the system was solvated using the TIP3 water model. Neutralize system charge by adding ions. The energy minimization of the whole system is achieved using the OPLS3e force field (all-atom type force field). The geometry of water molecules, bond lengths and bond angles of heavy atoms are all constrained by the SHAKE algorithm. The continuous system is simulated by applying periodic boundary conditions and long-range electrostatics are maintained by the particle mesh Ewald method. The system was equilibrated using an NPT ensemble with a temperature of 300 k and pressure of 1.0 bar. The Berendsen coupling algorithm was used for the coupling of temperature-pressure parameters. Before performing the kinetic simulation, the model was first optimized: (1) limiting the solvent for energy minimization; (2) no limitation for energy minimization; (3) temperature 10 K, Berendsen thermodynamic method, running for 12 ps, resampling every 1 ps; (4) temperature 10 K, pressure 1.01 × 105 Pa, running for 12 ps; (5) NPT ensemble, temperature 300 K, pressure 1.01 × 105 Pa, running for 24 ps; (6) NPT ensemble, temperature 300 K, pressure 1.01 × 105 Pa, normal pressure release coefficient, running for 24 ps. (7) MD finished product: After the end of the upfront preparation, a 100 ns run was performed at a time step of 1.2 fs, and trajectory recording was performed every 50 ps for a total of 20,00 frames. The RMSD (root mean square deviation) of the main chain atoms was calculated and graphically analyzed to understand the nature of protein-ligand interactions.
Results
Screening Candidate Compounds and Potential “Drug-Disease” Targets
By searching the TCMSP database, we obtained a total of 23 eligible YYH compounds and 20 eligible BJT compounds (Table 1). After assessing the UniProt database and eliminating repeated drug targets, we finally obtained 207 CKZ potential targets.
Basic Information of Some Active Components of CKZ.
We obtained 321 predicted COVID-19 targets from the NCBI gene database and 4585 targets from the GeenCards database. After removing the duplicates, we finally obtained 4681 potential COVID-19 targets. Taking the union of the drug and disease target sets by the Venn Diagram, we obtained 74 overlapping targets, which were considered as potential therapy gene targets.
Construction of Protein-Protein Interaction (PPI) Networks
Based on 74 potential targets, a PPI network was established by importing the 74 official gene symbols of the targets into the STRING 11.0 database (https://string-db.org). The PPI network contained a total of 74 nodes and 1012 protein interaction edges (Figure 1). We then analysed the PPI network via Cytoscape 3.8.0 software. According to the analyzing results, the average degree was 54.70, and the average betweenness centrality was 0.0094. We identified 18 key targets (Table 2) and the top 8 key compounds (Table 3) in the following analysis, suggesting the main therapeutic targets and directions.

Potential target protein-protein interaction (PPI) network.
The partial information of 18 key gene targets.
The top 8 key compounds of CKZ.
GO and KEGG Pathway Enrichment Analysis
Multiple biological processes (BPs), cellular components (CCs), and molecular functions (MFs) were identified in CKZ treatment of COVID-19 based on 74 key targets. We obtained 217 dominant enriched GO terms. The top 20 GO terms are listed in Figures 2–4. (1) According to the analysis of biological processes, the core targets were connected mainly to cytokine response, negative regulation of cell population proliferation, and positive regulation of cellular component movement. (2) The related cellular components included mainly the side of the membrane, serine protease inhibitor complex, and perinuclear region of cytoplasm. (3) At the same time, the main molecular functions consisted of cytokine receptor binding, ubiquitin-like protein ligase binding, protein homodimerization activity, and nitric-oxide synthase regulator activity (Table 4).

Go enrichment analysis. Biological process (fig 2), cellular component (fig 3). molecular function (fig 4).
The top 18 KEGG pathway analyses are shown in Figure 5. The main KEGG pathways included the AGE-RAGE signalling pathway in diabetic complications, rheumatoid arthritis, lipid and atherosclerosis, MAPK signaling pathway, etc.

KEGG pathway enrichment analysis.
The Functional Analysis of Identified key Targets by Gene Ontology Analysis.
The docking results for target with different compound (kcal/mol).
Note: Binding energy fuction [38]:.
Component-Target Molecular Docking
Firstly, screened eight active compounds were molecularly docked with SARS-CoV-2 related proteins—ACE2, Spike1, and 3CLh target proteins, respectively. The full molecular docking results obtained by the SP method are shown in Tables 5–7. 38 Molecular docking results showed that the eight compounds had a good binding effect and high matching with the three target proteins (binding energies were less than − 5 kcal/mol). The complexes formed by compounds and proteins after docking were visualized using PyMOL 2.1 software to determine the binding mode of compounds and proteins, and the amino acid residues bound by compounds and protein pockets could be clearly seen according to the binding mode (Figures 6–8). The amino acid residues of quercetin that interact with the active site of the Spike1 protein are TRP-163, ASP-165, LYS-103, LYS-142, and ILE-144. MOL000098 is a flavonoid containing two benzene rings that can form π-π conjugated interactions with the active site amino acid (RP-163) and play an important role in stabilizing small molecules in the protein cavity. In addition, the compound can also form hydrogen bond interactions with amino acids ASP-165, LYS-103, LYS-142, and ILE-144 at the protein site, with short hydrogen bond distances and strong binding ability, which play an important role in small molecules in the anchor protein cavity. The amino acids of MOL004373 that interact with ACE2 protein include mainly ASN-210, GLU-208, GLY-205, LYS-94, etc, which indicates that the entry of compounds into the active pocket of the protein is driven by hydrogen bonds, which can effectively promote the formation of stable complexes between small molecules and proteins. MOL003542 matches well with the 3CLh target and is able to form multiple hydrogen bond interactions with amino acids (ASP-187, GLU-166, LEU-141) at the protein site, which plays an important role in anchoring small molecules in the protein pocket and indicates that small molecules bind well to the protein pocket. In summary, MOL000098, MOL004373, and MOL003542 have good performance in docking scoring and binding mode with ACE2, Spike 1, and 3CLh target proteins, which can form a stable complex with the protein, and have a correlation with the targets. Furthermore, we focused on top 11 potential targets and made docking with some strongly related compounds according to the previous analysis (Table 8). Among those, quercetin showed good binding effects with all top 11 targets, indicating its notable value for against COVID-19.

The binding mode of Spike1 protein with MOL000098. (A) The 3D structure of complex. (B) The electrostatic surface of MOL000098 with Spike1 protein. (C) The detail binding mode of MOL000098 with Spike1 protein. Yellow dash represents hyd1rogen bond distance or π-stacking.

The binding mode of ACE2 protein with MOL004373. (A) The 3D structure of complex. (B) The electrostatic surface of MOL004373 with ACE2 protein. (C) The detail binding mode of MOL004373 with ACE2 protein. Yellow dash represents hyd1rogen bond distance or π-stacking.

The binding mode of 3CLh protein with MOL003542. (A) The 3D structure of complex. (B) The electrostatic surface of MOL003542 with 3CLh protein. (C) The detail binding mode of MOL003542 with 3CLh protein. Yellow dash represents hyd1rogen bond distance or π-stacking.
The RMSD (Å) Results for Compound Before and After Docking.
The cavity (Å3) results for compound in target protein.
The docking results for target with different compound (kcal/mol).
Molecular Dynamics Simulations
To further investigate the interaction of the complex of protein with small molecules in a solvent environment, MD simulations of the complex were performed for 100 ns using molecular dynamics. The RMSD is calculated as the sum of all atomic position deviations of the conformation at a certain moment from the initial conformation. The larger the RMSD value, the more unstable the complex is. The wider the distribution of RMSD values, the more abundant the changes in molecular conformation. As can be seen from Figures 9–10, the average RMSD of the 3CLh protein complex with MOL003542 was less than 1.5 Å, and the complex reached dynamic equilibrium in a short time (5 ns), which indicated that MOL003542 matched the target of 3CLh well and was able to form a stable complex. In addition, the conformational change of the complex did not present a significant tomographic problem, which also indicates that the two can bind well. The average RMSD of ACE2 protein with MOL004373 compound was less than 1.8 Å, and the complex reached dynamic equilibrium in a short time (10 ns), which indicated that MOL004373 matched the ACE2 target well and was able to form a stable complex (Figures 11–12). Comparing with the two complexes, the two complexes have good stability and will not be separated from the active pocket of ACE, indicating that the compounds have good binding efficiency and inhibitory characteristics. The average RMSD of Spike1 protein with MOL000098 compound is less than 1.4 Å, which indicates that the complex of Spike1-MOL000098 is stable and that small molecules are able to form a stable match with the protein (Figures 13–14). In addition, we can also find that only a small part of the amino acids has greater conformational changes according to the RMSF diagram (Figures 10,12,14), which is mainly due to the fact that the amino acids of this part are located in the hinge region of the protein, themselves have greater flexibility, and their conformations are prone to some changes during the simulation. However, most of the amino acid conformational changes are small, which also reflects the stability of the complex and is consistent with the RMSD discussion.

RMSD plot during molecular dynamics simulations for 3CLh with MOL003542.

RMSF plot during molecular dynamics simulations for 3CLh with MOL003542.

RMSD plot during molecular dynamics simulations for ACE2 with MOL004373.

RMSF plot during molecular dynamics simulations for ACE2 with MOL004373.

RMSD plot during molecular dynamics simulations for Spike1 with MOL000098.

RMSF plot during molecular dynamics simulations for Spike1 with MOL000098.
Discussion
The 2019-nCoV has been identified to belong to the beta-BAT-SARS-CoV-2 lineage and is located in the Golgi apparatus, whose motifs are considered to also match those of human and mammals, such as zinc fingers, DNA-binding domains, and basic helix-loop-helix factors.
39
Over the past 2 years, COVID-19 has proven to be a large burden to many countries and a threat to life. Effective therapy against COVID-19 still needs to be explored, even if the COVID-19 vaccine and some antivirals have been widely used.
40
At the beginning of the key invasion of the virus, the SARS-COV-2 virus’s protein could bind to the host cell receptor ACE2; SARS-COV-2 3CLpro could generate some nonstructural proteins (NSPs) for virus replication. Thus, we focused mainly on these three proteins to identify the potential therapies.
41
Some studies have pointed out that compounds derived from medicinal plants such as
Some potential natural compounds extracted from plants for against COVID-19 infection.
We screened several key compounds of CKZ via the TCMSP database, and the top 8 components were considered to be closely related against COVID-19. Quercetin was the most active ingredient of CKZ and was related to 68 gene targets. Kaempferol was closely related to the regulation of 28 gene targets. In previous silico studies have proved that potential of ACE2, Spike protein S1, and SARS-CoV-2-3CL pro to act as treating targets in suppressing SARS-CoV-2 replication.
44
Our molecular docking results suggested that the top eight compounds all had good binding affinity with three target proteins. Quercetin, anhydroicaritin and 8-Iisopentenyl-kaempferol showed better binding effects among those 8 compounds. Selected active compounds also performed well binding ability with top 11 key drug-diseases targets. Based on the results of molecular docking, we also performed molecular simulation studies to further validate the well binding of quercetin, anhydroicaritin and 8-Isopentenyl-kaempferol with SARS-CoV-2 proteins. In summary, these compounds can form stable complexes with target proteins and were potential active molecules. It was worth noting that most key compounds are all classified as flavonoids. Flavonoids exist mainly in vegetable and fruit seeds and have shown significant treatment potential in anti-inflammatory, antioxidant and antiviral functions.
45
Almas Jabeen
In addition, we obtained 74 overlapping genes that were considered candidate genes in treatment. From the following PPI analysis, cytokines such as AKT1, TNF, IL6, VEGFA, IL1B, TP53, JUN, and CASP3 may be the hub genes among those. By analysing lung samples of COVID-19 patients, the expression of AKT1 was increased when compared with the expression of AKT1 in control groups.
48
IL-6 and IL-1β can promote inflammatory progression and playimportant roles in some key signalling pathways in severe COVID-19.
49
David M. Smadja
We also performed GO function and KEGG pathway enrichment analyses. Positive regulation of cellular component movement was the most key biological process. Membrane raft and cytokine receptor binding were the most significant cellular component and molecular functions, respectively. According to the KEGG analysis, the JAK-STAT signalling pathway and AMPK signalling pathway were considered key pathways in the treatment of COVID-19 by CKZI. The initiation of JAK-STAT signalling pathways is related mainly to high-affinity interactions between extracellular signalling cytokines and cognate receptors; the JAK-STAT signalling pathway is one of the significant promoters of CRS in COVID-19.
52
Bahjat
Conclusion
Network pharmacology, molecular docking and dynamic simulation approaches were used to identify potential targets and therapeutic mechanisms of CKZI in the treatment of COVID-19. Our results suggest that CKZI may be an effective and safe drug for COVID-19 treatment. Overall, these factors still need some work for validation.
Footnotes
Acknowledgments
Not applicable
Authors’ Contributions
Availability of Data and Materials
The datasets supporting the results of this article are included within the article. Further inquiries can be directed to the corresponding author.
Ethics Approval and Consent to Participate
Not applicable.
Competing Interests
The authors confirm that there are no conflicts of interest.
Declaration of Conflicting Interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by the Shanghai Changhai Hospital Scientific Research Fund, Natural Science Foundation of Shanghai, National Natural Science Foundation of China, (grant number 2019SLZ002, 2019YXK018, 21ZR1479200, 82170033).
