Abstract Despite the existence of differences in gene expression across numerous genes between males and females having been known for a long time, these have been mostly ignored in many studies, including drug development and its therapeutic use. In fact, the consequences of such differences over the disease mechanisms or the drug action mechanisms are completely unknown. Here we applied mechanistic mathematical models of signaling activity to reveal the ultimate functional consequences that gender-specific gene expression activities have over cell functionality and fate. Moreover, we also used the mechanistic modeling framework to simulate the drug interventions and unravel how drug action mechanisms are affected by gender-specific differential gene expression. Interestingly, some cancers have many biological processes significantly affected by these gender-specific differences (e.g., bladder or head and neck carcinomas), while others (e.g., glioblastoma or rectum cancer) are almost insensitive to them. We found that many of these gender-specific differences affect cancer-specific pathways or in physiological signaling pathways, also involved in cancer origin and development. Finally, mechanistic models have the potential to be used for finding alternative therapeutic interventions on the pathways targeted by the drug, which lead to similar results compensating the downstream consequences of gender-specific differences in gene expression. Keywords: mechanistic models, gene expression, signaling pathways, signal transduction, cancer therapies, drug mechanism of action, gender bias 1. Introduction It has long been known that males and females present important differences that may influence the interpretation of traits, such as disease phenotypes [[30]1,[31]2] and their treatment [[32]3]. In fact, since the introduction of microarrays, which allowed a systematic screening of the molecular differences between sexes, the existence of a large degree of sex-biased gene expression that could explain the molecular basis of such phenotypic differences became apparent [[33]4]. However, in most studies, sex is ignored or is not properly taken into account despite a vast majority of common diseases displaying clear sex differences in symptoms or prevalence [[34]5]. Reviews of studies based on animal models reveal an over-representation of experiments based exclusively on males [[35]6]. Moreover, in many experiments including male and female animals, the results were not analyzed by sex [[36]7,[37]8]. In spite of this, it has been suggested that just adding sex as a variable could lead to conceptual and empirical errors in research unless differences between human men and women are properly modeled [[38]9]. Thus, understanding the molecular basis of these differences is of utmost importance to identify the functional mechanisms behind them and being able to distinguish real sex-dependent cell activities from those ones due to confounding variables. Accordingly, in a recent study that revealed a strong gender-specific bias gene expression in osteoarthritis, conventional pathway enrichment analysis showed that female specific miRNAs were estrogen responsive and targeted genes in toll-like receptor signaling pathways, suggesting mechanistic links between inflammation and osteoarthritis [[39]10]. In addition, recently, the discovery of differences in a brain signaling pathway involved in reward learning and motivation that make male mice more vulnerable to autism seems to provide a mechanistic explanation on why autism spectrum disorders are more common in males [[40]11]. Therefore, a proper interpretation of the effect that differences in gene expression have over phenotypes, such as drug response or disease progression, involves understanding the mechanisms of the disease or the mode of action of drugs, which can be interpreted through mechanistic models of cell signaling [[41]12] or cell metabolism [[42]13]. Mechanistic models have helped to understand the disease mechanisms behind different cancers [[43]14,[44]15], including neuroblastoma [[45]16,[46]17], breast cancer [[47]18], rare diseases [[48]19], complex diseases [[49]20], the mechanisms of action of drugs [[50]21,[51]22], and other biologically interesting scenarios such as the molecular mechanisms that explain how stress-induced activation of brown adipose tissue prevents obesity [[52]23] or the molecular mechanisms of death and the post-mortem ischemia of a tissue [[53]24]. Among the few available proposals of mechanistic modeling algorithms that model different aspects of signaling pathway activity, Hipathia has demonstrated having superior sensitivity and specificity [[54]12]. Here, we propose the use of mechanistic models [[55]13,[56]14] of signaling activity related with cancer hallmarks [[57]25], other cancer-related signaling pathways, and some extra relevant cellular functions to understand the functional consequences of the gender bias in gene expression. Such mechanistic models use gene expression data to produce an estimation of profiles of signaling or metabolic circuit activity within pathways [[58]13,[59]14]. An interesting property of mechanistic models is that they can be used not only to understand molecular mechanisms of disease or of drug action but also to predict the potential consequences of gene perturbations over the circuit activity in a given condition [[60]26]. Actually, in a recent work, our group has successfully predicted therapeutic targets in cancer cell lines with a precision over 60% [[61]15]. Therefore, we will use this mechanistic framework to understand what is the molecular basis of the different effects of cancer drugs by directly simulating their effect in the patients. This approach has recently been used by us to understand the generation of resistances in cancer at the single cell level in glioblastoma [[62]27]. Therefore, circuit activity, which can easily be linked to specific cell functionalities, has been used here to discover the different molecular mechanisms triggered by the biased gene expression between human males and females in cancers and, what is even more interesting, in their differential response to treatments. 2. Materials and Methods 2.1. Data Source, Selection Criteria, and Data Preprocessing Gene expression data from patients belonging to The Cancer Anatomy Genome Project (TCGA) were downloaded from the International Cancer Genome Consortium (ICGC) data portal ([63]https://dcc.icgc.org/). To create datasets containing males and females with features as homogeneous and unbiased as possible, the Propensity Score Matching (PSM) technique [[64]28] (MatchIt R package) was used. This methodology allows selecting samples that are different in gender but as similar as possible in the rest of the relevant features. To achieve so, first, a logistic regression model of gender (male/female) was created and regressed on the following covariates: age at initial pathologic diagnosis, histological type, pathologic stage, neoplasm histologic grade, race, tobacco smoking history, and tumor purity. All the covariates were taken from the ICGC data portal, except the tumor purity that was available at Synapse ([65]https://www.synapse.org/#!Synapse:syn3242754). Next, the unbiased samples which have matching covariate weight profiles were selected. The trimmed mean of M-values (TMM) method (18) was used for the normalization of gene expression data originally obtained as gene read counts of samples. Normalized samples were log-transformed and truncation by quantile 0.99 was applied. Batch effect was corrected with COMBAT [[66]29]. Finally, the values were normalized between 0 and 1, as required by the signaling circuit activity mechanistic model [[67]14]. 2.2. Differential Gene Expression To obtain differentially expressed genes (DEG) between conditions compared (normal versus cancer or male versus female), a negative binomial generalized log-linear model was used after gene expression normalization. The p-values were adjusted using the False Discovery Rate (FDR) method [[68]30]. The edgeR package [[69]31] was used for this purpose. 2.3. Rationale of the Signaling Circuit Activity Mechanistic Model Circuit activities are modelled as described in [[70]14]. Pathways in the Kyoto Encyclopedia of genes and Genomes (KEGG) repository [[71]32] are used to define circuits that connect any possible receptor protein to specific effector proteins that are ultimately responsible for triggering cell activities. A total of 98 KEGG pathways involving a total of 3057 genes that form part of 4726 protein nodes were used to define a total of 1287 signaling circuits. Normalized gene expression values are used as proxies of protein activity [[72]33,[73]34,[74]35]. The intensity value of signal transduced to the effector is estimated by starting with an initial signal with an arbitrary value of 1 in the receptor, which is propagated along the nodes of the signaling circuits according to the following recursive equation: [MATH: Sn=υ n·1−< mstyle mathsize="140%" displaystyle="true">∏sa ∈A1−sa⋅∏si ∈I1−si :MATH] (1) where S[n] is the signal intensity for the current node n, v[n] is its normalized gene expression value, A is the set of activation signals (s[a]), arriving to the current node from activation edges, and I is the set of inhibitory signals (s[i]) arriving to the node from inhibition edges [[75]14]. Here the Hipathia R/Bioconductor package ([76]https://doi.org/doi:10.18129/B9.bioc.hipathia), which implements the Hipathia model, is used. Additionally, a web server implementation is also freely available at: [77]http://hipathia.babelomics.org/. 2.4. Cell Functional Output Triggered by the Signaling Circuit The effector nodes at the end of the circuits trigger cell functionalities. The functionality of the circuit has been annotated as the function that the effector performs. Such functionalities have been taken from the Uniprot [[78]36] annotations. In the case of ambiguity (e.g., the general term of apoptosis can refer to its activation or repression), the Uniprot annotations were refined by manual curation using more detailed Gene Ontology [[79]37] annotations or Gene Cards [[80]38] information on gene functionality. 2.5. Association of Signaling Circuits Activities to Cancer Hallmarks As explained below, each effector is known to be associated with one or several cell functions. Since these effector genes have been related specifically with one or several cancer hallmarks [[81]25] in the scientific literature, the CHAT tool [[82]39], a text mining based application to organize and evaluate scientific literature on cancer, has been used to link gene names with cancer hallmarks. 2.6. Estimation of the Differential Signaling Activity The Hipathia R/Bioconductor package was used to test for differential signaling activity between male and female samples. Gene expression profiles are normalized as described above in [83]Section 2.1 and uploaded in the Hipathia package. Then, these are transformed into the corresponding signaling circuit activity profiles, as explained above in [84]Section 2.2. Finally, Hipathia applies a Wilcoxon test to check for significant differences in the activity of the circuits. The p-values are corrected for multiple testing using FDR [[85]30]. 2.7. Drug Effect Simulation The effect that a drug with known targets has over the different signaling circuits is simulated using the PathAct [[86]26] strategy. Briefly, the original gene expression profiles of the patients are taken as reference set and a simulated set of pseudo-treated patient gene expression profiles is generated by substituting the gene expression value(s) of the gene(s) targeted by the drug by a very low value (0.001) that simulates the inhibition of the drug. That is, the gene, even if it is expressed, is substituted by an “almost no expressed” gene (equivalent to an inhibited gene product. The reason for simulating the inhibition with an arbitrarily low value and not with a 0 is because in this way the simulation is more realistic (probably it is never an absolute inhibition) and, on the other hand, it preserves some basal low activity value in the circuit contributed by the rest of the genes, that it is useful for testing purposes (see [[87]26] for details). Then, the HiPahtia R/Bioconductor package is used to generate the corresponding signaling profiles for the reference patient set and the pseudo-treated patient set, that are further compared and tested for differences with a Wilcoxon test. The p-values are corrected for multiple testing using FDR [[88]30]. [89]Table S1 contains the drugs that are used for each cancer in the simulation [[90]40]. 2.8. Differential Drug Effect between Male and Female Patients In this case, for each cancer type and each drug, we will have two paired datasets of male and female patients untreated and with the simulation of the treatment. For each paired dataset, we can test whether the effect of the drug significantly affects any signaling circuit or not. However, we are looking for differences at circuit level when we compare the male versus the female datasets. Then, each paired comparison untreated versus simulated drug treatment produces a distribution of fold changes for each circuit in each patient. To check for gender-specific differences in drug treatments, we simply compare the mean fold change values obtained for male and female patients. 3. Results 3.1. Data Processing Gene expression matrices were downloaded from the ICGC data portal ([91]https://dcc.icgc.org/) and processed as described in Methods. After the application of the PSM method to these data, a total of 3327 tumor samples corresponding to 13 different cancer types, containing samples of both genders with males and females with similar covariates, were used in the study (see [92]Table 1). Table 1. Cancer types used in this study. Cancer Code Cancer Type Female Male Sample Size Proportion (Male/Female) BLCA Bladder urothelial carcinoma 57 202 259 3.54 COAD Colon adenocarcinoma 113 207 320 1.83 GBM Brain Glioblastoma Multiforme 37 89 126 2.41 HNSC Head and Neck squamous cell carcinoma 97 328 425 3.38 KIRC Kidney renal clear cell carcinoma 124 314 438 2.53 KIRP Kidney renal papillary cell carcinoma 38 108 146 2.84 LGG Brain Lower Grade Glioma 104 205 309 1.97 LIHC Liver hepatocellular carcinoma 44 118 162 2.68 LUAD Lung adenocarcinoma 131 213 344 1.63 LUSC Lung squamous cell carcinoma 81 299 380 3.69 PAAD Pancreatic Cancer 30 77 107 2.57 READ Rectum adenocarcinoma 41 77 118 1.88 THCA Thyroid Carcinoma 66 127 193 1.92 Total 963 2364 3327 [93]Open in a new tab Profiles of normalized values of gene expression were then transformed into the corresponding profiles of signaling circuit activities upon the application of the Hipathia method [[94]14] that can be used to detect Gender-Specific Differential Signaling Activity (GS-DSA) by testing significant differential signaling activity between males and females in each cancer type. DEG between cancer and normal samples were also estimated for the cancers as described in Methods. 3.2. Gender-Specific Functional Differences in Cancer While all cancer types contain signaling circuits with gender-specific differential behavior, the distribution in the number of these circuits is remarkably asymmetric ([95]Figure 1). Specifically, READ, THCA, COAD, and GBM with 22, 52, 42, and 43 circuits, respectively, are cancers with a relatively small number of circuits with differential gender-specific activity, whereas, on the other extreme of the range, cancers like LUSC, KIRC, or HNSC with 224, 239, and 202 circuits, respectively ([96]Table 2). Although for most cancers the number of circuits displaying a significant GS-DSA is similar among males and females, in three cancers HNSC, LUAD, and LIHC, and to a lesser extent also in THCA and KIRP, the number of circuits displaying significant signal activity differences is much higher in females than in males when the effect of the drug is simulated. Figure 1. [97]Figure 1 [98]Open in a new tab Number of signaling circuits with significant gender-specific differential signaling activity (GS-DSA) in the different cancer types studied. (A) the number of circuits with significant GS-DSA in each cancer, decomposed into those in which the activity of the signaling circuit is higher in males than in females and vice versa. (B) after simulation of the drug treatment the number of circuits showing significant GS-DSA increases. Table 2. Gender-specific differential gene expression and signaling circuit activation across cancers. Cancer Type Cancer Cancer (M > F) Cancer (F > M) Drug Simulation Drug (M > F) Drug (F > M) Drug Diff. Cancer GBM 43 21 22 50 24 26 14 READ 22 5 17 34 13 21 19 PAAD 31 15 16 48 24 24 22 LGG 59 22 37 75 31 44 26 THCA 52 12 40 61 21 40 26 COAD 42 18 24 61 13 48 34 KIRP 145 57 88 180 78 102 77 HNSC 202 66 136 242 78 164 80 BLCA 104 56 48 161 92 69 89 LUAD 203 50 153 242 53 189 96 LIHC 168 36 132 212 35 177 100 LUSC 224 105 119 238 112 126 105 KIRC 239 98 141 301 141 160 107 [99]Open in a new tab While the number of signaling circuits showing differential gender-specific activity is proportional to the number of genes showing gender-specific differential expression ([100]Figure 2A), this gender-specific differential signaling activity seems to be only slightly related to the level of differential expression between cancer and the normal tissue ([101]Figure 2B) but completely unrelated to other relevant cancer parameters such as the mutational burden ([102]Figure 2C). Figure 2. [103]Figure 2 [104]Open in a new tab Relationships between. (A) gender-specific differential expressed genes (GS-DEG) and gender-specific differential signaling activity (GS-DSA); (B) differentially expressed genes between cancer and normal (DEG) and GS-DSA, and (C) average mutations (mutation burden) and GS-DSA. 3.3. Potential Differences in Drug Effects Due to Gender-Specific Functional Differences In order to understand the molecular mechanisms behind gender-specific differential effect of drugs, their effect was simulated individually in each patient as described in Methods. Each cancer type studied ([105]Table 1) was treated in the simulation with the specific drug(s) indicated (see [106]Table S1). The simulation produced a new set of profiles of signaling activity corresponding to the simulated treatment for the patients. For both male and female patients, the simulated treatment sets were compared to the corresponding reference patient sets. [107]Table S2 contains the circuits affected by the action of the different indicated drugs used in each cancer, both in male and female patients. Then, we were interested in circuits showing a significantly different effect of the drug between both sexes. The signaling circuits showing gender-specific differential behavior most pervasively across cancers occur only in a maximum of six cancer types simultaneously, which suggests a high heterogeneity in signaling programs across cancers. The cell functionalities triggered by the most pervasive gender-specific signaling circuits (presenting GS-DSA in at least four cancer types) are summarized in [108]Table 3. Most of the affected circuits belong to cancer-specific pathways, such as renal cell carcinoma, pancreatic cancer, prostate cancer, glioma, etc. There are also some physiological pathways such as ErbB (KEGG: hsa04012), p53 (KEGG: hsa04115), Apoptosis (KEGG: hsa04210), or VEGF (KEGG: hsa04370) signaling pathways. The functionalities affected can easily be mapped to cancer hallmarks [[109]25], such as angiogenesis, DNA recombination, Cell cycle, apoptosis, etc. [110]Figure 3 depicts the distribution of the most pervasive GS-DSA circuits across cancer types. [111]Figure 4 shows the cancer hallmarks most affected by the gender differences in gene expression and their consequences on signaling and ultimately in cell functionality. [112]Table S2 contains a comprehensive list of all the circuits, with details on the cancer hallmarks affected. The distribution of circuits showing GS-DSA is uneven across cancer types, with cancers with many circuits affected, such as BCLA, HNSC of KIRP, and others with only a few circuits affected by the gender activity bias, such as GMB or READ. Finally, [113]Figure 5 depicts a comprehensive map of relationships among cancers, signaling circuits, functions, and cancer hallmarks. Table 3. Circuits showing gender-specific differential signaling circuit activation in four or more cancers simultaneously. Effector Circuit Uniprot Annotation of Effector Circuits Cancers with GS-DSA Renal cell carcinoma: VEGFA * Angiogenesis BLCA, COAD, HNSC, KIRP, LIHC, LUAD Fanconi anemia pathway: RAD51 DNA recombination BLCA, HNSC, KIRP, LUAD, LUSC Fanconi anemia pathway: RAD51C DNA recombination BLCA, HNSC, KIRP, LUAD, LUSC Fanconi anemia pathway: BRCA1 DNA recombination; COAD, HNSC, KIRP, LUAD, LUSC Pathways in cancer: PTCH1 * Tumor suppressor BLCA, HNSC, KIRC, LUAD, LUSC Pancreatic cancer: E2F1 Apoptosis; Cell cycle BLCA, KIRP, LIHC, LUAD, LUSC Prostate cancer: RB1 Cell cycle BLCA, COAD, HNSC, KIRP, LUSC ErbB signaling pathway: RPS6KB1 Translation regulation HNSC, KIRC, LUSC, THCA ErbB signaling pathway: ELK1 Transcription; Transcription regulation BLCA, KIRC, KIRP, LIHC ErbB signaling pathway: STAT5A * Transcription; Transcription regulation KIRP, LIHC, LUAD, LUSC ErbB signaling pathway: ELK1 * Transcription; Transcription regulation BLCA, HNSC, LUAD, LUSC ErbB signaling pathway: CBLC Ubl conjugation pathway BLCA, LIHC, LUAD, LUSC ErbB signaling pathway: ERBB3 ERBB3 Cell differentiation BLCA, KIRC, KIRP, LUSC p53 signaling pathway: IGFBP3 Apoptosis KIRC, LGG, LIHC, THCA Apoptosis: BBC3 Apoptosis LGG, LIHC, PAAD, THCA Axon guidance: ILK Cell growth, Metastasis KIRC, KIRP, LUAD, READ VEGF signaling pathway: PTK2 Angiogenesis KIRP, LGG, LUAD, LUSC Oxytocin signaling pathway: CDKN1A Cell cycle BLCA, KIRC, LUSC, THCA Pathways in cancer: FIGF Angiogenesis BLCA, KIRP, LIHC, LUAD Pathways in cancer: FIGF * Angiogenesis BLCA, KIRC, LIHC, LUSC Proteoglycans in cancer: CCND1 Cell division; DNA damage KIRC, LIHC, LUAD, PAAD Proteoglycans in cancer: CDKN1A Cell cycle COAD, HNSC, KIRC, LUAD Proteoglycans in cancer: VEGFA * Angiogenesis HNSC, KIRP, LUAD, PAAD Proteoglycans in cancer: KDR ** Angiogenesis BLCA, HNSC, KIRP, LUAD Colorectal cancer: MAPK8 Biological rhythms GBM, KIRC, LIHC, LUSC Pancreatic cancer: MAPK8 Biological rhythms BLCA, COAD, LIHC, READ Glioma: E2F1 Apoptosis; Cell cycle BLCA, KIRP, LIHC, LUSC Glioma: E2F1 * Apoptosis; Cell cycle BLCA, HNSC, KIRP, LUSC Bladder cancer: RB1 Cell cycle BLCA, HNSC, KIRP, LUSC Acute myeloid leukemia: PIM1 Apoptosis; Cell cycle BLCA, LUAD, LUSC, THCA Small cell lung cancer: RB1 Cell cycle BLCA, HNSC, KIRP, THCA [114]Open in a new tab * and ** are used for disambiguation, it refers to effector genes occurring more than once in the same KEGG pathway. Figure 3. [115]Figure 3 [116]Open in a new tab Distribution of the most pervasive GS-DSA circuits across cancer types. * and ** are used for disambiguation, it refers to effector genes occurring more than once in the same KEGG pathway. Figure 4. [117]Figure 4 [118]Open in a new tab Cancer hallmarks affected by GS-DSA circuits across cancer types. Figure 5. [119]Figure 5 [120]Open in a new tab Relationships among cancers, signaling circuits, functions, and cancer hallmarks. * is used for disambiguation, it refers to effector genes occurring more than once in the same KEGG pathway. 3.4. Validation So far, the results depict the effects caused in signaling activity by the observed gender-biased differences in gene expression. However, the phenotypic consequences of these changes can be diverse in relevance and nature. In order to detect changes associated with drug effect, an exhaustive search in the literature has been done and for the following drugs a different activity in males and females was experimentally demonstrated: bevacizumab [[121]41], cabozantinib [[122]42], gefitinib [[123]43], lapatinib [[124]44], nilotinib [[125]45], ruxolitinib [[126]46], sorafenib [[127]47], sunitinib [[128]48], and trametinib [[129]49]. In addition, for vemurafenib [[130]50] and sonidegib [[131]51], a low gender effect was also demonstrated, although not enough for different dose indications (see [132]Table S1). [133]Table 4 lists the circuits that display GS-DSA when the effect of the drug has been simulated. The simulation has been made with drugs known to have a different effect for males and females, and it is important to note that [134]Table 4 only reports those circuits that were only differentially regulated in a drug with different activity between males and females and never in drugs which do not show this differential activity (see the whole list of drugs tested in [135]Table S3). [136]Figure S1 presents a comprehensive picture of the number of circuits affected and those that are relevant in the context of drug action. It is interesting to see how some drugs have an extensive gender-specific effect across many pathways, and, within them, across many signaling circuits, like ruxolitinib, while others seem to be circuit-specific, like bevacizumab or sorafenib. It is interesting to note that cancer-related pathways seem to be more pervasively affected by gender-specific differential activity in the drug simulation than physiologic pathways. Table 4. Simulation of the effect that drugs, with described gender bias, have over signaling circuits (described as pathway and the final effector of the circuit). Circuits for which a GS-DSA is detected after the simulation of the drug are marked with “Y”. Pathway Effector Bevacizumab Cabozantinib Gefitinib Lapatinib Nilotinib Ruxolitinib Sorafenib Sunitinib Trametinib Vemurafenib Sonidegib Ras signaling pathway BRAP . . . . . . . . . Y . cGMP-PKG signaling pathway MAPK1 . . . . . . Y . . . . cAMP signaling pathway MYL9, PTCH1, HHIP, ACOX1, F2R AMH, ORAI1, BAD, NFKBIA NFKB1, RYR2, GRIN3A, GRIA1, CFTR, SLC9A1, ATP2B1, CACNA1C, PDE3A, ATP1B4 FXYD1, RHOA, [137]C00165, [138]C01245, PAK1, MLLT4, [139]C00416, MAPK8, HCN4 . . . . . Y . . . . . Chemokine signaling pathway STAT1 . . . . . Y . . . . . Wnt signaling pathway JUN . . Y . . . . . . . . Hedgehog signaling pathway PTCH1, SMO, PTCH1, GLI1, HHIP, CCND1, BCL2, PRKACA, GLI1 SUFU, . . . . . . . . . . Y Axon guidance ILK . . . . . . . . . . Y VEGF signaling pathway: NOS3 Y . . . . . . . . . . Osteoclast differentiation: MAPK1 . . . . . . . Y . . . Osteoclast differentiation: NFKB1 . . . . Y . . . . . . Signaling pathways regulating pluripotency of stem cells HNF1A . . Y . . . . . . . . Jak-STAT signaling pathway BCL2, BCL2L1, MYC, AOX1, GFAP, MCL1, PIM1, CCND1 . . . . . Y . . . . . Natural killer cell mediated cytotoxicity TNF . . . . . . . . Y . . TNF signaling pathway CASP7, JUN, CEBPB . . . . . . . . Y . . Leukocyte transendothelial migration MAPK14 . . . . Y . . . . . . Inflammatory mediator regulation of TRP channels: TRPM8, TRPV4 . . . . . Y . . . . . Ovarian steroidogenesis STAR, HSD3B1, PLA2G4B, ACOT2, CYP19A1, HSD17B2, CYP19A1 . . . . . Y . . . . . Melanogenesis MITF . . . . . Y . . . . . Thyroid hormone synthesis TG . . . . . Y . . . . . Thyroid hormone signaling pathway STAT1, ESR1, THRB . . . . . . . . Y . . Adipocytokine signaling pathway AGRP, NPY, POMC, PPARGC1A, PTPN11 . . . . . Y . . . . . Regulation of lipolysis in adipocytes PLIN1, LIPE . . . . . Y . . . . . Aldosterone synthesis and secretion CYP11B2 . . . . . Y . . . . . AGE-RAGE signaling pathway in diabetic complications FOXO1, CCND1, NFATC1 . . . . Y . . . . . . Pathways in cancer CCND1 . . . . . . . . . Y . Pathways in cancer FIGF Y . . . . . . . . . . Pathways in cancer CCNA1, CSF3R, CSF2RA, CSF1R . . . . . . . . . . . Pathways in cancer CSF1R . . . . . . . Y . . . Pathways in cancer BMP2, GLI1, HHIP, PTCH1 . . . . . . . . . . Y Proteoglycans in cancer HSPB2 . . . . Y . . . . . . Proteoglycans in cancer: AKT3 . Y Y . . . . . . . . Proteoglycans in cancer: PRKCA . . . . . . Y . . . Colorectal cancer: MAPK8 . . . . . . . . Y . . Renal cell carcinoma VEGFA Y . . . . . . . . . . Renal cell carcinoma RAP1A . Y . . . . . . . . Renal cell carcinoma AKT3 . Y Y . . . . . . . . Pancreatic cancer RAC1 . . . . . . . . Y . . Pancreatic cancer [140]C00416 . . . . . . . . Y . . Basal cell carcinoma PTCH1 . . . . . . . . . Y Acute myeloid leukemia CCNA1, SPI1 . . . . . . . . . . . Non-small cell lung cancer FOXO3 . . . Y . . . . . . . [141]Open in a new tab 4. Discussion The differences associated with gender have been previously assessed in pan-cancer studies, most of them using TCGA cancer datasets, resulting in divergent patterns for sex bias in gene expression or immune features across multiple cancer types have been revealed [[142]40,[143]52]. Nevertheless, the functional consequences at the level of cell behavior or fate of gender bias in gene expression have remained mainly unknown. To our knowledge, this is the first time that such gender specific differences in gene expression are evaluated in the context of perturbation response, taking into consideration cell mechanisms as a whole, an approach that has successfully been used to explain different cancer molecular mechanisms [[144]14,[145]15,[146]17,[147]20,[148]53]. Differences in cancer epidemiology, susceptibility, and prognostics have been widely described, but exactly why this occurs at a molecular level has been poorly understood. Many cancers show dissimilarities in incidence and mortality rates associated to sex-specific disparities; some can be the result of different hormone levels, especially estrogen, or sexual chromosome dose [[149]54,[150]55]. However, other differences, such as chemotherapy [[151]56] or targeted therapy response [[152]57,[153]58], are the result of more complex cell processes that need to be evaluated in its cell mechanism context in order to be able to detect patterns. When evaluating individual circuits and pathways, most of them are indeed related with several cancer processes, apoptosis, and proliferation. Signaling circuits of the Fanconi Anemia pathway, involved in DNA repair [[154]59], and therefore the genome instability and mutation cancer hallmark, show the highest values of GS-DSA. Other relevant signaling circuits belong to the proteoglycans in cancer pathways. Proteoglycans abundance, their metabolism and their relationship with genomic instability is clearly gender-related [[155]60,[156]61]. The ErbB signaling pathway, whose regulation is highly associated with estrogen and androgen levels [[157]62,[158]63] and Oxytocin signaling pathway, associated with vasopressin, a known sex dependent pathway [[159]64], also contain signaling circuits showing significant GS-DSA. In order to assess the gender-specific differences in global mechanisms across cancers, we grouped the circuits by cancer hallmark ([160]Figure 4). Most of the hallmarks showed a gender-specific perturbation response pattern, but above all we find sustaining proliferative signaling, resisting cell death and evading growth suppressors hallmarks, mainly composed of circuits belonging to cell proliferation and apoptosis pathways. Indeed, gender-related differences in cell proliferation and differentiation pathways have already been described for several tissues in humans [[161]65,[162]66,[163]67,[164]68]. Genome instability hallmark shows a considerable number of GS-DSA circuits across cancers as well, which is concordant with previous studies, showing gender-associated differences in expression of genes involved in the aforementioned DNA repair [[165]69,[166]70], a different pattern of copy-number aberrations [[167]71] or oxidative stress [[168]72]. Moreover, available data suggest that sex influences measures of age-associated genomic instability, which increases in both males and females with age. However, how sex affects genome instability is less clear, as tissue studied, genetic background, and the method selected can influence results immensely, as well as environmental factors that are difficult to address [[169]73]. Another cell process showing GS-DSA is lipid metabolism. It is well known that fat pad shows a different pattern in females and males, and these differences are the result of differences in metabolism at a molecular level [[170]74,[171]75,[172]76]. Since some cancers present a high involvement of lipid metabolism in tumor initiation and progression, considering the intrinsic gender differences seems logical [[173]77,[174]78,[175]79]. Besides the general deregulation of hallmarks, some of them are of special relevance in certain cancers. Particularly important is the role of angiogenesis, tumor promoting inflammation, and metastasis, which shows a clear pattern of GS-DSA in PAAD, KIRC, and LUAD, where some gender-specific events have already been described, as mutations [[176]80,[177]81], and in general, risk [[178]82,[179]83]. KIRC, together with HSNC and LUSC, are the cancers with the highest number of gender-specific signaling circuits, all of them showing sex-dependent differences in prognosis, mortality, and treatment response, as well as in molecular characteristics associated with them [[180]71,[181]84,[182]85,[183]86]. It is interesting to note that some cancers, such as PAAD, LUAD, and LUSC, are highly influenced by environmental factors, and therefore the gender differences might not be of physiological origin but rather could be determined by gender-specific lifestyles. In order to evaluate the implications of gender in cancer, there are many factors that need to be addressed beyond the scope of this work, such as the possible implication of sex-biased transcription factors, miRNAs expression [[184]87], methylation pattern [[185]88,[186]89], or even innate immune response [[187]90,[188]91]. However, no one can argue that biological intrinsic differences between females and males exist, and these differences are influencing all kind of cell behavior, thus the identification of the processes underlying these differences will facilitate the exploration of sex-biased disease susceptibility and therapy. Many of the proliferation-related circuits targeted by drugs presenting sex-bias are indeed regulated by estrogen in a direct or indirect manner, such as Ras, cGMP-PKG, or cAMP signaling and all the circuits related to MAPK proteins, as osteoclast differentiation [[189]92,[190]93]. Therefore, these circuits can show a different response depending on estrogen and other hormones levels, which are highly variable between sexes. Moreover, as aforementioned, lipid metabolism and proteoglycans may be influenced by hormone levels and by sex. Interestingly, melanogenesis, the synthesis of melanin and responsible of pigmentation, has been demonstrated to be regulated by hormones in a sex-specific manner in model organisms [[191]94], and some studies suggest that these gender-associated pigmentation differences also occur in humans [[192]95,[193]96,[194]97]. Moreover, skin hyperpigmentation is indeed more frequent in women and may be linked to sexual hormones [[195]98,[196]99]. The results presented here highlight the fact that gender needs to be considered when choosing the appropriate treatment in cancer. The approach presented here, based on mechanistic modeling of cell signaling pathways, shows its potential in evaluating the gender-specific differences in certain mechanisms of action of several drugs, and, therefore, in predicting potential non-responders and resistances [[197]58,[198]100]. It has also been shown that mechanistic models can be an excellent tool for the simulation of the effect of drugs, as we have recently demonstrated in cancer [[199]27]. In particular, [200]Table 4 shows circuits that display GS-DSA after the simulation of drugs for which a gender-specific activity has been reported, but never display GS-DSA in simulations of drugs with similar activity in both genders. In this way, the modeling framework used here provides the mechanistic link between the effect of the drug at a molecular level and at a phenotypic level. Like in any other study based on gene expression, it must be considered that any port-transcriptional modification, which can be relevant in cancer, is not primarily captured in the data. However, if such modification has an effect on the behavior of the neoplastic cell, it will be better detected, even indirectly, by its impact in the global signaling pattern of the cell, rather than by a conventional gene-centric analysis. In any case, mechanistic modeling can also be applied to proteomic or phosphoproteomic data, which would better account for the effect of post-transcriptional modifications in protein activity. However, given the difficulty of obtaining direct measurements of protein levels, an extensively used proxy for protein presence is the observation of the corresponding mRNA within the context of the module [[201]34,[202]53,[203]101,[204]102]. 5. Conclusions The use of mechanistic models that quantify cell behavioral outcomes provides a unique opportunity to understand the molecular mechanisms of cancer development and progression [[205]103], and ultimately paves the way to suggest highly specific, individualized therapeutic interventions [[206]26,[207]104]. Here we demonstrate how mechanistic models are suitable for uncovering the functional consequences that the gender-biased gene expression triggers downstream signaling circuits. Mechanistic models offer an opportunity to reconsider alternative targets on the pathways relevant for the therapeutic interventions that lead to similar results compensating the downstream consequences of the gender-specific differences in gene expression. Acknowledgments