COVID-19 Research Paper Volume 13, Issue 21 pp 23913—23935
New tale on LianHuaQingWen: IL6R/IL6/IL6ST complex is a potential target for COVID-19 treatment
- 1 Department of Pharmacology, College of Basic Medical Sciences, Jilin University, Changchun, Jilin Province 130021, People’s Republic of China
How to Cite
LianHuaQingWen (LHQW) improves clinical symptoms and alleviates the severity of COVID-19, but the mechanism is unclear. This study aimed to investigate the potential molecular targets and mechanisms of LHQW in treating COVID-19 using a network pharmacology-based approach and molecular docking analysis. The main active ingredients, therapeutic targets of LHQW, and the pathogenic targets of COVID-19 were screened using the TCMSP, UniProt, STRING, and GeneCards databases. According to the “Drug-Ingredients-Targets-Disease” network, Interleukin 6 (IL6) was identified as the core target, and quercetin, luteolin, and wogonin as the active ingredients of LHQW associated with IL6. The response to lipopolysaccharide was the most significant biological process identified by gene ontology enrichment analysis, and AGE-RAGE signaling pathway activation was prominent based on the interaction between LHQW and COVID-19. Protein-protein docking analysis showed that IL6 receptor (IL6R)/IL6/IL6 receptor subunit beta (IL6ST) and Spike protein were mainly bound via conventional hydrogen bonds. Furthermore, protein-small molecule docking showed that all three active ingredients could bind stably in the binding model of IL6R/IL6 and IL6ST. Our findings suggest that LHQW may inhibit the lipopolysaccharide-mediated inflammatory response and regulate the AGE-RAGE signaling pathway through IL6. In addition, the N-terminal domain of the S protein of COVID-19 has a good binding activity to IL6ST, and quercetin and wogonin in LHQW may affect IL6ST-mediated IL6 signal transduction and a large number of signaling pathways downstream to other cytokines by directly affecting protein-protein interaction. These findings suggest the potential molecular mechanism by which LHQW inhibits COVID-19 through the regulation of IL6R/IL6/IL6ST.
IntroductionCOVID-19 has spread across the globe since its outbreak at the end of 2019 due to its high susceptibility and transmissibility . Patients often display symptoms of fever, dry cough, fatigue, and even lung inflammation and infiltration, as well as systemic cytokine storms . All countries worldwide are trying to develop drugs to treat this disease. In the absence of a specific cure for COVID-19, traditional Chinese medicine (TCM) plays an important role in inhibiting inflammation, preventing the development of disease, and shortening the treatment cycle because it has multiple components and targets, broad-spectrum antibacterial and antiviral activity, and can improve immunity. LianHuaQingWen (LHQW), a type of TCM used to treat respiratory diseases that are identified during public health events, was developed based on the experience of well-known Chinese doctors in Han, Ming, and Qing dynasties for treating cold, influenza, and other exogenous febrile diseases combined with modern pharmacological research. According to the Pharmacopeia of the People’s Republic of China 2015 Edition , LHQW is a compounded medication (Supplementary Table 1) composed of 13 components, including primary, secondary, and adjuvant components, and the toxicity of TCM is generally controllable. It can clear away the “lung heat,” inhibit bacterial infection, enhance immunity, reduce inflammation and fever, relieve cough, reduce sputum, reduce high fever, and alleviate aversion to cold, muscle pain, headache, sore throat, and other symptoms of heat toxicity. LHQW has a history of over 2,000 years of being used as a TCM medication during epidemic diseases. Previous studies have shown that LHQW had a broad-spectrum of antiviral effects targeting many viruses, including SARS coronavirus (SARS-CoV), Middle East respiratory syndrome coronavirus (MERS-CoV), H1N1, H3N2, and H7N9. It also exhibits antibacterial, anti-inflammatory, antipyretic, cough-relieving, expectorant, immunity-regulating, and other systematically intervening properties [4–5]. Since the COVID-19 outbreak, LHQW has been listed as the recommended drug during the observation period of COVID-19 in more than 20 national and local COVID-19 diagnostic and treatment plans . Owing to the participation of the Chinese medical team in local COVID-19 treatments in Italy and other countries, LHQW has obtained registration approval in six countries, including Brazil, Romania, Thailand, Ecuador, Singapore, and China, as well as Hong Kong and Macao SAR. On April 14, 2020, during the Press Conference of the Joint Prevention and Control Mechanism of the State Council, it was announced that LHQW could significantly relieve cough, fever, and fatigue caused by COVID-19 and could effectively reduce the rate of mild-to-severe symptoms; thus, it could be used for the routine treatment of COVID-19. In addition, with the approval of the National Medical Products Administration (NMPA), “In the routine treatment of COVID-19, it can be used for fever, cough, and fatigue caused by the mild and common type” was added in “Functions and indications” in the drug instructions, and “The course of treatment for mild and common COVID-19 is 7–10 d” was added in “Usage and dosage” . A recent study by Hu et al.  showed that the recovery rate for COVID-19 patients treated with LHQW was significantly higher than that in the control group because it improved the clinical symptoms and reduced the severity of the disease. It is suggested that the drug has clinical significance for the treatment of COVID-19, but the mechanism of LHQW in treating COVID-19 has not yet been identified. TCM prescriptions are known for being multi-component and multi-target, and involved in systemic regulation. Therefore, network pharmacology is an effective method for studying the mechanisms of action of these prescriptions. Yan et al.  and Wang et al.  have used network pharmacology and molecular docking technology to show that LHQW capsule could provide improved efficacy as a cure for COVID-19 through its multi-herb, multi-target, multi-signaling pathway and multi-biological function. The Spike (S) protein, including the S1 and S2 subunits, of COVID-19 is the key protein that mediates the entry of the virus into the host cell. The S protein promotes the fusion of the viral envelope with the target cell membrane or endosome membrane that allows the viral nucleic acid to enter the cell [10–11]. However, it is unclear whether LHQW can interfere with this process. This study aimed to investigate the mechanism of LHQW protective action during COVID-19 infection using network pharmacology-based approach and molecular docking analysis.
Results of network pharmacology analysis of LHQW protective action during COVID-19 infection
LHQW ingredients-targets network and protein-protein interaction network including COVID-19 pathogenic targets, LHQW therapeutic targets, junction targets of COVID-19, and LHQWGeneCards, TCMSP, and UniProt databases were used to screen 277 COVID-19 pathogenic targets. In addition, a network of 80 LHQW active ingredients and 234 therapeutic targets was constructed after screening (Supplementary Figure 1A and Supplementary Table 2). To analyze the association between targets, 277 COVID-19 pathogenic targets and 234 therapeutic targets were imported into the STRING database respectively. Protein-protein interaction (PPI) networks of 266 COVID-19 pathogenic targets and 233 LHQW therapeutic targets were obtained after removing duplicate values separately. Cytoscape was used to optimize the two networks and compute the values of degrees of network connections. Based on the maximum degree value, TNF and AKT1 were respectively identified as the core targets in the network of COVID-19 pathogenic targets (Supplementary Figure 2) and the network of LHQW therapeutic targets (Supplementary Figure 1B). In addition, secondary and tertiary targets were identified based on the values of their degree (in decreasing order). The darkness of the color and size of the circle are positively correlated with the role of the target in the network. To investigate the reflective relationship between COVID-19 and LHQW, the potential targets of COVID-19 were mapped with those of LHQW to obtain 45 junction targets. Using these targets, the PPI network of LHQW in treating COVID-19 was constructed using STRING and Cytoscape. Interleukin 6 (IL6) was detected as the core target according to a degree ≥41, and seven secondary targets, including TNF, MAPK8, CXCL8, TP53, MAPK3, IL10, and CASP3, were detected based on 38 ≤ degree <41. In terms of 33 ≤ degree < 38, 12 tertiary targets were further identified, including MAPK1, CCL2, IL1B, IFNG, PTGS2, ALB, IL4, ICAM1, FOS, IL2, MAPK14, and RELA. (See Figure 1).
Gene-Concept network of gene ontology enrichment analysis and Kyoto Encyclopedia of Genes and Genomes pathway analysis of junction targets of COVID-19 and LHQWGene ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis identified 2,274 biological processes (P ≤ 0.05), 50 cellular components (P ≤ 0.05), 137 molecular functions (P ≤ 0.05), and 152 KEGG pathways (P ≤ 0.05) in which junction targets of COVID-19 and LHQW were involved, separately. In terms of P value (in increasing order), the five most significant biological processes (Figure 2A), cellular components (Figure 2B), molecular functions (Figure 2C), and KEGG pathways (Figure 3) were screened to build the Gene-Concept network respectively.
Figure 2. Gene-Concept network of GO enrichment analysis of junction targets of COVID-19 and LHQW. (A) Biological processes. (B) Cellular components. (C) Molecular functions.
Results of core target IL6 for LHQW in treating COVID-19
IL6 was involved in the most significant GO enrichment and KEGG pathwaysThe analysis of IL6 as the core target of LHQW in treating COVID-19 showed that there were five biological processes that included IL6 in the top five biological processes in which the junction targets were involved (P ≤ 0.05). Response to lipopolysaccharide (LPS) was identified as the biological process in which IL6 was most actively involved (P ≤ 2.23 × 10−25). There were no cellular components with IL6 involvement among the first five cellular components in which the junction targets of COVID-19 and LHQW were included. There were three molecular functions with IL6 involvement among the first five molecular functions in which the junction targets of COVID-19 and LHQW were involved (P ≤ 0.05). According to P ≤ 4.15 × 10−17, the most significant molecular function with IL6 involvement was cytokine receptor binding. There were 5 KEGG pathways with IL6 involvement in the first 5 KEGG pathways in which the junction targets of COVID-19 and LHQW were involved (P ≤ 0.05). IL6 involvement was most active in the AGE-RAGE signaling pathway (P ≤ 2.36 × 10−23) (Figure 4A).
AGE-RAGE signaling pathway in diabetic complicationsThe AGE-RAGE signaling pathway in diabetic complications contained 19 targets. In addition to the core target IL6 (red), there were five secondary targets (pink), including TNF, CXCL8, MAPK8, MAPK3, and CASP3, and six tertiary targets (seashell), including IL1B, CCL2, MAPK1, ICAM1, MAPK14, and RELA, as well as seven quaternary targets (seashell), including STAT1, NOS3, SERPINE1, IL1A, BAX, BCL2, and PRKCA, indicating that the effect of LHQW on the prevention and control of COVID-19 was closely related to the regulation of various targets (Supplementary Figure 3).
Effective active ingredients of IL6The results of this study suggest that the effective active ingredients, including luteolin (Lonicerae Japonicae Flos, Ephedra Herba, Forsythiae Fructus), quercetin (Lonicerae Japonicae Flos, Ephedra Herba, Forsythiae Fructus, licorice, Houttuyniae Herba, Pogostemon Cablin (Blanco) Benth.), and wogonin (Forsythiae Fructus) are all connected with the core target IL6, indicating that the effect of LHQW on the prevention and treatment of COVID-19 is closely related to these three active ingredients (Figure 4B).
Verification of molecular docking for LHQW in treating COVID-19 by the results of network pharmacology
Protein-protein docking analysis of IL6R/IL6/IL6ST and COVID-19-Spike interactionsProtein-protein docking between COVID-19-Spike and IL6 receptor (IL6R)/IL6/IL6 receptor subunit beta (IL6ST) was simulated to verify whether IL6, the core target identified through our network pharmacology-based approach, was a pathogenic target. The results showed that IL6R/IL6/IL6ST (PDB ID: 1P9M) had high level of CDOCKER interaction energy with the Spike protein of COVID-19 (PDB ID: 6VXX). The maximum ZDOCK score was 1955.699.
The binding model of IL6R/IL6/IL6ST and SpikeTo simulate the process by which COVID-19 attaches to the host cell receptor, Pymol was used to analyze the binding model of IL6R/IL6/IL6ST and Spike. This analysis revealed that conventional hydrogen bonds occurred between the IL6R/IL6/IL6ST amino acid residues Arg276 (IL6ST), Trp291 (IL6ST), Trp288(IL6ST), Lys234 (IL6ST), Trp287 (IL6ST), Ser111 (IL6ST), Cys112 (IL6ST), Asn109 (IL6ST), and the Spike amino acid residues Asp198, Asp198, Thr167, Cys166, Tyr170, Tyr170, Tyr170, and His207, respectively (Figure 5A).
The binding model between IL6R/IL6 and IL6STDue to its biological characteristics, IL6R binds to IL6 and then binds to IL6ST with high affinity, thus inducing signal transduction process. This study analyzed the binding model of IL6R/IL6 and IL6ST using Pymol. The results indicated that IL6ST amino acid residues Ser229, Tyr168, and His145 formed conventional hydrogen bonds with IL6R/IL6 amino acid residues Asp34 (IL6), Gln28 (IL6), and Gln124 (IL6), respectively (Figure 6A).
Protein-small molecule docking analysisThe activity of the three active ingredients obtained through our network pharmacology-based approach was validated through protein-small molecule docking analysis.
Protein-small molecule docking analysis based on the binding model of IL6R/IL6 and IL6STProtein-small molecule docking, based on the binding model of IL6R/IL6 and IL6ST, was performed to determine if the three active ingredients affected binding between IL6R/IL6 and IL6ST. We found that the CDOCKER interaction energy values of quercetin, luteolin, and wogonin docking with the interaction center of IL6R/IL6 and IL6ST (1P9M-6VXX complex) were −26.1801 kcal mol−1, −24.3014 kcal mol−1, and −24.9362 kcal mol−1, respectively. The binding model analyzed by Pymol showed that there was a conventional hydrogen bond between quercetin and Gln124 (IL6R/IL6) of the 1P9M-6VXX complex, and conventional hydrogen bond and Pi-Pi stacking could be found between quercetin and His145 (IL6ST) of the 1P9M-6VXX complex. It should also be noted that Gln124 (IL6R/IL6) and His145 (IL6ST) are the key amino acids required for stable binding between IL6R/IL6 and IL6ST (Figure 6B–6D).
Protein-small molecule docking analysis based on the binding model of IL6R/IL6/IL6ST and SpikeTo investigate whether the three active ingredients affected COVID-19 attachment to the host cell receptors, protein-small molecule docking analysis was used based on the binding model of IL6R/IL6/IL6ST and Spike. This analysis showed that the CDOCKER interaction energy levels of the 1P9M-6VXX complex with quercetin, luteolin, and wogonin based on this binding model were −34.5920 kcal mol−1, −32.8415 kcal mol−1, and −30.5083 kcal mol−1, respectively. The analysis results of the binding model showed that although quercetin had the highest binding energy, it did not directly interact with the key amino acids that formed bonds between IL6R/IL6/IL6ST and Spike. In contrast, Van der Waals interaction between wogonin and the key amino acids Trp291 (IL6ST) and Asp198 (Spike) which formed stable bonds between IL6R/IL6/IL6ST and Spike. (See Figure 5B–5D)
DiscussionBased on previous studies, the present study improved the analysis method of network pharmacology and molecular docking, and selected the core target in the network based on the maximum degree value, which indicated the role of nodes in network pharmacology. The clusterProfiler package included in the R language was used to obtain the latest Gene Ontology and KEGG pathway for the core target to predict the mechanism of the inhibitory effect of LHQW on the COVID-19 inflammatory stress response. The core target and its receptor complex obtained by network pharmacology were first attached to the S protein by protein-protein docking, and then to the active ingredients connected with the core target by protein-small molecule docking, in order to simulate the process of drug intervention after COVID-19 invades the host cell receptor. These analyses validate the results of network pharmacology and provide a molecular basis for investigating the mechanism of the inhibition of the COVID-19 inflammatory stress response by LHQW. COVID-19 infection causes a type of lung collateral disease that can be treated using TCM. Its clinical and pathological features include immune system disorders, damage to the deep airway and alveoli caused by an inflammatory reaction, and severe pulmonary congestion. According to the literature [12–14], patients with COVID-19 show symptoms of an excessive inflammatory response, especially in severe cases. Excess release of cytokines and chemokines (IL6, TNF-α, IL2, IL7, CCL2, CXCL10, MCP-1) cause cytokine storm (CS) in these patients. As a result, there is a diffuse damage caused to the lung endothelium and alveolar epithelial cells, leading to acute respiratory distress syndrome (ARDS) and multiple organ dysfunction syndrome (MODS). In the advanced stage of viral infection, CS is considered one of the most important factors of the death in COVID-19 patients. CS, also known as Cytokine Release Syndrome, is an overactive immune response caused by an external stimuli, as a result of the rapid release of cytokines such as Interleukin, Interferon (IFN), TNF, and Colony-Stimulating Factor (CSF) in large quantities [15–17]. By analyzing 99 confirmed cases of COVID-19, Chen et al. proposed that the virus could infect other cells through the respiratory mucosa and induce a CS, causing severe immune injury to the lungs and other organs . In our study, we found that the core pathogenesis target of COVID-19 was TNF, while the secondary pathogenesis targets included IL6, IL10, and IL2 cytokines. The analysis results were consistent with the above clinical findings and indicated that inflammatory cytokines played crucial roles in COVID-19 infections. Therefore, inhibition of CS is critical for the treatment of COVID-19. Moreover, we identified the core target IL6 and other secondary targets, such as MAPK8, CXCL8, TP53, MAPK3, IL10, and CASP3, by mapping 45 junction targets, although the core therapeutic target of LHQW AKT1 did not map with the core target TNF of COVID-19. This conclusion is consistent with recent clinical findings. Recent clinical studies have shown a significant increase in the levels of IL6, TNF, and many other inflammatory factors, especially in patients with severe pneumonia. Disease severity was positively correlated with IL6 levels. Once the disease is under control, IL6 levels also decrease rapidly. Notably, if IL6 levels do not decrease rapidly after treatment, the patient prognosis is often poor [12, 18, 19]. Therefore, early detection of IL6 may provide an early indication of severe infection. In addition, IL6 has been identified as a clinical early warning indicator for severe and critically ill patients in the diagnosis and treatment protocol for COVID-19 (Trial Version 8) . Therefore, it is suggested that LHQW may regulate CS by inhibiting IL6, controlling the immunological stress caused by COVID-19, and relieving inflammatory symptoms caused by pulmonary infection. The biological process involving response to LPS was highlighted in GO enrichment analysis. LPS can cause damage to alveolar epithelial cells and capillary endothelial cells and then change the intercellular space and permeability, causing inflammatory cells to release IL6. It has been reported that IL6 levels in the LPS group were significantly higher than those in the control group (P < 0.01), and LHQW could reduce IL6 levels in mouse lung tissue infected by the influenza virus, suggesting that LHQW may alleviate lung injury by inhibiting LPS and decreasing the level of IL6 . The AGE-RAGE signaling pathway based on IL6 was significantly involved in the mechanisms of LHQW treatment of COVID-19 in KEGG pathway analysis. The AGE-RAGE signaling pathway is often activated in diabetic complications, which mediate the downstream activation of a large number of inflammatory cytokines. It has been reported that advanced glycation end product (AGE)/receptor AGE (RAGE)/nuclear transcription factor-κB (NF-κB) signaling pathway played an important role in oxidative stress and excessive inflammatory responses. RAGE is expressed at a low level under normal or physiological conditions, but AGE levels increase during inflammation or trauma and further lead to the upregulation of RAGE expression. Overexpression of RAGE induces the expression and release of a large number of proinflammatory cytokines, such as adhesion molecules, growth factors, IL6, IL8, and TNF-α, subsequently causing tissue injury . A recent study showed that diabetic complications and COVID-19 infection could both activate the AGE-RAGE signaling pathway and promote the progression of COVID-19 . Therefore, LHQW may have a good therapeutic effect on COVID-19 patients with a history of diabetes mellitus. It has been reported that the invasion of COVID-19 into human cells may be due to ACE2 binding to the receptor-binding domain (RBD) region of the S protein, resulting in membrane fusion and entry of the virus [22, 23]. Single-cell sequencing data, however, indicated that the overall ACE2 expression level was low in various human tissues, especially in pulmonary and bronchial tissues [24, 25]. In addition, recent studies have revealed several neutralizing antibodies, such as mAbs (4A8, COVA1-03, COVA1-21) that bound to the S protein, but the binding site was not in RBD [26, 27]. These results indicate the possibility of other receptors or co-receptors binding to the different domains of the S protein to promote the COVID-19 entry into the cell. According to our protein-protein docking results, IL6R/IL6/IL6ST and Spike were bound stably via conventional hydrogen bonds between the N-terminal domain (NTD) region of the Spike protein S1 subunit and IL6ST. NTD in the S1 subunit plays an important role in recognizing and binding receptors on the surface of the target cell membrane [10, 11]. Therefore, IL6ST is involved in the process by which COVID-19 enters the cell. Interleukin 6 (IL6, the core target) is a multifunctional cytokine that can carry out signal transduction and plays its biological function only by forming the IL6/IL6R/IL6ST hexameric complex with its receptor. However, under physiological conditions, IL6ST cannot bind to IL6 directly. Only after IL6R binds to IL6 it can bind to IL6ST with high affinity to activate downstream signaling. Thus, IL6 expression can be regulated at the receptor level. Any measure and methods that can affect the binding of IL6 and IL6R or the signal transduction mediated by the binding of IL6/IL6R and IL6ST may be applied in the treatment of IL6-related diseases. Targeting the IL6 receptor to develop and screen drugs is a potential approach for treating IL6-related diseases. Although some artificially synthesized antagonists of IL6 and IL6R have been developed, there are only a few reports on the screening of antagonists from natural products . Conversely, IL6ST, also known as membrane glycoprotein 130, is a shared signal-transducing receptor for a family of more than 10 different four-helix-bundle cytokines, including IL6, leukemia inhibitory factor, ciliary neurotrophic factor (CNTF), oncostatin-M (OSM), and others [29–31]. Therefore, because of the biological characteristics of IL6R/IL6/IL6ST and Spike, this study designed a protein-small molecule docking analysis based on two different binding models (IL6R/IL6 and IL6ST; IL6R/IL6/IL6ST and Spike). The results showed that all three active ingredients (luteolin, quercetin, and wogonin) had perfect CDOCKER interaction energy with the 1P9M-6VXX complex (1P9M: IL6R/IL6/ IL6ST; 6VXX: Spike) based on these two different binding models, with the highest CDOCKER interaction energy found for quercetin. Meanwhile, the analysis results of the IL6R/IL6 and IL6ST binding models indicate that quercetin may affect IL6 signal transduction mediated by the binding of IL6R/IL6 and IL6ST to Gln124 (IL6R/IL6) and His145 (IL6ST). Nevertheless, wogonin may directly influence IL6ST-mediated cytokine signal transduction by competitively binding to Trp291 (IL6ST) and Asp198 (Spike) according to the binding model of IL6R/IL6/IL6ST and Spike. In this study, we screened for potential targets using network pharmacology, constructing PPI network, and exploring potential pathways by an enrichment analysis, to study the mechanism of protective action of LHQW against COVID-19 symptoms. Furthermore, verification based on molecular docking, which was developed to understand the mechanism and provide explanation for the effectiveness of LHQW for COVID-19 treatment, confirmed the importance of LHQW in the prevention and control of COVID-19 infection. Moreover, because of the use of computer simulation in this study, the potential risk of biosafety was avoided, and the efficiency of the research on sudden infectious diseases and TCM modernization was improved, thereby providing more scientific evidence for the clinical application of LHQW and direction for in-depth research. In addition, this study provides insight into prescription selection for the clinical treatment of COVID-19 and a reference for reducing the risk of the associated clinical CS. Our study is also expected to provide lead compounds for the development of new drugs for the treatment of COVID-19. In conclusion, this study revealed that the NTD region of the COVID-19 Spike protein had a high CDOCKER interaction energy with IL6ST using network pharmacology and molecular docking analysis. Quercetin in LHQW may directly affect the interaction between IL6R/IL6 and IL6ST by competitively binding to Gln124 (IL6R/IL6) and His145 (IL6ST) to mediate IL6R/IL6/IL6ST signal transduction. Wogonin may directly influence the interaction between IL6R/IL6/IL6ST and Spike via competitive binding to Trp291 (IL6ST) and Asp198 (Spike), affect IL6ST-mediated cytokine signal transduction, and regulate the excessive inflammatory stress caused by CS. LHQW may inhibit the LPS-mediated inflammatory response and regulate the AGE-RAGE signaling pathway by regulating the function of IL6, which in turn can influence oxidative stress, excessive inflammatory response, and immune function to control the inflammatory stress response during COVID-19 infection. Due to the limitations of the computational methods of chemistry and biology, the results of this study need to be verified by follow-up experiments to provide a basis for the treatment of COVID-19 using TCM.
Database and softwareThe following databases and software packages were used in this study: GeneCards  (https://www.genecards.org/), Traditional Chinese Medicine Systems Pharmacology Database  (TCMSP, http://tcmspw.com/tcmsp.php), UniProt  (https://www.uniprot.org/), STRING  (https://STRING-db.org/cgi/input.pl), PubChem  (https://pubchem.ncbi.nlm.nih.gov/), RCSB Protein Data Bank  (http://www.rcsb.org/), R(version R3.6.1) , Cytoscape (3.7.2) , ZDOCK Server 3.0.2 , Discovery Studio 2019 Client, Pymol (version 2.4.1).
Network pharmacology and molecular docking analysis of the mechanism of LHQW in treating COVID-19The detailed information is available in the Flow Charts 1 and 2.
Network pharmacology of the mechanism of LHQW in treating COVID-19
Screening of COVID-19 pathogenic targets, LHQW active ingredients and therapeutic targetsA search was conducted on the GeneCards database using the keywords “coronavirus disease 2019,” “coronavirus pneumonia,” “coronavirus” and “novel coronavirus 2019.” The results were exported to Microsoft Excel. For the purposes of this study, those with a score ≥1 were defined as pathogenic targets of COVID-19. Information has been added according to the literature. From LHQW components (Supplementary Table 1), in the TCMSP, the main active ingredients were screened according to oral bioavailability (OB) ≥30% and drug-likeness (DL) ≥0.18, and the therapeutic targets of the main active ingredients were screened through Targets information of TCMSP. Finally, gene symbol conversion was performed on the therapeutic targets screened in TCMSP according to UniProt. Additional information has been added to the literature.
Protein-Protein Interaction (PPI) network of COVID-19 pathogenic targets, LHQW therapeutic targets, junction targets of COVID-19 and LHQW and LHQW ingredients-targets networkCOVID-19 pathogenic targets, LHQW therapeutic targets, and junction targets of COVID-19 and LHQW through mapping were imported into the STRING database and optimized using Cytoscape, respectively. The degree value was calculated using a network analyzer. The degree value is an indicator of the importance of network nodes. The higher the value is, the stronger is the correlation between the corresponding node and other nodes. Subsequently, PPI-related information was visualized through a network diagram, and the core target of COVID-19 pathogenic targets, LHQW therapeutic targets, and junction targets were extracted based on the maximum value of degree, separately. Cytoscape was used to introduce the active ingredients and potential therapeutic targets of LHQW, and the ingredients-targets network was constructed.
Gene-Concept network of gene ontology enrichment analysis and KEGG pathway analysis of the junction targets of COVID-19 and LHQWGO includes three following aspects: biological process, cellular component, and molecular function. GO enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses were conducted for junction targets of COVID-19 and LHQW using the R and clusterProfiler  package in this study, respectively. In addition, the ggplot2  package was used to plot the Gene-Concept network to identify the biological process, the location of the reaction in the cell, the molecular function involved, and relevant KEGG pathway of LHQW in treating COVID-19, separately.
Network of active ingredients, GO and KEGG pathways corresponding to the core targetThe core target, corresponding GO and KEGG pathways, and active ingredients were imported into Cytoscape, respectively, to build the network.
Verification of molecular docking for LHQW in treating COVID-19 by the results of network pharmacologyAccording to the results of network pharmacology, IL6 is the core target of LHQW in treating COVID-19, and the active ingredients of LHQW, luteolin, quercetin, and wogonin are all connected to the core target IL6. Therefore, this study will carry out further molecular docking verification based on these results.
Protein-protein docking analysis between IL6R/IL6/IL6ST complex and COVID-19 spikeRCSB PDB databases were used to download the 3D structure PDB files of IL6R/IL6/IL6ST complex (PDB ID: 1P9M) and Spike protein (PDB ID: 6VXX), respectively, and then imported into Pymol to remove solvent molecules and ligands. After that, they were imported into ZDOCK SERVER 3.0.2 (Fast Fourier transform-based protein docking program) for protein-protein docking, and all possible binding models in the translation and rotation space between the two proteins were searched. The ZDOCK score was calculated using root-mean-square deviation (RMSD). Based on the maximum ZDOCK score, the optimal attitude of the two bindings was screened, and the output was the 1P9M-6VXX (IL6R/IL6/IL6ST-Spike) complex. The amino acid residues and hydrogen bonding forces of protein-protein binding were predicted using Pymol.
Protein-small molecule docking analysisCDOKER  is a high-accuracy molecular docking method based on the CHARMm force field. In this method, high-temperature kinetics are used to search for the flexible conformation of the ligand molecules. It is believed that the lower the CDOCKER interaction energy, the more stable the conformation of the ligand binding to the receptor. In addition, ≤−5.0 kcal mol−1 indicates that the protein and small molecule can bind, and ≤−7.0 kcal mol−1 indicates strong binding ability. The 1P9M-6VXX complex was introduced into the Discovery Studio 2019 Client, and polar hydrogens were added to the complex via Add Polar. According to the results of the network pharmacology and the characteristics of IL6 and Spike in physiological conditions, the active sites between IL6R/IL6 and IL6ST, between IL6R/IL6/IL6ST and Spike protein were identified from Receptor Cavities, and the Apply Forcefield was used to force field to the complex. The PubChem database was used to obtain the SDF files of the active ingredients to the core target IL6 as small molecules, which were imported into the Discovery Studio 2019 Client. The Apply Forcefield was used to generate the force field to the small molecules. Finally, the Dock Ligands (CDOCKER) was used for protein-small molecule docking, adjusting the Top Hits-Pose Cluster Radius to 0.5, and the other parameters were kept at default.
Availability of data and materialsPublicly available datasets were analyzed in this study. These data can be found here: [GeneCards] at (https://www.genecards.org/), [Traditional Chinese Medicine Systems Pharmacology Database (TCMSP)] at (http://tcmspw.com/tcmsp.php), [UniProt] at (https://www.uniprot.org/), [STRING] at (https://STRING-db.org/cgi/input.pl), [PubChem] at (https://pubchem.ncbi.nlm.nih.gov/), [RCSB Protein Data Bank] at (http://www.rcsb.org/).
3D: 3-dimension; ACE2: angiotensin I converting enzyme 2; AGE: Advanced Glycation End Product; AKT1: AKT serine/threonine kinase 1; ALB: Albumin; ARDS: Acute Respiratory Distress Syndrome; Arg: arginine; Asn: asparagine; Asp: aspartic acid; BAX: BCL2 associated X, apoptosis regulator; BCL2: BCL2 apoptosis regulator; BP: Biological Process; CASP3: Caspase 3; CC: Cellular Component; CCL2: C-C motif chemokine ligand 2; CHARMm: Chemistry at HARvard Macromolecular Mechanics; CNTF: ciliary neurotrophic factor; COVID-19: Corona Virus Disease 2019; CS: Cytokine Storm; CSF: Colony Stimulating Factor; CXCL10: C-X-C motif chemokine ligand 10; CXCL8: C-X-C motif chemokine ligand 8; Cys: cysteine; DL: Drug-Likeness; FOS: Fos proto-oncogene, AP-1 transcription factor subunit; Gln: glutamine; GO: Gene Ontology; His: histidine; ICAM1: Intercellular Adhesion Molecule 1; IFN: Interferon; IFNG: Interferon Gamma; IL10: Interleukin 10; IL1A: Interleukin 1 alpha; IL1B: Interleukin 1 beta; IL2: Interleukin 2; IL4: Interleukin 4; IL6: Interleukin 6; IL6R: Interleukin 6 receptor; IL6ST: Interleukin 6 receptor subunit beta; IL7: Interleukin 7; IL8: Interleukin 8; KEGG: Kyoto Encyclopedia of Genes and Genomes; LHQW: LianHuaQingWen; LIF: leukemia inhibitory factor; LPS: Lipopolysaccharide; Lys: lysine; mAb: monoclonal antibody; MAPK1: Mitogen-activated protein kinase 1; MAPK14: Mitogen-activated protein kinase 14; MAPK3: Mitogen-activated protein kinase 3; MAPK8: Mitogen-activated protein kinase 8; MCP-1: Monocyte Chemoattractant Protein-1; MERS-CoV: Middle East respiratory syndrome coronavirus; MF: Molecular Function; MODS: Multiple Organ Dysfunction Syndrome; NF-κB: nuclear transcription factor-κB; NMPA: National Medical Products Administration; NOS3: Nitric Oxide Synthase 3; NTD: N-terminus domain; OB: Oral Bioavailability; OSM: oncostatin-M; PDB: Protein Data Bank; PPI: protein-protein interaction; PRKCA: Protein kinase C alpha; PTGS2: Prostaglandin-endoperoxide synthase 2; RAGE: receptor for Advanced Glycation End Product; RBD: receptor-binding domain; RCSB: Research Collaboratory for Structural Bioinformatics; RELA: RELA proto-oncogene, NF-kB subunit; RMSD: root-mean-square deviation; S: Spike; SARS-CoV: SARS coronavirus; Ser: serine; SERPINE1: Serpin Family E Member 1; STAT1: Signal Transducer and Activator of Transcription 1; TCM: traditional Chinese medicine; TCMSP: Traditional Chinese Medicine Systems Pharmacology Database and Analysis Platform; Thr: threonine; TNF: Tumor Necrosis Factors; TP53: Tumor Protein p53; Trp: tryptophan; Tyr: tyrosine; UniProt: Universal Protein Resource.
Methodology, formal analysis, visualization, and writing—original draft preparation, Zhao TY; data curation, Cui XL, Wang YR, and Yue FL; software, Zhang M and He K; writing—editing, Chen Li and He K; conceptualization and supervision, Li J. All authors have read and agreed to the published version of the manuscript.
Conflicts of Interest
The authors declare no conflicts of interest. The funders had no role in the study design; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.
This research was supported by the Science and Technology Research Department of Jilin Province and the Education Department of Jilin Province (to Jing Li, No. JJKH20211228KJ).
- 1. Wu JT, Leung K, Leung GM. Nowcasting and forecasting the potential domestic and international spread of the 2019-nCoV outbreak originating in Wuhan, China: a modelling study. Lancet. 2020; 395:689–97. https://doi.org/10.1016/S0140-6736(20)30260-9 [PubMed]
- 2. Hu K, Guan WJ, Bi Y, Zhang W, Li L, Zhang B, Liu Q, Song Y, Li X, Duan Z, Zheng Q, Yang Z, Liang J, et al. Efficacy and Safety of Lianhuaqingwen Capsules, a repurposed Chinese Herb, in Patients with Coronavirus disease 2019: A multicenter, prospective, randomized controlled trial [Phytomedicine 85 (2021) 153242]. Phytomedicine. 2021; 93:153775. https://doi.org/10.1016/j.phymed.2021.153775 [PubMed]
- 3. Chinese Pharmacopeia Commission. Pharmacopeia of the People's Republic of China (version 1). Beijing: China Medical Science Press. 2015: 945–6.
- 4. Duan ZP, Jia ZH, Zhang J, Liu S, Chen Y, Liang LC, Zhang CQ, Zhang Z, Sun Y, Zhang SQ, Wang YY, Wu YL. Natural herbal medicine Lianhuaqingwen capsule anti-influenza A (H1N1) trial: a randomized, double blind, positive controlled clinical trial. Chin Med J (Engl). 2011; 124:2925–33. [PubMed]
- 5. Ding Y, Zeng L, Li R, Chen Q, Zhou B, Chen Q, Cheng PL, Yutao W, Zheng J, Yang Z, Zhang F. The Chinese prescription lianhuaqingwen capsule exerts anti-influenza activity through the inhibition of viral propagation and impacts immune function. BMC Complement Altern Med. 2017; 17:130. https://doi.org/10.1186/s12906-017-1585-7 [PubMed]
- 6. National Health Commission of the People’s Republic of China Diagnosis and treatment protocol for COVID-19 (Trial Version 8). National Health Commission of the People’s Republic of China: Beijing, China. 2020. http://www.nhc.gov.cn/yzygj/s7653p/202104/7de0b3837c8b4606a0594aeb0105232b.shtml.
- 7. The aim of the press conference of the joint prevention and control mechanism of the State Council is to promote the progress of scientific research on drug and vaccine development for COVID-19. National Medical Products Administration. [2020-04-14]. https://www.nmpa.gov.cn/zhuanti/yqyjzxd/yqyjxd/20200414205301912.html.
- 8. Yan H, Zou C. [Mechanism and material basis of Lianhua Qingwen capsule for improving clinical cure rate of COVID-19: a study based on network pharmacology and molecular docking technology]. Nan Fang Yi Ke Da Xue Xue Bao. 2021; 41:20–30. https://doi.org/10.12122/j.issn.1673-4254.2021.01.03 [PubMed]
- 9. Wang X, Zhang W, Wang M, Zhou Z, Zhu B, Li J, Zhang R, Zhang X, Li Q, Weng W. Study on the mechanism of Lianhua Qingwen capsule in the treatment of COVID-19 through network pharmacology and chemical components. Modernization of Traditional Chinese Medicine and Materia Media-World Science and Technology. 2020; 22:3169–77.
- 10. Ou X, Liu Y, Lei X, Li P, Mi D, Ren L, Guo L, Guo R, Chen T, Hu J, Xiang Z, Mu Z, Chen X, et al. Characterization of spike glycoprotein of SARS-CoV-2 on virus entry and its immune cross-reactivity with SARS-CoV. Nat Commun. 2020; 11:1620. https://doi.org/10.1038/s41467-020-15562-9 [PubMed]
- 11. Walls AC, Park YJ, Tortorici MA, Wall A, McGuire AT, Veesler D. Structure, Function, and Antigenicity of the SARS-CoV-2 Spike Glycoprotein. Cell. 2020; 183:1735. https://doi.org/10.1016/j.cell.2020.11.032 [PubMed]
- 12. Huang C, Wang Y, Li X, Ren L, Zhao J, Hu Y, Zhang L, Fan G, Xu J, Gu X, Cheng Z, Yu T, Xia J, et al. Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China. Lancet. 2020; 395:497–506. https://doi.org/10.1016/S0140-6736(20)30183-5 [PubMed]
- 13. Mehta P, McAuley DF, Brown M, Sanchez E, Tattersall RS, Manson JJ, and HLH Across Speciality Collaboration, UK. COVID-19: consider cytokine storm syndromes and immunosuppression. Lancet. 2020; 395:1033–34. https://doi.org/10.1016/S0140-6736(20)30628-0 [PubMed]
- 14. Chu H, Chan JF, Wang Y, Yuen TT, Chai Y, Hou Y, Shuai H, Yang D, Hu B, Huang X, Zhang X, Cai JP, Zhou J, et al. Comparative Replication and Immune Activation Profiles of SARS-CoV-2 and SARS-CoV in Human Lungs: An Ex Vivo Study With Implications for the Pathogenesis of COVID-19. Clin Infect Dis. 2020; 71:1400–09. https://doi.org/10.1093/cid/ciaa410 [PubMed]
- 15. Channappanavar R, Perlman S. Pathogenic human coronavirus infections: causes and consequences of cytokine storm and immunopathology. Semin Immunopathol. 2017; 39:529–39. https://doi.org/10.1007/s00281-017-0629-x [PubMed]
- 16. Xu Z, Li S, Tian S, Li H, Kong LQ. Full spectrum of COVID-19 severity still being depicted. Lancet. 2020; 395:947–8. https://doi.org/10.1016/S0140-6736(20)30308-1 [PubMed]
- 17. Chen N, Zhou M, Dong X, Qu J, Gong F, Han Y, Qiu Y, Wang J, Liu Y, Wei Y, Xia J, Yu T, Zhang X, Zhang L. Epidemiological and clinical characteristics of 99 cases of 2019 novel coronavirus pneumonia in Wuhan, China: a descriptive study. Lancet. 2020; 395:507–13. https://doi.org/10.1016/S0140-6736(20)30211-7 [PubMed]
- 18. Wang D, Hu B, Hu C, Zhu F, Liu X, Zhang J, Wang B, Xiang H, Cheng Z, Xiong Y, Zhao Y, Li Y, Wang X, Peng Z. Clinical Characteristics of 138 Hospitalized Patients With 2019 Novel Coronavirus-Infected Pneumonia in Wuhan, China. JAMA. 2020; 323:1061–9. https://doi.org/10.1001/jama.2020.1585 [PubMed]
- 19. Fang Y, Zhang H, Xu Y, Xie J, Pang P, Ji W. CT Manifestations of Two Cases of 2019 Novel Coronavirus (2019-nCoV) Pneumonia. Radiology. 2020; 295:208–9. https://doi.org/10.1148/radiol.2020200280 [PubMed]
- 20. Cui WW. Protective effect of Lianhua Qingwen capsule on lipopolysaccharide-induced acute lung injury in mice. Heilongjiang University of Chinese Medicine: Harbin, China. 2015.
- 21. Li J, Hu X, Liang F, Liu J, Zhou H, Liu J, Wang H, Tang H. Therapeutic effects of moxibustion simultaneously targeting Nrf2 and NF-κB in diabetic peripheral neuropathy. Appl Biochem Biotechnol. 2019; 189:1167–82. https://doi.org/10.1007/s12010-019-03052-8 [PubMed]
- 22. Zhou P, Yang XL, Wang XG, Hu B, Zhang L, Zhang W, Si HR, Zhu Y, Li B, Huang CL, Chen HD, Chen J, Luo Y, et al. Addendum: A pneumonia outbreak associated with a new coronavirus of probable bat origin. Nature. 2020; 588:E6. https://doi.org/10.1038/s41586-020-2951-z [PubMed]
- 23. Hoffmann M, Kleine-Weber H, Schroeder S, Krüger N, Herrler T, Erichsen S, Schiergens TS, Herrler G, Wu NH, Nitsche A, Müller MA, Drosten C, Pöhlmann S. SARS-CoV-2 Cell Entry Depends on ACE2 and TMPRSS2 and Is Blocked by a Clinically Proven Protease Inhibitor. Cell. 2020; 181:271–80.e8. https://doi.org/10.1016/j.cell.2020.02.052 [PubMed]
- 24. Lukassen S, Chua RL, Trefzer T, Kahn NC, Schneider MA, Muley T, Winter H, Meister M, Veith C, Boots AW, Hennig BP, Kreuter M, Conrad C, Eils R. SARS-CoV-2 receptor ACE2 and TMPRSS2 are primarily expressed in bronchial transient secretory cells. EMBO J. 2020; 39:e105114. https://doi.org/10.15252/embj.20105114 [PubMed]
- 25. Sungnak W, Huang N, Bécavin C, Berg M, Queen R, Litvinukova M, Talavera-López C, Maatz H, Reichart D, Sampaziotis F, Worlock KB, Yoshida M, Barnes JL, and HCA Lung Biological Network. SARS-CoV-2 entry factors are highly expressed in nasal epithelial cells together with innate immune genes. Nat Med. 2020; 26:681–7. https://doi.org/10.1038/s41591-020-0868-6 [PubMed]
- 26. Brouwer PJM, Caniels TG, van der Straten K, Snitselaar JL, Aldon Y, Bangaru S, Torres JL, Okba NMA, Claireaux M, Kerster G, Bentlage AEH, van Haaren MM, Guerra D, et al. Potent neutralizing antibodies from COVID-19 patients define multiple targets of vulnerability. Science. 2020; 369:643–50. https://doi.org/10.1126/science.abc5902 [PubMed]
- 27. Chi X, Yan R, Zhang J, Zhang G, Zhang Y, Hao M, Zhang Z, Fan P, Dong Y, Yang Y, Chen Z, Guo Y, Zhang J, et al. A neutralizing human antibody binds to the N-terminal domain of the Spike protein of SARS-CoV-2. Science. 2020; 369:650–5. https://doi.org/10.1126/science.abc6952 [PubMed]
- 28. Kang S, Tanaka T, Narazaki M, Kishimoto T. Targeting Interleukin-6 Signaling in Clinic. Immunity. 2019; 50:1007–23. https://doi.org/10.1016/j.immuni.2019.03.026 [PubMed]
- 29. Taga T, Kishimoto T. Gp130 and the interleukin-6 family of cytokines. Annu Rev Immunol. 1997; 15:797–819. https://doi.org/10.1146/annurev.immunol.15.1.797 [PubMed]
- 30. Grötzinger J, Kurapkat G, Wollmer A, Kalai M, Rose-John S. The family of the IL-6-type cytokines: specificity and promiscuity of the receptor complexes. Proteins. 1997; 27:96–109. https://doi.org/10.1002/(SICI)1097-0134(199701)27:1<96::AID-PROT10>3.0.CO;2-D. [PubMed]
- 31. Benigni F, Fantuzzi G, Sacco S, Sironi M, Pozzi P, Dinarello CA, Sipe JD, Poli V, Cappelletti M, Paonessa G, Pennica D, Panayotatos N, Ghezzi P. Six different cytokines that share GP130 as a receptor subunit, induce serum amyloid A and potentiate the induction of interleukin-6 and the activation of the hypothalamus-pituitary-adrenal axis by interleukin-1. Blood. 1996; 87:1851–4. https://doi.org/10.1182/blood.V87.5.1851.1851. [PubMed]
- 32. Stelzer G, Rosen N, Plaschkes I, Zimmerman S, Twik M, Fishilevich S, Stein TI, Nudel R, Lieder I, Mazor Y, Kaplan S, Dahary D, Warshawsky D, et al. The GeneCards Suite: From Gene Data Mining to Disease Genome Sequence Analyses. Curr Protoc Bioinformatics. 2016; 54:1.30.1–1.30.33. https://doi.org/10.1002/cpbi.5 [PubMed]
- 33. Ru J, Li P, Wang J, Zhou W, Li B, Huang C, Li P, Guo Z, Tao W, Yang Y, Xu X, Li Y, Wang Y, Yang L. TCMSP: a database of systems pharmacology for drug discovery from herbal medicines. J Cheminform. 2014; 6:13. https://doi.org/10.1186/1758-2946-6-13 [PubMed]
- 34. The UniProt Consortium. UniProt: the universal protein knowledgebase. Nucleic Acids Res. 2017; 45:D158–69. https://doi.org/10.1093/nar/gkw1099 [PubMed]
- 35. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, Simonovic M, Doncheva NT, Morris JH, Bork P, Jensen LJ, Mering CV. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019; 47:D607–13. https://doi.org/10.1093/nar/gky1131 [PubMed]
- 36. Kim S, Chen J, Cheng T, Gindulyte A, He J, He S, Li Q, Shoemaker BA, Thiessen PA, Yu B, Zaslavsky L, Zhang J, Bolton EE. PubChem 2019 update: improved access to chemical data. Nucleic Acids Res. 2019; 47:D1102–09. https://doi.org/10.1093/nar/gky1033 [PubMed]
- 37. Burley SK, Bhikadiya C, Bi C, Bittrich S, Chen L, Crichlow GV, Christie CH, Dalenberg K, Di Costanzo L, Duarte JM, Dutta S, Feng Z, Ganesan S, et al. RCSB Protein Data Bank: powerful new tools for exploring 3D structures of biological macromolecules for basic and applied research and education in fundamental biology, biomedicine, biotechnology, bioengineering and energy sciences. Nucleic Acids Res. 2021; 49:D437–51. https://doi.org/10.1093/nar/gkaa1038 [PubMed]
- 38. R Core Team. R: Language and Environment for Statistical Computing. R Foundation for Statistical Computing: Vienna, Austria; 2019. https://www.R-project.org/.
- 39. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, Schwikowski B, Ideker T. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003; 13:2498–504. https://doi.org/10.1101/gr.1239303 [PubMed]
- 40. Pierce BG, Wiehe K, Hwang H, Kim BH, Vreven T, Weng Z. ZDOCK server: interactive docking prediction of protein-protein complexes and symmetric multimers. Bioinformatics. 2014; 30:1771–3. https://doi.org/10.1093/bioinformatics/btu097 [PubMed]
- 41. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012; 16:284–7. https://doi.org/10.1089/omi.2011.0118 [PubMed]
- 42. Wickham H. ggplot2: Elegant graphics for data analysis. Springer-Verlag: New York. 2016.
Ding X, Wu Y, Wang Y, Vilseck JZ, Brooks CL
3rd. Accelerated CDOCKER with GPUs, Parallel Simulated Annealing, and Fast Fourier Transforms. J Chem Theory Comput. 2020; 16:3910–9. https://doi.org/10.1021/acs.jctc.0c00145 [PubMed]