\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[greek,ngerman,polish,english]{babel}



\usepackage{float}








\begin{document}

\title{Exploring the Therapeutic Efficacy: Unveiling the Active Compounds of
Huashi Baidu Granules Against COVID-19}



\author[1]{Chuanxi Tian}%
\author[2]{Jinyue Zhao}%
\author[3]{Qian Wang}%
\author[4]{Keke Luo}%
\author[2]{Shuang Zhao}%
\author[2]{Jiarui Li}%
\author[1]{Li Wan}%
\author[1]{Kaile Ma}%
\author[5]{Chensi Yao}%
\author[4]{Yanyan Zhou}%
\author[1]{Min Li}%
\affil[1]{China Academy of Chinese Medical Sciences Guang'anmen Hospital}%
\affil[2]{Changchun University of Chinese Medicine}%
\affil[3]{Peking University School of Pharmaceutical Sciences}%
\affil[4]{China Academy of Chinese Medical Sciences Institute of Chinese Materia Medica}%
\affil[5]{Zhejiang University School of Medicine Second Affiliated Hospital Linping Hospital}%


\vspace{-1em}



  
  \date{June 10, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
\textbf{Background}: Huashi Baidu granule (HSBD), an approved herbal
formula for treating COVID-19, demonstrates safety and efficacy. Despite
its market approval, the detailed methodology and identification of its
active components remain unexplored, leaving its bioactive constituents
and action mechanisms unclear. \textbf{Methods}: This study investigated
the potential mechanisms of HSBD's active ingredient in treating
COVID-19. Our approach integrated various techniques, including the GEO
database, network pharmacology, surface plasmon resonance, molecular
docking and molecular dynamics simulations, to formulate a comprehensive
research strategy. \textbf{Results}: The UHPLC-QqQ-MS/MS method employed
for HSBD analysis proved stable, reliable, and reproducible. We
identified 25 principal components in HSBD, with seven compounds
detected in plasma, namely pogostone, p-hydroxybenzoic acid,
Paeoniflorin, Rhein, Emodin, ephedrine hydrochloride, and
pseudoephedrine hydrochloride. Protein-Protein Interaction (PPI) network
analysis identified MMP9 as a pivotal target, while KEGG pathway
analysis indicated that HSBD mediates its therapeutic effects on
COVID-19 primarily through the TNF-\selectlanguage{greek}α \selectlanguage{english}and PI3K-AKT pathways. Surface
plasmon resonance analysis revealed that Paeoniflorin and Rhein exert
their antiviral effects by interacting with RBD and ACE2, respectively.
In contrast, Emodin's antiviral mechanism predominantly involves binding
to MMP9. Molecular docking results indicated strong binding affinities
of Rhein and Paeoniflorin to the hACE2 protein, and high binding
affinities of Paeoniflorin, Rhein, and Emodin to the MMP9 protein, all
of which were corroborated by molecular dynamics simulations.
\textbf{Conclusion:} We investigated the methodology and identified the
active components of HSBD, focusing on those absorbed into the plasma,
to elucidate the effective material basis of HSBD in the treatment of
COVID-19, our research offer insightful exploration into its mechanisms
of action against COVID-19.%
\end{abstract}%



\sloppy


\textbf{1 Introduction}

COVID-19, caused by the SARS-CoV-2 virus, has been a global health
crisis since December 2019, with over 760 million cases and 6.9 million
fatalities as reported by the World Health Organization. In this
context, Huashi Baidu granule (HSBD), also known as Q-14, comprises 14
Chinese materia medica and has significantly contributed to COVID-19
treatment in China. Following its approval by the State Drug
Administration on March 2, 2021, HSBD has demonstrated considerable
efficacy in reducing clinical symptoms, improving inflammatory markers
and oxygenation levels, and enhancing viral clearance and lung recovery,
thus confirming its safety and effectiveness for severe COVID-19
cases\textsuperscript{1}. However, despite its clinical success, the
specifics of HSBD's active components and their mechanisms remain
inadequately explored.

This research aims to delineate the methodology and profile the 25
active compounds in HSBD, including the seven absorbed into plasma,
employing the the Ultra-High-Performance Liquid Chromatography-Triple
Quadrupole Mass Spectrometry (UHPLC-QqQ-MS/MS) technique for their
concurrent determination. The study adopts an integrated research
approach, combining the GEO database, network pharmacology, surface
plasmon resonance, molecular docking and molecular dynamics simulation,
to elucidate the underlying mechanisms of the active ingredients of HSBD
in combating COVID-19.

\textbf{2. Materials and methods}

\textbf{2.1 Materials and chemicals}

Ephedrine Hydrochloride, Pseudoephedrine Hydrochloride, Amygdalin,
Glycyrrhizin, Paeoniflorin, and Astragaloside were acquired from the
China Academy of Food and Drug Testing and Research. Additional
compounds such as Isorhamnetin-3-O-beta-D-Glucoside, Hyperoside, Sinapic
acid, Pogostone, Pachymic acid, Magnolol,
Calycosin-7-O-beta-D-glucoside, Honokiol, p-Hydroxybenzoic acid,
Quercetin-3-O-beta-D-glucose-7-O-beta-D-gentiobioside, Gentianic acid,
Glycyrrhetinic acid, Atractylon, and Hederagenin were sourced from
Beijing Saibaike Technology Co. AloeEmodin, Rhein, and Physcion were
procured from Chengdu Ruifenshi Biotechnology Co. Meanwhile,
Chrysophanic acid and Emodin were obtained from Shanghai Yuanye
Biotechnology Co., and CinobufaginIS was purchased from the China
National Institute for the Control of Pharmaceutical and Biological
Products. Glibenclamide IS was acquired from Shanghai Aladdin
Biochemical Technology Co.

HSBD batches were supplied by Guangdong E Fang Pharmaceutical Co.
High-quality reagents, including methanol, formic acid, acetonitrile
from Thermo Fisher Scientific (China) Co. Ltd, distilled water from
Watson's, and other analytically pure reagents, were utilized to ensure
the integrity of the research.

\textbf{2.2 Chromatographic Conditions}

Chromatographic conditions were specified using an ACQUITY UHPLC HSS T3
column (2.1 mm \selectlanguage{ngerman}× 100 mm, 1.8 \selectlanguage{greek}µ\selectlanguage{english}m). The mobile phase consisted of
acetonitrile-water with 0.1\% formic acid. The gradient elution profile
was as follows: from 0 to 3 minutes, 1\% acetonitrile; 3 to 5 minutes,
1\% to 30\% acetonitrile; 5 to 7 minutes, 30\% to 35\% acetonitrile; 7
to 12 minutes, 35\% to 95\% acetonitrile; 12 to 13 minutes, maintained
at 95\% acetonitrile; and 13 to 16 minutes, 95\% back to 1\%
acetonitrile. The flow rate was set at 0.25 mL/min, with a column
temperature of 30\selectlanguage{ngerman}°C and an injection volume of 2 \selectlanguage{greek}µ\selectlanguage{english}L.

For mass spectrometry, an electrospray ionization ion source (ESI) was
used in multiple reaction monitoring (MRM) mode with detection
capabilities in both positive and negative ion modes. The spray voltage
(ion spray) was set at 4500 V, the ion source temperature at 550\selectlanguage{ngerman}°C, and
the collision gas (CAD) flow rate at 9 psi (approximately 6.895 kPa).
The curtain gas (CUR) was set at 20 psi, and the nebulized gas (N2) at
14 psi. Specific conditions for ion pair detection are detailed in Table
S1.

\textbf{2.3 Standard solution preparation}

Accurately measure the required quantities of each standard and dissolve
them in chromatographic grade methanol to prepare individual standard
stock solutions. Subsequently, combine and dilute these stock solutions
sequentially to obtain a mixed standard solution. For the internal
standard preparation, progressively dilute the glibenclamide and
huazenin toxic base stock solutions to final concentrations of 1.5 \selectlanguage{greek}µ\selectlanguage{english}g/mL
and 2.5 \selectlanguage{greek}µ\selectlanguage{english}g/mL, respectively. Store all prepared solutions at -20\selectlanguage{ngerman}°C until
analysis.

\textbf{2.4 Solution Preparation}

\textbf{2.4.1 Test material}

Dissolve 0.5 g of the sample powder in 25 mL of 70\% methanol, followed
by ultrasonic extraction for 30 minutes and subsequent filtration.
Aliquot 980 \selectlanguage{greek}µ\selectlanguage{english}L of the resulting supernatant and enrich it with 10 \selectlanguage{greek}µ\selectlanguage{english}L of
glibenclamide internal standard solution and 10 \selectlanguage{greek}µ\selectlanguage{english}L of huazenin toxic
base internal standard solution. Thoroughly vortex the mixture to ensure
homogeneity, then transfer it into a liquid phase vial to facilitate
quantitative analysis.

\textbf{2.4.2 Plasma sample preparation}

Transfer 200 \selectlanguage{greek}μ\selectlanguage{english}L of plasma into a 1.5 mL centrifuge tube, then add 600 \selectlanguage{greek}μ\selectlanguage{english}L
of acetonitrile. Vortex mix the contents for 3 minutes, followed by
centrifugation at 12,000 rpm for 15 min at 4 . Carefully collect the
supernatant and place it into a liquid phase vial for subsequent
analysis.

\textbf{2.5 Differential gene acquisition in COVID-19 severe patients}

The GSE164805 dataset, sourced from the GEO database, comprises 15
samples, including 5 normal and 5 severe samples, which were selected
for this analysis. Differential analysis was conducted using the R
package ''limma,'' with thresholds for differential genes established at
\textbar{}logFC\textbar{} \textgreater{} 2 and p-value \textless{} 0.05.
The genes identified under these criteria were deemed to exhibit
significant differential expression and were considered potential drug
targets in severe COVID-19 patients.

\textbf{2.6 Target construction components absorbed into the plasma of
HSBD}

The Traditional Chinese Medicine Systems Pharmacology Database (TCMSP)
was queried for seven components: Pogostone, P-Hydroxybenzoic acid,
Paeoniflorin, Rhein, Emodin, Ephedrine Hydrochloride, and
Pseudoephedrine Hydrochloride. These components were screened to
identify those with an oral bioavailability (OB) of at least 30\% and a
drug-likeness (DL) value of 0.18 or higher. Subsequently, information on
the targets associated with these components was extracted. The protein
names of these targets were then converted to gene names utilizing the
UniProt database.

\textbf{2.7 Construction of protein-protein interaction networks and
identification of hub genes}

The STRING platform was utilized to analyze the protein-protein
interaction (PPI) network for the intersecting targets of the seven
components absorbed into the plasma from HSBD and associated with severe
COVID-19. An interaction threshold of 0.4 was established to exclude
isolated nodes, with other settings at default values. The resulting
data were then imported into Cytoscape 3.8.2 for comprehensive network
analysis, where nodes represent distinct targets, and edges illustrate
the interactions between these targets.

\textbf{2.8 Gene Ontology and Kyoto Gene Encyclopedia and Genome
Enrichment Analysis}

The intersecting targets of the seven HSBD components and severe
COVID-19 were inputted into the DAVID database, a web-based functional
annotation tool. Settings were adjusted to ''official gene symbol'' for
the select identifier and ''genelist'' for both list type and species,
specifying ''Homo sapiens'' as the species. This setup facilitated the
Gene Ontology (GO) enrichment analysis across three categories:
biological process, cellular component, and molecular function.
Additionally, a KEGG pathway enrichment analysis was conducted. The
results of these analyses were subsequently visualized using bubble
mapping on the hiplot platform.

\textbf{2.9 Surface plasmon resonance}

The binding affinities of Paeoniflorin, Emodin, and Rhein to hACE2,
SARS-CoV-2 S-RBD, and MMP9 were assessed using the Biacore 8K+
instrument (Cytiva, America). Proteins hACE2, SARS-CoV-2 S-RBD, and
MMP9, all sourced from Sino Biological, Beijing, China, were immobilized
on a CM5 sensor chip. The immobilization levels for these proteins were
approximately 10,000 response units (RU). A series of compound
concentrations, each containing 5\% DMSO, were sequentially introduced
into the channel to ascertain binding affinity. The dissociation
constants (K\textsubscript{D}) for the compounds were calculated by
fitting the data to a steady-state affinity model using the Biacore 8K+
Evaluation Software.

In the competitive inhibition assays, hACE2 or S-RBD proteins were
immobilized on the CM5 sensor chip via amine-coupling, with
immobilization levels around 3,000 RU. Paeoniflorin and Rhein were
initially injected for pre-incubation with hACE2 or S-RBD proteins,
respectively. Subsequent injections of 5 nM S-RBD or 30 nM hACE2 were
administered to evaluate their interaction with hACE2 or S-RBD in the
presence of Paeoniflorin or Rhein. The efficacy of inhibition was
determined by comparing response units with and without the active
compounds.

\textbf{2.10 Enzymatic assay}

The enzymatic activities of SARS-CoV-2 3CL\textsuperscript{pro},
PL\textsuperscript{pro}, and TMPRSS2 were evaluated using the respective
peptides Thr-Ser-Ala-Val-Leu-Gln-pNA, Z-Arg-Leu-Arg-Gly-Gly-AMC
(supplied by Bachem Bioscience), and Boc-Gln-Ala-Arg-AMC (provided by
CUSABIO). The reaction buffer for 3CL\textsuperscript{pro} contained 40
mM PBS (pH 7.3), 100 mM NaCl, 1 mM EDTA, and 0.1\% Triton X-100. For the
PL\textsuperscript{pro} assay, the buffer was composed of 50 mM HEPES
(pH 7.4), 100 mM NaCl, and 0.01\% Triton X-100 (v/v). The buffer for the
TMPRSS2 assay included 50 mM Tris-HCl (pH 8.0), 100 mM NaCl, and 0.01\%
Triton X-100 (v/v). Measurements of absorbance or fluorescence were
conducted using a Biotek Synergy Neo2 instrument.

\textbf{2.11 Molecular docking}

The molecular structure files for MMP9 (PDB codes: 8D35, 2OVX), human
angiotensin-converting enzyme 2 (hACE2) (PDB code: 3D0G), and the
SARS-CoV-2 spike receptor-binding domain (RBD) (PDB code: 7CH5) were
retrieved from the Protein Data Bank (PDB) and AlphaFold databases.
Subsequent preprocessing of these target proteins involved the removal
of water molecules and proto-ligands, utilizing PyMOL software version
2.3.0. Additionally, the small molecular structures of Paeoniflorin
(Compound CID: 442534), Rhein (Compound CID: 10168), and Emodin
(Compound CID: 3220) were acquired from the PubChem database. Following
the download, these small molecules were optimized for their
conformations through molecular mechanics calculations using Chem3D
software (version 2020), resulting in structures with minimized energy.

The prepared protein molecules were then hydrogenated using AutoDock
Tools version 1.5.6. Similarly, the small molecules, post-optimization,
were hydrogenated and their torsional bonds were identified. Docking
simulations between the target proteins and the small molecules were
performed using AutoDock Vina software version 1.2.0. The docking
employed a Lamarckian Genetic Algorithm (LGA) in a semi-flexible docking
mode, with the exhaustiveness parameter set to 128, to predict the
interaction conformations.

\textbf{2.12 Molecular dynamics simulation}

Molecular dynamics simulations were conducted using Gromacs software
version 2022.3. The small molecules underwent preprocessing with
AmberTools22 to incorporate the GAFF force field, while Gaussian 16W
facilitated their hydrogenation and the calculation of RESP potentials.
Simulation conditions were set at a constant temperature of 300 K and
atmospheric pressure (1 bar), employing the Amber99sb-ildn force field
and the Tip3p water model as the solvent. The system's total charge was
neutralized by adding an appropriate quantity of Na+ ions.

Initially, the molecular dynamics simulation system underwent energy
minimization using the steepest descent method. This was followed by
equilibration phases: an isothermal-isochoric (NVT) and an
isothermal-isobaric (NPT) equilibration, conducted over 100,000 steps
with a coupling constant of 0.1 ps and a duration of 100 ps,
respectively. Subsequently, an unrestrained molecular dynamics
simulation was executed over 5,000,000 steps with a time step of 2 fs,
totaling a simulation time of 100 ns.

Upon completing the simulations, trajectory analysis was performed using
the software's integrated tools to calculate the root mean square
deviation (RMSD), root mean square fluctuation (RMSF), and the radius of
gyration for each amino acid within the trajectories.

\textbf{3 Results}

\textbf{3.1 Methodological Examination}

\textbf{3.1.1 Linear relationship}

The internal standard method was employed to quantify 25 compounds in
HSBD samples. Following the protocol delineated in section 2.3, each
standard solution was serially diluted. For analysis, 980 \selectlanguage{greek}µ\selectlanguage{english}L of each
diluted solution was combined with 10 \selectlanguage{greek}µ\selectlanguage{english}L of glibenclamide and 10 \selectlanguage{greek}µ\selectlanguage{english}L of
huazenin toxicity internal standard solutions. The peak areas were
measured under conditions specified in section 2.2. The ratio of peak
area of each compound to that of the internal standard (Yi/Ys) was
calculated to conduct linear regression, facilitating the construction
of standard curves and derivation of regression equations and
correlation coefficients, which are tabulated in Table S2.

Upon plotting these standard curves, the 25 compounds exhibited strong
linear relationships, with correlation coefficients (r values) exceeding
0.99, indicating a high degree of linearity. These findings underscore
the effectiveness of the method in quantifying the compounds, as
reflected by the substantial correlation coefficients.

\textbf{3.1.2 Precision, stability and reproducibility study}

Pipette 960 \selectlanguage{greek}μ\selectlanguage{english}L of the standard mixed solution, then add 10 \selectlanguage{greek}µ\selectlanguage{english}L of
glibenclamide and 10 \selectlanguage{greek}µ\selectlanguage{english}L of huazufenin as internal standard solutions.
Following this, perform continuous injections six times within one day
under specified conditions (Section 2.2), and subsequently record the
peak areas for the 25 active components alongside the internal
standards. Utilize the standard curve to determine the concentrations,
revealing that the relative standard deviation (RSD) ranges from 0.94\%
to 14.24\%. This range suggests satisfactory precision of the
instrument.

For reproducibility assessment, precisely weigh six samples of HSBD,
prepare the test solutions according to Section 2.4, and measure the
peak areas following Section 2.2. The calculation of the RSD for the 25
active ingredients indicates a range of 1.49\% to 9.04\%, confirming the
method's reproducibility

To prepare the high-sensitivity bio-detection test solution, adhere to
the guidelines in Section 2.4. Conduct measurements at 0, 2, 4, 8, 12,
and 24 hours as per Section 2.2, recording the peak areas of the 25
active components and internal standards. Calculation of the
concentrations using the standard curve indicates an RSD range from
1.08\% to 5.97\%, demonstrating that the sample maintains stability over
a 24-hour period (refer to Table S3).

\textbf{3.1.3 Sample recovery test}

Accurately weigh six samples from the same batch of High Sensitivity
Bio-Detection, each weighing 0.5 g, and add a control solution
equivalent to the sample content. Prepare the test solution as per the
protocol outlined in Section 2.4.1. Following the chromatographic
conditions specified in Section 2.2, the average recovery rates of the
samples were determined to be between 93.22\% and 98.29\%. The relative
standard deviations (RSDs) ranged from 0.8\% to 3.8\%. These findings
suggest that the method is both accurate and reliable (refer to Table
S4).

\textbf{3.1.4 Determination of sample content components}

The UHPLC-QqQ-MS/MS technique was employed to quantify 25 primary active
compounds in the HSBD granule sample, utilizing Multiple Reaction
Monitoring (MRM) mode. A 0.5 g sample of HSBD was precisely weighed, and
the test solution was prepared following the protocol in Section 2.4.1.
Chromatographic conditions were set according to Section 2.2, and the
content of each compound was calculated using the corresponding standard
curve. The findings are presented in Table S5.

Analysis revealed that seven constituents, including Patchouli ketone,
p-hydroxybenzoic acid (grasshopper ketone), Paeoniflorin, rhodacanthic
acid, rhubarbine, ephedrine hydrochloride, and pseudoephedrine
hydrochloride, were detectable in the plasma(Table S6). Figure 1A
displays the MRM profiles for these 25 compounds and the internal
standard in the Resolving Dampness and Defeating Toxin Granules.

\textbf{3.2 COVID-19 differential gene acquisition in patients with
severe COVID-19}

The dataset GSE164805 was retrieved from the Gene Expression Omnibus
(GEO) database. The original expression profile dataset contained
multiple duplicate symbols. To address this, the median expression for
each gene was calculated, and the symbol with the higher median
expression was retained (see Appendix: GSE\_Matrix.txt). The R package
'limma' was utilized to conduct a differential expression analysis on
the refined dataset, comparing severe to normal conditions (see
Attachment: GSE\_DEGlist.txt). The criteria for differential gene
identification were set at an absolute log fold change
(\textbar{}logFC\textbar{}) greater than 2 and a p-value less than 0.05.
This analysis identified 2740 differentially expressed genes, with 1128
being up-regulated and 1612 down-regulated.

To provide a comprehensive overview of the distribution of these
differential genes, a volcano plot was generated, as shown in Figure 1B.
Additionally, the expression patterns of these genes are visualized in a
heatmap presented in Figure 1C.

\textbf{3.3 PPI network of targets of action of HSBD plasma-absorbed
components and severe COVID-19 intersecting targets}

From the TCMSP database, we identified 429 targets associated with the
seven components absorbed into the plasma from HSBD. These targets were
cross-referenced with genes related to severe COVID-19, totaling 2,740,
resulting in 36 common targets between severe COVID-19 and HSBD
plasma-absorbed components\textsuperscript{2} (Figure 2A).

Utilizing the STRING platform, we constructed a Protein-Protein
Interaction (PPI) network for the intersecting targets between severe
COVID-19 and the HSBD-absorbed components. This network distinguished
up-regulated major targets (outer circle) from down-regulated core
targets (inner circle). By leveraging data from up-regulated
differentially expressed genes (DEGs), we mapped out signaling networks
connecting these genes with their corresponding transcription factors.
This analysis highlighted key genes and underscored the potential
therapeutic implications of active components of HSBD in treating
COVID-19. Within this network, the top five down-regulated genes in
terms of network connectivity were MMP9, PTGS2, MAOB, CCNB1, and MAOA,
whereas the top five up-regulated genes were MAOA, BCL2, GSK3B, ABCB1,
and ABCG2 (Figure 2B). Notably, MMP9 emerged as a significant gene, as
indicated by its highest connectivity degree, suggesting its crucial
role in the anti-COVID-19 response according to the network analysis and
corroborative literature.

\textbf{3.4 GO analysis of core targets}

The analysis of down-regulated genes revealed involvement in various
biological processes (BPs) such as regulation of chloride transport,
drug transport across the blood-brain barrier, and transepithelial
transport. Molecular functions (MFs) were associated with structures
like myelin sheath and plasma membrane, and cellular components (CCs)
included ubiquitin protein ligase binding and xenobiotic-transporting
ATPase activity (Figure 2C).

For up-regulated genes, identified BPs include response to xenobiotic
stimulus, positive regulation of cell migration, and lung development.
MFs predominantly involve the integral component of the plasma membrane,
while CCs feature activities such as endopeptidase and monoamine oxidase
activity, highlighted in Figure 2D. The visualization employs a color
and size scheme where redder bubble colors indicate smaller P-values and
hence greater significance, while larger bubble sizes denote a higher
number of genes enriched in the corresponding Gene Ontology (GO)
entries.

\textbf{3.5 KEGG pathway enrichment analysis}

The KEGG pathway enrichment analysis of upregulated and downregulated
differentially expressed genes (DEGs) revealed distinct pathways
associated with each gene set. Specifically, the downregulated DEGs were
primarily involved in the PI3K-AKT and IL-17 signaling pathways (Figure
2E), while the upregulated DEGs predominantly influenced the TNF
signaling pathway (Figure 2F).

\textbf{4 Mechanism of Action of the Active Components of HSBD}

Due to market unavailability, ephedrine hydrochloride and
pseudoephedrine hydrochloride were excluded from the validation process.
P-hydroxybenzoic acid, commonly used as an antiseptic, and pogostone,
known for its antimicrobial, anti-inflammatory, antioxidant, and
insecticidal properties, were also not included in subsequent validation
experiments. Based on the literature and our related studies,
Paeoniflorin, Rhein, and Emodin, compounds known to enter the
bloodstream, were identified as potentially key active ingredients in
HSBD for impacting COVID-19. To elucidate the action mechanisms of these
active components against COVID-19, we focused on six antiviral targets:
human angiotensin-converting enzyme 2 (hACE2), SARS-CoV-2 spike
receptor-binding domain (RBD), SARS-CoV-2 3C-like protease
(3CL\textsuperscript{pro}), SARS-CoV-2 Papain-Like Protease
(PL\textsuperscript{pro}), human transmembrane protease serine 2
(TMPRSS2), and Matrix Metalloprotein 9 (MMP9).

We initially investigated the binding affinities of Paeoniflorin, Rhein,
and Emodin with hACE2, RBD, and MMP9 using Surface Plasmon Resonance
(SPR) technology. The results, shown in Figures 3A-3D, and S1, indicated
that Paeoniflorin, Rhein, and Emodin all bind to hACE2 and RBD with the
dissociation constant (K\textsubscript{D}) values for hACE2 being 1.9,
8.3, and 31.9 \selectlanguage{greek}µ\selectlanguage{english}M and for RBD being 3.4, 13, and 29.1 \selectlanguage{greek}µ\selectlanguage{english}M, respectively.
Subsequent competition assays were performed to further characterize the
binding activity, particularly for Paeoniflorin and Rhein. Notably, at
concentrations of 10 and 100 \selectlanguage{greek}µ\selectlanguage{english}M, Rhein significantly inhibited the
binding to both hACE2 and RBD, while Paeoniflorin exhibited a
comparatively weaker inhibition (Figures 3E and S2). This reduced
inhibitory capacity of Paeoniflorin could be attributed to its weaker
binding affinity to RBD. Furthermore, SPR assays demonstrated that
Paeoniflorin, Rhein, and Emodin could bind to MMP9 with
K\textsubscript{D} values of 21.7, 12.7, and 9.3 \selectlanguage{greek}µ\selectlanguage{english}M, respectively
(Figures 3F and S3). However, enzyme activity assays revealed that at 50
\selectlanguage{greek}µ\selectlanguage{english}M, none of these three compounds significantly inhibited the activities
of 3CL\textsuperscript{Pro}, PL\textsuperscript{Pro}, and TMPRSS2. These
findings suggest that the antiviral effects of Paeoniflorin and Rhein
are primarily mediated through their interactions with RBD and ACE2,
whereas Emodin's antiviral mechanism is mainly associated with its
binding to MMP9.

\textbf{5 Molecular docking and dynamics simulation}

To elucidate the binding sites of Paeoniflorin, Rhein, and Emodin,
molecular docking and molecular dynamics simulations were employed to
assess their interactions with RBD, hACE2, and MMP9. Molecular docking
results revealed that Rhein and Paeoniflorin have binding energies with
hACE2 at -8.673 kcal/mol and -7.891 kcal/mol respectively, and with RBD
at -7.079 kcal/mol and -6.727 kcal/mol respectively (Figures 4A-4B, S4).
Specifically, Rhein formed hydrogen bonds with hACE2 at residues GLN98,
ASP206, and LYS562. Additionally, its anthraquinone moiety engaged in
conjugated interactions with hACE2 at residues LEU95, GLU208, and
VAL209. Regarding its interaction with RBD, Rhein established hydrogen
bonds with TYR365, LEU390, and CYS391, and exhibited hydrophobic
interactions with ALA363 and CYS391. Emodin presented a binding energy
of -8.702 kcal/mol with MMP9, while Rhein and Paeoniflorin showed
binding energies of -7.869 kcal/mol and -7.542 kcal/mol,
respectively(Figures 4C and S5). Emodin established hydrogen bond
interactions with MMP9 at residues TYR48, PRO97, CYS99, and GLU402.
Additionally, it demonstrated hydrophobic interactions with residues
CYS99, VAL398, PRO421, and TYR423. Concurrently, Emodin engaged in a
Pi-lone pair interaction with CYS99, and a Pi-sulfur interaction is
observed at the same residue. These findings elucidate the multifaceted
interaction dynamics of Rhein with hACE2 and RBD, and Emodin with MMP9,
providing valuable insights into its potential inhibitory mechanisms at
a molecular level.The binding affinity of Rhein and Paeoniflorin to
hACE2, as well as their strong interactions with MMP9 alongside Emodin,
were further validated through root mean square deviation (RMSD) and
radius of gyration (Rg) analyses. These analyses indicated that the
binding of Paeoniflorin to hACE2 stabilized after 30 ns of simulation
time, and the binding of Rhein stabilized after 55 ns, suggesting stable
and reliable interactions (Figures 5A and 5B). The radius of gyration
results suggested that the binding of these molecules to hACE2 did not
significantly alter the overall structural compactness of the protein
(Figure 5C).

Similarly, the RMSD values for the interactions of Paeoniflorin, Rhein,
and Emodin with MMP9 were analyzed, showing stabilization times of 35
ns, 20 ns, and 52 ns, respectively (Figures S6A-S6C). The decreasing
trend in the radius of gyration for all three compounds upon binding to
MMP9 suggests a more compact overall protein folding (Figures S6D-S6F).
These findings provide valuable insights for further research on
protein-small molecule interactions and the development of therapeutic
agents.

\textbf{6 Discussion}

We employed liquid chromatography-mass spectrometry to validate our
methodology and quantify the active ingredients in HSBD, confirming the
stability, reliability, and reproducibility of the methods. Among the 25
active constituents of HSBD, seven, including Pogostone,
p-Hydroxybenzoic acid, Paeoniflorin, Rhein, Emodin, Ephedrine
Hydrochloride, and Pseudoephedrine Hydrochloride, were identified as
being absorbed into the bloodstream.

Traditional Chinese Medicine (TCM) syndromes and herbal treatments play
a pivotal role in the foundation and evolution of network
pharmacology\textsuperscript{3}. While previous studies on HSBD formula
in network pharmacology have focused on the collective targets of the 14
herbs within HSBD and gene acquisition for COVID-19 from
GeneCards\textsuperscript{4-5}, our study introduces a novel approach.
We identified the components of HSBD that are absorbed into the
bloodstream using UHPLC-QqQ-MS/MS. Moreover, we utilized differentially
expressed genes from severe COVID-19 patients compared to normal
subjects, derived from the GEO database, as disease targets. The dataset
GSE164805, related to Chinese patient data, demonstrated HSBD granules'
effectiveness in treating COVID-19. Our analysis becomes more meaningful
by intersecting the seven identified blood-absorbed HSBD components with
severe COVID-19 disease targets for protein-protein interaction (PPI)
network analysis. This intersection suggests that HSBD could modulate
the severity of the disease by downregulating genes like MMP9, PTGS2,
MAOB, CCNB1, MAOA, and upregulating genes such as MAOA, BCL2, GSK3B,
ABCB1, and ABCG2. Literature supports that MMP9, notably elevated in
COVID-19 patient lung tissues and integral in cytokine recruitment,
plays a crucial role in combating the virus\textsuperscript{6}. Its
upregulation is associated with the disease's severity, positioning MMP9
as a potential biomarker for COVID-19 severity assessment.

We conducted an in-depth analysis to uncover the potential mechanisms by
which the active components of HSBD affect COVID-19, utilizing Gene
Ontology (GO) enrichment and KEGG signaling pathway analyses. Our
findings indicate that the impact of HSBD on COVID-19 primarily engages
the TNF, PI3K-AKT, and IL-17 signaling pathways, all of which are
pivotal in the pathogenesis of the disease. TNF, an inflammatory
cytokine produced by macrophages/monocytes, plays a critical role in
mediating various intracellular signaling processes, including the
expression of growth factors, cytokines, and transcription factors. The
receptors for TNF are instrumental in modulating the immune response,
particularly affecting T-cell activation, proliferation, and
survival\textsuperscript{7}. Furthermore, TNF has been linked to the
induction of cytokine storms, a phenomenon closely associated with the
severity of COVID-19 due to the resulting immunopathological
damage\textsuperscript{8}. Specifically, the upregulation of IL-17A, and
potentially IL-17F, within these cytokine storms, contributes
significantly to the immunopathology of both COVID-19 and acute
respiratory distress syndrome\textsuperscript{9}. The PI3K-Akt pathway,
another crucial cellular signaling route, when activated, can mitigate
oxidative harm and inflammatory responses in pulmonary tissue by
inhibiting apoptosis, autophagy, or oxidative
stress\textsuperscript{10}. Thus, the bioactive compounds in HSBD might
potentiate TNF signaling, fostering appropriate immune reactions by
upregulating the TNF pathway. This upregulation could enhance TNF
expression or its receptor binding, deactivate negative regulators, or
activate positive ones, thereby offering therapeutic benefits.
Concurrently, HSBD may curb excessive inflammatory responses and prevent
tissue damage by downregulating the IL-17 and PI3K-AKT pathways,
ensuring a balance between the immune and inflammatory systems.

To elucidate the action mechanisms of the active constituents in HSBD
against COVID-19, we focused on several viral and human targets: hACE2,
RBD, 3CL\textsuperscript{pro}, PL\textsuperscript{pro}, TMPRSS2, and
MMP9, employing SPR for our investigations. Our findings revealed that
Paeoniflorin and Rhein primarily exhibited antiviral properties through
their binding to RBD and ACE2. Conversely, EmodinThe antiviral mechanism
of Emodin predominantly involved interaction with MMP9. Molecular
docking results further substantiated that Rhein and Paeoniflorin
exhibited strong binding affinities to the hACE2 protein, while
Paeoniflorin, Rhein, and Emodin all demonstrated significant binding
strengths to the MMP9 protein. These interactions were also corroborated
by molecular dynamics simulations, reinforcing the potential antiviral
efficacy of these compounds.

\textbf{Conclusion}

In this study, we investigated the methodology and identified the active
components of HSBD, focusing on those absorbed into the plasma, to
elucidate the effective material basis of HSBD in treating COVID-19.
This research aims to shed light on the therapeutic constituents of HSBD
and offer insightful exploration into its mechanisms of action against
COVID-19.

\textbf{Authors' Contributions}

Chuanxi Tian: Methodology, writing---original draft. Jinyue Zhao: Formal
analysis, data Curation. Qian Wang: Methodology, investigation, Data
Curation, Writing - Review \& Editing. Keke Luo: Investigation,
resources. Shuang Zhao: Writing - Review \& Editing. Jiarui Li: Writing
- Review \& Editing. Li Wan: Writing - Review \& Editing. Kaile Ma:
Writing - Review \& Editing. Chensi Yao: Project administration, Writing
- Review \& Editing. Yanyan Zhou: Investigation, resources. Min Li:
Supervision, project administration, funding acquisition.

\textbf{Acknowledgements}

This work was supported by National Key Research and Development Program
of China (2020YFE0205100) and the National Natural Science Foundation of
China (Grants 22107004).

\textbf{Conflict of interests}

The authors declare no competing interests.

\textbf{Data availability}

The data that support the findings of this study are available from the
corresponding author on reasonable request.

\textbf{Reference}

1. Wang, Y.; Lu, C.; Li, H.; Qi, W.; Ruan, L.; Bian, Y.; Shi, H.; Song,
H.; Tu, S.; Zhang, Y.; Bai, T.; Cao, R.; Hong, K.; Li, H.; Liu, L.; Lu,
S.; Rong, N.; Liu, Y.; Fang, J.; Shi, J.; Yang, W.; Zhao, B.; Yang, Y.;
Zhao, Y.; Li, S.; Fan, T.; Rong, P.; Huang, L., Efficacy and safety
assessment of severe COVID-19 patients with Chinese medicine: A
retrospective case series study at early stage of the COVID-19 epidemic
in Wuhan, China. \emph{Journal of ethnopharmacology}
\textbf{2021,}\emph{277} , 113888.

2. Bardou, P.; Mariette, J.; Escudi\selectlanguage{ngerman}é, F.; Djemiel, C.; Klopp, C., jvenn:
an interactive Venn diagram viewer. \emph{BMC
bioinformatics}\textbf{2014,} \emph{15} (1), 293.

3. Wang, X.; Wang, Z. Y.; Zheng, J. H.; Li, S., TCM network
pharmacology: A new trend towards combining computational, experimental
and clinical approaches. \emph{Chinese journal of natural
medicines}\textbf{2021,} \emph{19} (1), 1-11.

4. Zhu, Y. W.; Yan, X. F.; Ye, T. J.; Hu, J.; Wang, X. L.; Qiu, F. J.;
Liu, C. H.; Hu, X. D., Analyzing the potential therapeutic mechanism of
Huashi Baidu Decoction on severe COVID-19 through integrating network
pharmacological methods. \emph{Journal of traditional and complementary
medicine} \textbf{2021,} \emph{11} (2), 180-187.

5. Zhang, F.; Guo, F.; Zhang, Y.; Xu, H.; Liu, Y.; Lin, L.; Li, H.;
Yang, H.; Huang, L., Huashibaidu formula attenuates sepsis-induced acute
lung injury via suppressing cytokine storm: Implications for treatment
of COVID-19. \emph{Phytomedicine : international journal of phytotherapy
and phytopharmacology} \textbf{2023,} \emph{109} , 154549.

6. Hazra, S.; Chaudhuri, A. G.; Tiwary, B. K.; Chakrabarti, N., Matrix
metallopeptidase 9 as a host protein target of chloroquine and melatonin
for immunoregulation in COVID-19: A network-based
meta-analysis.\emph{Life sciences} \textbf{2020,} \emph{257} , 118096.

7. Ward-Kavanagh, L. K.; Lin, W. W.; \selectlanguage{polish}Š\selectlanguage{english}ed\selectlanguage{ngerman}ý, J. R.; Ware, C. F., The TNF
Receptor Superfamily in Co-stimulating and Co-inhibitory
Responses.\emph{Immunity} \textbf{2016,} \emph{44} (5), 1005-19.

8. Ramasamy, S.; Subbian, S., Critical Determinants of Cytokine Storm
and Type I Interferon Response in COVID-19 Pathogenesis. \emph{Clinical
microbiology reviews} \textbf{2021,} \emph{34} (3).

9. Shibabaw, T., Inflammatory Cytokine: IL-17A Signaling Pathway in
Patients Present with COVID-19 and Current Treatment
Strategy.\emph{Journal of inflammation research} \textbf{2020,}
\emph{13} , 673-680.

10. Matthay, M. A.; Zimmerman, G. A., Acute lung injury and the acute
respiratory distress syndrome: four decades of inquiry into pathogenesis
and rational management. \emph{American journal of respiratory cell and
molecular biology} \textbf{2005,} \emph{33} (4), 319-27.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\textbf{Figure 1. A. MRM profiles of 25 compounds and internal standards
in HSBD.} 1. Quercetin-3-O-beta-D-glucose-7-O-beta-D-gentiobioside;
2.Isorhamnetin-3-O-beta-D-Glucoside; 3.Pseudoephedrine Hydrochloride;
4.Ephedrine Hydrochloride; 5.Amygdalin; 6. p-Hydroxybenzoic acid;
7.Gentianic acid; 8.Calycosin-7-O-beta-D-glucoside; 9.Hyperoside;
10.Glycyrrhizin; 11.Sinapic acid; 12.Astragaloside; 13.Glycyrrhetinic
acid; 14.AloeEmodin; 15.Rhein; 16.Paeoniflorin; 17.Magnolol; 18.Emodin;
19.Honokiol; 20.Pogostone; 21.Hederagenin; 22.Chrysophanic acid;
23.Physcion; 24.Pachymic acid; 25.Atractylon; S1. Cinobufagin; S2.
Glibenclamide. \textbf{B. Differential gene volcano map for COVID-19
severe patients.} In this plot, genes are denoted as points, with red
indicating up-regulation, blue indicating down-regulation, and gray
highlighting no significant difference in expression between severe and
normal samples.\textbf{C. Differential Gene Heat Map for COVID-19 severe
patients.} Each small square represents each gene, and its color
indicates the expression level of the gene, the higher the expression
level, the darker the color (red is high expression, blue is low
expression). The first row indicates the sample grouping, gray indicates
Normal samples and dark blue indicates Severe samples. Each row
indicates the expression of each gene in different samples, and each
column indicates the expression of all differential genes in each
sample.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Figure 2. Network pharmacological analysis was conducted to
investigate the plasma-absorbed components of HSBD for treating severe
COVID-19.} A. Identification of shared targets between severe COVID-19
and HSBD plasma-absorbed components. B. PPI network analysis of the
action targets for HSBD plasma-absorbed components and their
intersection with severe COVID-19 targets. C. GO annotation for
down-regulated DEGs. D. GO annotation specifically for up-regulated
DEGs. E. KEGG pathway enrichment analysis for own-regulated DEGs. F.
KEGG pathway enrichment analysis for up-regulated DEGs.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\textbf{Figure 3. Affinity activity of Paeoniflorin, Rhein, Emodin with
hACE2, RBD and MMP9 by SPR technique.} Both Paeoniflorin and Rhein bound
to hACE2 and RBD, with K\textsubscript{D} values of 1.9 (A) and 8.3 (C)
\selectlanguage{greek}µ\selectlanguage{english}M for hACE2 and 3.4 (B) and 13 (D) \selectlanguage{greek}µ\selectlanguage{english}M for RBD. \textbf{A.} Competition
experiments characterize the binding epitopes of Rhein. It can be seen
that Rhein could significantly inhibit the binding of hACE2 and
RBD.\textbf{B.} Emodin bound to MMP9 with K\textsubscript{D} value of
9.3 \selectlanguage{greek}µ\selectlanguage{english}M.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

\textbf{Figure4. The binding interactions of Rehin or Emodin to ACE2,
RBD and MMP9.} A. Molecular docking studies revealed the interaction
modes of Rhein with hACE2. B. The docking results of Rhein to RBD. C.
Molecular docking provided insights into the binding mechanism of Emodin
with MMP9 proteins. The dark green dashed lines denote hydrogen bond
interactions, whereas the light green, pink, and purple dashed lines
depict interactions related to Pi bonds, encompassing both polar and
nonpolar interactions.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

\textbf{Figure 5. Molecular dynamics (MD) analysis was performed to
assess the interaction between Rhein, Paeoniflorin and hACE2, focusing
on root mean square deviation (RMSD) and gyration.} A. Molecular
dynamics simulations of Rhein and hACE2 (RMSD Merge). B. Molecular
dynamics simulations of Paeoniflorin and hACE2 (RMSD Merge). C.
Molecular dynamics simulations of Rhein, Peoniflorin and hACE2 (Gyrate).
Gyrate is an indicator of the overall structural compactness of the
protein molecule.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-1/Figure-1}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-2/Figure-2}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-3/Figure-3}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-4/Figure-4}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-5/Figure-5}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-S1/Figure-S1}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-S2/Figure-S2}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-S3/Figure-S3}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-S4/Figure-S4}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-S5/Figure-S5}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-S6/Figure-S6}

\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\end{document}

