Safety assessment of drug combinations used in COVID-19 treatment: in silico toxicogenomic data-mining approach.

Improvement of COVID-19 clinical condition was seen in studies where combination of antiretroviral drugs, lopinavir and ritonavir, as well as immunomodulant antimalaric, chloroquine/hydroxychloroquine together with the macrolide-type antibiotic, azithromycin, was used for patient's treatment. Although these drugs are "old", their pharmacological and toxicological profile in SARS-CoV-2 - infected patients are still unknown. Thus, by using in silico toxicogenomic data-mining approach, we aimed to assess both risks and benefits of the COVID-19 treatment with the most promising candidate drugs combinations: lopinavir/ritonavir and chloroquine/hydroxychloroquine + azithromycin. The Comparative Toxicogenomics Database (CTD; http://CTD.mdibl.org), Cytoscape software (https://cytoscape.org) and ToppGene Suite portal (https://toppgene.cchmc.org) served as a foundation in our research. Our results have demonstrated that lopinavir/ritonavir increased the expression of the genes involved in immune response and lipid metabolism (IL6, ICAM1, CCL2, TNF, APOA1, etc.). Chloroquine/hydroxychloroquine + azithromycin interacted with 6 genes (CCL2, CTSB, CXCL8, IL1B, IL6 and TNF), whereas chloroquine and azithromycin affected two additional genes (BCL2L1 and CYP3A4), which might be a reason behind a greater number of consequential diseases. In contrast to lopinavir/ritonavir, chloroquine/hydroxychloroquine + azithromycin downregulated the expression of TNF and IL6. As expected, inflammation, cardiotoxicity, and dyslipidaemias were revealed as the main risks of lopinavir/ritonavir treatment, while chloroquine/hydroxychloroquine + azithromycin therapy was additionally linked to gastrointestinal and skin diseases. According to our results, these drug combinations should be administrated with caution to patients suffering from cardiovascular problems, autoimmune diseases, or acquired and hereditary lipid disorders.

Abreviations

Introduction

During the last four months, the world has fought an unprecedented health hazard – an infection caused by an unknown species of coronavirus (SARS-CoV-2) first reported in December 2019 in China and named COVID-19 (Wujtewicz et al., 2020). On the 11th of March 2020, The World Health Organization has announced a pandemic of COVID-19 disease. At that time, there were 118,319 reported cases and 4292 deaths due to the infection with SARS-CoV-2 (WHO, 2020a). Thus far, there are 10,021,401 confirmed cases and 499,913 deaths caused by the new coronavirus in 213 countries and these numbers are continuously growing (WHO, 2020b, Situtation Report, June 29th). The clinical spectrum of COVID-19 ranges from the asymptomatic infection and mild upper respiratory tract illness, to severe pneumonia, progressive respiratory failure and death (Velavan and Meyer, 2020). Key steps in the COVID-19 treatment include the prevention of the virus spread, which is achieved by home isolation of suspected cases and those with mild illnesses, along with the strict infection control measures at hospitals that include contact and droplet precautions. Mild illness should be treated at home, managed by maintaining hydration and nutrition, with the avoidance of non-adequate use of antibiotics and antiviral drugs. In severe, hospitalized patients, provision of oxygen through nasal prongs, face mask, high flow nasal cannula or non-invasive ventilation are indicated as well as mechanical ventilation, extra corporeal membrane oxygen support and renal replacement therapy (Singhal, 2020). Today, there are different recommendations for COVID-19 treatment, including the use of “off-label” drugs, traditional Chinese medicine approach, and intravenous application of immunoglobulins, interferons and plasma of patients recovered from COVID-19 (Abedi et al., 2020; Cao et al., 2020; Ling et al., 2020; Saavedra, 2020). The most common drugs used in order to ameliorate patient's condition are oseltamivir, ganciclovir, remdesivir, arbidol, as well as lopinavir, ritonavir, chloroquine and azithromycin (Cao et al., 2020; Gautret et al., 2020a; Ye et al., 2020). It has been noticed that most of them are aiming better therapy perspective when used in combinations. However, those drug combinations might carry some risks in the terms of safety and potential adverse effects (Cao et al., 2020).

Lopinavir, antiretroviral drug of the protease inhibitor class, was identified to have in vitro inhibitory activity against SARS-CoV, the virus that causes severe acute respiratory syndrome (SARS) in humans (Cao et al., 2020). Ritonavir, another antiretroviral drug, is usually combined with lopinavir to increase its plasma half-life by inhibiting cytochrome P450 enzyme (Cao et al., 2020). A combination of these two drugs significantly improved the clinical condition of SARS patients (Wu et al., 2020). Early clinical trial performed in Chinese COVID-19 patients showed that chloroquine had a positive effect on clinical outcome, as well as viral clearance (Gao et al., 2020; Gautret et al., 2020a). The results of the same trial have demonstrated that azithromycin added to hydroxychloroquine to prevent bacterial super-infection significantly alleviated virus elimination, suggesting a possible synergistic effect of the combination (Gautret et al., 2020a). However, clinical trials have shown that for lopinavir/ritonavir combination most frequent adverse effects were gastrointestinal, while serious adverse effects included severe anaemia, shock, acute kidney injury, acute hearth failure, ect., consequently leading to the early treatment cessation (Cao et al., 2020). As for chloroquine, cardiovascular were the most serious adverse effects, and it is generally recommended to avoid concurrent administration of other drugs known to prolong the QT interval (Borba et al., 2020; Cortegiani et al., 2020).

One way to further identify potential side effects of drug combinations and assess the risks and benefits of the multi-drug therapy is by using bioinformatics. Considering that toxicogenomics incorporates bioinformatics and conventional toxicology to investigate the interactions between genes and environmental diseases (Boverhof and Zacharewski, 2006; Tung et al., 2020; Waters and Fostel, 2004), it can be used to predict gene function and explore molecular pathways in the onset of disease (Van Breda et al., 2014). It can also be applied in the mixture assessment process, since it provides strategies to explore all the possible genomic interactions that may be of crucial importance in assessing their mechanism of toxicity (Boverhof and Zacharewski, 2006).

Thus, by using in silico toxicogenomic data-mining approach, we evaluated both risks and benefits of proposed anti-COVID-19 therapy. Such systems biology approach enables understanding of the molecular pathways and biological processes behind drugs-induced side effects which may be crucial for their management and prevention.

Methods

In silico data-mining sources

The Comparative Toxicogenomics Database (CTD; http://CTD.mdibl.org), Cytoscape software package version 3.8.0 (https://cytoscape.org) and ToppGene Suite portal (https://toppgene.cchmc.org) were used as the main in silico data-mining tools in our research.

The Comparative Toxicogenomic Database is a publicly available database that processes and integrates data describing relationships between chemicals/drugs, genes, and diseases (Davis et al., 2019). This database is updated on a monthly basis, thus allowing the accuracy, cohesiveness and accessibility of data (Davis et al., 2019, Davis et al., 2008). It primarily focuses on chemicals present in the environment and outcomes relevant to human health, but also provides data from animal studies if the gene of interest is also relevant to humans (Davis et al., 2008; Meng et al., 2013).

The database contains various tools that can be used for the extraction of toxicogenomic data. For our in silico analysis we used SetAnalyzer and VennViewer CTD tools. SetAnalyzer CTD tool performs set-based enrichment for collections of chemicals or genes, and pathway generation for collections of genes while VennViewer CTD tool forms a diagram to compare datasets of the observed drug/chemical, gene, and disease relationships.

Cytoscape is a free software package used for visualizing, modelling and analysing molecular and genetic interaction networks. It can be used in functional genomics and proteomics experiments for obtaining an interaction network for genes of interest (Cline et al., 2007). It can also be combined with large databases of protein-protein, protein-DNA, and genetic interactions and numerous plug-ins which provide a way to implement new algorithms, additional network analyses, and biological semantics (Cline et al., 2007; Shannon et al., 2003).

ToppGene Suite portal (https://toppgene.cchmc.org) is a freely available online tool which aids in identification of genes sets and their prioritization based on functional annotations and interactions with proteins (Chen et al., 2009). ToppFun ToppGene tool (https://toppgene.cchmc.org/enrichment.jsp) was applied in this research. This tool can be used for transcriptome, ontology, phenotype, proteome, and pharmacome annotations based gene list functional enrichment analysis (Chen et al., 2009).

Data extraction

The flow chart for the different steps of our analyses is shown in the Fig. 1 . In this study, information about the genes affected by the investigated drugs was obtained from the “Genes” data-tabs on the CTD database for each investigated drug. Lopinavir and ritonavir were found in CTD separately and as fixed lopinavir/ritonavir combination under 3 different names: Aluvia, Kaletra, and Lopimune. Thus, we evaluated the lopinavir/ritonavir – gene – disease interactions in two ways. Firstly, we examined fixed lopinavir/ritonavir combination and secondly, we compared it with the results obtained when lopinavir and ritonavir were regarded separately. Since chloroquine/hydroxychloroquine, and azithromycin are present in CTD database only as separate compounds, we investigated their combined effects. Extracted data was further analysed by using freely-available in silico tools: (1) Cytoscape (version 3.8.0) plug-ins: GeneMANIA (http://genemania.org/plug-in/), CytoHubba (http://apps.cytoscape.org/apps/cytohubba), and ClueGO version 2.5.6 (http://apps.cytoscape.org/apps/cluego) + CluePedia version 1.5.6 (http://apps.cytoscape.org/apps/cluepedia); (2) ToppGene Suite portal (https://toppgene.cchmc.org) and its tool ToppFun function; and (3) analytical tools from the CTD database – VennViewer (http://ctdbase.org/tools/vennViewer.go) and SetAnalyzer (http://ctdbase.org/tools/analyzer.go).

Fig. 1

Flow chart explaining the steps of our analysis.

Formation of gene interaction network and identification of hub genes

In this research, CTD VenViewer tool was used to obtain a list of genes common to the investigated drugs. Next, to generate a tight network of genes/proteins associated with the investigated drug combinations and to obtain the lists of related genes, extracted common genes were inserted into the GeneMANIA Cytoscape plug-in, while Homo sapiens was selected as target organism. Since constructing a disease specific network and using its topological parameters to rank the hub genes is one of the widely used approaches in biomarker prioritization studies (Ghatge et al., 2018), we used CytoHubba plug-in to calculate the top 15 hub genes out of the constructed GeneMANIA network and ranked them by the Maximal Clique Centrality (MCC) score. MCC, found to be the most accurate of all the methods available in the CytoHubba, was used considering its ability to generate precise predictions of essential proteins (Chin et al., 2014).

Gene ontology (GO) enrichment

In order to understand the biological importance of the generated hubs and to elucidate biological networks and pathways most relevant to the investigated drug combinations, we performed the gene ontology (GO) enrichment analysis.

GO: Biological processes and molecular functions enrichment

ToppGene ToppFun function was used to obtain the lists of the 10 most significant molecular functions and biological processes which may elucidate the connection between the use of the investigated drug combinations and disease development. The default settings of the ToppGene Suite portal were selected for this analysis, meaning that the p-value was set to 0.05 and FDR corrected.

GO: Pathway enrichment analysis

In our research, pathway analysis was performed by Cytoscape ClueGO plug-in version 2.5.6. KEGG, Reactome, and WikiPathways databases were selected to extract the list of pathways. The two-sided hypergeometric test was used for the enrichment with a Bonferroni step down correction and a κ score of 0.3 to link the terms.

Gene-diseases relationships

SetAnalyzer CTD tool was used for the identification of diseases associated with our set of 15 hub genes and thus, the use of investigated drugs. By default, p-value was set at 0.01. It is important to acknowledge that this tool predicts the link based on interacting hub genes without referring to the specific role of the genes in the found disease. Furthermore, we used VenViewer CTD tool to predict additional diseases associated with the data sets connected to the investigated drugs.

Results

Lopinavir/ritonavir combination

Common genes and gene ontology

Top ten interacting genes reported in CTD for fixed antiviral combination of ritonavir and lopinavir were namely: GLB1, TP53, CCL2, CDKN1A, CXCL8, IL6, APOA1, APOB, APOC3, and LMNA. It was noted that this drug mixture increases the expression, activity or simply secretion of all the mentioned genes/proteins which may explain the development of unwanted effects of its use. This drug combination acts on 171 gene ontology processes and 44 molecular pathways.

On the other hand, when regarded as separate compounds, lopinavir interacted with 16, while ritonavir affected the expression of 114 genes among which 13 were in common (Fig. 2A): ABCB11, ADIPOQ, CCL2, CCL3, CYP3A4, HMOX1, ICAM1, IL6, LEP, LMNA, RB1, TNF, and ZMPSTE24. They influence the expression of genes in the same manner, with the exception of CCL3 and RB1, on which ritonavir and lopinavir show dual effects. In particular, expression of the TNF, LMNA, IL6, ICAM1, HMOX1, CCL2 was increased by the use of ritonavir and lopinavir, while downregulation was seen in the case of ABCB11, ADIPOQ, CYP3A4, LEP, and ZMPSTE24. The number of common enriched biological processes was as high as 581 (Fig. 2B), while 11 molecular functions (Fig. 2C) were found to be mutual for ritonavir and lopinavir. The total of 655 and 1109 biological processes, and 13 and 107 molecular pathways were detected for lopinavir and ritonavir, respectively.

Fig. 2

Venn diagram showing the number of genes associated with lopinavir and ritonavir, single and combined (A) Venn diagram showing the enriched Biological Processes (B) and Molecular functions (C) associated with lopinavir and ritonavir, single and combined (VennViewer CTD tool).

Notably, only 4 genes were described as interacting genes in both analyses (CCL2, IL6, LMNA, and TNF) and thus, we further examined the safety of the fixed lopinavir/ritonavir combination separately from the combination of ritonavir and lopinavir, present in CTD as isolated compounds.

Gene network

Further step of our analysis included constructing a tight networks of lopinavir/ritonavir interacting genes, along with 20 related genes. In the case of fixed antiviral combination, the results have shown that more than half of these genes (52.12%) were in physical interactions, which means that there is an interaction between the protein products of these genes while 10.23% of them were in co-expression, which means that their expression levels are similar under the defined conditions in the gene expression study. 8.31% belonged to the same molecular pathway, 5.95% of the interactions were predicted by the server, 5.70% interactions were genetic, while 0.13% were shared protein domains (Fig. 3A). Gene network between 13 lopinavir/ritonavir interacting genes along with 20 related genes has shown that the majority of these genes (77.27%) were in co-expression. 12.47% of the investigated genes were in co-localization, 8.82% were predicted by the server, 1.28% were shared protein domains, while 0.16% were genetic interactions (Fig. 3B).

Fig. 3

Tight network between genes (A) connected with lopinavir/ritonavir fixed combination and (B) common to lopinavir and ritonavir when regarded separately, together with 20 related genes (GeneMANIA Cytoscape plug-in). Color legend: co-expression-purple; co-localization-blue; physical interaction-red; pathway-light blue; predicted-orange. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Top 15 hub genes extracted from the constructed networks are presented in Fig. 4A (for fixed lopinavir/ritonavir combination) and Fig. 4B (when lopinavir and ritonavir are regarded separately).

Fig. 4

15 hub genes (ranked by the Maximal Clique Centrality (MCC) score) for lopinavir/ritonavir fixed combination (A) and lopinavir/ritonavir combination (B). Color of the nodes represents the intensity of correlation with yellow being the lowest and red the highest (Cytohubba plug-in + MCC method). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

In the two sets, 6 genes were in common: CCL2, CXCL1, CXCL2, CXCL8, ICAM1, and IL6, all highly involved in the regulation of immune system related biological processes. Moreover, even the intensity of correlation for mentioned genes was similar in two sets, with an exception of CXCL1. Fixed lopinavir/ritonavir combination increased its expression to a higher level when compared to lopinavir plus ritonavir treatment.

Gene enrichment analysis

Gene ontology (GO) describes gene products in terms of their associated biological processes, cellular components and molecular functions. The top 10 molecular functions and biological processes for the generated gene sets are shown in the Table 1 . Molecular functions of the hub genes related to fixed lopinavir/ritonavir combination were shown to be highly associated to chemokine activity and lipid metabolism, such as cholesterol and sterol binding, or phosphatidylcholine-sterol O-acyltransferase activity. Similarly, the most important biological processes were related to chylomicron remodeling, lipid metabolism, as well as interleukin signaling and inflammation, as shown below. On the other hand, we noticed that majority of the hub genes linked to lopinavir/ritonavir combination show cytokine related activity and consequently, regulate inflammatory response. Additionally, under the top 10 biological processes, we detected cellular responses to lipopolysaccharides and lipids (Table 1).

Table 1

Top 10 molecular functions potentially linked to the use of lopinavir/ritonavir as fixed combination, or when regarded as separate substances in the CTD database (ToppGene ToppFun function, p-value: 0.05, FDR corrected).

	Lopinavir/Ritonavir fixed combination				Lopinavir + Ritonavir (VennViewer + SetAnalyzer)
GO	ID	Name	p-value	Genes (number)	ID	Name	p-value	Genes(number)
Molecular functions	GO:0008009	chemokine activity	2.840E-10	5	GO:0005125	cytokine activity	1.640E-14	9
	GO:0070653	high-density lipoprotein particle receptor binding	3.835E-10	3	GO:0005126	cytokine receptor binding	3.174E-13	9
	GO:0015485	cholesterol binding	5.595E-10	5	GO:0008009	chemokine activity	1.109E-12	6
	GO:0032934	sterol binding	9.369E-10	5	GO:0042379	chemokine receptor binding	9.914E-12	6
	GO:0042379	chemokine receptor binding	1.729E-9	5	GO:0048018	receptor ligand activity	2.077E-11	9
	GO:0120020	cholesterol transfer activity	1.734E-9	4	GO:0030546	signaling receptor activator activity	2.271E-11	9
	GO:0120015	sterol transfer activity	2.098E-9	4	GO:0030545	receptor regulator activity	4.883E-11	9
	GO:0005102	signaling receptor binding	4.933E-9	11	GO:0001664	G protein-coupled receptor binding	7.102E-8	6
	GO:0060228	phosphatidylcholine-sterol O-acyltransferase activator activity	7.654E-9	3	GO:0005102	signaling receptor binding	1.060E-7	10
	GO:0098772	molecular function regulator	8.160E-9	11	GO:0098772	molecular function regulator	2.578E-6	9
Biological processes	GO:0034371	chylomicron remodeling	2.732E-14	5	GO:0071222	cellular response to lipopolysaccharide	2.832E-19	11
	GO:0034378	chylomicron assembly	2.732E-14	5	GO:0032496	response to lipopolysaccharide	3.346E-19	12
	GO:0034370	triglyceride-rich lipoprotein particle remodeling	2.167E-13	5	GO:0071219	cellular response to molecule of bacterial origin	3.824E-19	11
	GO:0006954	inflammatory response	4.636E-13	11	GO:0002237	response to molecule of bacterial origin	5.520E-19	12
	GO:0032496	response to lipopolysaccharide	6.927E-13	9	GO:0071216	cellular response to biotic stimulus	1.345E-18	11
	GO:0002237	response to molecule of bacterial origin	1.004E-12	9	GO:0009617	response to bacterium	3.439E-15	12
	GO:1901700	response to oxygen-containing compound	2.589E-12	13	GO:0006954	inflammatory response	6.326E-15	12
	GO:0006952	defense response	2.795E-12	13	GO:0019221	cytokine-mediated signaling pathway	7.120E-15	12
	GO:0033993	response to lipid	4.803E-12	11	GO:0034097	response to cytokine	2.594E-14	13
	GO:0034368	protein-lipid complex remodeling	1.532E-11	5	GO:0071396	cellular response to lipid	4.713E-14	11

Pathway enrichment analysis

Our pathway enrichment analysis has shown that the hub genes detected for lopinavir/ritonavir fixed combination could be grouped into 2 main molecular pathways, IL-10 signaling and chylomicron assembly,. When exploring the lopinavir/ritonavir interacting hub genes and corresponding pathways, we detected 5 distinct clusters: IL-10 signaling, viral protein interaction with cytokine and cytokine receptor, Legionellosis, RANKL/RANK signaling, and mRNAs involvement in the immune response in sepsis.

Associated diseases analysis

Finally, we detected associated diseases (marker/mechanism and/or therapeutic) for both, fixed combination and lopinavir/ritonavir mixture (Table 2 ), among which were cardiovascular, respiratory diseases, immune system related disorders, skin diseases, to name a few. While apolipoprotein-coding genes were critical for development of distinct cardiovascular diseases, inflammatory genes have been reported as mechanistically related to, cardiovascular, respiratory, and immune system diseases. COVID-19 was seen in the list of hub genes – related diseases in both analyzes.

Table 2

Diseases potentially linked to the use of lopinavir/ritonavir as fixed combination, or when regarded as separate substances in the CTD database (SetAnalyzer CTD tool, p-value: 0.01).

	Disease	Disease Categories	Correctedp-value	Genes (Number)	Annotated Genes
Lopinavir/Ritonavir fixed combination	Cardiovascular Diseases	Cardiovascular disease	3.95E-15	13	APOA1,APOA4,APOB,APOC3,CCL2,CDKN1A,CXCL1,CXCL2,CXCL3, CXCL8,ICAM1,IL6,PTGS2
	Lung Diseases	Respiratory tract disease	1.24E-13	11	APOA1,APOC3,CCL2,CDKN1A|CXCL1,CXCL2,CXCL8,ICAM1,IL6,PTGS2,PTX3
	Reperfusion Injury	Cardiovascular diseasePathology (process)	1.44E-13	8	CCL2,CDKN1A,CXCL1,CXCL2,CXCL8,ICAM1,IL6,PTGS2
	Myocardial Ischemia	Cardiovascular disease	1.59E-12	9	APOA1,APOB,APOC3,CCL2,CXCL2,CXCL8,ICAM1,IL6,PTGS2
	Arteriosclerosis	Cardiovascular disease	4.92E-12	7	APOA1,APOA4,APOC3,CCL2,ICAM1,IL6,PTGS2
	Hypersensitivity	Immune system disease	2.58E-11	8	APOA1,APOC3,CCL2,CXCL1,CXCL2,ICAM1,IL6,PTGS2
	Inflammation	Pathology (process)	5.72E-11	7	APOA1,CCL2,CXCL2,CXCL8,ICAM1,IL6,PTGS2
	Stroke	Cardiovascular diseaseNervous system disease	2.93E-10	6	APOA1,CCL2,CXCL2,ICAM1,IL6,PTGS2
	Hyperlipidemias	Metabolic disease	6.92E-09	5	APOA2,APOA4,APOB,APOC3,ICAM1
	Hyperlipoproteinemia	Metabolic disease	3.29E-08	4	APOA2,APOA4,APOB,APOC3
	⁎COVID19	Respiratory tract diseaseViral disease	0.00606	2	CCL2,IL6
Lopinavir + Ritonavir (SetAnalyzer)	Inflammation	Pathology (process)	1.74E-18	10	CCL2,CCL3,CCL4,CXCL2,CXCL8,ICAM1,IL1B,IL6,TNF,TNFAIP3
	Brain Ischemia	Cardiovascular diseaseNervous system disease	8.69E-18	9	CCL2,CCL3,CCL4,CXCL2,ICAM1,IL1B,IL6,NFKBIA,TNF
	Pneumonia	Respiratory tract disease	1.31E-16	8	CCL2,CCL3,CCL4,CXCL1,CXCL8,IL1B,IL6,TNF
	Immune System Diseases	Immune system disease	2.73E-16	13	CCL2,CCL3,CCL4,CXCL1,CXCL2,CXCL8,ICAM1,IL1B,IL6,NFKB2,NFKBIA,TNF,TNFAIP3
	Reperfusion Injury	Cardiovascular diseasePathology (process)	4.67E-16	9	CCL2,CCL3,CXCL1,CXCL2,CXCL8,ICAM1,IL1B,IL6,TNF
	Vascular Diseases	Cardiovascular disease	2.19E-15	12	CCL2,CCL3,CCL4,CXCL1,CXCL2,CXCL8,ICAM1,IL1B,IL6,NFKBIA,TNF,TNFAIP3
	Myocardial Ischemia	Cardiovascular disease	1.01E-14	10	CCL2,CCL3,CCL4,CXCL2,CXCL8|ICAM1,IL1B,IL6,NFKBIA,TNF
	Shock, Hemorrhagic	Pathology (process)	3.34E-14	5	CXCL1,CXCL2,ICAM1,IL6,TNF
	Hypersensitivity	Immune system disease	1.62E-13	9	CCL2,CCL3,CCL4,CXCL1,CXCL2,ICAM1,IL1B,IL6,TNF
	Glomerulonephritis	Urogenital disease (male and female)	1.70E-12	7	CCL2,CCL3,CCL4,CXCL8,IL1B,IL6,TNF
	COVID19⁎	Respiratory tract diseaseViral disease	3.18e-12	5	CCL2,CCL3,IL1B,IL6,TNF

was not among the top 10 related diseases.

Our further CTD analysis predicted the total of 225 diseases linked to the fixed combination associated hub genes and 295 for the other set of highly connected genes. Among them 191 were in common (58.1%).

Moreover, we constructed a list of diseases to which development might contribute the use of a lopinavir and ritonavir combination. There were 9 such diseases: bipolar disorder, chemical and drug induced liver injury, diabetes insipidus, nephrogenic, Fanconi syndrome, HIV-associated lipodystrophy syndrome, metabolic syndrome, neutropenia, renal insufficiency, status epilepticus.

Chloroquine/hydroxychloroquine + azithromycin combination

Obtained results have shown that chloroquine interacted with 2296, hydroxychloroquine with 45 and azithromycin with 14 genes. Among these genes, chloroquine and azithromycin had a common effect on 8 (BCL2L1, CCL2, CTSB, CXCL8, CYP3A4, IL1B, IL6, TNF) (Fig. 5A). Both drugs increased the activity of CTSB gene and decreased the activity of CYP3A4 gene. Furthermore, they both decreased the expression of IL6 and TNF proteins. On the other hand, hydroxychloroquine and azithromycin had common effect on 6 genes (CCL2, CTSB, CXCL8, IL1B, IL6, TNF) (Fig. 5B). They both decreased the expression of IL6 and TNF protein as well.

Fig. 5

Venn diagram showing the number of genes associated with (A) azithromycin and chloroquine, single and combined; and (B) azithromycin and hydroxychloroquine, single and combined (VennViewer CTD tool).

Furthermore, our results have shown that chloroquine and azithromycin are involved in 238 common biological processes (Fig. 6A) and 2 molecular functions (cytokine activity and cytokine receptor binding) (Fig. 6B), while hydroxychloroquine and azithromycin were involved in 179 common biological (Fig. 6C) processes and 2 molecular functions (cytokine activity and cytokine receptor binding) (Fig. 6D).

Fig. 6

Venn diagram showing the enriched (A) Biological Processes and (B) Molecular functions associated with azithromycin and chloroquine, single and combined; and (C) Biological Processes and (D) Molecular functions associated with azithromycin and hydroxychloroquine, single and combined (VennViewer CTD tool).

A tight network between 8 genes common to chloroquine and azithromycin (Fig. 7A) or tight network between 6 genes common to hydroxychloroquine and azithromycin (Fig. 7B), along with 20 related genes was constructed. In the case of chloroquine and azithromycin combination, the results have shown that more than half of these genes (82.95%) were in co-expression, while 10.82% were in co-localization, which means that they are expressed in the same tissue or their gene products are identified in the same cellular location (Fig. 7A). The rest of the 6.23% of the interactions between genes were predicted by the server, while none of the genes were in physical interactions.

Fig. 7

Tight network between genes common to (A) chloroquine and azithromycin; or (B) hydroxychloroquine and azithromycin, together with 20 related genes (GeneMANIA Cytoscape plug-in). Color legend: co-expression-purple; co-localization-blue; physical interaction-red; pathway-light blue; predicted-orange. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

In the case of hydroxychloroquine and azithromycin combination, the results have shown that more than half of these genes (81.21%) were in co-expression, while 10.98% were in co-localization (Fig. 7B). Similarly, to the chloroquine/azithromycin combination, the rest of the 7.81% of the interactions between genes were predicted by the server, while none of the genes were in physical interactions.

The top 15 hub genes for chloroquine/azithromycin and hydroxychloroquine/azithromycin are presented in Fig. 8A and B, respectively. Obtained hub genes were identical for both combinations, while the intensity of correlation for some of them was different when compared among the two set

Fig. 8

15 hub genes (ranked by the Maximal Clique Centrality (MCC) score) obtained for (A) chloroquine/azithromycin combination; and (B) hydroxychloroquine/azithromycin combination. Color of the nodes represents the intensity of correlation with yellow being the lowest and red the highest (Cytohubba plug-in + MCC method). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Gene ontology connected with our set of 15 hub genes and, thus, potentially linked to the use of the two drug combinations was analysed. Top 10 obtained biological processes and molecular functions are shown in the Table 3 .

Table 3

Top 10 molecular functions and biological processes potentially linked to the use of chloroquine/hydroxychloroquine + azithromycin combination (ToppGene ToppFun function, p-value: 0.05, FDR corrected).

	ID	Name	p-value	Genes from Input	Genes in Annotation
Molecular functions	GO:0005125	cytokine activity	1.116E-16	10	225
	GO:0005126	cytokine receptor binding	3.042E-15	10	312
	GO:0008009	chemokine activity	3.269E-15	7	50
	GO:0042379	chemokine receptor binding	4.319E-14	7	71
	GO:0048018	receptor ligand activity	3.232E-13	10	497
	GO:0030546	signaling receptor activator activity	3.571E-13	10	502
	GO:0001664	G protein-coupled receptor binding	1.450E-9	7	308
	GO:0005102	signaling receptor binding	4.933E-9	11	1817
	GO:0098772	molecular function regulator	1.665E-7	10	1905
	GO:0048020	CCR chemokine receptor binding	6.088E-6	3	47
Biological processes	GO:0071222	cellular response to lipopolysaccharide	9.939E-22	12	221
	GO:0032496	response to lipopolysaccharide	1.340E-21	13	356
	GO:0071219	cellular response to molecule of bacterial origin	1.381E-21	12	227
	GO:0002237	response to molecule of bacterial origin	2.309E-21	13	371
	GO:0071216	cellular response to biotic stimulus	5.471E-21	12	254
	GO:0006954	inflammatory response	3.488E-19	14	806
	GO:0019221	cytokine-mediated signaling pathway	4.008E-19	14	814
	GO:0009617	response to bacterium	3.078E-17	13	766
	GO:0071345	cellular response to cytokine stimulus	8.112E-17	14	1188
	GO:0034097	response to cytokine	2.490E-16	14	1287
	GO:0071396	cellular response to lipid	5.165E-16	12	654

In order to understand the biological importance of the hubs, we performed the pathway analysis. The results have shown that chloroquine/hydroxychloroquine + azithromycin combination interacting hub genes can be grouped into 8 pathways: NOD-like receptor signaling pathway, viral protein interaction with cytokine and cytokine receptor, NF-kappa B signaling pathway, interleukin 10 signaling, Legionellosis, rheumatoid arthritis, IL-17 signaling pathway and IL-18 signaling pathway.

Associated diseases analysis

Top 10 diseases (marker/mechanism and/or therapeutic) associated with our set of 15 hub genes are presented in the Table 4 , sorted by the p-value. The listed diseases include both those for which the investigated drug combinations could be used as a therapy and those which could be induced by the said combinations.

Table 4

Diseases potentially linked to the use of chloroquine/hydroxychloroquine + azithromycin combination (SetAnalyzer CTD tool, p-value: 0.01).

Disease	Disease Categories	Corrected p-value	Number of genes	Annotated Genes
Inflammation	Pathology (process)	3.24e-21	11	CCL2, CCL3, CCL4, CXCL2, CXCL8, ICAM1, IL1B, IL6, PTGS2, TNF, TNFAIP3
Brain Ischemia	Cardiovascular diseaseNervous system disease	1.28e-20	10	CCL2, CCL3, CCL4, CXCL2, ICAM1, IL1B, IL6, NFKBIA, PTGS2, TNF
Reperfusion Injury	Cardiovascular diseasePathology (process)	1.09e-18	10	CCL2, CCL3, CXCL1, CXCL2, CXCL8, ICAM1, IL1B, IL6, PTGS2, TNF
Cardiovascular Diseases	Cardiovascular disease	1.15e-17	13	CCL2, CCL3, CCL4, CXCL1, CXCL2, CXCL3, CXCL8, ICAM1, IL1B, IL6, NFKBIA, PTGS2, TNF, TNFAIP3,
Cerebrovascular Disorders	Cardiovascular diseaseNervous system disease	2.12e-17	10	CCL2, CCL3, CCL4, CXCL2, ICAM1, IL1B, IL6, NFKBIA, PTGS2, TNF
Myocardial Ischemia	Cardiovascular disease	4.58e-17	11	CCL2, CCL3, CCL4, CXCL2, CXCL8, ICAM1, IL1B, IL6, NFKBIA, PTGS2, TNF
Skin and Connective Tissue Diseases	Skin disease	1.03e-16	14	CCL2, CCL4, CXCL1, CXCL2, CXCL3, CXCL8, HCAR3, ICAM1, IL1B, IL6, NFKBIA, PTGS2, TNF, TNFAIP3
Immune System Diseases	Immune system disease	2.93e-16	13	CCL2, CCL3, CCL4, CXCL1, CXCL2, CXCL8, ICAM1, IL1B, IL6, NFKBIA, PTGS2, TNF, TNFAIP3
Hypersensitivity	Immune system disease	7.40e-16	10	CCL2, CCL3, CCL4, CXCL1, CXCL2, ICAM1, IL1B, IL6, PTGS2, TNF
Skin Diseases	Skin disease	1.19e-15	13	CCL2, CCL4, CXCL2, CXCL3, CXCL8, HCAR3, ICAM1, IL1B, IL6, NFKBIA, PTGS2, TNF, TNFAIP3
COVID19a	Respiratory tract diseaseViral disease	3.41e-12	5	CCL2, CCL3, IL1B, IL6, TNF

Was not among the top 10 related diseases.

However, results obtained by VenViewer CTD tool have shown that the therapeutic effect of both chloroquine/azithromycin and hydroxychloroquine/azithromycin is connected with their effect on 4 genes (CCL2, IL6, TNF and IL1B1) directly associated with COVID-19 disease. Both combinations reduced the expression of IL-6 and TNF proteins. These genes were also listed in inflammation and immune system diseases.

Further VennViewer CTD analysis was conducted to obtain a list of diseases to which development might contribute the use of a combination of chloroquine and azithromycin or hydroxychloroquine and azithromycin. There were 22 such diseases for chloroquine/azithromycin combination: acute kidney injury, arthralgia, delirium, depressive disorder, dizziness, gastrointestinal diseases, headache, hearing disorders, hearing loss, heart arrest, heart block, lipidoses, long qt syndrome, muscular diseases, myasthenia gravis, nausea, pain, pruritus, psychoses (substance-induced), tachycardia (ventricular), torsades de pointes, vasculitis (leukocytoclastic, cutaneous). While, for hydroxychloroquine/azithromycin combination there were 11 common diseases, including dizziness, exanthema, gastrointestinal diseases, headache, hearing loss (sensorineural), heart arrest, heart block, long qt syndrome, muscular diseases, nephritis (interstitial), Stevens-Johnson syndrome and torsades de pointes.

Results of the CTD analysis also reveal that acute renal failure, depression, headache, nausea, pain, and tachycardia (ventricular) could be linked to simultaneous decreased expression of TNF and IL-6 proteins caused by the investigated drug combinations.

Discussion

Generally, it has been shown that drug combinations are more effective against certain viruses, such as HIV and SARS-CoV-2, in comparison with the individual ones (Muralidharan et al., 2020). Lopinavir and ritonavir combination has been proposed to control the virulence to a great extent in the COVID-19 affected patients within 48 h. It has even been suggested that these two drugs could be accompanied with other antiviral agents to enhance antiviral effects and improve clinical outcomes (Cao et al., 2020). For example, molecular docking and molecular dynamics simulation studies have shown that combination of lopinavir, oseltamivir and ritonavir were found to be highly effective against SARSCoV-2 protease (Muralidharan et al., 2020). Furthermore, the combination of hydroxychloroquine and azithromycin on SARS-CoV-2 infected cells was tested in vitro, showing that there was a considerable synergy of these two substances applied at doses correlating the concentrations likely to be obtained in humans (Gautret et al., 2020b). Epidemiological studies have also shown that hydroxychloroquine and azithromycin therapy was far more successful compared to untreated patients or treatment with hydroxychloroquine alone (Liu et al., 2020a, Liu et al., 2020b). The benefits of azithromycin in the treatment of viral infections can be supported by a limited number of data indicating inhibitory activity for Zika virus proliferation (Retallack et al., 2016), as well as a decrease in inflammation in a murine model of viral bronchiolitis after azithromycin administration (Beigelman et al., 2010).

However, since balance between benefit and negative effects of drugs is sometimes hard to achieve, we considered both risks and benefits of proposed anti-COVID-19 therapy in order to understand the molecular pathways and biological processes behind drugs-induced side effects, which may be crucial for their management and prevention.

We attempted to understand the potential of lopinavir/ritonavir combination to cause side effects or complications of already existing diseases by applying system biology approaches and comprehensive data-mining. Our results obtained from the CTD show that fixed ritonavir/lopinavir combination interacts with different set of genes in comparison to the mixture of these two drugs signalizing the limitation of the performed research. Bearing in mind that the analysis is heavily dependent on the underlying online sources from which the annotations are retrieved, interaction data might be incomplete and unreliable (Chen et al., 2009). Thus, higher number of chemical-gene interactions could be obtained when antiviral drugs were regarded as separate compounds than in the case of the less investigated fixed combination, for which only key interacting genes were reported. However, both analysed options showed common effects on IL6, ICAM1, TNF, and CCL2, genes involved in the process of inflammation. While IL6 and TNF code proinflammatory molecules, CCL2 stimulates the migration of inflammatory cells to the place of inflammation where they enter the tissue with the help of the endothelial and leukocyte-associated transmembrane protein, coded by ICAM1 (Wang et al., 2012). In the CTD database, CCL2, CCL3, IL1B, TNF and IL6 were recognized as lopinavir/ritonavir – modulated genes involved in the COVID19 treatment, leading to a suspicion that these proteins may be critical for some of the therapeutic properties of the given drugs. Increased expression of the pro-inflammatory cytokines, CCL2 and IL-6, triggers the activation of downstream pathways such as NOD-like receptor signaling pathway, involved in the detection of various pathogens and generation of innate immune responses (Davis et al., 2011; Philpott et al., 2014). Furthermore, IL-17, found in the molecular pathways triggered by the investigated drug combination, is known to connect innate and adaptive immunity. It activates the induction of IL-6 and TNF-alpha, while CCL2 is among its signature genes. On the other hand, IL-10 family cytokine activates many of the same innate inflammatory genes as IL-17 (Gaffen, 2008). In other words, these lopinavir/ritonavir – IL-17/IL-10/CCL2/IL6 interactions appear as important stimulators of the host immunity and pyroptosis and, thus, viral clearance (Hachim et al., 2020). However, their over-activation might cause inflammation – mediated side effects or induce AGE-RAGE signaling pathway in diabetic complications and might elicit the activation of multiple intracellular signaling pathways, including NF-kappaB activity (Nguyen et al., 2012). NF-kappa B further promotes the expression of pro-inflammatory cytokines and a variety of atherosclerosis-related genes, such as VCAM-1, tissue factor, VEGF, and RAGE (Santoro et al., 2003). Thus, it seems that antiretroviral drugs used in SARS-CoV-2 infection may exert both beneficial and adverse effects by interacting with the set of common genes making it harder to distinguish between their pharmacological characteristics and possible toxic effects.

Hence, the increased expression of all the mentioned genes by lopinavir/ritonacir combination, could subsequently stimulate inflammatory response in the treated patients and cause serious side effects shown in the Table 2, such as brain ischemia, reperfusion injury, myocardial ischemia or hemorrhagic shock, particularly in those patients who have other risk factors and secondary diseases. Furthermore, we could even speculate that the proposed therapy could increases the risk of severe complications in SARS-CoV-2-infected patients, the excessive production of proinflammatory cytokines called ‘cytokine storm’, which might be associated with the severity and outcome of COVID-19 disease (Conti et al., 2020; Tufan and Matucci-Cerinic, 2020).

In the next step of our analysis, after generating gene networks for both treatment options, we predicted hub genes that were then clustered based on their biological functionality into groups. In both data sets, first group of genes was involved in regulation of IL-10 signaling pathways and its sub-pathways such as NF-kB signaling, IL-17 signaling, and IL-4 and IL-13 signaling, to name a few. While fixed lopinavir/ritonavir combination modified their activity through IL6, CXCL8, CXCL2, CCL2, PTX3, APOA1, ICAM1, CXCL1, CXCL3, CDKN1A, APOB and PTGS2, combined use of these drugs may interact with NFKBIA, CXCL1, ICAM1, IL1B, CXCL2, TNF, CCL3, CCL4, CXCL8 and IL-6 to cause the change in IL-10 signaling.

As explained, IL-10 pathways represent the link between the investigated drugs and possible development of immune-related side effects. Indeed, our pathway enrichment analyzes predicted inflammation, hypersensitivity, glomerulonephritis, and other immune system related diseases among 10 most probable ritonavir/lopinavir – induced diseases. In other words, this drug combination may cause higher risk than benefit in patients suffering from autoimmune diseases such as rheumatoid arthritis, asthma or similar, possibly through identified biological processes such as cytokine-mediated signaling pathway and response to cytokines.

On the other hand, there were some differences between the data reported for fixed ritonavir/lopinavir combination and their mixture. Among the top 10 interacting genes, fixed antiviral combination increased the expression of 3 apolipoprotein coding genes: APOA1, APOB, and APOC3, while VennViewer CTD tool has shown that ritonavir and lopinavir commonly interact with 2 genes involved in lipid metabolism, ADIPOQ and LEP, by decreasing their activity. Thus, in the first examined set, enriched pathways were grouped under chylomicron assembly category, while the second one clustered around viral protein interaction with cytokine and cytokine receptor, Legionellosis, RANKL/RANK signaling, and mRNAs involvement in the immune response in sepsis. Subsequently, this led to difference in the top 10 predicted drug-induced diseases.

Fixed lopinavir/ritonavir combination was associated with the risk of developing lipid metabolism disorders such as hyperlipidemia, dyslipidemia or hyperlipoproteinemia, most likely by influencing biological processes such as chylomicron remodeling, chylomicron assembly, triglyceride-rich lipoprotein particle remodeling, and forming of protein-lipid complex. This could also be seen through the impaired molecular functions such as changes in binding of high-density lipid molecules or cholesterol to their receptors, or phosphatidylcholine-sterol O-acyltransferase activator activity.

Similarly, Liang et al. (2001) have reported that the human immunodeficiency virus (HIV) protease inhibitors can be associated with hyperlipidemia through inhibition of proteasomal degradation of nascent apolipoprotein B, the principal protein component of triglyceride and cholesterol-rich plasma lipoproteins. Dysregulation of lipid metabolism is a known risk factor for cardiovascular events (Maggi et al., 2017), which correlates with the results obtained in this study. In summary, as seen in HIV-infected patients, the use of lopinavir/ritonavir, alone and especially in combination, may contribute to the development of a lipid profile known to increase the risk of coronary heart disease (Friis-Møller et al., 2003). Therefore, based on result reported here, this combination should be used with caution in SARS-CoV-2-infected patients with hereditary or acquired metabolic diseases, as well as in patients with cardiovascular problems.

Chloroquine/hydroxychloroquine + azithromycin combination

The experimental studies have demonstrated that both chloroquine phosphate and hydroxychloroquine sulfate inhibit SARS-CoV-2 virus (Gautret et al., 2020b). Apart from the demonstrated ability of chloroquine to inhibit growth of SARS-CoV-2, an early clinical trial conducted in COVID-19 Chinese patients showed that this drug had a significant effect on clinical outcome (Gautret et al., 2020a), by reducing fever, improving CT imaging, and delaying disease progression (Gautret et al., 2020b).

However, hydroxychloroquine is considered to have lower toxicity and greater safety than chloroquine, especially when administrated at higher doses (Chatre et al., 2018; J. Liu et al., 2020).

Our results have demonstrated that both chloroquine/azithromycin and hydroxychloroquine/azithromycin combinations interacted with 6 genes (CCL2, CTSB, CXCL8, IL1B, IL6, and TNF), while chloroquine and azithromycin interacted with two additional genes (BCL2L1 and CYP3A4). Further genes taken into account in our investigations were the hub genes obtained by the Cytoscape Cytohubba plug-in. We noted that the same hub genes were generated for both drug combinations, which is understandable having in mind that chloroquine and hydroxychloroquine share similar chemical structures and mechanisms of acting as a weak bases and immunomodulators (J. Liu et al., 2020; Yao et al., 2020). However, the intensity of correlation for some of them was different when compared among the two sets. For chloroquine/azithromycin combination the correlation was the highest with ICAM1, IL1B, IL6, and HCAR3, while for hydroxychloroquine/azithromycin it was the highest with TNF and PTGS2. This difference might be behind the less pronounced adverse effects observed in the case of hydroxychloroquine treatment compared to chloroquine. However, our investigation has shown that, in the case of chloroquine/hydroxychloroquine + azithromycin combinations, it is difficult to draw a line between adverse and beneficial effects, since some of the genes (CCL2, IL6, TNF, and IL1B1) are involved not only in the ethiology of their toxic effects, but also in therapeutic action against COVID-19.

This could also be seen in the list of top 10 diseases connected with the use of chloroquine/azithromycin and hydroxychloroquine/azithromycin combinations conducted by our CTD analysis. Inflammation was listed in the first place, having in mind that both chloroquine and hydroxychloroquine are anti-inflammatory agents that, other than their direct antiviral activity, can significantly decrease the production of cytokines, in particular pro-inflammatory factors (J. Liu et al., 2020). Indeed, cytokine and chemokine activity and binding to combined receptors were seen as the top molecular functions impacted by the investigated drug combinations. Consequently, cellular response to cytokine stimulus and cytokine signaling pathways might be dysregulated, which could be regarded as mechanistic basis of their possible adverse reactions. Hypersensitivity in the immune system category was among the top diseases, which is expected having in mind that these two drugs are also used to treat various immune system diseases, such as rheumatoid arthritis and lupus erythematosus (Savarino et al., 2003; Wozniacka et al., 2006). Hence, rheumatoid arthritis was also present among our enriched pathways list. However, all of the other top 10 diseases listed in our study included cardiovascular diseases, which might be linked to the potential adverse effects of the said combinations, with 8 genes present in all of the listed cardiovascular diseases (CCL2, CCL3, CXCL2, ICAM1, IL1B, IL6, PTGS2, TNF).

Out of these genes, it should be noted that both CCL2 and ICAM1 are considered important in the pathogenesis of atherosclerosis (Khyzha et al., 2019; Motawi et al., 2012). CCL2 could be linked to unstable symptomatic human atherosclerotic plaques, with particularly important role both in post-ischemic and post-reperfusion inflammation. This gene modulates fibroblasts and endothelial cells functions and can be connected to myocardial infarction, thrombotic occlusion and myocardial ischemia (Khyzha et al., 2019). On the other hand, ICAM1 plays a crucial role in lymphocyte migration and activation, interacts with leukocyte integrins and promotes the atherosclerotic process at the surface of endothelial cells (Motawi et al., 2012).

Our data-mining results indicate that both chloroquine and hydroxychloroquine increase the expression of CCL2 mRNA. However, azithromycin acts in the opposite way, decreasing the expression of CCL2 protein, which might indicate that the use of this drug along some of the other two might act beneficially in the light of possible cardiovascular adverse effects. Yet, there was not enough data about ICAM1 to draw a conclusion about the possible synergistic or antagonistic effects of the drugs from the investigated combinations. This is also the case with CCL3, another gene closely associated with myocardial ischemia. Nevertheless, in their toxicogenomic data-mining study, Chen et al. (2019) have found that this gene was enriched Toll-like receptor signaling pathway, which includes immune response through circulating blood cells, increasing infarct size and influencing ventricular remodeling.

In the process of CTD analysis, we also generated a list of other diseases to which development might contribute the treatment with chloroquine/hydroxychloroquine + azithromycin combination. It was noted that fewer diseases were present in the case of chloroquine/azithromycin combination than hydroxychloroquine/azithromycin. Diseases common to all the investigated substances included dizziness, gastrointestinal diseases, headache, hearing loss, heart arrest, heart block, long qt syndrome, muscular diseases and torsades de pointes. Additional diseases connected with chloroquine/azithromycin, but not hydroxychloroquine/azithromycin, were acute kidney injury, arthralgia, delirium, depressive disorder, lipidoses, myasthenia gravis, pruritus, psychoses (substance-induced), tachycardia, (ventricular) and vasculitis (leukocytoclastic, cutaneous).

Predictably, our further GO enrichment analysis of the generated hub genes set has demonstrated that the majority of their molecular processes were linked to cytokine and chemokine activation process.

The pathway enrichment analysis has also demonstrated that chloroquine/hydroxychloroquine + azithromycin combination hub genes could be grouped into several pathways, including NOD-like receptor signaling pathway, viral protein interaction with cytokine and cytokine receptor, NF-kappa B signaling pathway, IL-10 signaling pathway, IL-17 signaling pathway and IL-18 signaling pathway.

All of these pathways are mutually linked and connected to the immune system response. For example, upon recognition of pathogen-associated molecular patterns, NOD-like receptors activate NF-κB or MAP kinases to induce the production of inflammatory cytokines such as IL-1β and IL-18 (Kumar et al., 2011). Results of a study in which cytokines were measured in blood plasma of 12 COVID-19 patients have shown that 2019-nCoV viral load was highly positively associated with the plasma levels of both IL-10 and IL-17 and could be used as markers to predict COVID-19 disease severity (Y. Liu et al., 2020). IL-17 was among the cytokines found abundant in COVID-19 patients and associated with the severe lung inflammation. This cytokine promotes the induction of other cytokines, such as IL-1β, IL-6 and TNF-α, as well as growth factors, G-CSF; chemokines; and matrix metalloproteinases (Casillo et al., 2020; Tufan and Matucci-Cerinic, 2020). Chemokines upregulated by IL-17 also include CXCL1, CXCL2, CXCL5, CXCL8 (IL-8), CXCL9, CXCL10, CCL2, and CCL20 (Xie et al., 2011), many of them listed as hub genes in our study.

In general, studies have shown that immunomodulatory effect of chloroquine and hydroxychloroquine may be useful in controlling the cytokine storm that occurs in the late-phase of critically ill SARS-CoV-2 infected patients (Yao et al., 2020). Both of the investigated drug combinations decreased the expression of IL-6 and of TNF proteins and thus, suppressed generated immune response. According to the CTD base, this could be further linked to the therapeutic action against COVID-19. However, it should also be considered that both IL-6 and TNF are components of a complex system which consists of many positive and negative regulators. IL-6 has both pro-inflammatory and anti-inflammatory properties. It is produced not only by immune cells, as well as monocytes and macrophages, but also by cardiovascular components, like endothelial cells, vascular smooth-muscle cells, and ischemic myocytes. Thus, IL-6 is involved not only in inflammation but also in the regulation of cardiac metabolism (Kanda and Akahashi, 2004). TNF is also considered a double-natured cytokine which may express protective and adaptive effects in cases of hemodynamic overload and ischemia-reperfusion (Bellisarii et al., 2001), since TNF signaling pathway has exceedingly been associated with cardiac remodeling following myocardial infarction (Chen et al., 2019). Accordingly, decreasing IL6 and TNF protective activity might also lead to the possible adverse cardiovascular effects present in the case of chloroquine/hydroxychloroquine + azithromycin. Furthermore, PTSH2 gene, enriched in TNF signaling pathway, was also listed as one of the genes that might influence the occurrence of possible cardiotoxic effects of the investigated combinations, especially since it has been found to be associated with the decreasing risk of stroke (Chen et al., 2019).

Additionally, our results have shown that both chloroquine and azithromycin decreased the activity of CYP3A4 gene. CYP3A4 is the most abundant metabolic enzyme in human liver, regulating more than 50% drug metabolism (Wang et al., 2019). CYP3A4 and CYP3A5 play important roles in chloroquine metabolism and are the major enzymes responsible for the chloroquine N-demethylation into N-desethylchloroquine in human liver microsomes (Kim et al., 2003). Inhibitory effect on the activity of this enzyme would mean the inhibition of metabolism and elimination of these drugs from the organism, thus increasing the risk of their possible toxic effects. However, CYP3A4 was not listed among the interactions for hydroxychloroquine. This could also be a reason why hydroxychloroquine/azithromycin combination could be regarded safer than chloroquine/azithromycin.

Another adverse event observed in our investigation included skin changes and muscular diseases. However, these effects are associated with the long-term use of chloroquine/hydroxychloroquine therapy (Chatre et al., 2018) and, again, may be explained through over-activation of cytokine activity and immune response.

Limitations

Our toxicogenomic data mining approach has proven useful for exploring and gaining better understanding of the safety profile of the selected COVID19 drug combinations. However, freely available tools such as CTD, various Cytoscape plugins and ToppGene Suite portal, have some boundaries which confirm that in silico toxicogenomic investigations could be used only as an addition to other in vitro and in vivo methods in the overall toxicity testing process. Having in mind that the data in such functional annotation-based prioritization processes is extracted by drawing statistical associations between genes and drugs of interest, it is not possible to directly address the dose-response relationship (Harris et al., 2020). Furthermore, the analysis deeply relies on the online sources the annotations were obtained from (in this case, CTD database), as well as the quality of interaction data present in them (Chen et al., 2009). It is often necessary to omit certain data due to the insufficient number of gene-drug annotations, and sometimes the analysis is not even possible, which is the case with some other drugs used in COVID-19 treatment, such as remdesivir and interferon beta. It should also be recognized that not only dose, but many other factors might influence the manifestation of drug toxic effects, such as route of administration, duration, individual drug metabolism rates, etc. (Davis et al., 2008), which should be explored in further throughput investigations.

Conclusion

In this study we used toxicogenomic data-mining approach to understand the safety profile of two selected COVID-19 treatment strategies. The use of available databases and software helped us predict the molecular mechanisms of adverse reactions which may be useful when proposing the side effects management in the next step of our research. Chosen drug combinations are seen to interact with genes (TNF, IL6, CCL2, and IL1B1) reported to be potential therapeutic targets in COVID-19 disease. However, while ritonavir/lopinavir increased their expression, chloroquine/hydroxychloroquine + azithromycin downregulated their activity, which was confirmed by the anti-inflammatory characteristics of the second drug combination. Obtained results highlighted the development of cardiovascular events or immune-related reactions as the highest risks for SARS-CoV-2-infected patients treated with antiretroviral drugs (lopinavir/ritonavir) or chloroquine/hydroxychloroquine + azithromycin combination. Based on our data, the use of these medications might lead to the severe complications, particularly in patients with comorbidities such as dyslipidemias, unstable cardiovascular or inflammatory diseases. Moreover, detected hub genes in both analyze were directly and indirectly linked to brain or myocardial ischemia, reperfusion injury and different respiratory problems. While clinical results show benefit of the proposed treatment options in SARS-CoV-2-infected patients, our results imply special caution when administering the patients with chronic illnesses. However, further in vitro and in vivo investigations are needed to confirm or infirm these findings.

Authors' contributions

All authors contributed to the study conception and design.

Katarina Baralić: Conceptualization, Formal analysis and investigation, Methodology, Software, Writing – original draft preparation Dragica Jorgovanović: Conceptualization, Formal analysis and investigation, Methodology, Software, Writing – original draft preparation Katarina Živančević: Conceptualization, Formal analysis and investigation, Methodology, Software, Writing – original draft preparation Evica Antonijević: Toxicity data review Biljana Antonijević: Supervision, Funding Acquisition Aleksandra Buha-Djordjevic: Writing – review and editing, Supervision Marijana Ćurčić: Toxicity data review, Writing – review and editing Danijela Đukić-Ćosić: Writing – review and editing, Supervision.

All authors have read and approved the final manuscript.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

ABCB11	ATP binding cassette subfamily B member 11
ADIPOQ	adiponectin, C1Q and collagen domain containing
APOA1	apolipoprotein A1
APOB	apolipoprotein B
APOC3	apolipoprotein C3
BCL2L1	BCL2 like 1
CCL2	C-C motif chemokine ligand 2
CCL3	C-C motif chemokine ligand 3
CDKN1A	cyclin dependent kinase inhibitor 1A
COVID-19	Corona Virus Disease 2019
CTSB	cathepsin B
CXCL1	C-X-C motif chemokine ligand 1
CXCL2	C-X-C motif chemokine ligand 2
CXCL5	C-X-C motif chemokine ligand 5.
CXCL8	C-X-C motif chemokine ligand 8
CXCL9	C-X-C motif chemokine ligand 9.
CYP3A4	cytochrome P450 family 3 subfamily A member 4
CCL20	C—C motif chemokine Ligand 20.
G-CSF	granulocyte-colony stimulating factor.
GLB1	galactosidase beta 1
GM-CSF	granulocyte-macrophage colony-stimulating factor.
HMOX1	heme oxygenase 1
ICAM1	intercellular adhesion molecule 1
IL10	interleukin 10
IL1B	interleukin 1 beta
IL6	interleukin 6
IL-12	interleukin 12
IL-13	interleukin 13
IL-14	interleukin 14
IL-15	interleukin 15
IL-17	interleukin 17
IL-17A	interleukin 17A
IL-18	interleukin 18
IFNα/β	interferon alfa/beta
IFNγ	interferon gamma
LEP	leptin
LMNA	lamin A/C
NF-kB	Nuclear factor kappa-light-chain-enhancer of activated B cells
NOD-like receptors	Nucleotide-binding oligomerization domain-like receptors
MAP	mitogen-activated protein.
RB1	RB transcriptional corepressor 1
RANK	Receptor activator of nuclear factor κB
RANKL	Receptor activator of nuclear factor kappa-Β ligand
SARS-CoV-2	Severe acute respiratory syndrome coronavirus 2
TNF	tumor necrosis factor
TP53	tumor protein p53
ZMPSTE24	zinc metallopeptidase STE24