Mechanistic Models of Signaling Pathways Reveal the Drug Action Mechanisms behind Gender-Specific Gene Expression for Cancer Treatments

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.


Introduction
It has long been known that males and females present important differences that may influence the interpretation of traits, such as disease phenotypes [1,2] and their treatment [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 [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 [5]. Reviews of studies based on animal models reveal an over-representation of experiments based exclusively on males [6]. Moreover, in many experiments including male and female animals, the results were not analyzed by sex [7,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 [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 [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 [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 [12] or cell metabolism [13]. Mechanistic models have helped to understand the disease mechanisms behind different cancers [14,15], including neuroblastoma [16,17], breast cancer [18], rare diseases [19], complex diseases [20], the mechanisms of action of drugs [21,22], and other biologically interesting scenarios such as the molecular mechanisms that explain how stress-induced activation of brown adipose tissue prevents obesity [23] or the molecular mechanisms of death and the post-mortem ischemia of a tissue [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 [12].
Here, we propose the use of mechanistic models [13,14] of signaling activity related with cancer hallmarks [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 [13,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 [26]. Actually, in a recent work, our group has successfully predicted therapeutic targets in cancer cell lines with a precision over 60% [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 [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.

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 (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 [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 (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 [29]. Finally, the values were normalized between 0 and 1, as required by the signaling circuit activity mechanistic model [14].

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 [30]. The edgeR package [31] was used for this purpose.

Rationale of the Signaling Circuit Activity Mechanistic Model
Circuit activities are modelled as described in [14]. Pathways in the Kyoto Encyclopedia of genes and Genomes (KEGG) repository [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 [33][34][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: 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 [14].

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 [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 [37] annotations or Gene Cards [38] information on gene functionality.

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 [25] in the scientific literature, the CHAT tool [39], a text mining based application to organize and evaluate scientific literature on cancer, has been used to link gene names with cancer hallmarks.

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 Section 2.1 and uploaded in the Hipathia package. Then, these are transformed into the corresponding signaling circuit activity profiles, as explained above in 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 [30].

Drug Effect Simulation
The effect that a drug with known targets has over the different signaling circuits is simulated using the PathAct [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 [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 [30]. Table S1 contains the drugs that are used for each cancer in the simulation [40].

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.

Data Processing
Gene expression matrices were downloaded from the ICGC data portal (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 Table 1).
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 [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.

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 (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 ( 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.

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 (Table 1) was treated in the simulation with the specific drug(s) indicated (see 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. 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 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 [25], such as angiogenesis, DNA recombination, Cell cycle, apoptosis, etc. Figure 3 depicts the distribution of the most pervasive GS-DSA circuits across cancer types. 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. 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, 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.  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.

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 [41], cabozantinib [42], gefitinib [43], lapatinib [44], nilotinib [45], ruxolitinib [46], sorafenib [47], sunitinib [48], and trametinib [49]. In addition, for vemurafenib [50] and sonidegib [51], a low gender effect was also demonstrated, although not enough for different dose indications (see Table S1). 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 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 Table S3). 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".

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 [40,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 [14,15,17,20,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 [54,55]. However, other differences, such as chemotherapy [56] or targeted therapy response [57,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 [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 [60,61]. The ErbB signaling pathway, whose regulation is highly associated with estrogen and androgen levels [62,63] and Oxytocin signaling pathway, associated with vasopressin, a known sex dependent pathway [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 (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 [65][66][67][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 [69,70], a different pattern of copy-number aberrations [71] or oxidative stress [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 [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 [74][75][76]. Since some cancers present a high involvement of lipid metabolism in tumor initiation and progression, considering the intrinsic gender differences seems logical [77][78][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 [80,81], and in general, risk [82,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 [71,[84][85][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 [87], methylation pattern [88,89], or even innate immune response [90,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 [92,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 [94], and some studies suggest that these gender-associated pigmentation differences also occur in humans [95][96][97]. Moreover, skin hyperpigmentation is indeed more frequent in women and may be linked to sexual hormones [98,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 [58,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 [27]. In particular, 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 [34,53,101,102].

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 [103], and ultimately paves the way to suggest highly specific, individualized therapeutic interventions [26,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.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4409/9/7/1579/s1, Figure S1: Heatmap with the simulation of the effect of the drug on the circuits, Table S1: Drugs that are used for each cancer in the simulation; Table S2: Circuits affected by the action of the different indicated drugs used in each cancer. Table S3: Results of the simulations of drug effects over any individual circuit for the drugs with targets within the corresponding circuit. Funding: This work is supported by grants SAF2017-88908-R from the Spanish Ministry of Economy and Competitiveness and "Plataforma de Recursos Biomoleculares y Bioinformáticos" PT17/0009/0006 from the ISCIII, both co-funded with European Regional Development Funds (ERDF) as well as H2020 Programme of the European Union grants Marie Curie Innovative Training Network "Machine Learning Frontiers in Precision Medicine" (MLFPM) (GA 813533) and "ELIXIR-EXCELERATE fast-track ELIXIR implementation and drive early user exploitation across the life sciences" (GA 676559).