Anti-inflammatory or anti-SARS-CoV-2 ingredients in Huashi Baidu Decoction and their corresponding targets: Target screening and molecular docking study

2.1 Preparation of active ingredients in HSBD

The active ingredients of all herbs that make up the HSBD were retrieved from the Traditional Chinese Medicines database (TCMdb) and the online public database Traditional Chinese Medicine Systems Pharmacology (TCMSP) using the standard herb names as the search terms. Then, we imported all active ingredients with anti-inflammatory or antiviral effects into Discovery Studio 2020 (DS2020) software for ligands preparation and minimization.

2.2 Screening active ingredients with good druggability

In the early stage of drug development, targeted selection of active ingredients according to the Absorption, Distribution, Metabolism, Excretion, and Toxicity (ADMET) properties is necessary to improve the success rate of drug development and reduce the waste of funds (Arora et al., 2016). We calculated ADMET properties, including blood–brain barrier penetration (BBB) and human intestinal absorption (HIA). Toxicity Prediction by Komputer Assisted Technology (TOPKAT) is a fast and accurate prediction of the toxicity and environmental effects of active ingredients based on the 2D structure of the molecule (Khalaf et al., 2020). The toxicological properties to be considered in this study are Rodent Carcinogenicity, oral LD50 in rats, Mutagenicity, and Aerobic Biodegradability. Active ingredients obtained by ADMET prediction and TOPKAT were then filtered through Lipinski's and Veber's rules.

2.3 Screening the target proteins of HSBD for COVID-19

The reverse finding target process in DS2020 was to match the active ingredients with good druggability against the target proteins. We then converted the target protein names into standard gene names using the UniProtKB search function in the UniProt (https://www.uniprot.org/) database. COVID-19-associated targets were screened from the Human Gene Database (GeneCards; http://www.genecards.org/) using the keywords “SARS-CoV-2” and “COVID-19”. The overlaps between COVID-19 and HSBD targets were supposed to be the potential targets of HSBD for the treatment of COVID-19. Moreover, the efficacy of these potential targets was also identified and distinguished, such as anti-inflammatory targets or anti-SARS-CoV-2 targets.

2.4 Construction of “herb-active ingredient-target” network and protein–protein interaction (PPI) network

A “herb-active ingredient-target” network was constructed using Cytoscape 3.8.0 software to show the active ingredients in HSBD and their corresponding targets in the treatment of COVID-19. The search tool for the Retrieval of Interacting Genes database (STRING, https://string-db.org/) was performed to analyze the potential targets of HSBD for treating COVID-19. The research species was set to “human species” for the PPI network. The results could reflect the core targets among these potential targets.

2.5 Molecular docking verification

The 3D crystal structures of the key target proteins were obtained from the RCSB Protein Data Bank (PDB, https://www.rcsb.org/). Then, we need to prepare the receptor binding site for the target protein with the help of DS2020 software. Molecular docking calculations were performed using CDOCKER, which can simulate the hypotheses before the implementation of specific experiments (Gentile et al., 2020). We applied a consensus scoring process based on different scoring functions (-CDOCKER ENERGY and -CDOCKER INTERATION ENERGY) to determine the best candidate active ingredients for the molecular docking process.

2.6 Molecular dynamics simulation and binding free energy calculation

The ff14SB force field was used for the protein part. The force field for ligands was generated from the generalized Amber force field potential using the antechamber module of Amber 18. All the systems were solvated in a 10 Å truncated octahedral water box using a TIP3P solvent model under periodic boundary conditions. Firstly, a total of 5000 cycles of steepest descent and then 5000 cycles of conjugated gradient minimizations were carried out with restrictions on all solute molecules. Secondly, a total of 5000 cycles of steepest descent and then 5000 cycles of conjugated gradient minimizations were carried out again without any restriction. Thirdly, a 300 ps heating procedure followed by a 500 ps equilibrium simulation was conducted. Finally, a 100 ns MD simulation was performed for each system under NPT conditions to produce a trajectory without any constraints. The temperature was maintained at 310 K by coupling to a Langevin heat bath using a collision frequency of 1 ps-1. The step side was set to 2 fs in all the simulations.

The MMGBSA method was used to calculate the binding free energy for the systems. The calculation formulas are as follows:

In Equation (2), the EMM, G_sol, and TΔS represent the gas-phase molecular mechanics' energy, the solvation free energy, and the vibrational entropy term, respectively. EMM is the calculated as the sum of internal (E_int), van der Waals (E_vdw), and electrostatic (E_ele) terms. G_sol represents the solvation energy, which is composed of the electrostatic solvation free energy (G_GB) and the nonpolar solvation free energy (G_SA). The binding free energies were the average values of calculating 4000 snapshots sampling from the last 40 ns trajectory.

3.1 The active ingredients of HSBD

99 active ingredients in HSBD were retrieved from TCMdb and TCMSP, as detailed in Table 1. Since this study required consideration of receptor-ligand interactions, the ligand preparation process is crucial. Different protonated states, isomers and tautomers usually have different 3D geometric structures and binding properties. If the exact binding configuration was unknown, one approach was to enumerate several likely configurations before screening the active substance and running docking protocols. Ligand Preparation helped us to accomplish this task. Therefore, the number of active ingredients reached 298 after ligands preparation and minimization.

3.2 Screening active ingredients with good druggability

We used DS2020 to evaluate the ADMET properties of 298 active ingredients. The results are shown in Fig. 1. The blue ellipse represents the 99 % confidence interval of the absorption, and the green ellipse represents the 99 % confidence interval of the BBB. An important factor that promotes the efficacy of drugs in the body depends on the absorption of drugs in the intestinal tract (Song et al., 2014). The absorption outside the 95 % ellipse tends to drop pretty quickly. The BBB is a significant barrier that blocks the physical diffusion of most proteins or compounds from the blood to the brain (Quaegebeur et al., 2011). Notably; the screened compounds must have blood–brain barrier permeability to ensure their activity during systemic administration (Banks, 2016; Albarracin et al., 2012). If a molecule is outside the 99 % confidence interval of the BBB model; the prediction of the molecule is considered unreliable. 95 active ingredients in the intersection of blue and green ellipses were retained for further screening.

Close

Fig. 1

ADMET properties of active ingredients.

Export to PPT

Adopting the TOPKAT process saves time to market, reduces animal testing, and assesses health risks (Wang et al., 2021). Based on the TOPKAT results, we removed all active ingredients from our list that exceeded the four properties' (Rodent Carcinogenicity, oral LD50 in rats, Mutagenicity, and Aerobic Biodegradability) optimal prediction spaces (OPS). The number of active ingredients was down to 60. Lipinski's and Veber's rules can determine whether the compounds have good absorption, penetration, and oral bioavailability (Lipinski et al., 1997). Finally, the 60 active ingredients were further filtered by applying reverse finding target, reducing the number of active ingredients to 48.

3.3 The target proteins of HSBD for COVID-19

We found that the anti-inflammatory or anti-viral target proteins of the 48 active ingredients in HSBD were Antigen peptide transporter 1 (TAP1), Retinol-binding protein 4 (RBP4), Interleukin 1 receptor antagonist (IL1RN), Centromere-associated protein E (CENPE), Beta-secretase 1 (BACE1), Transthyretin (TTR), Tyrosine-protein kinase Fyn (FYN), Thyroid hormone receptor alpha (THRA), Pulmonary surfactant-associated protein D (SFTPD), Nuclear receptor subfamily 0 group B member 2 (NR0B2), Cellular tumor antigen p53 (P53), SRSF protein kinase 1 (SRPK1), RAC-alpha serine/threonine-protein kinase (AKT1) and Protein Mdm4 (MDM4). Next, we used “COVID-19” and “SARS-COV-2” as keywords to search for targets related to COVID-19 in the GeneCards database. Table 2 showed 8 overlaps between COVID-19 and HSBD targets, namely RBP4, IL1RN, TTR, FYN, SFTPD, P53, SRPK1, and AKT1.

3.4 Analysis of the “herb-active ingredient-target” network and PPI network

The 8 potential targets and related herbs and their active ingredients in HSBD were used to construct the “herb-active ingredient-target” network by Cytoscape-3.8.0 software. 9 essential herbs were Huangqi, Tinglizi, Xingren, Mahuang, Huoxiang, Gancao, Dahuang, Caoguo, and Mahuang. As shown in Fig. 2, 9 herbs of HSBD were marked in green; 48 active ingredients screened from 9 herbs with anti-inflammatory or anti-SARS-CoV-2 effects were marked in pink; 8 potential targets of HSBD are marked in yellow, and the effects of the targets are marked in red. The black line represented that a specific active ingredient comes from a specific herb; the gray line represented the interaction between the active ingredient and the target, and the red line represented a particular target against inflammation or SARS-CoV-2. Some active ingredients were common to several herbs. For example, Kaempferol1, Kaempferol2, and Kaempferol3 were common ingredients of Huangqi, Tinglizi, and Mahuang. In addition, each target is bound to two or more active ingredients, indicating that these targets could be affected by multiple active ingredients simultaneously to exert different effects. The result revealed the mystery of the “multi-component and multi-target” of TCM formulas.

Close

Fig. 2

“herb-active ingredient-target“ network.

Export to PPT

Most proteins do not function in isolation, and their interactions with other proteins define their cellular functions (Liu et al., 2018). Understanding PPI will help to identify pharmacological targets and design drugs (Chen et al., 2010). The PPI analysis of 8 target proteins was carried out using the STRING plug-in of Cytoscape 3.8.0 (Fig. 3). There were 6 nodes and 12 edges in the PPI network. 2 disconnected targets (IL1RN and SFTPD) were excluded from the network. The nodes represent targets, and the size of the nodes represents the degree value. The larger the shape of the nodes was, the larger the degree value and the more important the targets were. Each edge represented an interaction between targets; the thicker the edge, the higher the combined score. The minimum score was set as the highest confidence of 0.9 to ensure high cross-confidence information. The targets within the highest confidence are P53, AKT1, TTR, and RBP4. Based on Fig. 3, P53 and AKT1 with higher degree value and confidence scores played more important roles in the network, which were likely to be the core therapeutic targets of HSBD on COVID-19.

Close

Fig. 3

PPI network.

Export to PPT

3.5 Analysis of molecular docking

Molecular docking is a calculation method used to predict the binding affinity of two molecules, usually a small molecule ligand and a target protein (Banchi et al., 2020). It can be seen from Table 2 that the PDB IDs corresponding to P53 and AKT1 were 2X0U, 2XOV, 3CQU and 3CQW. We used 2X0U and 3CQW with good resolution and few mutations as target receptors for docking studies. There are 48 active ingredients that were docked into the protein active site in this study. The molecular docking results (top 5) evaluated by 3 different scoring functions were shown in Table 3, including -CDOCKER_ENERGY, -CDOCKER_INTERATION_ENERGY, and Consensus score. The structures of the top 5 active ingredients for 2X0U and 3CQW are presented in Fig. 4.

Close

Fig. 4

Structures of the docking compounds for 2X0U and 3CQW.

Export to PPT

The CDOCKER score was reported as the negative value (i.e., -CDOCKER_ENERGY), which was calculated based on the internal ligand strain energy and receptor-ligand interaction energy (Rampogu and Rampogu Lemuel, 2016). In addition, the -CDOCKER_INTERATION_ENERGY represented the energy of non-bonding interactions between proteins and ligands. Higher -CDOCKER energy and -CDOCKER_INTERATION_ENERGY value indicated more favorable binding. From Table 4, the -CDOCKER ENERGY (35.58) and -CDOCKER_INTERATION_ENERGY (36.11) of Kaempferol2 were all higher than -CDOCKER ENERGY (14.02) and -CDOCKER_INTERATION_ENERGY (25.82) of the control compound (2-amino substituted benzothiazole scaffold) (Basse et al., 2010). In addition, the -CDOCKER ENERGY (49.69) and -CDOCKER_INTERATION_ENERGY (57.62) of Kaempferol3 were all higher than -CDOCKER ENERGY (21.23) and -CDOCKER_INTERATION_ENERGY (29.60) of the control compound (Spiroindoline 13j) (Lippa et al., 2008). Moreover, some experiments have shown that kaempferol is able to protect cells from virus-induced cell death and may be a promising antiviral drug for SARS-CoV-2 (Khan et al., 2021). Taken together, Kaempferol2 and Kaempferol3 can show better effects in the treatment of COVID-19, both in computer simulation screening and in vitro inhibition experiments.

Selecting a high activity 2-amino substituted benzothiazole scaffold as a control, it can be found that the active binding site of 2X0U is close to the amino acid residues Thr155, Pro151, Thr230, Val147, Cys220, Pro222, Leu145, and Pro223. Table 4 showed the interaction of amino acids in the ligand–protein for the top 5 docking active ingredients. By comparing the interacting amino acids in the top 5 docking compounds, it was found that Val147, Leu145, Pro222, Pro223, and Thr230 are the crucial amino acid residues for 2X0U. Selecting a high activity Spiroindoline 13j as a control, it can be found that the active binding site of 3CQW is close to the amino acid residues Ala230, Ala177, Leu156, Met227, Val164, Asp292, Thr291, Asn279, Thr291, Glu234, Met281and Pro223. By comparing the interacting amino acids in the top 5 docking active ingredients, it was found that Asp292 and Asn279 are the crucial amino acid residues for 3CQW.

The docking pose of Kaempferol2 was shown in Fig. 5 (A), which displayed 6 key residues around Kaempferol2, consisting of Pro151, Pro222, Pro223, Thr150, Val147, and Tyr220. Val147, Pro222, and Pro223 are the crucial residues for Kaempferol2 based on the result of Table 4. The 2D docking interaction plots of Kaempferol2 with the 2X0U active sites are shown in Fig. 5 (B). The docking pose of Kaempferol2 showed 3 Pi-Alkyl interactions binding with Pro151, Val147, and Pro223; 1 Pi-Antion interaction binding with Tyr220; and 1 conventional hydrogen bond interaction binding with Thr150. The docking pose of Kaempferol3 was shown in Fig. 6 (A), which displayed 4 key residues around Kaempferol3, consisting of Lys179, Phe161, The160, and LYS276. The 2D docking interaction plots of Kaempferol3 with the 3CQW active sites are shown in Fig. 6 (B). The docking pose of Kaempferol3 showed 2 salt bridges keying with Lys276 and Lys179; 2 conventional hydrogen bond interactions binding with Thr160 and Phe161; and 1 Pi-Pi T-shaped interaction binding with Phe161.

Close

Fig. 5

(A) The docking pose of Kaempferol2. (B) 2D docking Interaction plots of Kaempferol2 at the 2X0U active site.

Export to PPT

Close

Fig. 6

(A) The docking pose of Kaempferol3. (B) 2D docking Interaction plots of Kaempferol3 at the 3CQW active site.

Export to PPT

3.6 Global fluctuation of each system and binding free energy analysis

In order to explore the rationality of conformational samplings, the time evolutions of root-mean-square deviation (RMSD) for heavy atoms referenced to the corresponding energy-minimized initial structure were calculated on each system. As shown in Fig. 7, all systems reached equilibrium at about 60 ns. Thus, the last 40 ns trajectories were used to calculate the binding free energy.

Close

Fig. 7

The RMSD curves for the complexes.

Export to PPT

In our work, MMGBSA was applied to calculate the binding free energy for each complex. As listed in Table 5, the binding free energy of 2X0U-Kaempferol2 was −12.9 kcal mol⁻¹, and that of 3CQW-Kaempferol3 was –22.5 kcal mol⁻¹. These values are all less than zero, implying that the binding of both complexes is energetically feasible.