Kinetic Modeling–Based Detection of Genetic Signatures That Provide Chemoresistance via the E2F1-p73/DNp73-miR-205 Network

Resistance to genotoxic drugs as the major cause of cancer therapy failure is a serious problem for oncologists and their patients that requires the understanding of the pivotal triggering events and its evolution when cancer progresses. Chemoresistance can be innate or acquired and may attributable to a single agent or a class of drugs. Potential mechanisms to counteract the therapeutic effects of DNA-damaging agents include the reduction of effective drug concentrations via enhanced efflux, detoxification enzymes or drug sequestration, modification of drug targets, changes or mutation in mitotic checkpoint signals, and hyperactivation of DNA repair mechanisms (1).

The cellular transcription factor E2F1 is a unique member of the E2F family of proteins as it regulates the tumor suppressor p53 and its homologue p73 thereby promoting apoptosis by the activation of a plethora of death pathways (2). One of the first signals recognized to induce E2F1′s apoptotic activity was DNA damage (3). Overexpression of E2F1 was shown to increase the sensitivity of malignant cells to apoptosis upon genotoxic treatment (4), which is critical for tumor growth reduction (5). In this regard, activation of p73 by E2F1 functions as backup once p53 is defective to ensure that damaged cells can undergo apoptosis (6,7). In such a scenario, high levels of endogenous E2F1 should increase the effectiveness of DNA-damaging agents and discriminate between tumors clinically responsive or resistant to anticancer therapy. This is in line with clinical surveys that associated E2F1 expression with improved survival in patients treated with adjuvant chemotherapy and underscores E2F1′s role as an endogenous chemosensitizer in patients with cancer and beyond, as a predictive factor for therapeutic success and clinical outcome (8).

This beneficial view on E2F1 was thwarted as other studies showed that overexpression or amplification of the E2F1 genomic locus in conjunction with RB1 loss frequently observed during cancer progression is associated with metastatic disease and chemoresistance (9). From this perspective, aberrantly elevated E2F1 levels in high-grade tumors can be considered as a marker for unfavorable patient survival prognosis. In support of these findings, we have recently shown that dysregulated E2F1 causes malignant progression in therapy-resistant metastatic melanoma xenografts in which depletion of endogenously high E2F1 levels abrogates tumor invasion and pulmonary metastasis (10). According to our results, the aggressive behavior of E2F1 in melanoma cells is partially mediated through the induction of the EGF receptor (EGFR) pathway.

It is well established that E2F1 stimulates the expression of the tumor suppressor p73 and its N-terminally truncated isoform (named DNp73) via direct transactivation of the TP73 gene (11). While full-length p73 inhibits cancer development by inducing cell-cycle arrest and apoptosis through its ability to bind p53 DNA target sites, DNp73 acts as antagonist of wild-type p53 family members by either directly interfering with DNA binding or forming inactive heteromeric complexes with transcriptionally active p73 (12). Strikingly, in melanoma metastases DNp73 isoforms are strongly upregulated in conjunction with full-length p73 compared with primary tumors (13). Because of the previous results, DNp73 exhibits antiapoptotic activity in human melanoma cells and specific suppression of individual isoforms enhances the sensitivity towards cytotoxic drugs (14). Investigating mechanisms responsible for dysregulated E2F1 losing its apoptotic function in human skin cancer, we identified miR-205 as a specific target of p73 and found that upon genotoxic stress its expression is sufficiently abrogated by endogenous DNp73 (15). Significantly, we showed that metastatic cells can be rescued from drug resistance by selective knockdown of DNp73 or overexpression of miR-205 in p73 depleted cells, leading to increased apoptosis and the reduction of tumor growth in vivo. These results suggest the E2F1-p73/DNp73-miR-205 network as a crucial mechanism for chemoresistance and thus as a target for prevention of metastasis.

The structure of the E2F1-centered biochemical network that mediates the potential resistance against anticancer therapies is a complex system involving signaling and transcriptional regulation. This network is enriched with network motifs, including feedback and feedforward loops. Furthermore, these network motifs are well interconnected, overlap, and cross-talk with other cancer-related signaling pathways (4, 16). Networks containing several of these motifs often show nonintuitive regulatory patterns, which require the use of mathematical modeling to understand their function and regulation (17,18).

In this work, we constructed a regulatory map that summarizes the current literature on E2F1 signaling dysregulation in cancer and complemented it with information about E2F1 and miR-205 targets, which have known or putative roles in drug response. On the basis of this map, we derived a mathematical model and used it to investigate under which conditions dysregulation of this network provides tumor cells with resistance mechanisms against several families of anticancer drugs. In connection to this, we further explored the consequences of tumor heterogeneity in the effectiveness of these therapies. Our results indicate that kinetic modeling under the systems biology paradigm can be an effective method to detect genetic signatures providing tumor cells with chemoresistance mechanisms and thereby to design therapeutic strategies to counteract them.

The human melanoma cell line SK-Mel-147 characterized by morphology and cytogenetics was obtained from M. Soengas (Department of Dermatology, University of Michigan, Comprehensive Cancer Center, Ann Arbor, MI). Tested and authenticated cells were cultured as described (15). Plasmid miRZIP-205 encoding antagomir-205 was purchased from BioCat. Transfections were conducted with Lipofectamine 2000 (Invitrogen). Adenoviral vector with shE2F1 was described previously (10). Expression of hsa-miR-205 and hsa-RNU6B (for normalization) was detected using specific TaqMan assays (Applied Biosystems). PCR, Western blot, and apoptosis assays were conducted as described. The wet lab data were obtained previously and are published in (15).

Information about the network components and their interactions was extracted from relevant published reports and databases (HPRD and STRING for protein–protein interactions; TRED for E2F1 transcriptional targets; TarBase 6.0, miRTarBase, miRecords and miR2Disease for validated targets of miR-205; and TransmiR for transcription factors regulating the expression of miR-205). Putative miR-205 targets and miR-205 transcription regulators were extracted from additional databases (see Supplementary Materials and Methods for a detailed description of the data retrieval process). As we were interested in the role of E2F1 and miR-205 in pro- and antiapoptotic processes as well as in drug responses, we filtered the obtained set of molecular interactions based on the corresponding Gene Ontology (GO) terms (GO:0043065, GO:0043066 and GO:0042493). Finally, all the information collected was integrated into a regulatory map in Cytoscape (available in Supplementary Materials and Methods).

Representative interactions in the constructed regulatory map were translated into a kinetic model based on ordinary differential equations. The model constructed accounts for the dynamics of miRNA, mRNA and protein expression levels, as well as cell populations in different cancer-relevant scenarios, including tumor proliferation and genotoxic or cytostatic drug administration. The model has the following general structure:

where Prot_i accounts for the expression level of the proteins involved in the network, mRNA_i for the messenger RNAs mediating their synthesis, miRNA_i for miRNAs repressing expression of given proteins, and C_i for the size of proliferative, arrested, or apoptotic cell populations. GS is an input signal that accounts for the effect of genotoxic or cytostatic drugs in the synthesis and regulation of proteins and mRNAs. D_k1, F_k1 and G_k2 account for the description of biochemical processes as rate equations that depend on the expression levels of proteins and mRNAs, or the time-dependent genotoxic stress signal GS. H_k3 are rate equations describing the dependence of cell population dynamics on the expression level of critical proteins in the network. Values of the model parameters characterizing the reaction rates were assigned following a hybrid strategy, composed of: (i) extracting parameter values from published information (e.g., protein, mRNA, and miRNA half-lives); (ii) estimating a subset of parameter values by tuning them to fit published quantitative and qualitative data; (iii) reducing some parameters by normalizing variables around the basal, nonstressed levels of mRNA and protein. Model equations and chosen parameter values are provided in Supplementary Materials and Methods.

To search for genetic signatures providing tumor cells with a chemoresistant phenotype, we produced in silico (kinetic model-based) cell lines by randomly perturbing values of selected model parameters using Latin hypercube sampling. Thereby, we generated 10⁵in silico cell lines, each one with a distinctive set of parameter values accounting for an individual genetic background.

For each in silico cell line, we conducted simulations in 3 cancer relevant scenarios: (i) nonstressed growth conditions; (ii) apoptosis under; and (iii) proliferation under cytostatic drug administration. We investigated genetic signatures providing the tumor cells with a chemoresistant phenotype, by randomly perturbing parameter values and simulating the response of the model to different cancer-relevant scenarios. We analyzed the effects of perturbations in the expression levels of critical network components by iteratively modifying their synthesis rate constants accounting for their induction and simulating the behavior of the network. We assessed the efficacy of conventional anticancer treatment in in silico cells with different E2F1-related genetic signatures by simulating periodic cycles of genotoxic or cytostatic drug injection. Model calibration, computational simulations, and data analysis were conducted using MATLAB running on a high-end workstation Fujitsu Celsius V840, 2x CPU AMD Opteron 2214, 2.2 GHz. ECC. Complete information about the construction of the regulatory map, as well as about the model structure, calibration, and analysis is included in Supplementary Materials and Methods.

Our recent results suggest that the E2F1-p73/DNp73-miR-205 pathway is involved in the regulation of pro- and antiapoptotic genes, which confers chemoresistance (15). We investigated whether this system is part of a wider and more complex network regulating the response of tumor cells to anticancer drugs. Towards this end, we constructed a regulatory map that contains the E2F1-p73/DNp73-miR-205 pathway as core module and includes genes targeted by E2F1 or miR-205, with known or putative relation to the drug response and drug-induced apoptosis (Fig. 1A). Overall, the map unravels a tightly interlocked regulatory system. Computational analysis of the network (81 nodes, 194 edges) shows a robust structure and a topology that is composed of many regulatory motifs, especially enriched in feedforward loops (for details, see Supplementary Materials and Methods). Interestingly, we detected several instances of network motifs in which a direct E2F1 target gene interacts with a miR-205–repressed target. Because E2F1 regulates miR-205 via p73, processes downstream of the interaction between both targets are regulated by E2F1 following 2 independent routes and can be considered as feedforward loop (Fig 1A). Several of those loops are associated with the efficient apoptosis initiation and involve the positive regulation of proapoptotic genes by E2F1 and the repression of antiapoptotic genes by miR-205. For example, the E2F1 target Hrk regulates apoptosis through interaction with the miR-205 targets BCL-2 and BCL-XL (Fig 1A; ref. 19). In addition, we detected at least one additional instance of this motif, in which E2F1 promotes EGFR expression and miR-205 represses ERBB3 (both belonging to the HER receptor family). In several tumors, both receptors are associated with abnormal proliferation and upon heterodimerization they influence sensitivity to cytostatic drugs like erlotinib (20). Furthermore, DNp73 can inhibit the expression of p73 transcriptional targets, including miR-205. The expression of these targets may depend on a tight trade-off between p73 and DNp73 expression and activity, which are both transcriptionally regulated by E2F1. These E2F1-p73/DNp73-target motifs can be considered as incoherent feedforward loops. Overall, the activation of miR-205 and therefore of its targets may depend on the interplay among E2F1, p73, DNp73, and other regulators of miR-205 like TWIST1 and TGFβ-1 signaling, as suggested by our regulatory map. Because of the network complexity, mathematical modeling is necessary to investigate nonintuitive regulation of this system in cancer-relevant scenarios.

We set up and analyzed a kinetic model using ordinary differential equations. To construct the model, we selected parts of the regulatory map that are representative for the network response to different anticancer drugs (Fig. 1B). The model is organized into four interconnected modules: the core regulatory module, two modules of transcriptional targets, and one for tumor cell population. The core regulatory module of the kinetic model is composed of E2F1, both isoforms of p73 (respectively, p73 and DNp73) and miR-205 (miR-205) and accounts for their temporal evolution under drug administration (15). We include here the transcriptional regulation of miR-205 expression by the known oncogene TGFβ-1 (TGFB1), which exemplifies the set of regulators of miR-205 expression external to the core module. The inclusion of this variable allows for the investigation of cross-talk, via miR-205, between E2F1 signaling and other oncogenic signals involved in the epithelial–mesenchymal transition (EMT). EMT has been associated with the emergence of chemoresistance (21).

In addition, the model includes a module for pro- and antiapoptotic transcriptional targets, which are regulated by the network components in response to genotoxic stress. Precisely: (i) for representing the dynamics of proapoptotic proteins whose expression is promoted by E2F1, we included a variable accounting for Harakiri expression (22); (ii) similarly, for proapoptotic proteins whose expression is positively regulated by p73 and negatively regulated by DNp73, a variable accounting for Bax expression is defined (23); and finally (iii) for antiapoptotic proteins whose expression is repressed by miR-205, we chose BCL2 as a representative and defined a corresponding variable (15). These pro- and antiapoptotic proteins are used here as an example and represent a wider set of proteins which undergo similar regulation. This kind of simplification has been successfully applied to reduce the complexity of models of biochemical networks, for example, in the study conducted by Aguda and colleagues (24). Furthermore, the model includes a module for receptors whose activity is regulated in response to cytostatic drugs. This module includes variables accounting for the dynamics of EGFR, whose expression is promoted by E2F1 (10), and ERBB3, which is repressed by miR-205 (25,26). In several tumors the heterodimerization of these receptors has been linked to abnormal cell proliferation and changes in the sensitivity of tumor cells to cytostatic drugs (20). Additional input variables account for the effect of genotoxic (GxD) and cytostatic (CyD) drugs triggering DNA damage-E2F1–mediated apoptosis and abolishment of tumor cell proliferation, respectively (4).

Finally, the last module includes a phenomenological kinetic equation, which connects the components of the network to the dynamics of the tumor cell population. In this equation, the representatives for pro- and antiapoptotic proteins regulate the cell death rate after genotoxic stress, whereas EGFR and ERBB3 control the tumor cell proliferation rate. Taken together, our model is able to describe the cascade of regulatory events after drug administration including stress signals, downstream intracellular responses, the regulation of pro- and antiapoptotic signals (for genotoxic drugs) or proliferative signals (for cytostatic drugs), and finally the global effect on tumor cell population (Fig. 1C; see Supplementary Materials and Methods for a complete description of the model).

In our analysis, we considered two families of drugs: genotoxic drugs, like doxorubicin and cisplatin, inducing apoptosis; and cytostatic drugs, like erlotinib and lapatinib, repressing tumor cell proliferation. In some scenarios, it has been observed that these drugs lose efficacy when some of the network components are dysregulated (15,27). With the help of the kinetic model, we analyzed whether specific genetic signatures of the core module can be linked to a phenotype response of resistance to these anticancer drugs. We define a genetic signature as a group of genes in a tumor cell whose combined expression pattern is linked to a specific phenotype (i.e., chemoresistant or chemosensitive). In line with this, we designed a nominal in silico chemosensitive genetic signature, inspired by that of SK-Mel-147 (15). This signature accounts for the expression of the network components for which an in silico cell line holds the phenotype of a chemosensitive tumor cell; that is, abnormal proliferation but responsiveness to genotoxic drugs (via triggering of apoptosis) and cytostatic drugs (via inhibition of proliferation, see Fig. 2B).

To detect genetic signatures providing chemoresistance, we generated a population of 104 model configurations obtained by random perturbation of the model parameter values of the nominal chemosensitive genetic signature. In the following, we define these distinct configurations as in silico cell lines. For each one of them, we simulate the model and obtained the expression levels of the network components and the tumor cell population in three cancer-relevant scenarios: (i) under nonstress conditions, at time zero, and after 120 hours; (ii) after genotoxic drug administration, at 48 hours; and (iii) after cytostatic drug administration, at 120 hours. The in silico cell lines were classified into the following groups: (i) cell lines that show abnormal proliferation in nonstress conditions, named as tumor cells; (ii) those resistant to genotoxic drugs; (iii) those resistant to cytostatic drugs; and (iv) those resistant to both genotoxic and cytostatic drugs. For each group, the data describing the expression levels of the network components at time zero were normalized with respect to the nominal chemosensitive values (Fig. 2C). For the computational definition of the simulation scenarios and classification of the resulting cell phenotypes, see Supplementary Materials and Methods.

Our results indicate that 30% of the in silico cell lines are resistant to genotoxic drugs (Fig. 2C). These cells display high basal levels (at time zero, nonstress conditions) of E2F1 and DNp73, whereas miR-205 is downregulated compared with the nominal chemosensitive tumor cell and the group of in silico tumor cells (Fig. 2A). In addition, approximately 1% of the cell lines display resistance to cytostatic drugs. These cells have a genetic signature composed of much higher basal E2F1, DNp73, and ERBB3 levels, whereas miR-205 appears strongly downregulated and TGFβ-1 is moderately upregulated. In both cases, EGFR is expressed at its maximum level, but our analysis indicates that this is a property already acquired by the group of tumor cells as a consequence of high levels of E2F1 (see Supplementary Materials and Methods). Finally, a small fraction of the cell lines display a phenotype of double resistance and their genetic signature is very similar to that of cytostatic drug–resistant cells.

Taken together, our results suggest a relevant role for E2F1, DNp73, and miR-205 in the regulation of resistance to several anticancer drugs. To validate this model-based prediction, we compared them with our recently published data on the interplay of E2F1 and miR-205 under genotoxic stress administration (15). To this end, we conducted simulations in which we iteratively perturbed the values of the parameter accounting for E2F1 and miR-205 induction and computed the percentage of apoptotic cells 48 hours after genotoxic drug administration. As shown in Fig. 3A, in silico cell lines with high levels of E2F1 display lower percentage of apoptotic cells in response to chemotherapy. Our previous experiments with melanoma SK-Mel-147 cells confirm the simulation results (Fig. 3B). SK-Mel-147 cells were treated with cisplatin and high percentage of cells became apoptotic (Scenario 1 in Fig. 3B). Inhibition of miR-205 with antagomir-205 reduces the percentage of apoptotic cells after treatment (Scenario 2), whereas selective knockdown of endogenously high E2F1 leads to a pronounced increase in apoptosis (Scenario 3). Finally, repression of both E2F1 and miR-205 reduces the apoptotic effect of cisplatin (Scenario 4). These experimental results confirm the model predictions displayed in Fig. 3A for the scenarios named with the same numbers.

To verify whether dysregulation of couples of network components affects chemosensitivity to genotoxic and cytostatic drugs, we conducted additional simulations in which we iteratively modified the values of the protein synthesis rates (Fig 4). Precisely: (i) because the p73/DNp73 ratio is known to play a major role in the aggressiveness of cancer cells (11, 28, 29), we iteratively modified the parameters accounting for p73 and DNp73 synthesis; ii) given that overexpression of E2F1 is associated with aggressiveness and chemoresistance (4,15), we also considered iterative modulation of E2F1 and DNp73 synthesis parameters to test whether their dysregulation suffices to induce resistance; (iii) considering that dysregulation of TGFβ-1 signaling promotes invasion and metastasis (30), we modified the E2F1 synthesis parameter and TGFβ-1 expression level to investigate whether they can synergize in the emergence of a resistance phenotype; and (iv) given that the expression of ERBB3 associates to the resistance to certain cytostatic drugs (20,31,32), we modified E2F1 and ERBB3 synthesis parameters to test whether their combined dysregulation can cause enhanced chemoresistance.

As shown in Fig. 4A, in silico cell lines with strong DNp73 and/or weak p73 synthesis display a very high percentage of cells resisting apoptosis after genotoxic drug administration, whereas strong p73 synthesis and reduced DNp73 induction provide efficient cell death. Furthermore, our kinetic simulations suggest that very strong E2F1 induction alone is sufficient to provoke chemoresistance to genotoxic drugs and can therefore compensate low DNp73 levels (Fig. 4B). In addition, high induction of TGFβ-1 reduces the level of E2F1 required to induce resistance. Taken together, strong E2F1 synthesis can synergize with strong induction of DNp73 and TGFβ-1 to provide genotoxic drug resistance (Fig. 4A). When analyzing the emergence of resistance to cytostatic drugs through the considered network, our analysis indicates that strong ERBB3 and E2F1 induction can synergize to achieve chemoresistance (Fig. 4B). Taken together, the results suggest that resistance to both genotoxic and cytostatic drugs can emerge with strong ERBB3 and E2F1 induction. We summarized our results in a qualitative manner in Fig. 4C, which lists possible genetic signatures conferring chemoresistance. We note that in our model, the regulation of antiapoptotic genes like BCL-2 and ERBB3 by E2F1, p73/DNp73, and TGFβ-1 is mediated via miR-205. Therefore, we hypothesize that the coregulation of miR-205 by its transcriptional regulation can foster chemoresistance (Fig. 4D).

It is a well-established fact that cell populations with different genetic signatures influencing cancer hallmarks coexist in malignant tumors (33). In line with this, we used our kinetic model to investigate how tumor heterogeneity (i.e., coexistence of several populations with different genetic signatures) affects the efficacy of conventional anticancer therapy (Fig 5). In our analysis, we focused on genotoxic drugs. Towards this end, we reformulated our model to consider two kinds of tumors: (i) tumors with homogeneous expression of E2F1 and (ii) tumors with heterogeneous expression of E2F1. In our simulations, heterogeneity in E2F1 expression is modeled by making 90% of cells with a defined E2F1 level, 5% cells with 0.5-fold downregulation, and 5% cells with 2.5-fold upregulation of E2F1. In accordance with our previous simulations, we considered two biologic scenarios: E2F1 overexpression (10¹, represented by E2F1▴) and extreme overexpression (10², represented by E2F1▴▴).

We used the model to simulate the effect of an anticancer treatment composed of periodic injections of a conventional genotoxic agent. In a tumor with homogenous E2F1 overexpression, our results indicate that the tumor is abolished after several injection cycles (Fig. 5 top left). However, in a tumor with homogenous E2F1 extreme overexpression, our simulation suggests that the tumor is resistant to the therapy and recovers fast growth afterwards (Fig. 5 top right). Interestingly, in a tumor with heterogeneous E2F1 overexpression, the tumor becomes virtually abolished after several cycles of drug injection, but the tumor reemerges after the drug injections are discontinued (Fig. 5 bottom left). Finally, for a tumor with heterogeneous E2F1 extreme overexpression, model simulations indicate that the tumor is fully drug resistant with continuous proliferation (Fig. 5 bottom right).

We further analyzed the scenario of heterogeneous E2F1 overexpression and found that in our simulations during the treatment, the composition of the tumor shifts from one of predominantly chemosensitive cells at time point zero (Fig. 6A, dark blue curve) towards a composition of entirely genotoxic-drug resistant cells at twenty weeks (Fig. 6A, red curve). This suggests that the genotoxic drug treatment can induce selective pressure on tumor cells. This is because chemosensitive cells undergo apoptosis during treatment, whereas cells with the resistant genotype withstand and continue proliferating. Therefore, the tumor reemerges after discontinuation of the treatment after 20 weeks (Fig. 6B). Our simulations indicate that during treatment the genetic signature of E2F1 and miR-205 in the tumor shift from that of chemosensitive cells to the one in chemoresistant cells (upregulation of E2F1 and downregulation of miR-205, Fig. 6C). Strikingly, we were able to confirm our prediction with experiments, in which we stimulated SK-Mel-147 cells with iteratively increasing levels of cisplatin (1–50 μmol/L) and measured E2F1 and miR-205 expression before and after the treatment with this genotoxic drug (respectively and in Fig. 6C; ref. 15).

We developed a kinetic model to investigate the emergence of chemoresistance in tumor cells, mediated by a network with a core module composed of E2F1, p73, DNp73, and miR-205. In previous efforts, mathematical modeling was applied to address the role of E2F1 and other network components in cancer and cell proliferation (24, 34–36). However, to our knowledge, this work is the first modeling-based analysis of E2F1-mediated chemoresistance.

Our results suggest that dysregulation of this network can provoke chemoresistance to multiple anticancer drugs. In line with this, our model unraveled several genetic signatures of tumor cells that can promote chemoresistance. Regarding genotoxic drugs, we found a signature of high E2F1 and low miR-205 expression that promotes resistance. Downregulation of miR-205 in this signature can be mediated by an imbalance in the p73/DNp73 ratio or by the dysregulation of other cancer-related pathways like TGFβ-1 signaling, whose components positively or negatively regulate miR-205 expression (37,38).

With respect to the imbalance in the p73/DNp73 ratio, our simulations predict that tumors with high E2F1 and normal/high DNp73 expression show a higher percentage of surviving cells after genotoxic stress than tumors expressing high E2F1 and low DNp73. This supports the idea that DNp73-mediated repression of miR-205 is a mechanism by which drug resistance is achieved. Because DNp73 is promoted by E2F1, high levels of E2F1 can induce chemoresistance via the E2F1-DNp73-miR-205 axis due to the negative regulation of miR-205 by DNp73 (Fig. 7A). However, DNp73 is upregulated in melanoma metastases with a concomitant rise of the full-length p73, a positive regulator of miR-205 (13). Taken together, this module has the structure of an incoherent feedforward loop. This suggests that the specific ratio between both p73 isoforms expression and activity determines the functional outcome of E2F1 and cell fate upon genotoxic treatment. This ratio is regulated in malignant melanoma by as yet unknown factors involved in alternative splicing of p73 pre-mRNA, possibly independent of E2F1 (39). In tumor cells with high E2F1 induction, our analysis indicates that downregulation of DNp73 and/or upregulation of p73 can restore chemosensitivity, which emphasizes the importance of the p73/DNp73 ratio.

When constructing the network, we found that some proteins involved in the epithelial-to-mesenchymal transition (EMT), like TGFβ-1, BMP, or TWIST1 can regulate miR-205 expression in some tumors (38,40). Furthermore, our simulations indicate that when E2F1 is upregulated, high levels of TGFβ-1 or other EMT-related signals regulating miR-205, external to the E2F1 network, can provoke genotoxic drug resistance in cancer cells. This suggests an intriguing hypothesis, in which E2F1 signaling cross talks and synergizes with some EMT-related signals, and their regulation of common miRNAs can mediate this synergy in the context of chemoresistance (Fig. 7B).

Moreover, our analysis gives additional insights into the resistance to cytostatic drugs that inhibit the activation of HER receptor family members. We found a genetic signature composed of high E2F1, low miR-205, and high ERBB3 that can render tumor cells insensitive to cytostatic drugs. As cited before, miR-205 posttranscriptionally represses ERBB3. Interestingly, recent results suggest that tumor cell lines sensitive to erlotinib show high levels of miR-205 and other miRNAs (41), which supports our predictions. According to our model analysis, tumor cells with E2F1 overexpression and miR-205 downregulation do not show cytostatic drug resistance, but synergy with high ERBB3 gene activation is required to acquire this feature. In addition, the simulations suggest that in most of the cases, this genetic signature can lead to double resistance, this means the tumor cells become insensitive to some genotoxic and cytostatic drugs at the same time, a feature that may have important consequences in the design of personalized cancer treatments. The underlying molecular mechanism providing this double resistance is to be further elucidated.

Taken together, our results indicate that differences in the expression of some of the proteins in this network are critical for tumor cells to become chemoresistant, and that this resistance is substantiated through regulatory loops including miR-205. We note that the predictions we made are based on qualitative modeling, similar to those made in the studies conducted by Aguda and colleagues and Yao and colleagues (28,39). To make the predictions of our model quantitative, the model should be further refined in iterative cycles of experimentation and model calibration (17,18).

In the model parameterization analyzed, the feedback loop between E2F1 and miR-205 has minor effects in the dynamics of the network. Our previous results of modeling miRNA regulation of cancer-related genes suggest that the strength of miRNA regulation is context-dependent and may be affected by the coregulation of other miRNAs or proteins (42). Thus, we hypothesize that this feedback loop could be negligible in the scenarios investigated in this article, but still crucial in other biologic contexts which have not been investigated in this paper. In this way, miR-205–mediated repression of some antiapoptotic targets in our network may be affected by similar coregulation. Our modeling results complement those of Aguda and colleagues (24) and others (43) who also used kinetic modeling to investigate the role of miRNA regulation on E2F1 in cancer. Aguda and colleagues (24) found that the coupling between the E2F/Myc positive feedback loop and the E2F/Myc/miR-17-92 negative feedback loop is crucial to regulate a bistable on–off switch in E2F/Myc protein levels, a feature that is critical in the dysregulation of E2F1 activity during cancer progression. We here claim that E2F1-mediated regulation of miRNAs may as well be critical to explain the emergence of chemoresistance.

Finally, our simulations indicate that tumor heterogeneity in terms of the genetic signature of the E2F1-p73/DNp73-miR-205 network can affect the efficacy of anticancer therapy. In addition, our simulations suggest that the therapy can induce selective pressure on tumor cells favoring those showing chemoresistance. Other authors have suggested that under conditions of tumor genetic instability, heterogeneous tumors contain populations of one or more resistant clones (44, 45). These resistant subpopulations can be residual compared with other clones in a pretreated tumor due to the phenotypic cost associated with chemoresistance, but could get favored when the anticancer therapy is applied (46). This hypothesis is supported by our results. In our analysis, we only considered the role of the E2F1-p73/DNp73-miR205 network in developing chemoresistance. However, these molecules interact with additional genes involved in the regulation of other phenotypic responses and therefore display pleiotropy (4,16). Thus, it is possible that some features of the resistant phenotype here described may have negative effect on other cancer-associated traits, which are out of the scope of this article. In a more realistic case, one could enrich the analysis by expanding the network with interaction partners representative of these other phenotypes and conducting a more extensive analysis.

Our analysis indicates that the network under investigation is enriched in feedforward loops. Most of these motifs are central to the efficient apoptosis initiation after genotoxic stress and therefore their dysregulation can cause chemoresistance (47,48). Previous publications have emphasized the important role of feedforward loops regulating the transient dynamics of biochemical pathways (49). We here support the hypothesis that these loops may further control important features of the long-term response in signaling and transcriptional systems. The paradigmatic case is the regulation of miR-205 by p73 (positive) and DNp73 (negative), two transcriptional targets of E2F1. This system is a canonical case of an incoherent feedforward loop. Several publications suggest that the regulation and activity of p73 and DNp73 is nonlinear and varies for different levels of E2F1 expression (28,50). Under these conditions, our simulations indicate that the incoherent feedforward loop allows a single protein (e.g., E2F1) to nonmonotonically regulate components downstream of the loop (e.g., miR-205; see Fig. 7A). Thus, for an interval of E2F1 expression levels the incoherent feedforward loop is dominated by the p73-mediated positive branch and promotes expression of miR-205. However, for higher expression of E2F1, the negative, DNp73-mediated branch predominates and the downstream target is repressed (Fig. 7A). This suggests that incoherent feedforward loops can be active structures of regulation when the different branches of the loops get distinctively regulated.

Furthermore, we identified several network motifs with the structure of incoherent feedforward loops, in which a component downstream the network is regulated by a direct E2F1 target gene and a miR-205–regulated target. In our study, we focus on apoptosis-related targets (e.g., BCL-2 and HRK, see Fig. 7B). As we mentioned before, miR-205 expression is regulated by other cancer-related signals, external to the E2F1 network (e.g., EMT-related signals). Therefore, changes in the activity of those external signals may alter the activation status of the antiapoptotic branch and as a consequence the balance between pro- and antiapoptotic signals is also altered (Fig. 7B). In addition, miR-205 regulates proteins downstream the EMT signaling, like E-cadherin, and some members of the BCL2 family are regulated by E2F1. This suggests the existence of additional feedforward loops which were not investigated in this work (25). Taken together, feedforward loops involved in chemosensitivity can be considered active regulatory structures, whose activation status may depend on the cross-talk with other cancer-related signals.

No potential conflicts of interest were disclosed.

Conception and design: J. Vera, D. Engelmann, B.M. Pützer

Development of methodology: J. Vera, X. Lai, F.M. Khan

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): U. Schmitz

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): J. Vera, U. Schmitz, X. Lai, D. Engelmann

Writing, review, and/or revision of the manuscript: J. Vera, U. Schmitz, X. Lai, D. Engelmann, F.M. Khan, O. Wolkenhauer, B.M. Pützer

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): U. Schmitz

Study supervision: B.M. Pützer

This work was supported by the German Federal Ministry of Education and Research (BMBF) as part of the projects eBio:miRSys (0316175A to J. Vera and O. Wolkenhauer), eBio:SysMet (0316171 to J. Vera, O. Wolkenhauer, and B.M. Pützer), GERONTOSYS-ROSAge (0315892A to J. Vera, O. Wolkenhauer), and Rostock University Medical Center for the project Systems Medicine of Cancer Invasion and Metastasis.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.