Utilizing the abundant data from biological and clinical studies on CQ and HCQ (Gasmi et al., 2021), we sought to characterize the drug-gene interaction network and explore the entire mechanisms of CQ and HCQ treatment (including targeted and off-target mechanisms). The overlapping genes between CQ, HCQ and SARS-CoV-2 were identified by a combination of PubChem/SwissTarget/STITICH and constructed using Cytoscape (Figure 2). As displayed in Figure 2, CQ/HCQ-SARS-CoV-2 target network involved 170 nodes (47 common targets and 2 corresponding chemical components) and 189 edges. The results indicated that the 47 nodes might act as potential gene targets for the treatment of COVID-19 using CQ or HCQ. A Venn Diagram visualizing overlap between the gene targets of CQ, HCQ, and SARS-CoV-2 genes identified from PubChem. 47 overlapped target genes were identified (Figure 2B) and then were used to construct a drug-gene interaction network on Cytoscape (Figure 2C). 22 out of 47 genes in the drug-gene interaction network were involved in TLR and IL-6 signaling pathways (highlighted in yellow, Figure 2C). Notably these 22 genes contributed to the common genes (9 out of 12 genes) between CQ and HCQ gene networks (Figure 2C).

Given the interesting correlation of common genes in CQ and HCQ network to TLR or IL-6 signaling pathway (Figure 1), we sought to further elucidate the mechanism of SARS-CoV-2 drug treatment at the system biology level. Thus, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of cellular components and pathways using Annotation, Visualization and Integrated Discovery (DAVID) on the 47 target genes in COVID-19 and CQ/HCQ drug-gene interaction network in Figure 3. The top 25 terms ranked by fold-enrichment and p < 0.05 for each analysis were represented using a bubble graph where the most enriched annotations were at the bottom and the number of genes involved were represented by the size of the circles (Figure 3). For cellular components, the functional annotations were associated with Cytosol (GO: 0005829) (Figure 3A). For pathway analysis, 47 gene targets in the COVID-19 and CQ/HCQ drug-gene network were mostly enriched in TLR signaling pathway (GO: 0002224) and nuclear factor kappa light chain enhancer of activated B cells (NF-κB) signaling pathway (GO: 0038061) (Figure 3B).

The drug-gene network characterization in Figure 2 and GO/pathway enrichment analysis in Figure 3 identified an interesting lead on the mechanism of CQ/HCQ in involving the gene targets from TLR signaling pathway and NF-κB signaling pathway. The characterization of the protein network would deepen the understanding on the mechanism targeting SARS-CoV-2. Therefore, we characterized the protein targets of CQ and HCQ that involved in SARS-CoV-2 network in Figure 4. A protein-protein interaction (PPI) network was constructed using a combination of SwissTarget (Gfeller et al., 2014), STITCH (Kuhn et al., 2008) and PubChem on Cytoscape (Smoot et al., 2011) to identify the direct and indirect drug targets of CQ and HCQ (Figure 4A). Following the overlap with SARS-CoV-2 proteins identified from PubChem, we identified 7 common protein targets (TLR9, IL-6, TNF, ACE2, GSTA2, HMGB1, GSTM1) of CQ and HCQ. Among these targets, four of them (TNF, GSTA2, HMGB1, GSTM1) were shown to be specific to CQ and three of them were the common targets (TLR9, IL-6, ACE2) of CQ and HCQ for SARS-CoV-2 (Figure 4B). These overlapped target proteins were then used to construct a drug-protein interaction network on Cytoscape (Figure 4B).

Thus, the results in Figure 4 suggested TLR9 and IL-6 were the common protein targets of CQ and HCQ in SARS-CoV-2 network. The three common protein targets (TLR9, IL-6, ACE2) were of significance in the network and thus used for the following molecular docking study in Figure 5. In order to verify the direct interactions between three protein targets (TLR9, IL-6, ACE2) in SARS-CoV-2 network and CQ or HCQ, we performed molecular docking analysis (Figure 5). The docking affinity values between the targets and CQ/HCQ were reported by AutoDock Vina. Meanwhile, the hydrophobic, electrostatic, cation-pi interactions and hydrogen bonds formed in the receptor-spike complex were analyzed using Ligplot+. There were six pairs (TLR9-CQ, ACE2-CQ, IL-6-CQ, TLR9-HCQ, ACE2-HCQ, IL-6-HCQ) delivered into the docking simulation, among which ACE2 was included as positive control due to its direct interaction with CQ and HCQ (Badraoui et al., 2021). The greater the absolute value of the docking affinity, the stronger binding ability between the compounds and the active site of the targets.

The binding affinity of positive control pairs ACE2-CQ docking (−6.5 kcal/mol) and ACE2-HCQ docking (−6.6 kcal/mol) validated this molecular docking system (Figure 5A). The results suggested TLR9 possessed strong binding affinity to CQ and HCQ based upon the analysis of TLR9-CQ docking (−5.9 kcal/mol) and TLR9-HCQ docking (−6.4 kcal/mol). In comparison, IL-6 exhibited weaker binding affinity to both the drugs based on the results of IL-6-CQ docking (−4.8 kcal/mol), IL-6-HCQ docking (−5.1 kcal/mol) (Figure 5A). The detailed interaction patterns were demonstrated in Figure 5B. Among all the pairs of ACE2, TLR9 and IL-6, CQ and HCQ fit into the interface pocket formed by interaction amino acid residues in all three proteins (Figure 5B). The results showed hydrogen bond formation involved in the bonding between CQ and residues in Gln779(A) at TLR9 and Asp350(A) at ACE2. These hydrogen bonds between CQ and TLR9/IL-6 were considered as strong interactions considering the distance of 3.18 Å (TLR9-CQ), 3.09 Å (ACE2-CQ) respectively. In comparison, the hydrogen bond formation between HCQ and residues involved in Leu37(A) at TLR9 and Arg393(A) at ACE2, and Gln28(A) at IL-6. The hydrogen bonds between HCQ and TLR9/IL-6 were also considered as strong interactions given the distance of 2.97 Å (TLR9-HCQ), 3.09 Å (ACE2-HCQ), 2.99 Å (IL-6-HCQ) (Figure 5C). The other essential residues interacted with CQ and HCQ through hydrophobic effects, van der Waals forces, etc. (Figure 5C).

In analyzing the stability of TLR9-ligand complexes we utilized RMSD as the main readout. A larger RMSD value indicates lower stability of the protein complex in the given simulation conditions (Das et al., 2021a; Das et al., 2021b). Comparing RMSD of TLR9 backbone solvated in water and the TLR9-CQ complex it is evident there is minimal differences between the stability of the complex compared to the protein itself with average RMSD values of 0.1402 and 0.1427 nm respectively (p > 0.05). Similarly, the average RMSD of TLR9-HCQ was 0.1400 nm showing similar stability to the CQ bound counterpart yet no significant difference exists in average RMSD compared to free TLR9 (Figure 6). In order to further validate any potential differences seen between free TLR9 and the respective ligands we conducted a longer MD simulation to resolve any differences in stability between the complexes. Upon performing a 10 ns simulation the RMSD values stabilized at distinct values. Average RMSD values of 3.226, 3.069, and 2.339 nm were observed for TLR9 solvated in water, TLR9-CQ, and TLR9-HCQ respectively (Figure 6). Coinciding with binding affinity the average RMSD value of TLR9-HCQ showed a significant difference (p < 0.05) compared to free TLR9 while the same comparison for TLR9-CQ yielded no significance. This suggests the deviation of atoms between TLR9 and HCQ is less than that observed between TLR9 and CQ.

Molecular docking results in Figures 5, 6 confirmed TLR9 as direct targets of CQ and HCQ. Then the next question is how we can use the tools of CQ and HCQ to explore the possibility of repurposed drug development to disrupt the SARS-CoV-2 network in human body. To address this question, we have established SARS-CoV-2-human-CQ/HCQ interactome in Figure 7 and investigated the connectivity between TLR9 and IL-6, their PPIs and SARS-CoV-2 proteins that play a critical role in SARS-CoV-2 pathogenesis. By merging direct interaction partners of SARS-CoV-2 proteins from BioGRID with the human-protein interactions from Integrated Interactions Database (IID), we generated the SARS-CoV-2-human-CQ/HCQ interactome map (Figure 7A). As constructed on Cytoscape, the innermost circle were SARS-CoV-2 proteins, followed by the orange layer of direct targets, and the blue outer circle of secondary proteins (Figure 7A).

Using Diffusion, a network propagation algorithm for Cytoscape (Oughtred et al., 2019), we constructed a close-up view showing the interactions between the CQ and HCQ drug targets (TLR9 and IL-6) were associated with human proteins (BRD4, TBK1, NPC2, ERC1, RIPK1) that were found to interact with SARS-CoV-2 proteins components of the replicase-transcriptase complex (Nsp12, Nsp13), accessory factors (ORF8) and structural protein (E) (Figure 7B). Among these interactions, SARS-CoV-2 replicase-transcriptase complex (Nsp12 and Nsp13) was targeted by TLR9. As targeting IL-6 is a potential therapeutic strategy for COVID-19 (Y and R, 2021), our results showed that IL-6 targeted replicase-transcriptase complex (Nsp13), accessory factor (ORF8) and structural protein (E) (Figure 7B). Specifically, the role of TLR9 on SARS-CoV-2 human interactome were achieved through the interaction of TLR9-TBK1-Nsp13; TLR9-ERC1-Nsp13; TLR9-RIPK1-Nsp12 (Figure 7B).

Drug-Protein targets of CQ and HCQ were identified through 3 different databases: The Search Tool for Interactions of Chemicals (STITCH) (http://stitch.embl.de/), SwissTargetPrediction (http://www.swisstargetprediction.ch/), and PubChem (https://pubchem.ncbi.nlm.nih.gov/). On STITCH, the minimum required interaction score was set to high (0.700) and the search was conducted within the Homo sapiens organism. On SwissTargetPrediction, the minimum probability of interaction was set at 0.50. Gene targets of hydroxychloroquine and chloroquine were identified through the PubChem database under “Chemical-Gene Interactions.” Proteins and genes involved in SARS-CoV-2 infection were identified under the PubChem database.

The network visualization tool Cytoscape 3.6.1 (http://cytoscape.org/, ver. 3.5.1) was adopted to obtain the PPI network map. Common targets between the compound-putative target network and SARS-CoV-2 target PPI network were identified as potential targets for the components of CQ and HCQ in COVID-19. In such a network, a compound or a gene/protein serves as an “edge” and reflects an association between nodes.

In order to clarify the functions of the potential targets and their involvement in signaling pathways, GO enrichment and KEGG pathway enrichment analyses of the targets in the compound-COVID-19 target network were performed using the Database for DAVID 6.8 (https://david.ncifcrf.gov/). The GO project elucidates the biological significance of a set of genes through three categories: cellular component, molecular function, and biological process. The 47 overlapping genes from the COVID-19-compound network were inputted into the DAVID database to reveal gene functions that were overrepresented among the set. KEGG pathway analysis which reveals pathway functional annotations of a gene set, was also performed using the 47 overlapping genes through the DAVID platform. The top 25 terms with the highest fold enrichment values and p-value < 0.05 from the GO and KEGG analyses were chosen and visualized using the “GOplot” package in R software.

The crystal structures of TLR9 (3WPF, 1.96 Å), ACE2 (1R42, 2.20 Å), IL-6 (4O9H, 2.42 Å) were retrieved from the Protein Data Bank (https://www.rcsb.org). Using AutoDock Tools 1.5.6, the complexes were prepared by removing the original ligand, water molecules, and other structures, leaving only the protein of interest. Next, polar hydrogen bonds were added, and the structures were saved in .pdbtq format.

The 3-D structures of the ligands CQ (CID_2719) and HCQ (CID_3652) were obtained from PubChem (https://pubchem.ncbi.nlm.nih.gov/) in .sdf format. The molecular geometries of the ligands were then optimized using Avogadro 1.2.0 (https://avogadro.cc/) with Force Field type MMFF94 and saved in .pdb format.

All ligand and receptor files were saved as. pdbqt format. Then we used Autodock Vina, a freely available open-source package, to evaluate and verify the binding affinity of compound-target relationship, and the prediction results from network pharmacology. A grid box was generated to cover the entire protein with the X, Y, Z coordinates set at the center of the macromolecule. Furthermore, the exhaustiveness was set at 24 as per the published protocol by Forli et al. The binding models were visualized by PyMol2.3.0 software and Ligplot+ 2.2. software.

The docked output structures from molecular docking of TLR9 and ACE2 were converted to pdb files through PyRx (https://pyrx.sourceforge.io/) & BIOVIA Discovery Studio Visualizer (https://discover.3ds.com/discovery-studio-visualizer-download) using the confirmation of either CQ or HCQ that had the lowest vina result value (Model 1) along with the protein without the ligand serving as controls. Charmm-gui’s solution builder (https://www.charmm-gui.org/?doc=input/solution) was then used to create the cubical water box, fitted to the protein complex, the force field, and protonation states. Charmm-gui was also used to add counter ions to make the system electrically neutral before each molecular dynamic simulation. NAMD software (https://www.ks.uiuc.edu/Research/namd/) was then used to run the equilibration step using canonical NVT and NPT ensembles at a stable temperature of 300 K and pressure of 1 bar. The production step followed using canonical parameters optimized through Charmm-gui. Once all files were made, VMD software (https://www.ks.uiuc.edu/Research/vmd/) was used to visualize the molecule and run the simulation by loading the input file of the protein or complex and loading the production dcd file into it. With VMD’s RMSD trajectory tool, RMSD graphs were produced with RMSD (nm) on the y-axis and time on the x-axis.

Human proteins that directly interact with SARS-CoV-2 proteins identified by Gordon et al. were retrieved from the BioGRID database (https://thebiogrid.org/220839/publication/a-sars-cov-2-protein-interaction-map-reveals-targets-for-drug-repurposing.html). In order to explore the relationship of these protein with the human PPI interactome, we expanded the network by inputting the human proteins into the Integrated Interactions Database (IID) (http://iid.ophid.utoronto.ca/), which is an extensive interaction network that includes PPIs from 9 curated databases, PPIs from orthologues, and high confidence PPIs from established computational methods. We then mapped the output PPI from IID onto Cytoscape 3.6.1 (https://cytoscape.org/, ver. 3.5.1) and merged it with the SARS-CoV-2 human PPIs. Using Cytoscape’s built-in diffusion algorithm, we traced the interaction of the proteins targeted by CQ and HCQ to the SARS-CoV-2 viral proteins.