key: cord-0300289-39730g3k
authors: Karunakaran, Kalyani B.; Ganapathiraju, Madhavi K.; Jain, Sanjeev; Brahmachari, Samir K.; Balakrishnan, N.
title: Drug Contraindications in Comorbid Diseases: a Protein Interactome Perspective
date: 2022-01-15
journal: bioRxiv
DOI: 10.1101/2022.01.11.475465
sha: 2fec2631686cd140ca6752084cc1af8e9aba8853
doc_id: 300289
cord_uid: 39730g3k

Adverse drug reactions (ADRs) are leading causes of death and drug withdrawals and frequently cooccur with comorbidities. However, systematic studies on the effects of drugs in comorbidities are lacking. Drug interactions with the cellular protein-protein interaction (PPI) network give rise to ADRs. We selected 6 comorbid disease pairs, identified the drugs used in the treatment of the individual diseases ‘A’ and ‘B’– 44 drugs in anxiety and depression, 128 in asthma and hypertension, 48 in chronic obstructive pulmonary disease and heart failure, 58 in type 2 diabetes and obesity, 58 in Parkinson’s disease and schizophrenia, and 84 in rheumatoid arthritis and osteoporosis – and categorized them based on whether they aggravate the comorbid condition. We constructed drug target networks (DTNs) and examined their enrichment among genes in disease A/B PPI networks, expressed across 53 tissues and involved in ~1000 pathways. To pinpoint the biological features characterizing the DTNs, we performed principal component analysis and computed the Euclidean distance between DTN component scores and feature loading values. DTNs of disease A drugs not contraindicated in B were affiliated with proteins common to A/B networks or uniquely found in the B network, similarly regulated common pathways, and disease-B specific pathways and tissues. DTNs of disease A drugs contraindicated in B were affiliated with common proteins or those uniquely found in the A network, differentially regulated common pathways, and disease A-specific pathways and tissues. Hence, DTN enrichment in pathways, tissues, and PPI networks of comorbid diseases will help identify drugs contraindications in comorbidities.

Comorbidity is the phenomenon in which one or more diseases co-exist with a primary disease in patients. Comorbidities are the norm rather than exceptions among chronic conditions and pose a significant threat to the physical and psychosocial wellbeing of patients [1] . Comorbidities increase with age, and the risk of mortality increases with the number of comorbidities. A longitudinal study (1992) (1993) (1994) (1995) (1996) (1997) (1998) (1999) (2000) (2001) (2002) (2003) (2004) (2005) (2006) has shown that the mortality risk increased by 25% in patients with 3-4 chronic comorbidities and by 80% in those with 5 or more comorbidities, both in comparison with individuals having no chronic conditions [2] . The prevalence of comorbidities increases from 10% in 0- 19 yearolds to 78% in individuals aged 80 or more [3] . The prevalence of comorbidity in women of age groups of 18-44 years, 45-64 years, and ≥ 65 years was 68%, 95%, and 99% and in men, it was 72%, 89%, and 97% [4] . As per the US National Comorbidity Survey Replication (NCS-R) survey, 73.8-98 .2% of the respondents reported having at least one comorbid condition along with a primary condition [1] . The most striking finding from this report was that the estimates of individual disease burden based on the respondents' perception of their health condition decreased substantially when adjusted for comorbidity [1] . This effect was particularly magnified for neurological disorders, chronic pain, anxiety disorders, major depressive disorder, and diabetes, all of which contribute immensely to the global disease burden [1] . For example, anxiety disorders collectively affect 284 million people (63% females, 2.5-7% variation by country) around the world, and are among the most prevalent mental health and neurodevelopmental disorders (WHO and IHME, 2017) [5] .

Disease comorbidity may increase the likelihood of experiencing adverse drug reactions [6] [7] [8] . Drugs that are beneficial in the treatment of one disease may aggravate or even cause comorbid conditions, giving rise to adverse drug reactions, e.g. beta-blockers that treat hypertension and heart disease may aggravate asthma [6] , trimethoprim and sulfamethoxazole to treat AIDS may increase the patient's susceptibility to Stevens-Johnson syndrome and toxic epidermal necrolysis [7] ; malaria patients with AIDS and osteoarthritis treated with artemisinin based combination antimalarial therapy were 3 times more likely to experience adverse side effects [8] . Serious adverse drug reactions constitute the 

The PPI networks of the proteins encoded by the disease-associated genes were assembled by extracting their protein interactors from the PPI repositories BioGRID [27] and HPRD [26] using BisoGenet [28] and the network building options specified before. The input nodes for the construction of each of the disease networks were the 100 top-ranking genes compiled from the DisGeNET database.

Matching node ratio (N M ) was measured as the ratio of the total number of common nodes shared between the two PPI networks of a comorbid pair and the total number of unique nodes in the two disease networks [30] . Matching link ratio (L M ) was measured as the ratio of the total number of common links (i.e. edges) shared between the two PPI networks of a comorbid pair and the total number of unique links in the two disease networks [30] . ‫ܮ‬ ெ ൌ ‫ת‬ ‫‬ (2) A l = Number of links in the PPI network of disease A B l = Number of links in the PPI network of disease B

The same formula shown above was also used to calculate the matching link ratio for common links of path length 2 and path length 3. Links of specific path lengths were retrieved using the Cytoscape application called NetworkAnalyzer [31, 32] . individuals who had applied for support from the U.S. Medicare program during 1990-1993 [33] .

Comorbidity data was available for five out of our six comorbid disease pairs and two out of the three non-comorbid pairs in HuDiNe. Specifically, data was not available for Anxiety -Depression and Multiple sclerosis -Peroxisomal disorders. Hence, N A , N B and N AB were extracted for seven out of the nine disease pairs. The diseases were specified in the form of their ICD-9 codes (at three digits level): asthma (ICD-9: 493), hypertension (ICD-9: 401), type 2 diabetes (ICD-9: 250), obesity (ICD-9: 278), chronic obstructive pulmonary disease (ICD-9: 496), heart failure (ICD-9: 428), Parkinson's disease (ICD-9: 332), schizophrenia (ICD-9: 295), rheumatoid arthritis (ICD-9: 714) and osteoporosis (ICD-9: 733). The population size N was considered to be 13,039,018, i.e. the total number of individuals represented in the HuDiNe dataset.

WebGestalt [34] was used to compute the distribution of genes involved in specific signalling pathways in the drug target networks, and compare it with the background distribution of genes belonging to this pathway among all the genes associated with any pathway in the selected database (Reactome) [35] . Statistical significance of the enrichment was computed using Fisher's exact test and corrected using the Benjamini-Hochberg method for multiple test adjustment.

The enrichment of the drug target networks in genes expressed in specific tissues was computed using RNA-sequencing data from 53 postnatal human tissues extracted from GTEx [36] (version 8). Genes with high or medium expression (transcripts per million (TPM) ≥ 9) in 53 tissues were included, provided that they were not housekeeping genes, i.e. genes detected in all the tissues with transcripts per million ≥ 1, as identified in the Human Protein Atlas [37] . TPM is a metric for quantifying gene expression; it directly measures the relative abundance of transcripts. The GMT files served as inputs for a gene over-representation analysis (GSEA) based on hypergeometric distribution. The following GWAS datasets were selected in TSEA-DB [38] for identification of disease-specific tissues (trait IDs are given in parentheses): anxiety (4679), depression (5315), chronic obstructive pulmonary disease (571), heart failure (5333), asthma (5259), hypertension (169), type 2 diabetes (4628), obesity (1031), Parkinson's disease (4607), schizophrenia (5215), rheumatoid arthritis (4614) and osteoporosis (746).

BaseSpace Correlation Engine (https://covid-19.ce.basespace.illumina.com/c/nextbio.nb) was used to identify the correlations between the gene expression profile induced by maprotiline in PC3 cells (Broad Connectivity Map (CMAP 2.0) [39] ), the expression profile associated with major depressive disorder and generalized anxiety disorder (GSE98793 [40] ) and the expression profile of adrenal cortex. The software uses a non-parametric rank-based approach to compute the extent of enrichment of a particular set of genes (or 'bioset') in another set of genes [41] .

Principal component analysis (PCA) was used to capture relationships between the drug target networks and the disease networks/biological pathways/tissues. For each disease pair, negative log-transformed p-values indicating the statistical enrichment of the disease networks/biological pathways/tissues in the 4 drug target networks were assembled into a data matrix containing disease networks/biological pathways/tissues as rows and drug target networks as columns; each cell in the matrix contained a -log 10 P value. Following the established approach [42] , log transformation was performed to reduce the influence of extreme values on the obtained PCs. PCA was performed with a web-based tool called ClustVis (https://biit.cs.ut.ee/clustvis/) [43] . The data matrix was pre-processed such that 70% missing values were allowed across the rows and columns. The -log 10 P values in the matrix were further centred using the unit variance scaling method, in which the values are divided by standard deviation so that each row or column has a variance of one; this ensures that they assume equal importance while finding the components. The method called singular value decomposition (SVD) with imputation was used to extract principal components. In this method, missing values are predicted and iteratively filled using neighbouring values during SVD computation, until the estimates of missing values converge. The factor/component loadings corresponding to the disease networks/pathways/tissues that contributed to the selected principal components were also extracted.

Component loadings are correlation coefficients between the variables in rows and the factors (i.e. PC1, PC2 etc.). The squared value of a component loading gives the percentage of the variance explained by a particular original variable, and essentially its contribution to the principal components. Finally, for each of the disease pairs, the Euclidean distance between the principal component scores of each of the drug target networks were computed for all the component loading values pertaining to the particular biological modality. This resulted in a list of the specific disease protein sets/pathways/tissues that may be closely related to each of the different drug target networks.

To identify potential mechanisms of adverse drug interactions within comorbid diseases, we systematically studied pairs of comorbid diseases ('disease A' and 'disease B') and their FDAapproved drugs. We separated the drugs into two groups, namely, disease A drugs that are (a) contraindicated and (b) not contraindicated in disease B, and disease B drugs that are (c) contraindicated and (d) not contraindicated in disease A We then constructed the interactomes of the proteins targeted by these drugs and examined these drug target interactomes in the context of three biological factors, namely, (i) proteins exclusive to interactomes of diseases A and B and those that are in their intersection, and (ii) biological pathways and (iii) tissues associated with these drug target interactomes.

Specifically, we selected three pairs of non-comorbid diseases as negative controls and six pairs of comorbid diseases for our analysis. The non-comorbid pairs were: (I) Multiple sclerosis -Peroxisomal disorders [44] , (II) Schizophrenia -Rheumatoid arthritis [45] [46] [47] , (III) Asthma -Schizophrenia [48] . The comorbid pairs were (IV) Anxiety -Depression [49] , (V) Asthma -Hypertension [50, 51] , (VI) Chronic obstructive pulmonary disorder (COPD) -Heart failure [52, 53] , (VII) Type 2 diabetes -Obesity [54, 55] , (VIII) Rheumatoid arthritis -Osteoporosis [56] and (IX) Parkinson's disease -Schizophrenia [57] .

The drugs indicated for use in each of the diseases were retrieved from Drug Bank (version 5.1.8) [23] . For each pair, we categorized the drugs into the four groups (a-d) mentioned earlier, based on their clinical activity in the diseases, collected from the TWOSIDES database (version 0.1) [24] , a compendium of drugs and their contraindications (see Additional File 2: Table S2 ). Drugs contributing to specific adverse effects were collected by manually selecting relevant 'condition concept names' (Additional File 1: Table S1 ). For example, to identify the anxiolytic drugs that may cause depression, the condition concept names, depression, major depression, depressive symptom, depression suicidal, depression postoperative, postpartum depression, depressive delusion, and agitated depression, were selected. The list of anxiolytic drugs was then compared with the list of drugs associated with these condition concept names. The matching drugs were compiled into groups 'a' and 'c', for example, "drugs effective in anxiety and contraindicated in depression". Similarly, groups 'b' and 'd' drugs were compiled. The proteins targeted by the drugs belonging to groups a and b were retrieved by querying the Drug Bank database through the DGIdb drug-genee interaction database) web portal [25] (see Additional File 3: Table S3 ). Finally, the protein-protein interaction (PPI) networks of the drug targets were assembled by extracting their protein interactors from the PPI repositories BioGRID [27] (version 4.3.194) and HPRD [26] (version 9) using a Cytoscape plugin called BisoGenet [28] .

The methodology of our study is illustrated in Fig. 1 . To characterize the 4 classes of drug target networks (DTNs), we examined 3 types of data that may reflect their biological profiles, namely (i) disease PPI networks, (ii) biological pathways and (iii,) tissue gene expression. Specifically, we conducted gene overrepresentation analyses based on hypergeometric distribution to check the enrichment of the DTNs among proteins that are unique to/shared between networks of disease A and disease B, genes showing high/moderate expression in 53 tissues across the human body, and proteins involved in ~1000 biological pathways. Overlaps computed in this manner with each of the 3 types of biological data were considered to be statistically significant at p-value < 0.05 after multiple test adjustments with the Benjamini-Hochberg method.

As a first step towards identifying the specific biological data modalities (disease subnetworks/pathways/tissues) that were relatively more 'closer' to each of the different types of DTNs in terms of Euclidean distance, we generated a data matrix of the DTNs (columns) versus the various members of the biological data modality (rows) (for example, for the data modality 'disease subnetwork', the members would be 'common to both the networks', 'unique to disease A network' and 'unique to disease B network' and for the data modality 'tissue', the members would be 'amygdala', 'aorta', 'lungs' etc.). Each cell contained the negative of log-transformed p-values. -log 10 transformed p-values have been used as inputs for PCA in previous studies [58, 59] . Following the established approach [42] , log transformation was performed to reduce the influence of extreme values on the obtained PCs. Single value decomposition (SVD) with imputation and unit variance scaling was applied to this matrix to extract principal components that explained the variance observed with each of the data modalities across the DTNs. Principal component analysis (PCA) has been applied to matrices containing gene-level association scores in several studies [59] . PCA is primarily used to capture systematic variations underlying datasets. All the principal components generated after this analysis were considered for our study, since they may together reveal underlying clustering patterns among the different DTNs. 

Relative risk is an experiential measure of comorbidity as it compares the observed prevalence of a pair of comorbid diseases in the population with the expected number, which is calculated based on the prevalence of the individual diseases in the population. We then explored whether this information was embedded in the disease networks, i.e., whether the relative risk of comorbidity of the disease pairs would be reflected in the similarity of the disease networks. For each of the comorbid pairs, we computed four established network similarity measures, namely, matching node ratio (N M ) for all the nodes shared between the two disease networks, and the matching link ratio (L M ) [30] for all the (i) shared links (i.e. edges), (ii) shared links of path length 2 (connecting two nodes via one intermediate node) and (iii) shared links of path length 3 (connecting two nodes via two intermediate nodes)

between the two disease networks.

We computed the relative risk for each of the disease pairs observed in hospital claims data of 13,039,018 U.S. individuals who had filed for support from the Medicare program during the period of 1990-1993, made available as the HuDiNe dataset [33] . The ICD-9 codes corresponding to pairs of diseases diagnosed as primary and secondary conditions, along with the number of individuals who were diagnosed with diseases A or B or both (N A , N B and N AB , respectively) were available (see Methods). Comorbidity data was available for five out of our six comorbid disease pairs (i.e. except for Anxiety -Depression) and two out of the three non-comorbid pairs (i.e. except for Multiple sclerosis -Peroxisomal disorders in HuDiNe.

For each of the diseases considered, the top 100 genes associated with the disease were curated from the DisGeNET database (version 7) [29] based on their gene-disease association (GDA) scores (see Additional File 4: Table S4 ). The GDA score ranges from 0 to 1 and is computed for a gene based on the number of publications supporting its association with the disease, and the number and types of database sources (levels of curation (expert-curated/computationally-predicted) and the model four measures of network similarity, namely, matching node ratio (green data points), matching link ratio of all shared edges (red data points), matching link ratio of all shared organisms in which the association was validated). The 100 top-ranking genes collected in this manner were used as starting points for the construction of disease networks. Here also, the network is assembled by extracting PPIs from BioGRID and HPRD using a Cytoscape plugin BisoGenet, similar to assembling of DTNs. Then, we systematically conducted network overlap analyses with each of the 9 disease pairs and identified the proteins (a) shared between the two disease networks, (b) unique to disease A and (c) unique to disease B ( Table 1) . The relative risk between diseases was proportional to the matching node and link ratios (Fig. 2) . The control disease pairs showed low relative risks and smaller disease network overlaps, whereas three out of five comorbid disease pairs showed high relative risks and larger network overlaps, namely, Asthma -Hypertension, COPD -Heart failure and Type 2 diabetes -Obesity. However, this trend was not seen in the comorbid pairs, Rheumatoid arthritis -Osteoporosis and Parkinson's disease -Schizophrenia. Specifically, their higher relative risks (compared with other comorbid pairs), were not accompanied by a corresponding increase in the network overlaps. ~85% of the human interactome awaits experimental discovery [60] . Hence, two factors may have led to the underestimation of the network overlaps. Firstly, the inherent incompleteness of these disease networks [60] . Secondly, the tendency of incomplete networks to exhibit small overlaps [60] .

Next, we tested the potential of each of the disease subnetworks to be acted upon by drugs or their susceptibility to pharmacological modulation (druggability), by examining their enrichment among a group of 4,463 proteins deemed to be druggable [61] , similar to the approach followed in a previous study [62] . These proteins are bound with high affinity at specific binding sites by drugs that follow the Lipinski's 'rule-of-five', i.e. orally bioavailable drugs with specific molecular characteristics that influence their pharmacokinetic ability to enter systemic circulation and act on their target sites ( We found that the proteins shared between the two diseases were the most significantly enriched for druggable targets in 5 out of the 6 tested comorbid pairs ( Table 2 ). In case of the sixth pair, namely anxiety and depression, the proteins that are exclusive to the depression network were found to be more enriched for druggable targets. In 2 out of the 5 disease pairs that shared many common drug targets, the drug target proteins were significantly enriched in G protein-coupled receptor activity (pvalue <0.05) ( Table 3) . Based on these observations and the finding in the previous section that relative risk varies in tandem with network similarity measures, we speculated that contraindications in comorbidities may arise from drug action on druggable proteins shared between the networks of comorbid diseases ( 

To test these corollaries, we systematically computed the overlaps between three groups of disease proteins, namely, proteins that are (a) common to disease A and disease B networks, (b) unique to Table 4) ; previous studies have examined the overlaps between the PPI networks of drug targets and disease-associated proteins [20, 64] . For each of the six disease pairs, we created a data matrix of DTNs (columns) versus disease subnetworks (rows), which contains -log(p-values)

indicating the statistical significance of these enrichments. This data matrix was used as the input for PCA. In order to identify the specific disease subnetworks that were the nearest to each of the DTNs, we calculated the Euclidean distance between the PC scores of each of the DTNs across all the extracted axes and the corresponding component loading values of all the disease subnetworks across these axes (following the methodology depicted in Fig. 1) . By counting the two disease subnetworks that were the closest to each of the different DTNs, we identified two predominant patterns.

In 10 out of the 12 cases, the DTNs of drugs used for a specific disease and not contraindicated in a comorbid condition were found to be closest/second closest to the proteins uniquely found in the network of the comorbid condition. Additionally, in 9 out of the 12 cases, they were closest/second closest to the proteins shared between the networks of both the diseases. In contrast, the DTNs of drugs used for a specific disease and contraindicated in a comorbid condition were found to be closest/second closest to the proteins uniquely found in the network of the disease for which these drugs were primarily used in 8 out of the 12 cases. Additionally, in 9 out of the 12 cases, they were closest/second closest to the proteins shared between the networks of both the diseases.

These observations led us to speculate two scenarios. Firstly, disease A drugs that are not contraindicated in disease B may target proteins unique to the disease B subnetwork involved in mechanisms that are either inconsequential/beneficial for disease B, but whose modulation is certainly beneficial for the treatment of disease A. Alternatively, they may target common mechanisms that are dysregulated in a similar manner in both the diseases and pharmacologically modulate them in a similar direction. Secondly, disease A drugs may become contraindicated in disease B when they target either (a) common mechanisms that are pharmacologically oppositely modulated in a manner that benefits disease A but aggravates disease B or (b) mechanisms unique to disease A that aggravate disease B. Additionally, we hypothesized that biological processes such as signalling pathways that function at a higher level than disease subnetworks could be regulating the action of drugs under comorbid conditions. *, ** and *** indicate low, medium and high levels of statistical significance. , and indicate non-significant overrepresentation, non-significant underrepresentation and significant underrepresentation respectively.

Common 

We identified the pathway associations of the DTNs using the gene set analysis toolkit called

WebGestalt [34] . WebGestalt computes statistical significance enrichment of a functional group (e.g., G alpha (12/13) signalling events' and 'muscarinic acetylcholine receptors' were identified among the top-10 pathways that were close to anxiety drugs not contraindicated in depression in our study (Fig.   3) . Adrenergic receptor signalling could be regulated via G α (12/13); Gα12 and Gα13 have been shown to mediate alpha-1 adrenergic receptor-induced JNK activation in rat cardiomyocytes [65] . The drug maprotiline was among our list of anxiety drugs without contraindications in depression (Fig.   4) . Corroborating this, clinical data suggested that the drug is effective in alleviating anxiety symptoms co-occurring with depression [66] . Maprotiline acts an inhibitor of SLC6A2 (sodiumdependent noradrenaline transporter) and inhibits noradrenaline reuptake in the brain. It also acts as an and/or depression (disease B) associated genes have been shown. Note that maprotiline (an anti-anxiety (disease A) drug not contraindicated in depression (disease B)) targets a higher number of proteins associated uniquely with depression in the adrenergic, serotonergic and cholinergic systems, which is in line with our observation that disease A drugs that are not contraindicated in disease B are closely associated with proteins uniquely found in the disease B network (i.e. depression in this specific example). Serotonin receptors were found to be associated in our analysis with depression drugs not contraindicated in anxiety; antagonistic activity on serotonin receptors is shown by two such drugs shown in the diagram (flupentixol and mirtazapine).

antagonist to alpha-1 adrenergic receptors (ADRA1A, ADRA1B and ADRA1D) and alpha-2 adrenergic autoreceptors and heteroreceptors (ADRA2A, ADRA2B and ADRA2C), and enhances central noradrenergic and serotonergic functions, which have been linked to alleviation of anxiety and depression [67] . Maprotiline also acts as a weak antagonist to muscarinic acetylcholine receptors (CHRM1, CHRM2, CHRM3, CHRM4 and CHRM5); enhanced cholinergic signaling has been linked to both anxiety and depression [68] . It is notable that maprotiline targets a higher number of proteins associated uniquely with depression (ADRA2A, HTR2C, SLC6A2 and CHRM2) in the adrenergic, serotonergic and cholinergic systems (Fig. 4) . It targets only one receptor associated with both anxiety and depression (HTR2A), and no gene uniquely associated with anxiety (Fig. 4) . These observations are in line with our findings with the overlap of DTNs with disease subnetworks, i.e. DTNs of disease A drugs that are not contraindicated in disease B (e.g. maprotiline) are closely associated with proteins uniquely found in the disease B network (i.e. depression in this specific example). Since maprotiline has been discontinued from usage since 2020 in U.S. [69] , note that we are citing this drug only as a demonstrative example. Drugs acting on the serotonergic system is known to be effective in both short-term and long-term treatment of patients with major depressive disorder and anxiety disorders [70] . 'Serotonin receptors' was identified among the top-10 pathways that were close to depression drugs not contraindicated in anxiety (Fig. 3) . This may suggest the broadspectrum efficacy of drugs acting on serotonin receptors in both the conditions. Two such drugs in our study displayed antagonistic activity on serotonin receptors -flupentixol [71] and mirtazapine [72] (both acting on HTR2A and HTR2C) -and have been used to treat depression accompanied by anxiety symptoms (Fig. 4) .

The pathway 'dopamine receptors' was found to be close to PD drugs contraindicated in SCZ (Fig. 5) ,

indicating that the enhancement in dopamine levels brought about by PD drugs may in fact induce Parkinson's disease and/or schizophrenia associated genes have been shown. Note that levodopa and ropinirole are used in the treatment of Parkinson's disease (disease A), but contraindicated in schizophrenia (disease B), and flupentixol is used in the treatment of schizophrenia, but contraindicated in Parkinson's disease. Note that levodopa and ropinirole target a higher number of dopamine receptors associated uniquely with Parkinson's disease, which supports our finding that disease A drugs that are contraindicated in disease B are closely associated with proteins uniquely found in the disease A network (i.e. Parkinson's disease in this specific example). SCZ, which has been linked to a hyperdopaminergic state [57] . The dopamine agonists belonging to this group of PD drugs have been shown to induce psychosis, namely, levodopa (acting on DRD1, DRD2, DRD3, DRD4 and DRD5) and ropinirole (DRD2, DRD3 and DRD4) (Fig. 6) [73, 74] . It is notable that levodopa and ropinirole target a higher number of dopamine receptors associated uniquely with Parkinson's disease (DRD1 and DRD2) (Fig. 6) . It targets only one dopamine receptor (DRD3) uniquely associated with schizophrenia (Fig. 6) . These observations are in line with our findings with the overlap of DTNs with subnetworks, i.e. DTNs of disease A drugs that are contraindicated in disease B (e.g. levodopa and ropinirole) are closely associated with proteins uniquely found in the disease A network (i.e. Parkinson's disease in this specific example). Figures S7-S13) . Following this, we employed the tissue-specific enrichment analysis database (TSEA-DB) [38] to retrieve the top-3 tissues that may be preferentially affiliated with the diseases in each of the pairs. TSEA-DB is a reference database for information on disease-associated tissues, specifically, the tissues in GTEx that show significant enrichment of genes harbouring diseasesassociated variants compiled from the GWAS catalog [38] . We checked whether the top-3 tissues associated with each of the diseases in a disease pair (according to TSEA-DB) appeared among the list of tissues identified to be nearest to each of the 4 DTNs pertaining to this disease pair in our analysis. Out of the 11 tissues identified to be closer to the target networks of drugs used for a specific disease and not contraindicated in a comorbid condition, 6 were found to be associated with the comorbid condition as per TSEA-DB, whereas 3 were associated with the specific disease for which the drugs were used and 2 were associated with the disease as well as the comorbid condition.

Conversely, out of the 9 tissues identified to be closer to the target networks of drugs used for a specific disease and contraindicated in a comorbid condition, 5 were found to be associated with the specific disease, whereas 3 were associated with the comorbid condition in which the drugs were contraindicated and one was associated with the disease as well as the comorbid condition.

These percentages obtained with a low number of tissues suggest cautious interpretation.

Nevertheless, these results seem to corroborate our previous findings with disease subnetworks and biological pathways. Specifically, the networks of disease A drugs that are not contraindicated for disease B seemed to be nearest to tissues preferentially affiliated with disease B. This could indicate that these tissues could be equally important to the pathophysiology of disease A and its therapeutic alleviation (as they might be to these same aspects of disease B), despite showing a high enrichment for genes harbouring disease B-associated variants. For example, the adrenal gland was detected as a tissue highly specific to depression by TSEA-DB. In our analysis, this tissue appeared to be nearest to DTNs of anxiety and depression along PC1 and PC2, which explain 90.7% and 6.5% of the total variance respectively. The tissues that were exclusively associated with each of the 4 DTNs among the top-ten tissues that were identified to be highly related to the DTNs, after computing the Euclidean distance between the component loading values and the component scores, are shown as square-shaped data points for the DTN of drugs effective in anxiety and not contraindicated in depression, diamond-shaped data points for the DTN of drugs effective in depression and not contraindicated in anxiety, triangle-shaped data points for the DTN of drugs effective in anxiety and contraindicated in depression and cross-mark-shaped data points for the DTN of drugs effective in depression and contraindicated in anxiety. The tissues shown in circular and rectangular boxes were also identified to be highly specific to anxiety and depression respectively by TSEA-DB (due to a significant enrichment of anxiety/depression-associated variants). Note that adrenal cortex, which was identified to be associated with anti-anxiety (disease A) drugs that are not contraindicated in depression (disease B), is a tissue enriched with depression (i.e. disease B) associated variants. This corroborates our finding that disease A drugs that are not contraindicated in disease B are affiliated with disease B-specific tissues.

the DTN of anxiety drugs that were not contraindicated in depression (Fig. 7) , indicating that targeting of the adrenal gland may be vital to treat anxiety without aggravating comorbid depressive symptoms. The adrenal gland is an organ in the endocrine system that secretes the cortisol hormone, following the activation of the hypothalamic-pituitary-adrenal (HPA) axis by psychological stressors [75] . Several studies support the role of the adrenal gland as a focal point for depression. The adrenal gland exhibits a 70% increase in its volume in depressed individuals before successful anti-depressant treatment as well as in comparison with their matched controls [76, 77] . The cortisol hormone secreted by the adrenal gland, upon stress-induced activation of the HPA axis, has been linked to depressive symptoms in humans and monkeys. Increased cortisol levels have been positively correlated with depressive behaviour in rhesus macaques [78] . Enhanced cortisol secretion has been observed in depressive individuals [79] , and has been proposed to (a) increase susceptibility to depression [80] and (b) be correlated with the stress experienced by depressed individuals [81] .

Hyperactivation of the HPA axis has been noted in generalized anxiety disorder [82] . Treatment with selective serotonin reuptake inhibitors (SSRIs) has been shown to reduce HPA hyperactivity in both depressed patients and patients with generalized anxiety disorder [83] [84] [85] . Therefore, it is possible that anti-anxiety drugs that do not aggravate depressive symptoms target the adrenal gland, which produces the cortisol hormone, an effector or 'endpoint' of the HPA axis that seems to be regulated in a similar manner in depression as well as anxiety. We performed comparative transcriptome analysis of disease-associated, tissue-associated and drug-induced gene expression profiles using the BaseSpace Correlation Engine to analyse this hypothesis. BaseSpace Correlation Engine software suite is a data analysis platform that is used to study the effect of diseases and drugs on publicly available gene expression data [86] .

As mentioned in the previous section, maprotiline was found among our list of anti-anxiety drugs that are not contraindicated in depression; clinical data supports its utility in the treatment of anxiety symptoms associated with depression [66] . The differential gene expression (DGE) profile induced by maprotiline (12.8 µM) in PC3 cells (Broad Connectivity Map (CMAP 2.0) [39] ) was negatively correlated with the profile identified in the blood samples of patients with major depressive disorder patients (MDD) with generalized anxiety disorder (GAD) versus MDD patients without GAD (GSE98793 [40] ) (Fig. 8a) . This negative correlation of maprotiline with MDD/GAD could illustrate the fact that drugs administered to treat diseases often revert the expression of perturbed diseaseassociated genes to their normal levels [87, 88] . Secondly, the MDD/GAD profile was negatively correlated with the expression profile of adrenal gland cortex (Fig. 8a) , indicating that this tissue could be critical to disease alleviation. Maprotiline-induced DGE profile was positively correlated with the profile of adrenal gland (Fig. 8a) , indicating that maprotiline-mediated MDD/GAD alleviation may be dependent on adrenal gland, i.e. the reversal of MDD/GAD-associated expression profile induced by maprotiline could occur in the adrenal cortex. We then asked whether the genes differentially expressed in each of these profiles converged on a common set of biological processes.

Specifically, we identified the top-10 Gene Ontology (GO) biological processes enriched among the genes differentially expressed in (i) MDD/GAD versus maprotiline (in different directions), (ii) MDD/GAD versus adrenal cortex (in different directions) and (iii) maprotiline versus adrenal cortex (in the same direction). We then used the web-based tool called NaviGO [89] to group these 30 enriched biological processes into functionally cohesive networks based on semantic similarity measures of GO terms. Two such functional networks not only had top-scoring edges between the GO terms, but also contained GO terms enriched among all the three differential expression profiles (Fig.   8b, c) . One network contained four GO terms associated with protein folding (Fig. 8b) , and another network contained eleven GO terms representing cell cycle events (Fig. 8c) . Interestingly, 'protein folding', 'cyclin D associated events in G1' and 'G1 phase' were independently retrieved among our top-10 Reactome pathways found to be nearest (in terms of Euclidean distance) to the DTN of antidepressants that are not contraindicated in anxiety. Together, these results suggest that adrenal cortex (a) Correlation of differential gene expression profiles associated with a comorbid condition (major depressive disorder and generalized anxiety disorder), a drug (maprotiline) and a tissue (adrenal cortex). -log 10 (p-values) indicating the overlap of the expression profiles have been shown; red and green colors indicate negative and positive correlations between the profiles respectively. Significant overlap was found among the genes that are upregulated in patients with both major depressive disorder (MDD) and generalized anxiety disorder (GAD) and downregulated on treating PC3 cells with maprotiline (p-value = 7.4E-06), among the genes that are upregulated in MDD/GAD patients and downregulated in adrenal cortex (p-value = 8.4E-28), and among the genes that are downregulated on treating PC3 cells with maprotiline and downregulated in adrenal cortex (p-value = 0.034). (b, c) The functional networks of the Gene Ontology biological processes related to (b) protein folding and (c) cell cycle events that were enriched in the three expression profiles. The GO terms associated with each of the expression profiles have been shown using different node colors. The thickness of the edges corresponds to the Resnik semantic similarity score for GO terms (greater the thickness of the edges, greater is the similarity between the linked GO terms). may be preferentially targeted by drugs such as maprotiline that produce beneficial effects in anxiety as well as in depression, and that their actions may converge on protein folding and cell cycle processes. Note that maprotiline has been discontinued from usage [69] and is only being cited here as a demonstrative example.

On the other hand, the networks of disease A that are contraindicated for disease B seemed to be nearest to tissues preferentially affiliated with disease A. This could indicate that these disease Aspecific tissues may play a role in producing beneficial effects in disease A, while producing deleterious effects in disease B. For example, spleen was detected as a tissue highly specific to rheumatoid arthritis. The primary functions of this lymphoid organ are blood filtration, recycling of iron from old blood cells and generation of adaptive immune responses against bacterial, fungal and viral infections [90] . However, spleen has also been shown to act as a reservoir of osteoclast precursor cells, which upon resorption into bones, differentiate into osteoclasts. [91] Splenomegaly (enlargement of the spleen) has been noted in 5-10% and 52% of rheumatoid arthritis patients in separate studies (based on physical examination and imaging studies respectively) [92] [93] [94] . Rheumatoid arthritis patients are also prone to developing spontaneous splenic ruptures [95] . In our analysis, spleen was identified to be nearest (in terms of Euclidean distance) to the DTN of rheumatoid arthritis drugs that were contraindicated in osteoporosis (Fig. 9) . This seemed to indicate that spleen mediated opposite effects in rheumatoid arthritis and osteoporosis. Anecdotal evidence seemed to support this conjecture. While splenectomy seemed to improve rheumatoid arthritis in a patient [96] , it seemed to inhibit (a) attenuation of osteoporosis in a rat model [97] and (b) fracture healing in patients [98] . Table 5 summarizes the general conclusions of our study. We discovered that the DTNs of disease A drugs that are not contraindicated for a disease B may be nearest (in terms of Euclidean distance) to Table 5) . On the other hand, disease A drugs that are contraindicated for a disease B may be nearest to (a) proteins that are either uniquely found in the PPI network of disease A or are shared between the PPI networks of disease A and disease B, (b) biological pathways that are associated with disease A or are commonly active in both the diseases, and are regulated in an opposing manner in both the diseases and (c) tissues showing a high enrichment of disease Aassociated variants and thereby preferential affiliation with the etiology of disease A, and mediating opposing effects in disease A and disease B ( Table 5) . Table 5 : Disease network, pathway and tissue-level characterization of drugs that are contraindicated/not contraindicated in comorbid conditions. A has been used to indicate the Component loading values of 39 tissues associated with the drug target networks (DTNs) of rheumatoid arthritis (RA) and osteoporosis to PC1 and PC2 have been plotted along the X and Y axes respectively. PCA was performed with the p-values of enrichment of the tissues significantly associated (p-value < 0.05) with the DTNs of RA and osteoporosis. These values were transformed to -log 10 P values, which were then assembled into a data matrix containing tissues as rows and DTNs as columns. Unit variance scaling was applied across this matrix. Single value decomposition (SVD) with imputation was used to extract the principal components (PCs). The component loading values shown in the figure correspond to component scores of 4 DTNs along PC1 and PC2 that explain 89% and 6% of the total variance respectively. The tissues that were exclusively associated with each of the 4 DTNs among the top-ten tissues that were identified to be highly related to the DTNs, after computing the Euclidean distance between the component loading values and the component scores, are shown as square-shaped data points for the DTN of drugs effective in RA and not contraindicated in osteoporosis, diamond-shaped data points for the DTN of drugs effective in osteoporosis and not contraindicated in RA, triangle-shaped data points for the DTN of drugs effective in RA and contraindicated in osteoporosis and cross-mark-shaped data points for the DTN of drugs effective in osteoporosis and contraindicated in RA. The tissues shown in circular and rectangular boxes were also identified to be highly specific to RA and osteoporosis respectively by TSEA-DB (due to a significant enrichment of RA/osteoporosis-associated variants). Note that spleen, which was identified to be associated with rheumatoid arthritis (disease A) drugs that are contraindicated in osteoporosis (disease B), is a tissue enriched with rheumatoid arthritis (i.e. disease A) associated variants. This corroborates our finding that disease A drugs that are contraindicated in disease B are affiliated with disease A-specific tissues.

close affiliation of a specific category of drug target network with specific disease protein sets, disease-associated pathways and tissues. 

Despite the increased prevalence of adverse drug reactions in comorbidities, knowledge on the mechanistic basis of drug contraindications in such conditions is limited. In our study, we attempted to characterize the biological profiles of the target networks of drugs used in specific diseases that are either contraindicated or not contraindicated in a comorbid disease. We sought to provide an integrated interactome, pathway and tissue level view of the drug target networks.

The first key finding in our study was that the relative risk of comorbidity between diseases was proportional to their network similarity measures (Fig. 2) . The four network similarity measures along with the relative risk were low in the case of our three negative control pairs, namely, Multiple sclerosis -Peroxisomal disorders, Schizophrenia -Rheumatoid arthritis, and Asthma -Schizophrenia. This confirmed that these were indeed non-comorbid pairs. The network similarity measures and relative risk were higher in the case of Anxiety -Depression, Asthma -Hypertension, Chronic obstructive pulmonary disorder -Heart failure, Type 2 diabetes -Obesity, Rheumatoid arthritis -Osteoporosis, and Parkinson's disease -Schizophrenia, confirming that they were comorbidities. However, these measures do not follow the same trend in the case of the comorbid pairs. The higher relative risks of Rheumatoid arthritis -Osteoporosis and Parkinson's disease -Schizophrenia (compared with the other comorbid pairs) were not accompanied by a corresponding increase in the network similarity measures. Several factors may explain these variations in our analysis. Firstly, it has been shown that relative risk overestimates the comorbid associations between rare diseases and underestimates the associations between highly prevalent diseases [43] . The number of cases in the HuDiNe database for Rheumatoid arthritis -Osteoporosis and Parkinson's disease -Schizophrenia are 24629 and 5439 respectively, which can be classified as rare occurrences when compared with the other comorbid pairs. Additional File 17: Figure S14 shows the relationship of the relative risks of the nine pairs of diseases with the individual prevalence of the diseases and the prevalence of the disease pairs as comorbidities. Secondly, the human interactome is a progressively developing network with ~85% remaining to be discovered. Therefore, the inherent incompleteness of the human PPI network, sampling biases introduced as a result of the selective discovery of PPIs, and the tendency of such incomplete networks to exhibit small overlaps [60] could have led to the underestimation of the network overlaps. Our second key finding was that druggable proteins were highly enriched among the proteins shared between the networks of two comorbid diseases ( Table 2) .

Based on these results, we speculated that drug action on targets shared between the two diseases may give rise to contraindications in comorbidities. Interestingly, this hypothesis was only partially supported in our study.

The major finding in this respect was that the target network of the drugs used in the treatment of a specific disease A and contraindicated in a comorbid disease B showed preferential affiliation to proteins shared between the PPI networks of both the diseases or proteins uniquely found in the PPI network of the disease A, pathways shared by the two diseases or pathways associated with the disease A and tissues specifically associated with disease A ( Table 5) . As explained before, this was contrary to our hypothesis that these target networks would be preferentially affiliated with common mechanisms underlying the two diseases. This hypothesis was based on the assumption that adverse events stem from drugs inducing opposite pharmacological effects in comorbid diseases by targeting effectors that are shared between the two diseases. However, our findings indicate that mechanisms underlying the pathology of disease A may contribute to contraindications in the comorbid disease B.

Although further studies are required to examine the basis of this finding, it seems to indicate that the possibility of contraindications may be high when disease A drugs are highly specific to disease A in terms of the targeted PPI network, pathway and tissue. Instead, rational drug development should take into account the causative and correlational influences of the other comorbid conditions (disease B) that co-exist with disease A.

The target network of the drugs used in the treatment of a specific disease A and not contraindicated in a comorbid disease B showed preferential affiliation to proteins shared between the PPI networks of both the diseases or proteins uniquely found in the PPI network of the comorbid disease B, pathways shared between the two diseases or pathways associated with the comorbid disease B and tissues specifically associated with the comorbid disease B ( Table 5) . This was contrary to our expectation that these target networks would be preferentially affiliated with biological modalities pertaining to disease A. This conjecture was based on the assumption that for a drug to be specifically active against a specific disease A without aggravating a comorbid disease B, it had to reverse the phenotypes specifically associated with disease A. In this model, phenotypes of disease B were considered as 'off-targets' in line with the principles of conventional pharmacology, in which unintended effects of the drugs were attributed to interaction with pathways that may not be consequential to the pathology of disease A (i.e. pathways relevant to disease B) [13] . Our findings on the contrary indicate that the mechanisms underlying the pathology of the comorbid disease B may contribute to the therapeutic alleviation of disease A. Although further investigations may be necessary to dissect the basis for this observation, it is possible that an etiological association between the two diseases may cause their emergence or development to be interdependent. Specifically, future studies should concentrate on 3 etiological models of comorbidity [99] , namely, the direct causation For example, the alterations in such genes would have led to pathway perturbations in specific tissues, which if counteracted by the drugs, may lead to alleviation of disease A.

Despite disease A drugs contraindicated in disease B and disease A drugs not contraindicated in disease B showing preferential affiliation with disease A and disease B respectively, it was clear, at least in the case of the drug target and disease network analysis, that both these categories also showed affiliation with proteins shared between the two diseases ( Table 5 ). This is in line with the speculation that both beneficial and adverse outcomes of drug treatment may arise from shared effectors and pathways, and that it may be difficult to delineate the separate mechanisms underlying the two outcomes [13] . Future analysis should focus on biological variables with the potential to differentially affect the functions of such shared proteins, specifically their cellular, pathway and tissue landscapes.

Our current approach has some limitations. Firstly, our study is based on 6 pairs of diseases that were selected based on literature survey. Ideally, future studies must be expanded to include all the known pairs of comorbid disorders. Secondly, our analysis did not take the overlaps among the drug target networks into account; this would have allowed us to identify the network configurations of disease A -disease B -disease A drug not contraindicated in disease B -disease B drug not contraindicated in disease A. Secondly, although we were able to support our findings by citing evidence based on the known clinical activity of specific drugs, further investigations with the six comorbid disease pairs are essential to confirm the validity of our findings. These should focus on large-scale analysis of patient treatment data collected from observational studies and functional assays in animal models of human comorbidities.

In summary, our findings suggest that studies driven by biological modalities that influence comorbidities, such as disease PPI networks, pathways and tissue-specificity, are essential for rational drug development and minimization of adverse events. The results from our study have therapeutic applications and may directly benefit future assessments of drug contraindications in individuals with comorbidities.

We observed that the target networks of disease A drugs that were not contraindicated in disease B were mostly affiliated with the disease B network, and pathways and tissues associated with disease B. On the other hand, the target networks of disease A drugs that were contraindicated in disease B were affiliated with the disease A network, and pathways and tissues associated with disease A. This could indicate that etiological associations between the two diseases could play an active role in their therapeutic alleviation. In summary, our findings suggest that the enrichment patterns of drug target networks in pathways, tissues and the PPI networks of comorbid diseases will help identify drugs with/without contraindications in comorbidities.

Funding: This work was partially supported by INSA Senior Scientist grant of Prof. N. Balakrishnan.

Comorbidity and disease burden in the National Comorbidity Survey Replication

Comorbid chronic diseases, discordant impact on mortality in older people: a 14-year longitudinal population study

Multimorbidity in general practice: prevalence, incidence, and determinants of co-occurring chronic and recurrent diseases

Prevalence of multimorbidity among adults seen in family practice

Abdelalim A: Global, regional, and national disability-adjusted life-years (DALYs) for 359 diseases and injuries and healthy life expectancy (HALE) for 195 countries and territories, 1990-2017: a systematic analysis for the Global Burden of Disease Study

Respiratory effect of betablockers in people with asthma and cardiovascular disease: population-based nested case control study

Incidence of toxic epidermal necrolysis and Stevens-Johnson Syndrome in an HIV cohort

Impact of comorbidity on adverse drug reaction profile in a cohort of patients treated with artemisinin combination therapies for uncomplicated malaria in Nigeria

When good drugs go bad

Network medicine: a network-based approach to human disease

The fourth paradigm: data-intensive scientific discovery

Introducing the medical bioinformatics in

The emerging paradigm of network medicine in the study of human disease

Relating drug-protein interaction network with drug side effects

Analysis of drug-induced effect patterns to link structure and side effects of medicines

Target essentiality and centrality characterize drug side effects

Drug target identification using side-effect similarity

Network neighbors of drug targets contribute to drug side-effect similarity

Structure of protein interaction networks and their implications on drug design

Network-based prediction of drug combinations

Tissue-specific genes as an underutilized resource in drug discovery

A comprehensive functional analysis of tissue specificity of human gene expression

DrugBank: a knowledgebase for drugs, drug actions and drug targets

Data-driven prediction of drug effects and interactions

DGIdb: mining the druggable genome

Human protein reference database-2009 update

BioGRID: a general repository for interaction datasets

BisoGenet: a new tool for gene network building, visualization and analysis

DisGeNET: a comprehensive platform integrating information on human disease-associated genes and variants

A Unique Graph Similarity Metric for Anomaly Detection

Computing topological parameters of biological networks

Cytoscape: a software environment for integrated models of biomolecular interaction networks

A dynamic network approach for the study of human phenotypes

WebGestalt 2019: gene set analysis toolkit with revamped UIs and APIs

WebGestalt 2019: gene set analysis toolkit with revamped UIs and APIs

The Genotype-Tissue Expression (GTEx) pilot analysis: Multitissue gene regulation in humans

Asplund A: Tissue-based map of the human proteome

TSEA-DB: a trait-tissue association map for human complex traits and diseases

A next generation connectivity map: L1000 platform and the first 1,000,000 profiles

Replicable and coupled changes in innate and adaptive immune gene expression in two case-control studies of blood microarrays in major depressive disorder

Ontology-based meta-analysis of global collections of high-throughput public data

RNA-Seq workflow: gene-level exploratory analysis and differential expression

ClustVis: a web tool for visualizing clustering of multivariate data using Principal Component Analysis and heatmap

Uncovering disease-disease relationships through the incomplete interactome

Negative association between schizophrenia and rheumatoid arthritis

At issue: schizophrenia and rheumatoid arthritis: the negative association revisited

A nationwide study on the risk of autoimmune diseases in individuals with a personal or a family history of schizophrenia and related psychosis

Influence of comorbid conditions on asthma

Comorbidity patterns of anxiety and depressive disorders in a large cohort study: the Netherlands Study of Depression and Anxiety (NESDA)

The relationship between age of asthma onset and cardiovascular disease in Canadians

Hypertension and asthma: a comorbid relationship

Unrecognized heart failure in elderly patients with stable chronic obstructive pulmonary disease

Chronic obstructive pulmonary disease in patients admitted with heart failure

Prevalence and co-prevalence of comorbidities among patients with type 2 diabetes mellitus. Current medical research and opinion

Prevalence and recognition of obesity and its associated comorbidities: cross-sectional analysis of electronic health record data from a large US integrated health system

Osteoporosis in Rheumatoid Arthritis: Dangerous Liaisons. Frontiers in Medicine 2020

Increased Risk of Parkinson's Disease in Patients With Schizophrenia Spectrum Disorders

Principal component analysis characterizes shared pathogenetics from genome-wide association studies

Detecting shared genetic architecture among multiple phenotypes by hierarchical clustering of gene-level association statistics

Protein-protein interaction networks and biology-what's the connection?

The druggable genome

A comparative study of disease genes and drug targets in the human protein interactome

Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings

An integrative network-based approach for drug target indication expansion

Gα12/13 Mediates α 1-Adrenergic Receptor-Induced Cardiac Hypertrophy

Drug information handbook: a comprehensive resource for all clinicians and healthcare professionals. Lexi-Comp Incorporated

Role of serotonergic and noradrenergic systems in the pathophysiology of depression and anxiety disorders

Cholinergic signaling in the hippocampus regulates social stress resilience and anxietyand depression-like behavior

Orange book: Approved drug products with therapeutic equivalence evaluations

The overlap between anxiety, depression, and obsessive-compulsive disorder

Depression-inducing and antidepressive effects of neuroleptics

A review of therapeutic uses of mirtazapine in psychiatric and medical conditions. The primary care companion for CNS disorders

Pathophysiology and treatment of psychosis in Parkinson's disease

An exploratory retrospective evaluation of ropinirole-associated psychotic symptoms in an outpatient population treated for restless legs syndrome or Parkinson's disease

The role of the hypothalamic-pituitary-adrenal axis in neuroendocrine responses to stress

Adrenal gland volume in major depression: increase during the depressive episode and decrease with successful treatment. Archives of general psychiatry

Chronic stress induces adrenal hyperplasia and hypertrophy in a subregion-specific manner

Prolonged secretion of cortisol as a possible mechanism underlying stress and depressive behaviour

Plasma cortisol in depressive illness

Cortisol hypersecretion predicts early depressive relapse after recovery with electroconvulsive therapy

The corticotropin-releasing factor (CRF) hypothesis of depression: new findings and new directions

Graeff FG: Anxiety, panic and the hypothalamic-pituitary-adrenal axis

Elevated cortisol in older adults with generalized anxiety disorder is reduced by treatment: a placebo-controlled evaluation of escitalopram

Effects of eicosapentaenoic acid and fluoxetine on plasma cortisol, serum interleukin-1beta and interleukin-6 concentrations in patients with major depressive disorder

Salivary cortisol in women with major depressive disorder under selective serotonin reuptake inhibitors therapy. Archives of women's mental health

Ontology-based meta-analysis of global collections of high-throughput public data

Drug repositioning beyond the low-hanging fruits. Current Opinion in Systems Biology

Drug repurposing: progress, challenges and recommendations

NaviGO: interactive tool for visualization and functional similarity and coherence analysis with gene ontology

The spleen in local and systemic regulation of immunity

Spleen serves as a reservoir of osteoclast precursors through vitamin D-induced IL-34 expression in osteopetrotic op/op mice

Enlarged spleen detected by abdominal ultrasonography in patients with RA

Splenectomy in Felty's syndrome

Splenomegaly in rheumatoid arthritis

Splenic involvement in rheumatic diseases

Improvement of rheumatoid arthritis following splenectomy for Felty syndrome

Enhancement of osteogenesis post splenectomy does not attenuate bone loss in ovariectomized rats

Effect of immune function changes after splenectomy on fracture healing in patients with compound injury

Defining comorbidity: implications for understanding health and health services

Author contributions: NBK conceived and designed the research. KBK designed and performed the analyses. NBK and MKG supervised the interactome-based analyses, and SKB and SJ provided scientific inputs on the biological aspects of the study. KBK prepared the manuscript and NBK, SKB, MKG and SJ edited the manuscript through extensive mutual consultations.

Data availability: The condition concept names from the TWOSIDES database that were used to categorize the drugs associated with each of the comorbid disease pairs, the drug lists, the proteins targeted by the drugs, and the genes associated with the comorbid disease pairs have been made available as Table S1, Table S2, Table S3 and Table S4 , respectively.