An in silico analysis of drugs potentially modulating the cytokine storm triggered by SARS-CoV-2 infection

The ongoing COVID-19 pandemic is one of the biggest health and societal challenges of the recent decades. Among the causes of mortality triggered by SARS-CoV-2 infection, the presence of an inflammatory "cytokine storm" (CS) at later stages of the disease has been found to play a determinant role. Here, we used available transcriptomic data from the bronchoalveolar lavage fluid (BALF) of COVID-19 patients suffering from a CS to obtain gene-signatures associated to this pathological process. Using these signatures, we interrogated the Connectivity MAP (CMap) dataset that contains the effects of over 5,000 small molecules on the transcriptome of human cancer cell lines, and looked for molecules which effects on transcription mimic or oppose those associated to the CS. Consistent with their medical use, this analysis found a significant enrichment of glucocorticoids or inhibitors of the Janus Kinases (JAK) as drugs that could revert the CS. On the other hand, molecules that potentiate the immune response such as PKC activators are predicted to worsen the CS. Besides these expected findings, our analysis also reveals a potential effect of the antibiotic doxycycline or MAPK/RAF/MEK inhibitors in reverting the CS, or of topoisomerase inhibitors and the anti-alcohol abuse drug disulfiram in potentiating its effects. Finally, our analyses support that the gender-related differences in the severity of COVID-19 are related to the anti-inflammatory properties of female hormones. While acknowledging that this is an analysis based on limited available data, we decided to share it as a resource that might help others in the selection of drugs that could be tested in the context of experimental models of CS triggered by viral infections.


INTRODUCTION
Since the first patient was hospitalized in Wuhan on December 12 th 2019 (1), the ongoing COVID-19 pandemic has caused 6,288,771 confirmed infections and 379,941 deaths by June 3 rd 2020 (2). The disease-causing pathogen is SARS-CoV-2, a novel positive sense single-stranded RNA coronavirus belonging to the lineage B of the betacoronavirus genus (3). Infections by different coronavirus strains have been historically present in humans, although they most often only cause a mild cold (e.g. infections by HKU, 229 and hCOV-OC43 coronaviruses) (4). However, several new coronavirus strains have emerged in the last 2 decades that cause an Acute Respiratory Distress Syndrome (ARDS) with severe health consequences including the patients´ death. These include SARS-CoV, which between 2002 and 2003 led to 8,000 cases with a mortality of 9.5% and MERS-CoV, which in 2012 caused 2,500 confirmed infections with a mortality rate of 36% (5). As of today, and while there is significant variability between countries and final data are still to be fully defined, the averaged worldwide mortality associated to SARS-CoV-2 infections is around 6% (6).
Due to the lack of a curative treatment for COVID-19, overwhelming efforts are being dedicated to the development therapies including vaccines (7), plasma transfusions from convalescent patients (8) and antiviral treatments limiting viral replication (9) or infection (10,11). Nevertheless, and while respiratory failure associated to ARDS is the leading cause of mortality in COVID-19 patients, accumulating evidence shows that the lethality in a subgroup of the severe patients occurs due to the late onset of an inflammatory "cytokine storm" (CS) (12). First defined in 1993 in the context of graft-versus-host disease (13), the CS has been observed in a wide range of inflammatory diseases such as multiple sclerosis or pancreatitis, and in the context of viral infections including SARS-CoV (14,15). In COVID-19, the CS is caused by an overproduction of proinflammatory cytokines such as IL-6, IL-1b and TNFa by lung infiltrating alveolar macrophages, which triggers vascular hyperpermeability, multiorgan failure and ultimately cell death (16). Potential treatments for the CS of COVID-19 patients include broad spectrum anti-inflammatory therapies such as corticosteroids or more targeted therapies like anti-IL-6 receptor antibodies (Tocilizumab) (17). Other therapies like Janus Kinase (JAK) or IL-1 inhibitors have also been proposed as treatments for the CS based on their potent effects in suppressing key inflammatory signaling pathways driven by IFNg or NFkB, respectively, and are being evaluated in clinical trials (18). Importantly, the use of antiinflammatory therapies such as corticosteroids is only recommended in the severe cases undergoing a CS and not in mild cases of COVID-19, as the use of antiinflammatory therapies in earlier stages of the disease would limit the efficacy of the immune system in fighting the infection.
Databases containing transcriptional signatures of human cancer cell lines exposed to a specific perturbation (perturbagen) such as drugs or genetic manipulations have become an excellent platform for data mining in biomedical research (19)(20)(21). For instance, the Connectivity Map (CMap) from the Broad Institute at MIT stores over 1,5M signatures from a wide range of human cancer cell lines exposed to around 5,000 drugs and 3,000 genetic perturbations (overexpression or shRNA-mediated depletion) (https://www.broadinstitute.org/connectivity-map-cmap) (22,23). These signature databases can be used, for instance, to identify the mechanism of action of a drug, by comparing the transcriptional signature induced by the compound to that induced by genetic perturbation of an individual gene (23). Another relevant use of CMap is to identify drugs that induce a transcriptional signature which negatively correlates with that associated to a given disease as a strategy for discovering potential new therapies (24). In the context of the COVID-19 pandemic, the CMap has been recently used to identify drugs capable of modifying the levels of ACE2 (25), the receptor that SARS-CoV-2 uses for viral entry (10). Here, and based on recently available transcriptomic data from COVID-19 patients, we used CMap to identify drugs that could potentially alleviate or aggravate the severity of the CS associated to late stages of this disease.

Biological pathways related to the hypercytokinemia found in COVID-19 patients
In order to define a transcriptional signature related to COVID-19, we used recently published transcriptomics data from bronchoalveolar lavage fluid (BALF) cells where COVID-19 patients (n=8) were compared to healthy individuals (n=20) and community-acquired pneumonia patients (n=146) (26). Consistent with clinical observations, this study found that COVID-19 patients present a distinct activation of an IFN-dependent cytokine response. We thus defined a transcriptional signature based on the cytokines and cytokine receptor genes that were differentially expressed in COVID-19 patients (COVID CS ) (signature available in Table S1). Using this signature and the clue.io tool from CMap (http://clue.io) we first looked into biological pathways (CMap classes) that have transcriptional signatures that correlate with the COVID CS . Consistent with their key roles in inflammation (27), these analyses identified activation of NFkB and PKC signaling as the CMap classes most significantly similar to the COVID CS ( Fig. 1A and Table S2). Of note, PKC activation triggers NFkB and TNFa-dependent inflammatory responses in human bronchial epithelial cells, suggesting that PKC could act early on in triggering the inflammatory response (28). In what regards for genes that, when downregulated, trigger a signature similar to the COVID CS the list was not as clearly restricted to immune factors as the one from overexpressed genes, although it included anti-inflammatory factors such as NFkB inhibitor alpha (NFKBIA) (Fig. 1C and Table S4). Interestingly, Gene Ontology (GO) analyses of genes whose transcriptional signature opposes that of the COVID CS detected an enrichment in pathways related to female sexual hormones (Fig. 1D). We find this of interest since even if male and females are equally infected with SARS-CoV-2 the fatality rate is significantly higher in men, which could be related to genderrelated differences in the intensity of inflammatory responses (29). In fact, several studies have indicated that estrogen levels influence the severity of diseases linked to inflammation, including cancer, due the modulation of NFkB signaling (30)(31)(32)(33).
Moreover, gender differences in COVID-19 mortality rates are highest in younger patients and decline progressively, further supporting a connection with female hormone levels (34). Noteworthy, estrogen treatments are already been explored in clinical trials for reducing the severity of COVID-19 in men and women 55 and older (as estrogen levels decline after menopause) (35).

Analysis of drugs triggering transcriptional signatures that mimic the COVID CS
Next, we interrogated CMap for the identification of compounds triggering a similar effect on transcription as the one observed in COVID-19 patients, as these compounds could in principle potentiate the severity of the hypercytokinemia. Before discussing the effect of specific compounds, we adapted the Gene Set Enrichment Analysis (GSEA) to conduct a "drug GSEA" in order to identify compound classes that mimic the COVID CS (36). Consistent with the analysis from genetic perturbations, PKC activating drugs were the most significantly enriched class, followed by "ATPase inhibitors", "Topoisomerase Inhibitors", "Insulin sensitizers" and "CDC inhibitors" ( Fig.   2A,B). Interestingly, and in relation to insulin, a recent study has found that the triglyceride and glucose index (TyG) correlates with the severity and mortality of COVID-19 patients, which builds upon epidemiological observations that reveal obesity as a risk factor among COVID-19-related fatalities (37). If we restrict our CMap analysis to medically approved drugs, only Topoisomerase inhibitors showed a significant enrichment as drugs mimicking the COVID CS , although the absence of other compound classes from this list can be influenced by low number of clinically available drugs in other pathways (Fig. 2C).
In what regards to specific drugs that could potentiate the cytokine response in COVID-19 patients, three of the top compounds were PKC activators: prostratin, phorbol-12-myristate-13-acetate (PMA) and ingenol ( Fig. 2D and Table S5). The hydroxychloroquine-related antimalarial drug (mefloquine) was also found among the top compounds and was in fact the first compound from the list when the analysis was restricted to drugs already approved for medical use ( Fig. 2E and Table S6). We find this of potential interest given the debate that has been around the usefulness of hydroxycloroquine for the treatment of COVID-19, which was questioned by a recent registry analysis of data from 671 hospitals that concluded that the drug increases mortality in COVID-19 patients (38) (although this very study is also now being questioned (39)).

Analysis of drugs triggering transcriptional signatures opposed to the COVID CS
As mentioned in the introduction, one of the initial ideas behind the development of CMap was to identify drugs that trigger transcriptional responses that negatively correlate with those associated to a given disease, in order to identify potential new therapies through drug repurposing (24). Based on this concept, we next interrogated CMap for drugs that trigger signatures opposed to the COVID CS . Drug GSEA analyses identified glucocorticoid receptor agonists as the most significantly compounds that trigger a signature opposed to the COVID CS , followed by Vitamin K antagonists and MEK inhibitors (Fig. 3A). Furthermore, and in addition to glucocorticoid receptor agonists, the analysis of biological pathways (CMap classes) with signatures negatively correlating to the COVID CS also identifies drugs in clinical trials like JAK inhibitors (17,18) (Fig. 3B and Table S7). Other significantly enriched CMap classes include loss-of-function of S100 calcium-binding proteins and inhibitors of Bromodomains, Leucine-rich repeat kinases (LRRK) or the MEK kinase.
The connection with COVID-19 of LRRKs, Bromodomains or S100 calcium-binding proteins is unclear at this point. However, we believe that the potential effects of MEK inhibitors in counteracting the cytokine response seem of interest for a number of reasons. First, because even if MEK inhibitors have been mainly developed as antineoplastic agents, numerous studies support that targeting MEK limits cytokine responses in several contexts like graft-vs-host disease, cerebral ischemia, activation of T cells by superantigens or influenza A infection (40)(41)(42)(43)(44)(45). Second, because MEK inhibition has been previously shown to inhibit replication of other coronaviruses such as the mouse hepatitis virus (MHV) (46). Third, because besides MEK, several compounds with the highest scores of opposing signatures to the COVID CS are inhibitors of other factors from the same signaling route (MEK, RAF, EGFR, MAPK) ( Figure 3C). And, finally, because as far as we can tell inhibition of this signaling route is currently not being explored as a potential therapy in the context of COVID-19 and its associated CS.
Besides compound classes, among specific drugs with associated signatures that negatively correlate with COVID CS we found compounds known to limit inflammatory responses including the IKK inhibitor TPCA-1, which inhibits NFkB signaling and reduces acute lung injury (47) (Fig. 3D and Table S8). If the analysis was restricted to medically approved compounds the list was led by compounds like the antidepressant phensuximide or the anticoagulant warfarin ( Fig. 3E and Table S9), for which the relationship with inflammation is unclear at this point.

Genetic perturbations related to the COVID-19 CS transcriptional signature
While the previous analyses were performed on a transcriptional signature defined from COVID-19 patients, this was not explicitly mentioned to be obtained from the severe patients undergoing a CS. In order to define a transcriptional signature specifically related to the CS, we used data from another manuscript, where the BALF transcriptome of a COVID-19 patient suffering from a CS was compared to that of healthy individuals or a milder case of the disease (signature available in Table S1; CS s ) (48). Consistent with our previous findings on the COVID CS , NFkB and PKC signaling were the CMap classes most significantly similar to the CS s ( Fig. 4A and Table S2). Furthermore, IL1 (IL1R1) or TNFa receptors (TNFRS10A and TNFRSF1A) and IFN-response genes (IRF2, IFR5, IFI16) were once again found among the top list of factors that, when overexpressed, trigger a transcriptional signature similar to the CS s (Fig. 4B and Table S3). GO analyses from the top hits among overexpressed factors correlating to the CS s further identified inflammation-and cytokine productionrelated pathways amongst the most significantly enriched ones (Fig. 4C).
Interestingly, pathways related to Vitamin D metabolism and angiotensin-activated signaling were also significantly enriched among those that mimic the CS s , both of which have been linked to the severity of COVID-19 (49,50). As for genes that, when downregulated, would trigger a transcriptional signature similar to the CS s , this list also showed a significant overlap with that resulting from the COVID CS analyses (Fig. 4D and Table S4, compare to Fig. 1C). Importantly, GO analyses found once again an enrichment of pathways related to female hormone signaling among those linked to genes with signatures that negatively correlate with the CS s (Fig. 4E), further reinforcing the idea that the anti-inflammatory properties of estrogens might be behind the gender differences found on the severity of COVID-19. We also note an enrichment in pathways related to cholesterol metabolism, which together with the fact that "negative regulation of fat cell proliferation" was found in the corresponding COVID CS analysis, could be related to the role of obesity as a comorbidity factor for COVID-19 (37).

Identification of drugs triggering transcriptional signatures mimicking the CS s
We next looked into compounds eliciting a transcriptional signature that positively correlates to the CS s . Despite the disparity of the datasets used, the results from these analyses were also remarkably concordant with our previous findings on the COVID CS .
Drug GSEA found PKC activating drugs followed by "ATPase inhibitors" and "Apoptosis Stimulants" as compound classes that significantly correlate with the CS s (Fig. 5A,B); whereas only Topoisomerase inhibitors showed a significant enrichment if medically approved drugs were exclusively analyzed (Fig. 5C). Regarding the potential aggravating effects of topoisomerase inhibitors on the COVID-19 CS, we believe that this is most likely related to Topoisomerase-II. In fact, Topoisomerase-I activity has been found to mediate the expression of inflammatory genes (51). In contrast, Topoisomerase-II inhibition triggers the expression of inflammatory cytokines through activation of cGAS-STING innate immune sensing of cytoplasmic DNA (52).
Regarding specific drugs that could potentially aggravate the CS, and besides PKC activators (ingenol, PMA and prostratin), the histone methyltransferase inhibitor chaetocin has also been shown to boost cytokine production and inflammation (53, 54) (Fig. 5D and Table S5). Interestingly, if we focused on drugs already approved for medical use, the compound which signature most resembled the CS s was the ACE inhibitor fosinopril, which is of relevance given the debate on the potential aggravating effects of antihypertensive ACE inhibitors in the context of COVID-19 ( Fig. 5E and Table S6) (55). Unexpectedly, several drugs used for the treatment of viral or bacterial infections such as nitazoxanide, amantadine or primaquine were also found among approved compounds that trigger a signature similar to the CS s . Finally, another drug whose signature mimics the CS s is disulfiram (Antabuse), which is used for the treatment of alcohol abuse disorders (currently used by more than 200,000 patients in the USA (https://www.addictioncenter.com/alcohol/disulfiram/)). While how disulfiram could be related to the CS is at this point unclear, we want to note that, among other targets, this compound is a potent in vitro inhibitor of Topoisomerases (56).

Identification of drugs triggering transcriptional signatures opposing the CS s
We next queried CMap for drugs causing opposing signatures to the CS s . Once again, the general results from these analyses were highly consistent with those from the COVID CS . For instance, drug GSEA analyses yielded glucocorticoid receptor agonists as significantly enriched (Fig. 6A,B). Besides glucocorticoid receptor agonists, the analysis of CMap classes negatively correlating to the CS s again identified inhibitors of Bromodomains, Leucine-rich repeat kinases (LRRK), JAK kinases, V Type ATPases or the MEK kinase ( Fig. 6C and . 6D and Table S7).
The analysis yielded more subtle and interesting differences with the previous COVID CS analysis when we looked into specific compounds with signatures opposite to the CS s (Fig. 6E and Table S8). Among the top of this list there were several compounds reported to limit inflammatory responses such as treprostinil, a prostacyclin analog that blocks NFkB nuclear translocation in human alveolar macrophages (59); the antidepressant bupropion, which was shown to lower TNFa and IFNg production in mice (60); the IKK inhibitor TPCA-1 which, as mentioned, was shown to inhibit NFkB signaling and to reduce acute lung injury (47); and the farnesyltransferase inhibitor tipifarnib, which presents anti-inflammatory properties (61). When the analysis was restricted to drugs approved for medical use (Fig. 6F and Table S9), some of the compounds like niacin (vitamin B3) have already been proposed for the treatment of the CS in COVID-19 patients (62), and others like the calcium channel blocker verapamil reduce lung injury in mice after LPS injections, which triggers a CS-related pathology (63). Noteworthy, a recent study of drugs that could be repositioned for the treatment of COVID-19 based on interactome analyses also identified several compounds that we identify here including verapamil or inhibitors of Bromodomains and RAF kinases (64).
To end, we tried to investigate the potential mechanism by which the antidepressant Remarkably, the transcriptional signatures associated to phensuximide and warfarin were also inversely correlated to inflammation, particularly in macrophages, helping to understand the potential effect of these drugs in counteracting a CS (Table S10).
Altogether these analyses have allowed for a better understanding of the transcriptional perturbations associated to COVID-19 and provided a hypothesis for compound classes or specific drugs that could be investigated in the context of the CSs triggered by viral infections.

A CMap analysis of the COVID-19 CS based on single-cell RNAseq data
Finally, and in order to complement our analyses on whole transcriptomic data, we included TNF receptors and IFN-response genes ( Fig. 7B and Table S3). In what regards to compounds, drug GSEA analyses found again PKC activators as a significantly enriched compound class that triggers a transcriptional signature mimicking the CS s [sc], together with proteasome, HDAC or HSP inhibitors (Fig. 7C,D).
In striking contrast to whole transcriptomic analyses, in this case glucocorticoid receptor agonists were found as potentially stimulating the CS s [sc]. However, we should note that glucocorticoids are known to have cell-type specific effects and specifically in macrophages, where this signature has been defined, they have been shown to promote cytokine release (66). Nevertheless, and in agreement with our previous analyses of the COVID CS and CS S , the top of the list for similarity scores of individual compounds included the same three PKC activators (prostratin, PMA and ingenol) and the histone methyltransferase inhibitor chaetocin, as well as several others ( Fig. 7E and Table S5). By restricting the analysis to medically approved compounds, the list is also remarkably similar to our findings made on whole BALF transcriptomics, and again included disulfiram, menadione, amsacrine and nilotinib, as well as multiple topoisomerase inhibitors (amsacrine, teniposide, irinotecan, SN-38) (Fig. 7F and Table S6).
Finally, we searched for compounds signatures that negatively correlate with the transcriptase, bacterial 30S ribosome, PKA and glycogen synthase kinase ( Fig. 8A and Table S7). Among specific compounds that elicit transcriptional signatures opposed to the CS s [sc] the top hits included compounds not identified in whole transcriptome analyses and that could potentially be related to counteracting inflammation ( Fig. 8B and Table S8). These included the estrogen estriol, relevant in the context of the above mentioned gender-related differences in COVID-19 mortality; the dopamine receptor agonist fenoldopam, previously shown to reduce the severity of LPS-induced lung-injury and cytokine release in mice (67); the CDK inhibitor indirubin, known to reduce the LPS-induced cytokine expression in macrophages (68,69) and to decreased influenza virus derived lung damage (70,71); the Cytochrome P450 inhibitor metyrapone, found to reduce cytokine production in the context of inflammation (72); and the Calcium channel blocker bepridril and Myosin light chain kinase inhibitor ML-7 which have also been reported to present anti-inflammatory roles (73,74). From this list of compounds only fenoldopan, bepridil, metyrapone and estriol are already approved for medical use (Fig. 8C and Table S9).
The top hits among clinically available drugs with a signature negatively correlating to the CS s [sc] also included the RAF inhibitor vemurafenib and the antibiotic doxycycline.
For the reasons mentioned above, we are particularly intrigued by the potential of MAPK/MEK/RAF inhibitors in the context of the COVID-19 CS, which became evident in all of our 3 analyses (40)(41)(42)(43)(44)(45)(46). As for doxycycline, and besides its usefulness as an antibiotic that might combat bacterial infections triggering pneumonia in COVID-19 patients, doxycycline is also an anti-inflammatory agent (75) and several studies have already suggested its use to combat the COVID-19 CS (76,77). Further supporting these conclusions, drug GSEA analysis revealed MEK inhibitors and "Bacterial cell wall synthesis inhibitors", together with PI3K inhibitors (Fig. 8D,E) (Table   S11). We want to end this manuscript with a clear statement in that by no means we are proposing novel clinical indications for any of these agents. We simply wanted to contribute in the context of the current health crisis and provide a resource that others might find of use in the selection of drugs that could be tested for efficacy in experimental systems of CS.

Calculation of similarity scores
Three independent transcriptional signatures defined from COVID-19 patients (COVID CS , CS s and CS s [sc]) were used as inputs to query CMap using its clue.io tool (https://clue.io), in order to obtain similarity scores for the signatures associated to all the perturbagens (compounds, CMap classes and the over-expression or knock-down of genes) available at CMap (23). Similarity scores were downloaded, and results sorted based on their scores and the type of perturbagen. For restricting the analysis to medically approved drugs, we used drug information available at CMap. Except for simple filtering we used R version 3.6.3 (2020-02-29) for our computational analyses.

Drug "Gene Set Enrichment Analysis" (GSEA)
We adapted Gene Set Enrichment Analysis (GSEA) to enable enrichment analyses of "drug classes" based in their mechanism of action (MOA), for which we used the GSEA method implemented in the R package fgsea (36).

DREIMT analyses
The DREIMT database query tool (http://dreimt.sing-group.org/) was used to search for similarities between the transcriptional signatures triggered by specific compounds and those associated to immunological processes. Similarity scores were downloaded, and results sorted and filtered based on their scores (scores >90 or <-90).